Fate of the non-Hermitian skin effect in many-body fermionic systems

We revisit the fate of the skin modes in many-body non-Hermitian fermionic systems. Contrary to the single-particle case, the many-body ground state cannot exhibit an exponential localization of all eigenstates due to the Pauli exclusion principle. However, asymmetry can still exist in the density profile, which can be quantified using the imbalance between the two halves of the system. Using the non-Hermitian Su-Schrieffer-Heeger (SSH) chain as an illustration, we show the existence of two distinct scaling regimes for the imbalance. In the first one, the imbalance grows linearly with the system size, as generically expected. In the second one, the imbalance saturates to a finite value. By combining high-precision exact diagonalization calculations and analytical arguments, we observe that the imbalance does not scale when the occupied bands can be deformed to their Hermitian limit. This suggests a direct connection between the corresponding bulk topological invariants and the skin effect in many-body systems. Importantly, this relation also holds for interacting systems.

Another unique feature of nH systems is the anomalous localization of all eigenstates dubbed the non-Hermitian skin effect 28,34-37 . All single-particle eigenstates become exponentially localized at the boundaries of the system. Generally, skin modes arise due to the breakdown of reciprocity, i.e., introducing asymmetric hoppings. The conventional bulk-boundary correspondence is then broken, obscuring the prediction of universal boundary phenomena from the periodic bulk properties [38][39][40][41][42][43][44][45][46][47][48][49] . Indeed, the phase diagram of the single particle Hamiltonian with periodic boundary conditions (PBC) can significantly differ from the one with open boundary conditions (OBC). The skin effect is directly linked to the structure of the singleparticle spectrum with PBC 33,50 . A non-zero winding of the PBC energy bands around any point in the complex plane guarantees that the related OBC modes are exponentially localized. The corresponding OBC eigenenergies are (nearly all) located within these non-zero winding regions. The existence of the skin effect was confirmed experimentally 2,4, 22,[51][52][53][54] . Recently, the notion of skin effect has been expanded by introducing different types of skin states 55 such as symmetry-enforced variants 50,56,57 , reciprocal skin effect 21 in two or higher dimensions, higher-order skin effect [58][59][60][61][62][63] , or the critical skin effect 64,65 .
So far, studies of the skin effect have largely focused on the single-particle Hamiltonian. Most of the experi-ments described by these effective nH Hamiltonians are either weakly interacting or non-interacting bosonic systems (e.g. quantum optics) or classical systems (e.g. the various types of circuits and active media). Moreover, in Hermitian non-interacting systems, the properties of the many-body states can be straightforwardly deduced from those of the single-particle Hamiltonian; this is not the case for nH models. Due to the non-orthogonality of the single-particle eigenmodes, the Pauli principle no longer acts trivially. The fermionic repulsion reshapes the occupied orbitals so that no more than one fermion (equivalently one hard-core boson) populates each physical site. Hence, the exponential localization of all fermions at a boundary becomes impossible 66,67 . Recently, several authors have investigated the fate of the skin modes in the many-body context 66,[68][69][70][71][72] . For instance, for interacting nH systems such as topological Mott insulators, the skin effect was observed in fermionic systems 67,73 , but not in the bosonic case 74,75 . Nonetheless, no general criterion for the survival of the skin effect has been derived so far.
In this Letter, we investigate the connection between topological properties of non-interacting nH systems and the skin effect in a many-body wave-function. The exponential localization of the skin modes naturally translates to an asymmetry of the many-body density profile that grows linearly with the system size. Through numerical simulations, using high-precision exact diagonalization (ED) and density matrix renormalization group (DMRG), we show that the many-body eigenstates of OBC systems exhibit a transition from a regime with an infinite asymmetry in the thermodynamic limit to one where it saturates. The critical anisotropy corresponds to a change in the topology of the periodic single-particle spectrum, even though the open system remains gapped at all times. In addition, we show that both scaling regimes survive disorder and interactions.

Conventions
where c ( †) are fermionic annihilation (creation) operators, j denotes the lattice position, A/B corresponds to the sublattice, and t 1/2 is the hopping amplitude. g ≥ 0 is a hopping anisotropy that breaks Hermiticity, see Fig. 1(a). We consider three boundary conditions: PBC, OBC (with 2L sites) or ABA-BC (OBC with 2L+1 sites). ABA-BC has an explicit analytical solution (see Appendix A and Ref. 77). For PBC, the single-particle gap closes for |t 1 e ±2g | = |t 2 |, with two phases adiabatically connected to the trivial and topological Hermitian phases, and an intermediate phase where the system has no line gap. In all cases, the single-particle spectrum has a non-trivial winding, and therefore the nH skin effect is always present. Indeed, for OBC and ABA-BC, all single-particle eigenstates are exponentially localized to the right side of the system. The model defined in Eq. (1) is related to the Hermitian SSH chain via the similarity transformation such that Note that the OBC and ABA-BC phase diagrams are therefore independent of g with a gap closing at t 1 = t 2 , and their spectrum is always real. In the following, we focus on the ground state properties of the OBC and ABA-BC systems. We define the ground state as the right eigenstate that minimizes the real part of the energy: Our results generalize straightforwardly to other right eigenstates and to the left eigenvectors. The observables are taken to be: This expectation value is relevant when considering nH Hamiltonians obtained from post-selection 78,79 . The quantum state described by such an approach remains a proper Hermitian density matrix. As H SSH {g} is non-interacting, observables can be obtained efficiently from the diagonalization of the single-particle Hamiltonian 80 (see Appendix B). We numerically diagonalize the Hermitian system with high precision before applying the similarity transformation to obtain the eigenstates Imbalance I as a function of the system size for t1/t2 = 2 and three representative values of g. At g < gc, I saturates to a finite value. At g > gc, the scaling becomes linear. It also exhibits a series of damped oscillations whose frequency diverges when g → g + c . The inset shows the corresponding PBC spectra in the complex plane.
of the nH models.
Imbalance as a marker of the skin effect: As can be seen from the transformation U g , all the single-particle orbitals are exponentially localized at the right boundary of the chain. We expect the asymmetry of the orbitals to transfer to the many-body density distribution n j g = n j,A + n j,B g . Due to the Pauli principle, only one fermion can be located at a given site and therefore the many-body version of the skin effect is not a sum of the exponential orbitals. To quantify the degree of asymmetry of the density distribution, we use the imbalance Working at fixed charge density N F = j n j g /2L ≤ 0.5 81 , the imbalance can at most scale linearly with system size. We consider this linear scaling to be the manifestation of the skin effect in many-body systems. In Fig. 2, we show the scaling of the imbalance in the limit t 1 = t 2 where the SSH model maps to the Hatano-Nelson model 68,82 . We observe a linear scaling of the imbalance at all fillings. At small g, the slope increases with the anisotropy g, but saturates at a value equal to N F . The large g limit corresponds to a situation where all particles are located in the right half of the system. We note that the imbalance can be generalized to the higher moments of the density with similar results.

Skin effect and band topology:
We now turn to the general case where t 1 = t 2 and focus on half-filled systems (L/(2L + 1) filling for ABA-BC). Fig. 1 summarizes our results for ABA-BC. Below and at the critical value g c = 0.5 | log t 2 /t 1 |, the imbalance saturates with system size. For g < g c , the saturating value is reached exponentially fast. At the critical point, the convergence becomes algebraic, with an exponent which depends on the t 1 /t 2 ratio and appears to peak at t 1 = t 2 (see Appendix C for more details). Above g c , the imbalance scales linearly with L, marking a fundamental anisotropy in the thermodynamic limit. We also observe slowly decaying oscillations around the linear regime. The relative amplitude of the oscillations, as well as their period, increases with decreasing g. The limit g → g + c appears to correspond to a divergence of the period of the oscillation, as well as a vanishing of the linear term. Due to the instability of the non-normal Hamiltonians, or alternatively the high amplitude of the coefficients in the similarity transformation U g , we stress the importance of using high-precision libraries to distinguish between finite-size 83 and finite-precision effects. The system sizes shown here require a precision of up to several hundred of digits, hence the use of the ABA-BC where we can bypass diagonalization of the single-particle Hamiltonian. The same results were obtained for OBC.
Remarkably, the critical value g c for the open system corresponds to the gap closing of the periodic system. The single-particle Hamiltonian with OBC or ABA-BC remains gapped when crossing g c 84 , and its eigenvectors vary smoothly. The abrupt change in the imbalance (discontinuous in the thermodynamic limit) directly arises from the Pauli principle.
To understand the source of the discontinuity, we summarize here how the many-body state is obtained from the single-particle orbitals. Let us denote {|R j } the right eigenvectors of the single-particle Hamiltonian, and c † R,j = c † |R j , the corresponding orbitals. The eigenvec-tors of the many-body Hamiltonian are given by Due to the non-orthonormality of the {|R j }'s, the orbitals do not respect the conventional fermionic anticommutation relation and the right side of Eq. (7) is not normalized. It is convenient to instead work with {|Q ij }an orthonormal family generating the occupied {|R ij } -and the corresponding orbitals {q † ij }. They can be obtained by Gram-Schmidt orthonormalization. We then have The q ij now obey standard anticommutation relations. Due to the orthonormalization, the smooth deformation of the original eigenvectors when g crosses g c is amplified as the norm of the intermediate states can be arbitrarily close to zero. We now focus on the saturating regime, g < g c . Note that the saturation of the imbalance depends on the fine structure of the orbitals. Taking the corresponding PBC Hermitian state instead leads to a linear scaling of the imbalance after applying U g . Generalized Brillouin zone (GBZ) approaches 34,46,85,86 give the same results. It is informative to briefly study a simple example of saturation. Consider L independent orbitals such that a j,n = αb j,n for 1 ≤ n ≤ L. By dimensional arguments, we obtain The many-body state with L occupied orbitals can therefore be rewritten as independently of the exact form of a j,n 87 . This also holds when the orbitals are all exponentially localized at a boundary; the imbalance is constant and equal to zero. As all orbitals have the same local structures, the Pauli principle forces electrons to spread out through the system.
The suppression of the scaling is observed only at halffilling in two-band models. If the low-energy band is only partially occupied (resp. the high-energy band is partially occupied), the imbalance grows linearly with system size. This can be immediately deduced from the toy model in Eq. (9): if we do not occupy the full translationinvariant subspace in Eq. (10), the imbalance survives. In the nH-SSH model, for g < g c , the periodic single-particle spectrum has a line gap, while for g > g c , it forms a single band with no line gap and therefore its many-body spectrum is gapless. Let {E PBC/OBC j } be the set of energy bands with PBC or OBC, respectively, and the bands are separated by line gaps. Following Ref. 50, the OBC and PBC bands can be deformed into each other (merging several OBC bands to form a single one if required), up to some small set of states. We propose the following conjecture: Conjecture 1 For non-interacting systems exhibiting single-particle skin effect, if the sets of orbitals occupied in a many-body state |Ψ can be mapped to fully occupied PBC bands (up to a number of orbitals of measure 0 in the thermodynamic limit), the many-body skin effect is suppressed. Otherwise, the imbalance in |Ψ will scale linearly with system size.
Beyond the SSH model: disorder, interactions and symmetries: Our conjecture implies that the many-body skin effect has a topological origin that can be predicted from considering the periodic system. It is generally not the case for topological properties of nH systems which can be obtained from a GBZ approach or from the bulk of the OBC system.
To further show the resilience of the skin effect, we investigate the effect of disorder 82,88-92 . For single-particle states, skin effect and disorder compete to localize states either at the boundaries or within the bulk of the system. When the disorder-induced localization length is larger than the one induced by the breaking of reciprocity (g −1 in our example), the single-particle states are Andersonlocalized. They then occupy arbitrary positions in the lattice which therefore limits the scaling of the imbalance. On the other hand, the saturation is highly sensitive to the structure of the orbitals and therefore its survival is unclear. In Fig. 3, we show the imbalance in a disordered nH-SSH model with a random chemical potential for several values of the disorder strength W . We consider two different disorder configurations: (a) µ B j = −µ A j , acting as a local chemical potential, or (b) µ B j and µ A j independent. Saturated and linear regimes are resilient to the presence of small disorder in both setups. For configuration (a), the disorder increases the critical value of g by opening a line gap between the bands of the corresponding PBC models (see Fig. 3(c, d)). At small W , similar results are obtained for configuration (b). However, at large W , the disorder closes the line gap at (and below) g c and the imbalance scales linearly accordingly.
In addition, we investigate the effect of interactions in the two scaling regimes. We use DMRG 93 to access larger system sizes with OBC before applying the similarity transformation which is a matrix product operator of bond dimension 1. We are mainly limited by the floating (1) the following interaction: For U real, the system remains gaugeable to a Hermitian model by the similarity transformation given in Eq.
(2). In the limit t 1 = t 2 (and g = 0), the model is equivalent to the XXZ chain without transverse field. In Fig. 4, we show the scaling of the imbalance for (a) attractive and (b) repulsive interactions. Both linear and saturated regimes survive the presence of small interactions. For attractive interactions (U < 0), fermions prefer to occupy neighbouring orbitals, and therefore the critical g c decreases. Conversely, for repulsive interactions (U > 0), the interactions favor the electrons spreading throughout the system, so g c increases. The notion of bands is no longer well-defined in the presence of interactions, but systems without a single-particle line gap have a gapless many-body spectrum. We verify that the variation of the critical g c matches a shift of the gap closing point. We use the response of the periodic system to flux insertion to detect the presence of a gap given the limited system sizes we have access to.
In the Appendices, we show that our results are also valid when the system is no longer gaugeable to a Hermitian limit or when we break particle-hole symmetry by introducing a third band.
Conclusions: In this Letter, we have investigated the skin effect in many-body wavefunctions. Generically, for fermionic systems, the single-particle skin effect translates into a linear scaling of the imbalance. However, we have numerically shown that if the orbitals occupied in an OBC many-body state correspond to fully occupied PBC bands, the skin effect is suppressed in the thermo-dynamic limit. Remarkably, this sharp transition occurs while the OBC spectrum remains completely gapped. We have demonstrated that the suppression of the skin effect, despite being sensitive to the details of the orbitals, survives both the presence of disorder and interactions. This strongly implies a topological origin of this phenomenon, similar to the topological nature of the nH skin effect in single-particle Hamiltonians. An analytical proof of this connection to topology is left for future work. Our results pave the way to a better understanding of the skin effect in higher dimensions, with the interplay between the dimensionality of the Fermi surface and other forms of skin effect.   The nH-SSH model with ABA-BC takes the form: Following Ref. 77, we can write explicitly its eigenvalues and eigenstates. Note that due to the choice of boundary conditions, there is always a zero-energy edge mode. The eigenvalues are given by: t 2 1 + t 2 2 + 2t 1 t 2 cos θ n , n = 1, . . . , L − t 2 1 + t 2 2 + 2t 1 t 2 cos θ n , n = L + 1, . . . , 2L 0, n = 2L + 1 (A2) and the corresponding eigenstates u and The zero energy eigenstate is then given by The many-body Hamiltonian H can be rewritten as With this expression, we see directly that its right eigenstates are proportional to Due to the non-orthogonality of the diagonalizing basis of nH operators, these states are not normalized. Moreover, contrary to the Hermitian case, observables cannot be directly computed from the expression above. Instead, one needs to orthonormalize the family of occupied orbitals. If {|φ j } is an orthonormal basis of Span({|ψ R ij }, and we define we see that the q ij 's are conventional fermionic operators. The eigenstate in the form is indeed normalized. |Ψ is a conventional fermionic Gaussian state, on which we can apply all the standard tricks.

Appendix C: Scaling analysis
As discussed in the main text, the imbalance as a function of the system size exhibits intriguing features depending on the value of g. To underline the universality of our results, we perform a scaling analysis for three g regimes. Below g c , the imbalance saturates to a finite and constant value. Away from the dimerized limits (t 1 = 0 or t 2 = 0) and from the gapless point t 1 = t 2 , the imbalance I can be well approximated by an exponential I(L) = I ∞ + λe −L/ξ , as shown in Fig. 5(a). At g = g c , the correlations become algebraic with an exponent that depends on t 1 /t 2 , as presented in Fig. 5(b).
Finally, in the g > g c regime, we observe oscillations of the imbalance around the linear scaling (see Fig. 6(a)). To analyze these oscillations, we study the derivative dI/dL, which clearly shows distinct peaks as demonstrated in Fig. 6(b). The amplitude of the oscillations slowly decays with L, while their period appears to diverge as (g − g c ) −1 : the saturation can be seen as an infinitely long plateau (see Fig. 6(c)).

Appendix D: Requirements for numerical precision
Apart from the errors originating from too small system sizes, i.e., finite size effects, numerical instabilities in   6. (a) Imbalance scaling for several g > gc. Each curve exhibits a sequence of plateaus, with a period depending on the exact value of g. To quantify this dependence, we compute the numerical derivative dI/dL and perform a spline interpolation to estimate the positions of the peaks, shown in (b). The distance between peaks T diverges as ∼ (g − gc) −1 .
nH matrices pose additional difficulties. 64-bits floating point precision becomes insufficient due to a combination of large system sizes and large values of the reciprocitybreaking terms. It is therefore necessary to systematically study how numerical accuracy affects the results, as it is known that nH matrices are exponentially unstable 43,98 . In Fig. 7, we compare the many-body ground state densities obtained from Gram-Schmidt orthonormalization at different numerical precisions. At low precision (up to 150 digits), the local density close to the edge exceeds 1, which is mathematically incorrect for normalized fermionic states. The minimal number of required decimal digits scales roughly as gL.
FIG. 7. Many-body densities |ΨRR| 2 obtained from the analytical solutions for ABA-BC with g = gc + 0.1 and L = 400 (801 sites) at different numerical accuracies in the Gram-Schmidt procedure. The numerical inaccuracies lead to incorrect density profiles at the edges of the system. For the given system size, at least 200 digits precision is required to faithfully represent the density. A standard double type variable uses 64 bits to store the value, which translates to around 16 decimal digits. (1) and (13) with PBC for (a) g < gc, (b) g = gc, and (c) g > gc. The computations were performed for U = −1, 0 and 1 and at fixed system size L = 12 (24 orbitals).

Appendix H: Many-body computation of the imbalance
In this Appendix, we briefly introduce some convenient formulas for computing the imbalance and its variants in many-body systems. The imbalance I can be expressed as: where |Ψ(0) is the Hermitian ground state. From this structure, it is mathematically more convenient to work withĨ This is a minor generalization of the imbalance introduced in the main text, with essentially the same properties. It allows us to write: ∂ g log Ψ(0)| exp(2gLĨ 1 )|Ψ(0) . (H3) I 1 is therefore proportional to the derivative of the cumulant generating function of LĨ 1 . We numerically check that in the non-interacting case, the probability distribution is which explains our results.