Statics and Dynamics of non-Hermitian Many-Body Localization

,


Introduction
The investigation of isolated quantum systems has led to a remarkable understanding of novel states of matter [1,2].However, naturally occurring and engineered quantum systems are typically coupled to an environment, even if the coupling is weak.The dynamics of such open systems, over a reduced set of trajectories, can often be described by an effective non-Hermitian Hamiltonian which breaks unitarity [3].These have been realized in pioneering experiments on photonic [4][5][6][7] and matter-light [8][9][10][11][12] systems.The study of non-Hermitian Hamiltonians allows an enriched classification scheme for describing quantum matter [13,14].For example, non-Hermitian systems with parity and time-reversal symmetry can have real eigenvalues, much like their Hermitian counterparts [15][16][17].Non-Hermitian perturbations can also lead to novel phases and phase transitions beyond the equilibrium paradigm [12,[18][19][20].For a review of the applications of non-Hermitian systems see [21].
An important class of isolated quantum systems are those that fail to equilibrate on long timescales.This includes many-body localized (MBL) phases which retain memory of their initial conditions, as observed in analog and digital quantum simulators [22][23][24][25][26][27][28].The fate of MBL in open systems has also been investigated in cold atom settings via a controlled coupling to the environment [29][30][31][32][33][34][35].This is of considerable interest in solid state devices due to their intrinsic coupling to other degrees of freedom [36,37].The stability of MBL has also been studied [38,39] including the effects of local dissipation in the Lindblad formalism [29][30][31][32][40][41][42][43], and by coupling to delocalized environments; see for example [35,44].Effective non-Hermitian models can also provide insight into the stability of quantum matter in the presence of coupling to an environment.In the context of non-Hermitian MBL, pioneering studies have focused on the link between the phase diagram and spectral properties [45], together with mobility edges [46,47].Extensions of these investigations have also considered wavepacket dynamics [48] and the effects of quasiperiodic potentials [49,50].
In this work we explore the phase diagram of the interacting Hatano-Nelson model as a function of non-Hermiticity and the interaction strength; see Fig. 1.Using exact diagonalization, we provide evidence for a two-step approach towards the localized regime, with a significant separation between eigenstate and spectral transitions.The latter yields a dynamical instability within the localized regime.We provide predictions for the relaxation time scales of the particle imbalance with a view towards cold atom experiments.

Model
We consider an interacting version of the one-dimensional Hatano-Nelson model [51][52][53][54][55] with L sites and periodic boundary conditions where J is the hopping strength, g parametrizes the hopping asymmetry, h i represents on-site disorder and U is the nearest-neighbor interaction strength.For simplicity we consider hard-core bosons bi and b † i where ni = b † i bi .The disorder is drawn from a uniform distribution h i ∈ [−h, h].The asymmetry in the hopping terms renders the model non-Hermitian when g ̸ = 0. Throughout this work we consider half-filling and set J = 1.
As discussed in Supplementary Note 1, the model ( 1) can be derived within the Lindblad formalism by assuming a specific form of the jump operators [56] and post-selecting on jump-free trajectories.It therefore provides a useful framework for investigating the stability of MBL when coupled to an environment.Although the experimental realization of such systems is hindered by the need for post-selection, this may be within reach for modest system sizes [57,58].In particular, non-reciprocal transport has been realized in cold atoms [59] and non-Hermitian effects have been simulated with trapped ions [60].Refs [56,61,62] discuss the experimental realizations of non-Hermitian models related to the model (1) using cold atoms.
Before embarking on a detailed study of the model (1) we first discuss some limiting cases.In the non-interacting limit with U = 0, all states are localized for g = 0 [63].However, for g ̸ = 0 it undergoes a delocalization transition [51][52][53][54][55].In the Hermitian case (g = 0) with U ̸ = 0, the model exhibits a many-body localization (MBL) transition [36,[64][65][66][67].In the clean limit (h = 0) the model exhibits PT -symmetry breaking transitions in both the ground state and excited state sectors with only the former surviving to the thermodynamic limit [68].In this paper we study the interplay between non-Hermiticity and MBL by considering both the mid-spectrum states of the model (1) and its dynamics.To obtain the mid-spectrum states we employ exact diagonalization (ED) and retain a total of N T = ⌈0.04N⌉ eigenstates (rounded up to the nearest integer) which are closest to the mid-point of the spectrum, Tr(H)/N ; here L is the dimension of the Hilbert space.

Symmetry
The non-Hermitian Hamiltonian (1) exhibits time-reversal symmetry which ensures that the eigenvalues are either real or occur in complex conjugate pairs [45,69].In the noninteracting limit with U = 0 it has been shown that the localized eigenstates correspond to purely real eigenvalues and the delocalized states correspond to complex conjugate pairs [51][52][53][54][55].At large disorder all states are localized for U = 0, corresponding to an entirely real spectrum.In contrast, in the interacting case it is only the fraction of complex eigenvalues that goes to zero.In fact, the number of complex eigenvalues remains non-zero and grows exponentially with increasing system size, as shown in Fig. 2. For weak disorder the number of complex eigenvalues N C approaches the total number of computed eigenvalues N T ∼ 2 L / √ L but the number of real eigenvalues N R remains non-zero, see Fig. 2(a).Conversely, in the strong disorder regime the number of real eigenvalues approaches the total number of eigenvalues but the number of complex eigenvalues remains non-zero, see Fig. 2(b).As we discuss more fully below, in the many-body problem localized eigenstates may correspond to complex eigenvalues.

Spectral Transition
Following Ref. [45] we first consider the behavior of the fraction of complex eigenvalues f C = N C /N T with increasing disorder strength, where the overbar denotes disorder averaging.As shown in Fig. 3(a) the spectrum of the Hamiltonian undergoes a transition at a critical disorder strength where f C goes from increasing to decreasing with increasing system size [45].However, the number of complex eigenvalues N C does not go to zero with increasing L in the strong disorder regime.This is in contrast to the non-interacting (U = 0) case [54].In Fig. 1 we plot the evolution of this boundary as a function of the non-Hermiticity g and the disorder strength h.It can be seen that at large non-Hermiticity the transition occurs for larger disorder strength.

Eigenstate Transition
Having established the locus of the spectral transition, we now turn our attention to the stability of the eigenstates.In the Hermitian case, one of the hallmarks of localized eigenstates is their stability to local perturbations [65].An extension of this to the non-Hermitian case was suggested in Ref. [45].First, we denote the left and right eigenstates of the non-Hermitian Hamiltonian (1) by Denoting the eigenstates which correspond to purely real eigenenergies by |E k ⟩ R and |E k ⟩ L we examine the quantity where VNH is a non-Hermitian perturbation and Here we take VNH = b † 1 b2 and the eigenstates are ordered with increasing E ′ k .In Fig. 3(b) we plot the evolution of the instability parameter G versus disorder strength for different system sizes.It can be seen that the eigenstates are stable (unstable) above (below) a critical disorder strength.However, the value of the critical disorder strength differs from that obtained from f C .In Fig. 1 we plot the locus of the stability line obtained from Eq. ( 2).It can be seen to be well-separated from the spectral transition in our finite-size simulations.That is to say, the instability of the real eigenstates occurs in advance of the spectral transition where f C → 0 with increasing L. This suggests the possibility of a two-step transition to the localized phase in open systems.This is consistent with results obtained from the study of mobility edges [47].

Role of Interactions
Having established the presence of two boundaries in the non-Hermitian problem we now turn our attention to the role of the interaction strength.In Fig. 4 we show the evolution of the boundaries with increasing U for a fixed value of the non-Hermiticity g.It can be seen that in the non-interacting limit with U = 0 the transitions coincide [52,54,70]; that is to say the localization transition coincides with the spectral transition from complex to fully real.However, in the presence of interactions the phase structure changes.The two transitions separate with increasing U .This suggests the possibility of an intermediate regime as one goes from weak disorder to strong, as found in Fig. 1.

Dynamics
To further characterize the localized and delocalized phases, we turn our attention to non-equilibrium dynamics.We focus on the particle imbalance as measured in experiments on isolated systems [22].Here n e (t) and n o (t) are the occupations of the even and odd lattice sites, respectively.In the case of an initial density wave |Ψ(0)⟩ = |0101 . ..⟩ the imbalance can be rewritten as where N = L/2 is the total number of particles.The state of the system evolves according to where the normalization is explicitly enforced; the operator exp(−i Ĥt) does not preserve the norm of |Ψ(0)⟩ when Ĥ is non-Hermitian.In the context of Lindbladian dynamics this corresponds to trajectories which are post-selected on the absence of quantum jumps [3,21].
To explore the impact of complex eigenvalues on the relaxation dynamics, we consider the formation of a steady state in I(t).In Fig. 5(a) we plot the time-evolution of I(t) in both the delocalized and localized regimes corresponding to h = 3 and h = 18, respectively.In each case we take a single realization of the disorder with at least one complex conjugate eigenvalue pair in the spectrum.It is readily seen that I(t) approaches a steady state after a time τ , as indicated by the red square markers.Heuristically, one expects that τ is governed by the largest imaginary eigenvalue Λ max , corresponding to the largest rate of amplification; the positive and negative imaginary parts of the complex eigenvalue pairs correspond to amplification and decay, respectively.This is borne out in Fig. 5(b) which shows the distribution of τ for different disorder realizations and two values of the disorder strength h; only realizations with complex eigenvalues in their spectra are considered for this analysis.The dashed line is a guide to the eye showing τ = Λ −1 max .A notable feature of Fig. 5(b) is that the distribution of timescales changes markedly in going from the weak to the strong disorder regime.In the localized phase the distribution of the scaled quantity τ Λ max shows a concentration of timescales in the vicinity of τ = αΛ −1 max , where log 10 α ≈ 1.3.However, in the weak disorder regime we see a broader distribution of timescales which are only bounded from below by Λ −1 max .Moreover, the distribution of timescales becomes elongated along the axis τ = αΛ −1 max in the vicinity of the transition region (region II), as depicted in Figs 1 and 4.
To explore this further we plot the distribution of τ Λ max for different values of the disorder strength.The distribution develops a sharp peak in the vicinity of the eigenstate transition at h ≈ 10, for the chosen parameters.The peak location at α = τ Λ max with log 10 α ≈ 1.3 coincides with the relationship observed in Fig. 5, and does not change with increasing disorder strength.To show this more clearly, in the inset of Fig. 6 we plot the evolution of τ Λ max at the overall peak of the distribution shown in the main panel.The peak location shifts to lower values with increasing h, becoming independent of h in the localized phase.This is consistent with a direct relationship between τ and Λ −1 max in the localized phase.This is in marked contrast with the thermal phase which shows a broader distribution without a sharp peak.

Discussion
The preceding analysis shows that the phase diagram in Figs 1 and 4 can be interpreted in terms of the dynamical behavior of physical observables.In particular, the presence of complex eigenvalues directly impacts on the memory lifetime of the initial state.For weak disorder (region I), this memory is lost due to the presence of complex eigenvalues, reflecting delocalized eigenstates and non-Hermiticity.For intermediate disorder (region II), the memory of generic initial states is still eroded due to complex eigenvalues, despite the onset of localization.Nonetheless, the presence of localized real eigenstates can lead to long-lived memory of the initial conditions.Only at strong disorder (region III), where the vast majority of the eigenvalues are localized, and the eigenvalues are real, do generic initial states have long lifetimes.However, even here, this lifetime is not infinite -there is always some residual decay due to the presence of complex eigenvalues.
Although we cannot exclude the possibility that the eigenstate and spectral transitions merge in the thermodynamic limit, the crossing points are well separated for a broad range of parameters and accessible system sizes.This suggests the possibility of an intermediate region for at least part of the parameter space.In contrast, in the single-particle limit of the model ( 1), delocalization is always accompanied by complex eigenvalue formation and vice versa.As such, an intermediate region is absent in the single-particle case.

Conclusion
In this work we have investigated the phase diagram of the interacting Hatano-Nelson model as a function of the interaction strength and the non-Hermiticity.We have mapped out two regimes for MBL via the mid-spectrum eigenstates.In particular, we have shown that the delocalization instability of the real eigenstates occurs at a weaker disorder strength than the transition to a predominantly real spectrum.We have also explored the non-equilibrium dynamics of this model and shown the appearance of a dynamical signature in the vicinity of the eigenstate transition.It would be interesting to explore the phase diagram of this problem in the framework of Lindbladian dynamics.Finally, we note that the apparent extrapolation of region II to infinitesimal non-Hermiticity in Fig. 1 is suggestive of it probing the separation between landmarks in the Hermitian MBL transition [42].To rigorously establish or exclude such a connection would be an interesting subject of future research.

Methods
The eigenvalues and eigenvectors of the model (1) are calculated via an equivalent spin model using ED.The hard-core bosons are mapped to spins using the transformation bi , where Ŝ± j = Ŝx j ± Ŝy j are raising and lowering operators and Ŝα j is the α-projection of the spin.The spin Hamiltonian is The matrix elements of Ĥ are computed in a basis of product states using QuSpin [71].Figs 1-4 are obtained from the midspectrum eigenstates and eigenvectors using the shift-invert ED routine in SciPy [72]; related algorithms have been used in the Hermitian case [73,74].The right eigenvectors are obtained directly, and the left eigenvectors are obtained from the right eigenvectors of the transposed matrix.The dynamics in Figs 5 and 6 is obtained by iterating Eq. ( 4) using small time steps δt.This avoids the growth of the norm due to complex eigenvalues.Explicitly, we decompose exp where N is the dimension of the Hilbert space.The largest time step δt is chosen so that the norm of a vector with components exp(−iE k δt) is kept below 10 −10 .The data in Fig. 5(a) is obtained by sampling on a linear grid until t = 10 −1 and logarithmically thereafter.

Declarations
Data availability.The simulation data that have been generated and analyzed during this study are deposited in the King's Open Research Data System (DOI:10.18742/25130666) [75].

Fig. 1 Fig. 2
Fig. 1 Phase diagram of the interacting Hatano-Nelson model as a function of the non-Hermiticity parameter g and the disorder strength h.We set the hopping strength to J = 1 and the interaction strength to U = 2.The results are obtained using exact diagonalization (ED) with L = 8, 10, 12, 14, 16 sites.The left boundary (blue line) corresponds to the instability G of real eigenstates to a non-Hermitian perturbation.The right boundary (black line) indicates the spectral feature where the fraction of complex energy eigenvalues f C goes from increasing to decreasing with increasing L. Region I corresponds to unstable real eigenstates and increasing f C .Region II corresponds to stable real eigenstates and increasing f C .Region III corresponds to stable real eigenstates and decreasing f C .The error bars indicate the variation of the finite-size crossing points for each data set.

Fig. 3
Fig.3Determination of the locations of the spectral and eigenstate transitions.a Variation of the fraction of complex eigenvalues f C as a function of the disorder strength h for different system sizes L, with the non-Hermiticity parameter g = 0.5 and the interaction strength U = 2.The crossing point at h ≈ 10.8 shows the separatrix between regions II and III in Fig.1.An eigenvalue is considered real if its absolute imaginary part is below 10 −13 and complex otherwise.b Variation of the eigenstate instability measure G as a function of the disorder strength h for different system sizes L, with g = 0.5 and U = 2.The crossing point at h ≈ 9.1 shows the separatrix between regions I and II in Fig.1.In our finite-size simulations the transitions in panels a and b are well-separated.The vertical shaded regions extend between the lowest and highest values of h where the data for different values of L cross.The vertical line (red) indicates the mid-point of the region.The mid-point and region width are plotted as the markers and the error bars, respectively, in Fig.1.In both panels the data points are computed as the mean over 10 4 disorder realizations.The error bars corresponding to the standard error of the mean are smaller than the markers.

Fig. 4 18 Fig. 5
Fig.4  Evolution of the phase diagram of the interacting Hatano-Nelson model as a function of interaction strength U and the disorder strength h.We set the non-Hermiticity parameter to g = 0.5 and the hopping strength to J = 1.The transition points and error bars are extracted using the method described in Fig.3using the same system sizes.Regions I-III are defined in Fig.1.The single-particle localization transition (red square) is computed from the crossing points of the fraction of complex eigenvalues f C for L = 100, 200, 300.Inset: evolution of the width ∆h of region II, corresponding to the horizontal separation between the data points in the main panel, as a function of U .

Fig. 6
Fig.6The distribution τ Λmax develops a sharp peak near region II of the phase diagram.Distribution of τ Λmax as a function of the disorder strength h.The non-Hermiticity parameter is set to g = 0.5 and the system size is L = 14, as used in Fig.5.The distribution exhibits a broad maximum which shifts towards lower values with increasing disorder strength h, as indicated by the black arrow.In the vicinity of the transition region (region II) shown in Figs1 and 4the distribution develops a sharp peak, as indicated by the red arrow.The location of this peak occurs at τ = αΛ −1 max , with log 10 α ≈ 1.3 in agreement with Fig.5.Inset: The location of the overall broad peak for L = 14 moves to lower values of α = τ Λmax with increasing disorder strength.The value of α becomes independent of h in the localized phase with log 10 α ≈ 1.3.