Entanglement in non-equilibrium steady states and many-body localization breakdown in a current driven system

TL;DR

This paper explores how entanglement properties and localization in a disordered 1D system change under non-Hermitian, current-driven conditions, revealing a transition from volume-law to area-law entanglement.

Contribution

It introduces a non-Hermitian model to study entanglement transitions in driven disordered systems, highlighting the impact of asymmetric hopping on many-body localization.

Findings

Identifies a transition from volume-law to area-law entanglement entropy.

Shows the existence of current-carrying non-equilibrium steady states.

Reveals rich entanglement structures in non-Hermitian eigenstates.

Abstract

We model a one-dimensional (1D) current-driven interacting disordered system through a non-Hermitian Hamiltonian with asymmetric hopping and study the entanglement properties of its eigenstates. In particular, we investigate whether a many-body localizable system undergoes a transition to a current-carrying non-equilibrium steady state under the drive and how the entanglement properties of the quantum states change across the transition. We also discuss the dynamics, entanglement growth, and long-time fate of a generic initial state under an appropriate time-evolution of the system governed by the non-Hermitian Hamiltonian. Our study reveals rich entanglement structures of the eigenstates of the non-Hermitian Hamiltonian. We find transition between current-carrying states with volume-law to area-law entanglement entropy, as a function of disorder and the strength of the non-Hermitian term.

Equations

$$\mathcal{H}$$

$$-J_z\sum _i^L{S}_i^z{S}_{i+1}^z.$$

$$\mathcal{H}$$

$$-V\sum _i^L{n}_i{n}_{i+1}.$$

$$\frac{d\rho }{dt}=-\mathfrak{i}[{\mathcal{H}}_h,\rho ]+\lambda \left(\{\mathcal{J},\rho \}-2\mathrm{Tr}(\rho \mathcal{J})\rho \right)\equiv \mathcal{L}\rho .$$

$$|\psi (t)⟩$$

$$|{\psi }_{\infty }⟩=\frac{⟨s_L|\psi ⟩}{|⟨s_L|\psi ⟩|}\frac{|s_R⟩}{\Vert |s_R⟩\Vert }$$

$$\frac{\partial \tilde{\rho }}{\partial t}$$

Full text

I Introduction

Classification of many-body quantum states in terms of their entanglement properties has become a major theme in condensed matter physics. These activities have revealed intriguing connections of entanglement with dynamics and thermalization of quantum systems Altman and Vosk (2015); Nandkishore and Huse (2015); Abanin et al. (2018). It has been understood that typical high-energy eigenstates of a generic isolated many-body system can be classified as either ergodic or non-ergodic depending on the subsytem entanglement entropy, e.g. $S_{EE}=-\mathrm{Tr}_{A}\rho_{A}\ln\rho_{A}$ of a subsytem $A$, where $\rho_{A}$ is the reduced density matrix of $A$ for the eigenstate. In particular, when $d$-dimensional system is divided into two subsytems with the smaller subsystem, say $A$, having length $L$, the entanglement entropy $S_{EE}\sim L^{d}$, i.e. scales with the volume of the subsystem. On the contrary, for many-body localized (MBL) states Gornyi et al. (2005); Basko et al. (2006); Oganesyan and Huse (2007); Pal and Huse (2010); Altman and Vosk (2015); Nandkishore and Huse (2015); Abanin et al. (2018), prominent examples of non-ergodic eigenstate, subsytem entanglement scales as $L^{d-1}$, i.e. with the area of the subsystem, the so-called area law scaling. Starting with a generic initial state, the MBL systems do not thermalize under the unitary dynamics, unlike the ergodic ones. Nevertheless, the MBL systems still give rise to a slow growth of entanglement entropy approaching a state with sub-thermal volume law scaling Žnidarič et al. (2008); Bardarson et al. (2012); Serbyn et al. (2013); Vosk and Altman (2013). One of the most natural realizations of MBL phases are found in strongly disordered systems where single-particle Anderson localization is stable to interaction Basko et al. (2006); Gornyi et al. (2005); Oganesyan and Huse (2007); Pal and Huse (2010), even at finite energy densities above the ground state. A many-body localizable system can undergo a volume-law (ergodic) to area-law (non-ergodic) transition Gornyi et al. (2005); Basko et al. (2006); Oganesyan and Huse (2007); Pal and Huse (2010); Altman and Vosk (2015); Nandkishore and Huse (2015); Abanin et al. (2018) as a function of the strength of the disorder or the interaction .

However, condensed matter systems are seldom isolated. In particular, some very important experimental setups require the system to be connected with external environment. A notable example is a system connected to leads and driven by current or a voltage bias or an electric field [see Fig.1(a)]. Generically, such systems are expected to attain an unique current-carrying non-equilibrium steady state (NESS). It is interesting to explore whether such NESSs could also be classified according to their entanglement content. For example, we could ask whether a many-body localizable system can undergo a dynamical transition when driven by a strong current or electric field, e.g. a transition between NESSs with distinct system-size scaling of current. Moreover, can such NESS transition be coincident with an entanglement transition?

Recently there have been some studies of entanglement properties Sharma and Rabani (2015); Dey et al. (2019) as well as entanglement transitions Gullans and Huse (2019, 2019) for the current carrying NESSs in non-interacting models using non-equilibrium Green’s functions or scattering states. However, computing entanglement properties of driven states of a disordered interacting system connected to infinite leads at the boundaries is an extremely challenging task. One possible way to access such boundary driven system is via Markovian Lindblad quantum master equation approximation Prosen and Žnidarič (2009); Žnidarič (2014); Mahajan et al. (2016); Zanoci and Swingle (2016). Nevertheless, generically, such Markovian evolution is destined to lead to a description of NESS in terms of mixed state having area law entanglement, e.g. quantified in terms of mutual information Wolf et al. (2008); Mahajan et al. (2016); Zanoci and Swingle (2016), and hence make the notion of entanglement transition obscure. Also, in non-interacting systems, there are examples of NESSs with higher than area-law mutual information Eisler and Zimborás (2014), or even volume-law entanglement Gullans and Huse (2019, 2019) . Such highly entangled NESSs can persist even in interacting systems with size smaller than electronic phase coherence length Gullans and Huse (2019). An important avenue, at present, is to try to mimic such NESS in less microscopic toy models, e.g. boundary driven random unitary circuit model Gullans and Huse (2019). Here we take a different route, and try to address driven states of interacting disordered systems through an effective model which incorporates the dissipation and the current drive via a non-Hermitian term.

The model studied here is an interacting version of the one dimensional Hatano-Nelson model Hatano and Nelson (1996, 1997, 1998) for non-Hermitian single-particle localization-delocalization transition. The model has a time ($\mathcal{T}$)-reversal symmetry, that could be spontaneously broken by the eigenstates for a large enough strength of the non-Hermitian term. In the non-interacting Hatano-Nelson model, $\mathcal{T}$-reversal symmetry breaking or real-to-complex eigenvalue transition Hatano and Nelson (1996, 1997, 1998) coincides with localization-delocalization transition. We study the interacting non-Hermitian model via exact diagonalization (ED). Our main results are the following.

(1) We show that a current-driven many-body localizable system can undergo a volume- to area-law entanglement transition in the eigenstates as a function of disorder and/or the strength of the non-Hermitian term for a fixed interaction strength.

(2) We find that entanglement transition has a direct correspondence with a transition, similar to MBL transition in Hermitian system Luitz et al. (2015); Macé et al. (2019), in terms of participation of the eigenstates in the many-body Hilbert space.

(3) We demonstrate that the entanglement and $\mathcal{T}$-reversal breaking transitions are distinct in the interacting case. We also find another distinct transition within the $\mathcal{T}$-reversal broken region in terms of the scaling of the current with system length.

(4) We show how the system approaches the NESS at long times, under an appropriately defined dynamics, in terms of time-evolution of the entanglement entropy and the memory of the initial state. In particular, as in the Hermitian MBL systems, starting with an initial unentangled (product) state the entanglement entropy grows linearly with time $t$ in the volume-law phase. In contrast, both in the $\mathcal{T}$-reversal broken and unbroken area-law phase, entanglement entropy grows as $\ln(t)$. However, the entanglement growth is followed by a decay at late times towards an unique NESS. Moreover, we find that the memory of the initial state always eventually gets lost as the NESS is attained at long times.

Overall, our study reveals a much richer eigenstate and dynamical phase diagram of the non-Hermitian model compared to the Hermitian systems.

The remainder of this paper is organized as follows. In Section II we describe the model and the dynamics governed by the non-Hermitian Hamiltonian. We give some physical motivations behind the model and its dynamics in Section III. Section IV discusses the results for the Hilbert-space localization-delocalization, entanglement and time-reversal symmetry breaking transitions in terms of finite-size scaling analysis. We also describe here the transient dynamics of the system during the approach to NESS and the properties of the NESS. In Sec. V we conclude with the implications and significance of our results. The appendix include some additional results.

II Model

In this section we introduce the non-Hermitian Hamiltonian and the dynamics. In the next section, we briefly review general motivations for the model, based on non-Hermitian approach to open systems, and discuss some specific physical realizations.

II.1 Non-Hermitian Hamiltonian

We study the following non-Hermitian one-dimensional ($1D$) $XXZ$ spin ($S=1/2$) model with an uniform random field $h_{i}\in[-W,W]$, $W$ being the disorder strength,

[TABLE]

Here $S_{i}^{\pm}=(S_{i}^{x}\pm iS_{i}^{y})$, $L$ is the number of sites, $\Psi$ is a real number (see below) and $J,~{}J_{z}$ are the spin exchanges; the latter controls the interaction strength. We apply periodic boundary condition. The above model can be rewritten as $\mathcal{H}=\mathcal{H}_{h}+\mathfrak{i}\lambda\mathcal{J}$, where $\mathcal{H}_{h}$ is the usual Hermitian random field $XXZ$ model with $J\to\tilde{J}=J\cosh{\Psi}$, $\lambda=J\sinh{\Psi}$ and $\mathcal{J}=-(\mathfrak{i}J/2)\sum_{i}(S_{i}^{+}S_{i+1}^{-}-S_{i}^{-}S_{i+1}^{+})$ is the sum of spin current across the system.

The model can be mapped, via Jordan-Wigner transformation, to a model of spin-less fermions hopping (with amplitude $t=J/2$) on the $1D$ lattice with random disorder potential ($\epsilon_{i}=h_{i}$) and nearest neighbor interaction ($V=J_{z}$), i.e.

[TABLE]

Here $c_{i},c_{i}^{\dagger}$ are fermion operators. In this case, $\mathfrak{i}\Psi$ is an imaginary vector potential corresponding to an imaginary flux $\mathfrak{i}\Psi L$ through a ring of circumference $L$ [Fig.1(b)] and the model is an interacting version of Hatano-Nelson model Hatano and Nelson (1996, 1997, 1998). The latter describes non-Hermitian single-particle localization-delocalization transition. The spinless fermion model is invariant under $\mathfrak{i}=\sqrt{-1}\to-\mathfrak{i}$, i.e. time ($\mathcal{T}$) reversal. However, as in the usual $\mathcal{PT}$-symmetric non-Hermitian models Bender and Boettcher (1998); Bender (2007), this $\mathcal{T}$- or pseudo $\mathcal{PT}$-symmetry can be broken by the eigenstates leading to complex eigenvalues. We refer to this real to complex transition as $\mathcal{T}$-reversal breaking for brevity even in the spin model [Eq.(1)]. Intriguingly, in the original Hatano-Nelson model the $\mathcal{T}$-symmetry breaking transition of the single-particle eigenstates coincides with the localization-delocalization transition and the delocalized states carry finite current. This leads to the question whether there is any localization-delocalization transition in the interacting model and whether the transition coincides with the $\mathcal{T}$-reversal breaking. Such congruence of the symmetry breaking and localization transition might lead to an exciting possibility of describing a “MBL transition” in terms of symmetry breaking.

II.2 Non-Hermitian dynamics

We further extend the model of Eq.(1) to describe the approach to the long-time NESS through the following dynamical equation for density matrix $\rho$ ($\hbar=1$),

[TABLE]

Here $\mathcal{L}$ defines a Liouvillian operator. As we discuss in the next section, the above dynamics for the density matrix can be obtained starting with the time-dependent non-Hermitian Schr$\ddot{\mathrm{o}}$dinger equation Sergi and Zloshchastiev G. (2013); Zloshchastiev and Sergi (2014). The dynamical evolution of Eq.(3) has been used previously in the context of $\mathcal{PT}$-symmetric quantum mechanics to describe system with gain and loss Brody and Graefe (2012). A similar model has been also used in the context of field-driven Mott transition Tripathi et al. (2016). Eq.(3) can describe the evolution of both pure and mixed states as discussed in Ref.Brody and Graefe, 2012 and keep $\mathrm{Tr}\rho(t)=1$. Unlike in a Lindblad master equation, an initial pure state remains pure during the time evolution, and, in this case, Eq.(3) reduces to a complex Schr$\ddot{\mathrm{o}}$dinger equation Brody and Graefe (2012). From Eq.(3), $\rho(t)$ can be formally written using the non-Hermitian Hamiltonian [Eq.(1)] as $\rho(t)=e^{-\mathfrak{i}\mathcal{H}t}\rho_{0}e^{\mathfrak{i}\mathcal{H}^{\dagger}t}/\mathrm{Tr}(e^{-\mathfrak{i}\mathcal{H}t}\rho_{0}e^{\mathfrak{i}\mathcal{H}^{\dagger}t})$, where $\rho_{0}$ is the initial density matrix. In particular, for an initial pure state $\rho_{0}=|\psi\rangle\langle\psi|$, the state at time $t$ is given by

[TABLE]

Here $\langle n_{L}|$ and $|n_{R}\rangle$ are the left and right eigenvectors of $\mathcal{H}$ with eigenvalue $E_{n}=\mathcal{E}_{n}+\mathfrak{i}\Lambda_{n}$; $\|.\|$ denotes the norm of a vector. As well known Faisal and Moloney (1981), the left and right eigenvectors form a bi-orthonormal basis, $\langle n_{L}|m_{R}\rangle=\delta_{nm}$ with the resolution of the identity $\sum_{n}|n_{R}\rangle\langle n_{L}|=\mathbb{I}$. The real part of the eigenvalue $\mathcal{E}_{n}=\langle n_{R}|\mathcal{H}_{h}|n_{R}\rangle/\langle n_{R}|n_{R}\rangle$, i.e. the expectation of the parent Hermitian Hamiltonian, and the imaginary part $\Lambda_{n}\propto\mathcal{J}_{n}=\langle n_{R}|\mathcal{J}|n_{R}\rangle/\langle n_{R}|n_{R}\rangle$. Here the bra vector $\langle n_{R}|=|n_{R}\rangle^{\dagger}$ is the Hermitian conjugate of the right eigenket, and similarly $|n_{L}\rangle=(\langle n_{L}|)^{\dagger}$. The identifications of real part and imaginary parts of the eigenvalue as the eigenstate expectation value of $\mathcal{H}_{h}$ and $\mathcal{J}$ are due to the fact that both the operators are Hermitian and hence their expectation values are real. It can be easily shown from Eq.(4) that, for eigen spectrum with non-zero imaginary part of the eigenvalues, the NESS $|\psi_{\infty}\rangle\equiv|\psi(t\to\infty)\rangle$ is given by

[TABLE]

where $\langle s_{L}|(|s_{R}\rangle)$ is the left (right) eigenvector with maximum imaginary part of the eigenvalue $\Lambda_{s}$.

III Motivations and physical realizations

III.1 Non-Hermitian approach to open quantum systems

Non-Hermitian Hamiltonians have been used in many past studies to describe decaying states (e.g. see Refs.Feshbach, 1958, 1962; Cohen-Tannoudji, 1968; Faisal and Moloney, 1981) and effects of dissipation and drive in open quantum systems (e.g. see Refs.Sergi and Zloshchastiev G., 2013; Zloshchastiev and Sergi, 2014; Rotter, 2009; Brody and Graefe, 2012; Tripathi et al., 2016; Fukui and Kawakami, 1998; Oka and Aoki, 2010). For the latter systems, the non-Hermitian approach has range of applicability different and complementary to more widely used Lindblad master equation approximation Weiss (2012). For instance, as already discussed, the non-Hermitian approach could describe both pure and mixed states and their dynamical evolutions Brody and Graefe (2012); Sergi and Zloshchastiev G. (2013); Zloshchastiev and Sergi (2014). The general theoretical framework, in principle, to rigorously derive such non-Hermitian Hamiltonian is rooted in the Feshbach projection operator technique Feshbach (1958, 1962); Cohen-Tannoudji (1968); Faisal and Moloney (1981); Rotter (2009). The latter can be used to construct non-Hermitian Hamiltonian to describe the effects of continuum of scattering states in the presence of an environment on a finite subsystem which has discrete states in isolation. However, here, as in earlier studies Tripathi et al. (2016); Fukui and Kawakami (1998); Oka and Aoki (2010), we use the model of Eq.(1) from a more phenomenological ground to describe NESS in a current driven interacting disordered open system. Apart from the connections to some of the possible physical realizations that we discuss below, the model of Eq.(1) could be thought of as an effective model to generate and study interesting current carrying pure states which are relatively simpler to access within exact numerical calculations, albeit for finite systems of modest sizes. Also, in the spirit of study of current carrying states in quantum field theories (see e.g. Ref.Cardy and Suranyi, 2000), the non-Hermitian term $i\lambda\mathcal{J}$ can also be interpreted as a Lagrange multiplier term that acts as a source field for inducing finite current in the eigenstates.

Once the non-Hermitian Hamiltonian of Eq.(1) is posited as an effective phenomenological model to describe current driven open system, the dynamical evolution [Eq.(3)] results naturally from the corresponding non-Hermitian time-dependent Schr$\ddot{\mathrm{o}}$dinger equation Faisal and Moloney (1981) as discussed in Ref.Sergi and Zloshchastiev G., 2013; Zloshchastiev and Sergi, 2014. We briefly review here the derivation of Eq.(3) for the sake of completeness.

For the non-Hermitian Hamiltonian the time-dependent Schr$\ddot{\mathrm{o}}$dinger equation, $\mathfrak{i}\partial|\tilde{\psi}(t)\rangle/\partial t=(\mathcal{H}_{h}+i\lambda\mathcal{J})|\tilde{\psi}(t)\rangle$ for a state $|\tilde{\psi}(t)\rangle$, leads to

[TABLE]

for a density matrix $\tilde{\rho}(t)=|\tilde{\psi}(t)\rangle\langle\tilde{\psi}(t)|$. As discussed in Refs.Sergi and Zloshchastiev G., 2013; Zloshchastiev and Sergi, 2014, the above evolution does not preserve $\mathrm{Tr}\tilde{\rho}(t)$. However, the expectation values of an operator can be obtained using the normalized density matrix $\rho(t)=\tilde{\rho}(t)/\mathrm{Tr}\tilde{\rho}(t)$. It is straightforward to obtain the time-evolution of Eq.(3) for the normalized density matrix $\rho$, and hence for the normalized pure state $|\psi(t)\rangle$ [Eq.(4)] with $\rho(t)=|\psi(t)\rangle\langle\psi(t)|$, from Eq.(6). Therefore, we use Eq.(3) to describe the dynamics in our model.

III.2 Physical realizations

1. Electric field-driven breakdown in interacting disordered systems: One of the motivations for the models of Eq.(1) and Eq.(2) comes from earlier works in Refs.Fukui and Kawakami, 1998; Oka and Aoki, 2010; Tripathi et al., 2016 on the electric field driven Mott transition. In this case, the imaginary gauge potential in Eq.(2) arises from a asymmetric hopping amplitudes, $t_{1L}=-te^{L\Psi}\neq t_{L1}=-te^{-L\Psi}$, between site 1 and $L$ on the ring [Fig.1(b)], while all other hoppings between nearest-neighbor sites are $-t$. It was argued in Ref.Fukui and Kawakami, 1998 that these asymmetric hoppings mimic the situation of an open chain, connected to leads at the two ends, where electron is supplied at site 1 and dissipated at site $L$ at the other boundary. Thus the asymmetric hoppings give rise to both dissipation and drive. The model of Eq.(2) can then be easily obtained via non-unitary gauge transformations, $c_{i}^{\dagger}\to e^{-i\Psi}c_{i}^{\dagger}$ and $c_{i}\to e^{i\Psi}c_{i}$ at site $i$, and $\mathfrak{i}\Psi L$ acts as an imaginary flux through the ring [Fig.1(b)].

In the case of electric field-driven Mott localization-delocalization (insulator-to-metal transition), the model with imaginary gauge potential has been shown Oka and Aoki (2010) to describe many-body Landau-Zener (LZ) Landau and Lifshitz (1958) quantum tunneling processes near field driven Mott transition within Dykhne Dykhne (1962) formalism. In the latter approach, a model with real gauge potential, such as for a constant electric field, is analytically continued to imaginary gauge potential to describe LZ transition.

Similar to the electric-field driven Mott transition in the clean fermionic Hubbard model Fukui and Kawakami (1998); Oka and Aoki (2010); Tripathi et al. (2016), the non-Hermitian model can be thought of as an effective model to describe possible electric-field driven transition in a MBL system through many-body Landau-Zener breakdown. As in the case of field-driven Mott transition Oka et al. (2003), the connection of the non-Hermitian model with field-driven MBL breakdown can be explored by considering interacting fermions on a ring [Fig.1(b)] under a constant electromotive force applied by a time-dependent (real) flux through the ring. It is indeed an interesting question whether such many-body LZ processes could provide an effective description for field or current driven transition in a MBL system. Here we do not attempt to address this issue, and, instead, use the non-Hermitian model [Eq.(1)] as an effective model to generate current driven pure states and describe possible entanglement transitions among them.

2. Anomalous boundary states of symmetry protected gapped topological (Hermitian) systems: Recently, new correspondence has been established between topological classification of $(d+1)$ dimensional gapped Hermitian systems and $d$ dimensional non-Hermitian systems Lee et al. (2019). For example, as discussed in Ref.Lee et al., 2019, the long-time steady state of non-interacting clean non-Hermitian model in Eq.(2) with $W=0$ and $V=0$ describes 1D chiral fermions, that are naturally realized at the edge of quantum Hall systems. Hence, the model of Eq.(2) could describe the effects of disorder and interaction on such chiral edge states.

3. Depinning transition in classical systems of interacting vortices: The original motivation behind the Hatano-Nelson model Hatano and Nelson (1996, 1997, 1998) was to study a classical systems of superconducting vortices with columnar defects and tilted magnetic field applied at an angle with columnar defects in one higher dimension than the quantum model. The latter can be mapped to the classical model via standard path integral approach, and the depinning transition is understood as the localization-delocalization transition in the quantum model. In the same spirit, the interacting model of Eq.(2) may describe the effects of inter-vortex interaction on the depinning transition Hamazaki et al. (2019).

IV Results

We discuss our results on the phases and phase diagram for the model of Eq.(1) in this section. To this end, we obtain the eigenvalues and the (left and right) eigenvectors of the non-Hermitian Hamiltonian [Eq.(1)] using numerical exact diagonalization for system sizes $L=10-16$ and sample over many disorder realizations (10000 for $L=10$, 6000 for $L=12$, 1440 for $L=14$, and 200 for $L=16$). The Hamiltonian [Eq.(1)] and the dynamics [Eq.(3)] conserves $S_{z}^{tot}=\sum_{i}S_{i}^{z}$. Hence we work in the $S_{z}^{tot}=0$ subspace. We expect the results to be similar for the other $S_{z}^{tot}$ sectors, however the effect of interaction is strongest for $S_{z}^{tot}=0$ subspace, which corresponds to half filling in the fermion sector. We take the interaction $J_{z}=J=1$ and vary $W$ and $\Psi$.

As in the studies of MBL systems Oganesyan and Huse (2007); Pal and Huse (2010); Altman and Vosk (2015); Nandkishore and Huse (2015); Abanin et al. (2018); Luitz et al. (2015); Macé et al. (2019), here we characterize the phases as a function of $W$ and $\Psi$ in terms of (a) properties of eigenstates and (b) from time evolution starting from a generic initial pure state. Moreover, to draw the eigenstate phase diagram we use the quantity, $\epsilon=(\mathcal{E}-\mathcal{E}_{0})/(\mathcal{E}_{M}-\mathcal{E}_{0})$, defined from the real part of the eigenvalues; $\mathcal{E}_{M}$ and $\mathcal{E}_{0}$ are the maximum and the minimum values of $\mathcal{E}$, respectively. For brevity, we refer to $\epsilon$ as energy density. The quantity $\epsilon$ is used as an index for the eigenstates. In the Hermitian limit $\Psi=0$, since energy is conserved under time evolution, the expectation value of energy density $\epsilon$ of the initial state can be used to characterize dynamics and long-time steady state ensemble. For example, in the ergodic phase of the Hermitian model, within eigenstate thermalization hypothesis (ETH)Deutsch (1991); Srednicki (1994, 1999), the long-time steady state can be related to a thermal Gibbs ensemble at a temperature $T$ with the thermal expectation value of energy density $\epsilon(T)=\epsilon$. Similarly, this relation can be used to define a temperature for an eigenstate using its energy density. As a result, an eigenstate phase transition at the many-body mobility edge ($\epsilon_{c}$) from ergodic to non-ergodic or MBL phase can be defined as a function of $\epsilon$ Basko et al. (2006); Gornyi et al. (2005); Oganesyan and Huse (2007); Pal and Huse (2010); Altman and Vosk (2015); Nandkishore and Huse (2015); Abanin et al. (2018); Luitz et al. (2015); Macé et al. (2019). In the non-Hermitian case the energy density $\epsilon$ is not conserved under the dynamics [Eq.(3)] and for evolution from a generic initial pure state the energy of the long-time steady state is uniquely determined by the energy density of the eigenstate, $|s_{R}\rangle$, as can be easily seen from Eq.(5). Hence $\epsilon$ cannot be related with a temperature. However, $\epsilon$ can still be used to denote the spectrum of states and eigenstate transitions, e.g. $\mathcal{T}$-reversal breaking and localization-delocalization transition or a many-body mobility edge, analogous to the single-particle mobility edge in the non-interacting Hatano-Nelson model Hatano and Nelson (1996, 1997, 1998). Furthermore, states other than $|s_{R}\rangle$ are accessed during the short and intermediate-time dynamics and they govern the approach towards the steady state as we discuss later.

Following the above considerations, we construct the eigenstate phase diagram in the $W-\epsilon$ plane for a few values of $\Psi$ based on – (a) participation entropy Luitz et al. (2015); Macé et al. (2019), which quantifies the extent of delocalization of an eigenstate in the Hilbert space for a chosen basis, as discussed in the next section, (b) entanglement entropy, and (c) fraction of eigenstate with non-zero imaginary part of eigenvalue. The last one serves as an order parameter for $\mathcal{T}$-reversal symmetry breaking in the eigenstates. In addition, we further quantify the $\mathcal{T}$-reversal breaking via scaling of current carried in the eigenstate with system size. As mentioned in the introduction, we find eigenstate transitions in terms of all these quantities. We assume the transitions to be continuous as in the numerical studies of MBL transition for the Hermitian case and corroborate the assumption via finite-size scaling collapse of the data of various quantities. However, we note that, even in the Hermitian case, the exact nature of MBL transition is currently debated Vosk et al. (2015); Potter et al. (2015); Dumitrescu et al. (2017); Thiery et al. (2018); Goremykina et al. (2019); Roy et al. (2019a, b), and cannot be ascertained within the modest system sizes accessed via ED studies.

To characterize the phases based on time evolution, we look into the long-time steady state and the time evolutions of entanglement entropy and Néel order parameter, which serves as the MBL order parameter Schreiber et al. (2015), starting from an initially un-entangled Néel state of staggered arrangement of spins. We describe our results in detail in the next sections.

IV.1 Ergodic to non-ergodic eigenstate transition in the Hilbert space

To characterize the eigenstates of the non-Hermitian Hamiltonian [Eq.(1)], we first obtain a phase diagram in terms of participation entropy, a diagnostic of ergodicity of the eigenstates. It’s defined as $S_{P}(|n_{R}\rangle)=-\sum_{\alpha}p_{\alpha}^{(n)}\ln{p_{\alpha}^{(n)}}$ for the eigenstates in the basis of the spin configurations, $|\alpha\rangle=|S_{1}^{z},S_{2}^{z},\dots,S_{N}^{z}\rangle$ Luitz et al. (2014, 2015) with $p_{\alpha}^{(n)}=|\langle\alpha|n_{R}\rangle|^{2}$.

In the delocalized phase, eigenstates are ergodic and has support over a finite fraction of sites in the Hilbert space, i.e. $\langle\alpha|n_{R}\rangle\sim 1/\sqrt{D_{H}}$ as required by normalization, and hence $S_{P}\simeq\ln{D_{H}}$, where $D_{H}$ is the dimension of the Hilbert space in the $S_{z}^{tot}=0$ subspace. On the other hand, in the MBL phase, the eigenstates are expected to exhibit a fractal character, i.e. a delocalized but non-ergodic behavior, with support over exponentially large number but vanishing fraction of sites in the Hilbert space, i.e. $\langle\alpha|n_{R}\rangle\sim 1/D_{H}^{a/2}$ or $S_{P}\simeq a\ln{D_{H}}$ with $a<1$ Macé et al. (2019); Logan and Welsh (2019). We obtain the disorder averaged $S_{P}(\epsilon)$ as a function $W$ and $\epsilon$ and plot the slope of $S_{P}$ with $L$ ($\sim\ln{D_{H}}$) in Fig.2 (a-c) for $\Psi=0,0.3$ and 0.6. To define $S_{P}(\epsilon)$, and other quantities discussed in the next sections, as a function of $\epsilon$, we have binned $\epsilon$ into intervals containing 100 or more eigenstates for a given disorder realization. Indeed, consistent with the Hermitian MBL case Luitz et al. (2015); Macé et al. (2019) in Fig.2 (a), we find a transition in the Hilbert space in terms of eigenstate participation for $\Psi=0.3,0.6$ Fig.2 (b),(c). The transition shifts to higher disorder for increasing $\Psi$. Hence, a non-zero current drive causes a breakdown of the MBL states of the parent Hermitian model over certain regime of disorder.

IV.2 Entanglement transition

We obtain the von Neumann entanglement entropy $S_{EE}$ for each the eigenstates, i.e. $S_{EE}(|n_{R}\rangle)=-\mathrm{Tr}(\rho_{A}\ln\rho_{A})$. Here $\rho_{A}=\mathrm{Tr}_{B}|n_{R}\rangle\langle n_{R}|$ is the reduced density matrix of the subsystem $A$, for the real-space bipartition of the system into left half, $A$, and right half, $B$. As already discussed in the introductory section, the MBL transition is defined from volume-law to area-law transition Luitz et al. (2015) in the Hermitian case, i.e. $\Psi=0$. A similar transition is observed for $\Psi=0.3$ with increasing disorder as shown in Fig.3(a), where disorder averaged $S_{EE}$ is plotted as a function of length $L$ at $\epsilon=0.5$, at the middle of the spectrum, for a range of disorder strength. As can be seen, for large values of disorder $S_{EE}$ tends to a constant value at larger $L$, whereas $S_{EE}$ increases linearly with $L$ for weaker disorder strength. The presence of a transition is evident from a clear crossing of $S_{EE}/L$ vs. $W$ curves for different $L$ in Fig.3(b). The dependence of $S_{EE}$ on $L$ is consistent with a volume-law scaling for $W<W_{c1}$, and an area-law scaling for $W>W_{c1}$. This crossing point is used to find out the critical disorder $W_{c1}\approx 3.6$, implying an entanglement transition, similar to the Hermitian case Luitz et al. (2015).

The transition is further corroborated by a reasonably good data collapse in Fig.4(a),(b) for $\Psi=0.3,~{}0.6$, obtained using the finite-size scaling ansatz $S_{EE}(L)=Lg(L^{1/\nu_{1}}(W-W_{c1}))$ Luitz et al. (2015), where $g(x)$ is the scaling function. The finite-size sclaing form assumes a volume-law scaling at the critical point, as in the case of Hermitian model Luitz et al. (2015). We find that the data for $S_{EE}/L$ at $\epsilon=0.5$ and $\Psi=0.3$ [Fig.4(a)] for various $W$ and $L$ can be collapsed quite well into two universal curves for $W<W_{c1}$ and $W>W_{c1}$ with a critical exponent $\nu_{1}\simeq 1$ and $W_{c1}=3.6\pm 0.1$. The two curves, one below and another above the transition, result from the sign of $W-W_{c1}$ in the argument of the scaling function $g(x)$. We have done similar analysis for the entire spectrum of $\epsilon$. The crossing of the curves and the data collapse are most prominent near the middle of the spectrum. The finite-size data collapse of $S_{EE}/L$ at $\epsilon=0.5$ for $\Psi=0.6$ is shown in Fig.4(b). For this larger value of $\Psi$, we find $W_{c1}=4.2\pm 0.1$, and $\nu_{1}=1.5$, i.e. the extracted critical exponent changes with $\Psi$.

Here it is important to note that, as in the exact diagonalization studies of Hermitian MBL systems Luitz et al. (2015), for $S_{EE}(L)$ and the other quantities discussed in the next sections, neither the scaling function $g(x)$ nor the critical exponents can be extracted very accurately from finite-size scaling limited to such small systems, and with only a few system sizes accessible via ED. However, the scaling function $g(x)$, as such, is not required to obtain the data collapse. The collapse is obtained by taking all the data as function of $W$ and $L$, e.g. from Fig.4(b) (inset), and plotting as a function of the $L|W-W_{c1}|^{\nu_{1}}$ (or $L^{1/\nu_{1}}|W-W_{c1}|$) with appropriate choice of $W_{c}$ and $\nu$ that generate a good scaling collapse. Additionally, we have also used a polynomial for $g(x)$ and fitted all the data as a function of $L^{1/\nu_{1}}(W-W_{c1})$ with the polynomial. This gives results consistent with that obtained using the simpler data collapse mentioned above. We do not perform a detailed error analysis for the scaling collapse. The main purpose of the finite-size scaling here is to bring out the existence of the transition, rather than estimating the critical exponents. The errorbars in $W_{c1}$ is estimated from the range of values of $W_{c1}$ for which the data collapse is reasonable while keeping the critical exponent $\nu_{1}$ fixed. The finite-size scaling analysis for the other quantities discussed below are also performed in a similar manner.

For the the phase boundaries $W_{c1}$ in the $W-\epsilon$ plane for $\Psi=0,0.3,0.6$ in Fig.2 we use the standard deviation, $\sigma_{E}$, of the entanglement entropy over disorder realizations. At the transition, $\sigma_{E}$ is expected to show a peak that diverges with $L$ Kjäll et al. (2014). In Fig.2, we plot $W_{c1}$ from the peak of $\sigma_{E}$ for system size $L=16$. The phase boundary is consistent with that obtained from the participation entropy for ergodic to non-ergodic transition. In Figs.5(a) and (b), we obtain scaling collapses of $\sigma_{E}/(L-c)$ vs. $W$ for $\Psi=0.3$ and $\Psi=0.6$ at $\epsilon=0.5$ with $c$ as a fitting parameter Luitz et al. (2015). Here the values of $\nu_{1}$ and $W_{c1}$ are fixed from the finite-size scaling of $S_{EE}$. Surprisingly we see an indication of two peaks in $\sigma_{E}$ for $L=16$ [Fig.5(b)]. The reason behind the two-peak structure is not clear at present and will be studied in a future work.

We note that the average $S_{EE}$ and $\sigma_{E}$ shown in Fig.4 and Fig.5 are obtained by averaging over the eigenstates with real eigenvalues. The averaging over the eigenstates with complex eigenvalues also gives similar results.

IV.3 $\mathcal{T}$-reversal breaking:

We now address one of the main questions raised in the introductory section, i.e., whether the entanglement or the localization transition coincides with the $\mathcal{T}$-reversal breaking, as in the non-interacting Hatano-Nelson model Hatano and Nelson (1996, 1997, 1998). We define a $\mathcal{T}$-reversal order parameter, $\phi(\epsilon)$, the fraction of imaginary eigenvalues at $\epsilon$. Numerically an eigenvalue is defined to have a non-zero imaginary part by setting an ad-hoc small cutoff. However, we have checked that our results are insensitive to the choice of cutoff. We find a clear crossing of the $\phi$ vs. $W$ curves for different $L$, as shown in Fig.6(a),(b) (insets) for $\Psi=0.3$ and $\Psi=0.6$ at $\epsilon=0.5$. As evident, the crossing point is at $W=W_{c2}\approx 4.75$ for $\Psi=0.3$, clearly larger than $W_{c1}$ for the entanglement transition in Fig.4(a). A good scaling collapse can again be obtained with an exponent $\nu_{2}\simeq 0.9$ and $W_{c2}=4.75\pm 0.1$ [Fig.6(a)]. The collapse of the data for $\phi$ vs. $W$ could not be obtained with $W_{c2}=W_{c1}$. This establishes the fact that $\mathcal{T}$-reversal breaking is distinct from the entanglement transition and occurs within the area-law phase. We find similar results for $\Psi=0.6$ [Fig.4(b) and Fig.6(b)], namely $W_{c2}=5.6\neq W_{c1}$. In this case, the critical exponent $\nu_{2}=1.5$. From the crossing points of $\phi(\epsilon)$ vs. $W$ curves we obtain the phase boundary $W_{c2}$ for the $\mathcal{T}$-reversal breaking in Figs.2.

To further characterize the $\mathcal{T}$-reversal breaking eigenstates, we compute the current $j(\epsilon)$, obtained by averaging over the magnitude of the currents, $|j_{n}|$, carried by the eigenstates with imaginary eigenvalues, and disorder realizations. Again, we find crossings, $W=W_{c3}\approx 2.3$ for $\Psi=0.3$ and $W=W_{c3}\approx 3.0$ for $\Psi=0.6$, in the $\mathcal{J}=Lj$ vs. $W$ plots for different $L$ [Figs.7(a),(b)] at $\epsilon=0.5$. The transition, which we refer to as current transition for brevity, is seemingly distinct from both entanglement and $\mathcal{T}$-reversal breaking transitions. We can obtain a scaling collapse for $j$ with $\nu_{3}=1$ and $W_{c3}=2.3$ for $\Psi=0.3$, and $\nu_{3}=1.5$ and $W_{c3}=3.0$ for $\Psi=0.6$, as shown in Figs.7(a),(b)(inset). We obtain a phase boundary for the current transition at other values of $\epsilon$ from the crossing points, as shown in Figs.2(b),(c).

The scaling of $j$ with $L$ for $W\ll W_{c3}$ is consistent with $j$ approaching a constant for $L\to\infty$ (not shown). The scale-invariant crossing point indicates a diffusive scaling of the current at the transition, namely $j\sim 1/L$. In fact, for W_{c3}<W\ {\raise-2.15277pt\hbox{\buildrel<\over{\sim}}}\ W_{c2}, the scaling of $j$ with $L$ could be consistent with $j\sim 1/L^{\gamma}$ with \gamma\ {\raise-2.15277pt\hbox{\buildrel>\over{\sim}}}\ 1, and we expect $j\sim e^{-L/\zeta}$ for $W\gg W_{c2}$, deep inside the $\mathcal{T}$-reversal unbroken localized phase; $\zeta$ is the characteristic localization length. However, this is hard to verify from the exact diagonalization numerics limited to such small system sizes. As discussed later [see Fig.11], the current, $j_{\infty}$, carried by the long-time NESS also exhibits similar transition. For $W>W_{c1}$, in the area-law phase, the scaling of the current and the transition at $W_{c2}$ could be studied by a matrix-product operator (MPO) based implementation Prosen and Žnidarič (2009); Žnidarič (2014) of the dynamical evolution in Eq.(3). It would be interesting to establish the existence of such current-carrying pure NESS with area-law entanglement in an interacting system, as discussed in ref.Gullans and Huse, 2019.

The phase boundaries in the $W-\epsilon$ plane for all the transitions, ergodic to non-ergodic, entanglement, time-reversal breaking and current transition, and their evolutions with $\Psi$ are summarized in Fig.2(b),(c). For comparison, we also show the ergodic to non-ergodic and entanglement transition for $\Psi=0$ in Fig.2(a). The three phase boundaries $W_{c1},~{}W_{c2}$ and $W_{c3}$ naturally shift to higher disorder with increasing $\Psi$. As already mentioned, we obtain the phase boundary, $W_{c1}(\epsilon)$, from the peak of the standard deviation of entanglement entropy $\sigma_{E}$ [Fig.5. Similar, albeit slightly higher, values for $W_{c1}$ are obtained from the crossing of $S_{EE}/L$ vs. $W$ curves [Fig.4]. However, we do not use the crossing point to plot phase boundary as a function of $\epsilon$, since clear crossings can only be detected over middle one-third and one-half of the spectra for $\Psi=0.3$ and $\Psi=0.6$, respectively. We note that the clear crossing point for $S_{EE}/L$ vs. $W$ curves could be obtained for the Hermitian case ($\Psi=0$) almost over the entire spectrum, except at the very edges.

IV.4 Growth and decay of entanglement entropy and the approach to NESS:

In this section we study the time evolution under the non-Hermitian evolution [Eq.(4)] starting with a simple un-entangled initial state. We choose the Néel state $|\psi_{0}\rangle=|\uparrow\downarrow\uparrow\downarrow\dots\rangle$ for this purpose. We compute the entanglement entropy, $S_{EE}(t)$, as function of time from the time-evolved state $|\psi(t)\rangle$ [Eq.(4)]. The results for $S_{EE}(t)$ are shown in Figs.8(a),(b) for two values of disorder strength, (a) $W=2.15$, in the volume-law phase, and (b) $W=4.3$, the $\mathcal{T}$-reversal broken area-law phase. The result for higher disorder, in the $\mathcal{T}$-reversal unbroken area-law phase, is similar to that for (b). In each of the cases, to bring out the crucial effect of interaction, we compare $S_{EE}(t)$ with that obtained for the non-interacting Hatano-Nelson model ($J_{z}=0$). The prominent features of $S_{EE}(t)$ are, (i) an initial growth, (ii) a broad peak followed by a decay or relaxation, and (iii) an eventual approach to a steady state value corresponding to the NESS. In the volume-law phase [Fig.8(a)(inset)], $S_{EE}(t)$ grows linearly with time, whereas, $S_{EE}(t)\propto\ln{t}$ initially in the area-law phase [Fig.8(b)], as in the Hermitian MBL case Bardarson et al. (2012).

Using Eq.(4), the growth and subsequent decay of entanglement entropy can be understood from the density matrix, $\rho(t)=|\psi(t)\rangle\langle\psi(t)|=e^{-\mathcal{L}t}\rho_{0}=\sum_{n,m}e^{-(\beta_{nm}+\mathfrak{i}\delta_{nm})t}C^{\psi}_{nm}(t)|n_{R}\rangle\langle m_{R}|$, where the coefficient $C^{\psi}_{nm}(t)$ is obtained from Eq.(4). Here $-(\beta_{nm}+\mathfrak{i}\delta_{nm})$ could be thought of as the eigenvalues of the Liouvillian operator $\mathcal{L}$; $\beta_{nm}=2\Lambda_{s}-\Lambda_{n}-\Lambda_{m}\geq 0$ and $\delta_{nm}=\mathcal{E}_{n}-\mathcal{E}_{m}$. The real part, $-\beta_{nm}$, of eigenvalue of $\mathcal{L}$ leads to relaxation and $\beta_{nm}=0$ corresponds to the long-time steady state. For a weak strength ($\Psi\ll 1$) of the non-Hermitian term, and for $t\ll e^{L}$ 111Since the many-body spectra for the real part of the eigenvalues are exponentially dense, the actual typical level spacing $\tilde{\delta}_{typ}\sim e^{-L}$. However this does not contribute to the time evolution in Figs.8(a), (b) for $t\ll e^{L}$, $\delta_{typ}\gg\beta_{typ}$, and the initial growth of the entanglement entropy appears over a time window, \delta_{typ}^{-1}\ {\raise-2.15277pt\hbox{\buildrel<\over{\sim}}}\ t\ {\raise-2.15277pt\hbox{\buildrel<\over{\sim}}}\ \beta_{typ}^{-1}, due to the dephasing from the exponential factor $e^{\mathfrak{i}\delta_{nm}t}$ in $\rho(t)$. In this time window, the factor $e^{-\beta_{nm}t}\approx 1$. Here $\delta_{typ}$ and $\beta_{typ}$ are the typical values of $\delta_{nm}$ and $\beta_{nm}$, respectively, that contributes to $\rho(t)$ for $t\ll e^{L}$.

This initial dephasing mechanism is similar to the one that leads to $\ln{t}$-growth of entanglement entropy in the MBL phase Serbyn et al. (2013). For the interacting system, left to itself, the dephasing would typically lead to a diagonal ensemble and long-time state with volume-law entanglement, as it does for the MBL phase Bardarson et al. (2012). Nevertheless, for the non-Hermitian case, the relaxation due to $\beta_{nm}$ kicks in for t\ {\raise-2.15277pt\hbox{\buildrel>\over{\sim}}}\ \beta_{typ}^{-1} and gives rise to the decay of $S_{EE}(t)$, as in Figs.8(a),(b), before the system could reach the diagonal ensemble. However, $\beta_{nm}$ has a gap Can et al. (2019), $\beta_{g}\sim\mathcal{O}(1)$, above its minimum value, $\beta_{nm}=0$ (Appendix A). On the contrary, the spectrum of $\delta_{nm}$ is gapless, $\mathrm{min}(\delta_{nm})\sim e^{-L}$. Hence, there could be interesting dephasing dynamics that goes on during the decay over \beta_{typ}^{-1}\ {\raise-2.15277pt\hbox{\buildrel<\over{\sim}}}\ t\ {\raise-2.15277pt\hbox{\buildrel<\over{\sim}}}\ \beta_{g}^{-1}, even though $\delta_{typ}\gg\beta_{typ}$. For $t>\beta_{g}^{-1}$, the entanglement entropy rapidly approaches the value for that of the NESS, dictated by the eigenstate with the maximum imaginary part for the eigenvalue, as discussed earlier.

IV.5 Time evolution of Néel order parameter

We also show the time-evolution of Néel order parameter $N(t)=(2/L)\sum_{i}(-1)^{i}\langle\psi(t)|S_{i}^{z}|\psi(t)\rangle$ in Fig.9. This acts like a MBL order parameter Schreiber et al. (2015) and characterizes the memory of the initial state. It approaches a finite value for $t\to\infty$ for the infinite system size in the MBL phase of the Hermitian model; $N(t)$ decays to zero with time in the ergodic phase. Here, for the non-Hermitian case, we find that $N_{\infty}=N(t\to\infty)$, the Néel order of the NESS, decreases with $L$, even for strong disorder, deep in the area-law phase. This is expected in a driven system, which loses its initially memory due to the drive. We also notice an interesting dynamical regime, presumably over the time window $\beta_{typ}^{-1}<t<\beta_{g}^{-1}$, in the decay of $N(t)$ [Fig.9(b), (c)]. The order parameter plateaus over a region as a function of $t$, as if trying to retain the initial memory.

In contrast, such regime is absent in the Hermitian model with $\Psi=0$. As shown in Fig.9(a) for $W=0.43<W_{c1}$, the long-time value $N_{\infty}$ goes to zero in the ergodic phase. This can be seen from the fact that $N_{\infty}$ decreases with $L$ tending to zero for $L\to\infty$. On the contrary, in the MBL phase, for $W=4.3,~{}7$, $N_{\infty}$ approaches a constant value with increasing $L$. Hence, the Néel order parameter in this case serves as the MBL order parameter, i.e. it diagnoses the persistence of the initial memory at arbitrary long times. In the non-Hermitian model [Eq.(1)], the NESS is approached at long times. This is true even deep in the MBL phase, $W\gg W_{c2}$, since there is always a finite number of eigenstates with complex eigenvalues, albeit with a vanishing fraction, for any finite-size system. Of course, for $W\gg W_{c2}$, it is expected that the imaginary part of the eigenvalues $\sim e^{-L/\zeta}$, and, hence, the NESS will ensue only after a very long time $\sim e^{L/\zeta}$ for large systems. As shown in Fig.9(b) for $\Psi=0.3$, and in Fig.9(c) for $\Psi=0.6$, the Néel order parameter for the NESS decreases rapidly with $L$, presumably approaching zero for $L\to\infty$.

IV.6 Properties of the NESS

We have so far discussed the entanglement entropy and current carried by the eigenstates over the entire spectrum of eigenstates. In this section we discuss the entanglement entropy and system-size scaling of current for the long-time NESS for time evolution from Néel state. The NESS is governed by the eigenstate with maximum imaginary part of the eigenvalue or the current [Eq.(5)].

IV.6.1 Entanglement entropy of NESS

We denote the entanglement entropy of the NESS by $S_{EE}^{\infty}$, i.e. $t\to\infty$ limit of $S_{EE}(t)$ [Fig.8]. We also analyze the maximum value, $S_{EE}^{peak}$, reached by $S_{EE}(t)$ during its time evolution towards NESS. We show $S_{EE}^{\infty}$ and $S_{EE}^{peak}$ in Figs.10(a), (b), as a function of $W$ and $L$. We find no volume-law to area-law transition in either $S_{EE}^{peak}$ or $S_{EE}^{\infty}$ with the disorder strength. In fact, $S_{EE}^{peak}$ and $S_{EE}^{\infty}$ neither scale as the volume nor the area, even though both of these increase with $L$. Hence, the entanglement for the NESSs in the non-Hermitian model of Eq.(1) obeys a system-size scaling intermediate between the area- and the volume-law. Similarly, the results for $S_{EE}^{peak}$ indicates that $S_{EE}(t)$ never attains the diagonal ensemble value during its time evolution, unlike that in the Hermitian case Bardarson et al. (2012). For a given $L$, unlike the eigenstate averaged $S_{EE}$ [Fig.3(a)], the maxima in $S_{EE}^{peak}$ and $S_{EE}^{\infty}$ appear at a finite disorder. Also, we do not see any system-size crossing for $S_{EE}^{\infty}/L$ ($S_{EE}^{peak}/L$) vs. $W$ in Figs.10(c), (d). We were also not able to obtain any finite-size data collapse for $S_{EE}^{\infty}$ and $S_{EE}^{peak}$.

IV.6.2 NESS current

There is also a current transition for the NESS, as in Fig.7 for the eigenstates. We show the system size scaling of current ($j_{\infty}$) carried by the NESS in Fig.11. The transition is shown in terms of a finite-size scaling collapse for $Lj_{\infty}$ vs. $W$ curves.

V Conclusions and discussions

To summarize, we have studied a non-Hermitian disordered model, the interacting version of Hatano-Nelson model Hatano and Nelson (1996). We propose that the model can be used to generate a rich variety of current-carrying states and study their entanglement properties. For example, we have found both ergodic and non-ergodic eigenstates with volume and area-law scalings, respectively. Furthermore, we show the existence of a time-reversal symmetry breaking phase transition within the area-law phase. We have established a detailed phase diagram as a function of disorder and strength of the non-Hermitian term based on the properties of the eigenstates and the long-time NESS, as well as, from the time evolution of entanglement entropy.

The models of Eqs.(1),(2) combine several non-trivial aspects of open quantum systems, namely disorder, interaction and coupling to environment via a non-Hermitian term that also induces a current drive. The disorder tries to localize the spin-flip excitations in the spin model [Eq.(1)] or the fermions in the model of Eq.(2). The non-Hermitian term, on the other hand, tries to induce a current and adds to the delocalizing tendency of the hopping term in Eq.(2). This can be easily understood by considering the limit $\Psi\rightarrow\infty$, where the eigenstates of the Hamiltonian are also the eigenstates of the current operator and hence the eigenstates carry finite current. Note that the hopping and the non-Hermitian term in Eq.(2) commute with each other. The competition of disorder and current drive and dissipation, already present in the non-interacting Hatano-Nelson model Hatano and Nelson (1996, 1997, 1998), leads to a localization-delocalization and $\mathcal{T}$-reversal breaking transition transition even in 1D. This is in contrast to the Hermitian 1D Anderson model which, as well known, has all the single-particle states localized and there is no transition.

A similar competition between disorder and non-Hermitian term modeling cooperative decay channel has been studied in other non-Hermitian Anderson models, in the context of collections of overlapping two level systems coupled to light, see, e.g., Refs.Celardo et al., 2014; Biella et al., 2013, which discuss a superradiance (delocalization) to localization transition as a function of disorder. However, the Hatano-Nelson model has the additional aspect of a possible $\mathcal{T}$-reversal breaking in the eigenstates.

As in the Hermitian system, the interaction term in Eq.(2), of course, adds to the delocalization tendency by generating matrix elements between exponentially large number of localized many-body product eigenstates of the non-interacting model Basko et al. (2006). In addition, interaction in the non-Hermitian model separates the localization-delocalization and the $\mathcal{T}$-reversal breaking transitions, thus leading to the existence of a novel $\mathcal{T}$-reversal broken area-law phase, unlike that in the non-interacting Hatano-Nelson model. The interaction is also crucial to give rise to intermediate time entanglement growth, absent in the non-interacting model, as shown in Fig.8. Overall, the non-Hermitian and the interaction terms makes the delocalized phase more robust to disorder.

During the completion of our work we became aware of an independent recent work Hamazaki et al. (2019) which has studied non-Hermitian MBL transition in the same model. Our focus, i.e. to model current driven systems and NESS, is entirely different from that of Ref.Hamazaki et al., 2019. Our results and conclusions are also substantially different from those in ref.Hamazaki et al., 2019. In particular, unlike Ref.Hamazaki et al., 2019, we find that the phase boundaries in $W-\epsilon$ plane for $\mathcal{T}$ reversal-symmetry breaking and entanglement transitions are distinct. The main reasons behind the difference can be attributed to the choice of parameters, the methodology of detecting phase transition, and, to some extent, interpretation of small-system ED results. The latter is a general problem even for well-studied Hermitian MBL case Luitz et al. (2015); Macé et al. (2019). Ref.Hamazaki et al., 2019 typically focuses on a smaller value of $\Psi=0.1$ and obtain the phase diagram based on level statistics, fraction of imaginary eigenvalues and entanglement entropy averaged over middle one-third of the spectrum, unlike the energy-resolved phase boundary in the $W-\epsilon$ plane in our case. Due to the smaller value of $\Psi$ the $\mathcal{T}$ reversal-symmetry breaking and entanglement transitions are much closer, even though different, in Ref.Hamazaki et al., 2019. Based on some heuristic arguments, the authors of Ref.Hamazaki et al., 2019 conjecture that the small difference between the two transitions will go to zero in the thermodynamic limit. On the contrary, our results with detailed finite-size scaling analysis in the $W-\epsilon$ plane shows that the two phase boundaries are clearly distinct for the particular choice of parameters ($\Psi=0.3,0.6$) in our work. We also note that our results for the ergodic to non-ergodic transition is consistent with a more recent work White et al. (2020) on the same model, studied in a different context. Moreover, we give detailed insight into the approach to NESS and identify new dynamical regime where memory of initial state can be harnessed over a relatively long intermediate time window during the relaxation to the NESS (Sections IV.5).

In future, it would be also interesting to understand whether the non-Hermitian model can indeed describe some features of a many-body localizable system under an actual current drive or voltage bias applied through leads.

Acknowledgement: We thank Subroto Mukerjee, Abhishek Dhar and Adhip Agarwala for useful discussions. We also thank Marko Znidaric and Nicolas Laflorencie for their valuable comments on our manuscript. SB acknowledges support from The Infosys Foundation (India), and SERB (DST, India) ECR award .

Appendix A Spectrum of the Liouvillian operator

As discussed in the main text, the time evolution of the density matrix is controlled by the eigenspectrum of the Liouvillian operator $\mathcal{L}$ in Eq.(3). The relaxation of the system to the NESS is controlled by the real part of the eigenvalues of $\mathcal{L}$, i.e. $-\beta_{nm}>0$, where $\beta_{nm}$ $=2\Lambda_{s}-\Lambda_{n}-\Lambda_{m}$. In Figs.12, we plot the disorder averaged distribution, $P(\bar{\beta})$, of the normalized quantity, $\bar{\beta}=(\beta-\beta_{0})/(\beta_{M}-\beta_{0})$ for $\Psi=0.3$. Here $\beta_{0}$ and $\beta_{M}$ are the minimum and the maximum values of $\beta_{nm}$, respectively. Since the imaginary eigenvalues of the non-Hermitian Hamiltonian in Eq.(1) appears in complex conjugate pairs, the probability distribution is symmetric about $\bar{\beta}=1/2$. Also, there is a peak at $\bar{\beta}=0$, which corresponds to the NESS, as mentioned in the main text. As shown in Figs.12, the peak is separated from most of the spectrum by a long tail. In the limit $L\to\infty$, the tail is expected to tend to a gap Can et al. (2019), $\tilde{\beta}_{g}$, that separates the peak from the relaxation modes with $\bar{\beta}>0$.