Liouvillian Skin Effect in an Exactly Solvable Model

The interplay between dissipation, topology and sensitivity to boundary conditions has recently attracted tremendous amounts of attention at the level of effective non-Hermitian descriptions. Here we exactly solve a quantum mechanical Lindblad master equation describing a dissipative topological Su-Schrieffer-Heeger (SSH) chain of fermions for both open boundary condition (OBC) and periodic boundary condition (PBC). We find that the extreme sensitivity on the boundary conditions associated with the non-Hermitian skin effect is directly reflected in the rapidities governing the time evolution of the density matrix giving rise to a Liouvillian skin effect. This leads to several intriguing phenomena including boundary sensitive damping behavior, steady state currents in finite periodic systems, and diverging relaxation times in the limit of large systems. We illuminate how the role of topology in these systems differs in the effective non-Hermitian Hamiltonian limit and the full master equation framework.


I. INTRODUCTION
Topological phenomena in the non-Hermitian (NH) realm has attracted ample interest during the past few years [1][2][3]. Compared to their conventional Hermitan counterparts [4][5][6], NH effective Hamiltonians exhibit an entirely different catalog of gapped and gapless topological phases [1][2][3][7][8][9][10][11][12][13][14]. The arguably most dramatic effect is caused by a macroscopic piling up of eigenstates at the boundaries of the system [15][16][17], leading to a spectral sensitivity that grows exponentially with system size when coupling the boundaries [18]. This phenomenology has been dubbed the NH skin effect [19] whereby the celebrated bulk-boundary correspondence in Hermitian systems is replaced by a dichotomy described by either non-Bloch band invariants [19] or in terms of a biorthogonal bulk boundary correspondence [18]. This has led to a blossoming field of research [3,.
Experiments displaying the NH bulk-boundary correspondence have so far been limited to classical systems including mechanical [40,41], electrical [42,43], and photonic [44,45] platforms. Moreover, the remarkable sensitivity to boundary conditions has recently been suggested to be harnessed in applications such as NH topological sensors [38], and a quantum input-output theory of such systems has been developed [46,47]. A fully consistent quantum mechanical description in terms of Lindbland master equations [48] appropriate for Markovian dissipative systems [49][50][51][52][53][54] has earlier been studied and fruitfully employed in the context of preparing or stabilizing Hermitian topological phases [55][56][57][58][59]. Recent pioneering work has highlighted that also in such fully quantum mechanically consistent descriptions, the sensitivity to boundary conditions remains [60][61][62][63][64][65][66][67], but it is fair to say that a comprehensive understanding is still lacking.
Here we advance the understanding of this problem by providing a complete analytical solution of a class of dissipative fermionic chains for both open and periodic boundary conditions (Fig. 1). Remarkably, we find that there is a NH SSH Hamiltonian H S (2.8) that fully diagonalizes the Liouvillian: All normal modes can be expressed in terms of the eigenstates H S , and the rapidities, β m , of the Liouvillian are simply related to the energy eigenvalues E m of H S as β m = const + iE m . (1.1) This directly implies that the topological properties and the skin effect of H S carry over to the quantum context, including fluctuations and quantum jumps, although the interpretation, as we illustrate, is somewhat altered and there is no skin effect in the steady state reached at sufficiently long times. shows the Hamiltonian HS in Eq. (2.8) whose eigenstates diagonalizes the full dissipative model (a), and (c) illustrates the effective short-time Hamiltonian H eff in Eq. (2.57) describing the short-time dynamics before any quantum jump has taken place. Intriguingly, both (b) and (c) are NH SSH models although their asymmetric hopping parameters γi and ηi in Eqs. (2.12) reflect different aspects of the microscopic dissipative processes in (a). γi evaluates the total strength of loss and gain dissipations on a given bond, while the imbalance between the two gives rise to ηi. All models exhibit skin effects as discussed in the text.
The effective Hamiltonian, H eff , of our dissipative model ignoring quantum jumps is also of the form of a NH SSH model in Eq. (2.57) though the effective parameters reflect different aspects of the underlying dissipation as illustrated in Fig. 2. Using recent insights into exactly solving the full set of eigenstates and eigenvalues of NH SSH models [34], this allows us to carry out a detailed study of dynamical phenomena in the Lindblad setting, comparing to the much more studied NH phenomenology and revealing several interesting features, including anomalous damping behavior and diverging relaxation times (in large systems).
The structure of the paper is organized as follows. In Section II, we introduce the fermionic bond-dissipative SSH model in a Majorana fermion representation. The exact solutions to the quadratic Lindbladian, which can be mapped to a generalized NH SSH chain up to a total rapidity shift, are constructed according to the periodic and open boundary conditions. We point out the differences and links between the full Liouvillian spectrum and the energy spectrum of the truncated effective Hamiltonian that neglects quantum jumps. In Section III, we analyze the configurations of the non-equilibrium steady state (NESS). When the boundary is closed, the Liouvillian gap vanishes in an intermediate window of hopping amplitudes and a persistent current flow is identified there. In Section IV, we discuss the Liouvillian skin effects in an open quantum system. The anomalous quantum dynamics are manifested in the relaxation of the current, the modulated damping behavior, and the lifetime dependence of an edge mode on the system size.

II. THE MODEL
Our model is built on a bond-dissipative SSH chain of spinless fermions, shown in Fig. 1. There are hopping terms along the chain with alternating amplitudes t 1 and t 2 : H = j t 1 a † j,A a j,B + t 2 a † j+1,A a j,B + H.c. In the presence of single-particle loss and gain bond dissipations, the set of Lindblad operators takes the general form: L l = (j,α) f l (j,α) a j,α and L g = (j,α) f g (j,α) a † j,α . a † j,α (a j,α ) creates (annihilates) a fermion in the jth unit cell that belongs to the sublattice α = A or B and they satisfy the fermionic anti-commutation relations: Under different boundary conditions, we choose a total of n = 2N (2N − 1) sites for PBC (OBC) with N an integer. Subjected to open boundary, the last unit cell is broken with an empty B site. The odd number of sites helps to stabilize a zeroenergy boundary mode and eventually leads to an exact Liouvillian spectrum.
The full dynamics of the generic model is captured by the Lindblad equation [48,54], where µ denotes the summation over all types of Lindblad dissipators.
It is important to note that with linear dissipators, the master Eq. (2.1) generates quadratic Lindbladians, from which the relaxation process of any observable can be studied by solving the equation of motion either in the Majorana fermion representation [55,56,[67][68][69][70][71][72][73] or in the original complex fermion basis [60,63]. While the two approaches are equivalent, we employ the first method due to its advantage of obtaining the exact Liouvillian spectrum instantaneously.

A. Liouvillians of non-interacting fermions
Before addressing the specifics of our model, we give a brief review of this general third quantization approach through which the quadratic Lindbladians are diagonalized in the canonical basis of Majorana fermions [68][69][70].
Starting with n complex fermions, the reduced density matrix of the system ρ = jk ρ jk |a j a k | lives in a Hilbert space of dimension 2 n × 2 n . One can always construct a pair of Majorana fermions out of one complex fermion, for instance, ω j = a † l + a l and ω k = −i(a † l − a l ) such that they become their own anti-particles ω † j = ω j and satisfy anti-commutation relations {ω j , ω k } = 2δ j,k . The Hamiltonian and the dissipators take the matrix form depending on the details of the mapping: H = j,k w j H j,k w k , L ν = j l ν j w j . The Hilbert space is now represented in the 2 2ndimensional Liouville space K expanded by the new set of Majorana operators: P α = w α1 1 w α2 2 · · · w α2n 2n with α j ∈ {0, 1}. We apply the notation x = (x 1 , x 2 , . . . , x k ) T to represent a vector of scalars or operators. Over the space K, it is convenient to define the adjoint creation and annihilation linear maps through ϕ fermions: The operator P F = (−1)N denotes the Fermi parity associated with the fermion numberN = 2n j=1 ϕ † j ϕ j = ϕ † · ϕ. Since [L, P F ] = 0 and (P F ) 2 = 1, the parity is conserved and takes the value ±1.
Since a physical observable must contain an even number of fermionic operators, one may focus on only a definite parity sector, e.g., the even parity sector (P F = 1) in this paper. Correspondingly, the dynamics of the system is governed bŷ . After proper diagonalization, one is able to express the Liouvillian in terms of rapidities β m and normal master modes (NMMs) b m , b m : with the band index m and NMMs satisfying the anticommutation relations {b m , b l } = δ m,l . Notably, the complex rapidity spectrum contains rich physics. An initial state ρ 0 of positive fermion parity approaches the NESS after a long-time evolution: ρ ss = eL + t ρ 0 t→∞ . While the imaginary part of the rapidity β m encodes the phase oscillation frequency, the real part of β m reveals the relaxation speed of the system to the steady state and Re(β m ) ≥ 0 is required naturally. To describe the asymptotic decay rate quantitively, a spectral Liouvillian gap [66,74] can be defined as (2.5) Extra insight comes from the upper triangular structure ofL + . It infers that the rapidity spectrum must coincide with the eigenvalues of the matrix X in the diagonal block, which is NH and also called the damping matrix. The off-diagonal block Y , on the other hand, shapes the configurations of NMMs and NESS [71,72].

B. Exactly solvable models
Next, we apply the third quantization approach to the bond-dissipative SSH chain in Fig. 1 and derive exact PBC and OBC Lindblad spectra. Compared with previous studies, our model is more general and include the special cases of Refs. [60,73]. We start from a generic set of linear bond dissipators acting on both t 1 and t 2 bonds, Under PBC, the index j runs over N unit cells. When the boundary opens up with the last B site taken away, the associated dissipators are curtailed simultaneously: It turns out that the damping matrix X in the Majorana representation can be transformed to the NH SSH Hamiltonian in Eq. (2.15): The strengths of asymmetric hopping terms γ 1 , γ 2 take the value 2γ i = |γ l i |+|γ g i |.
In the Appendix, we construct the exact spectrum of H S as a direct generalization of Ref. [18,34] to the new limit γ 2 = 0. It enables us to build an exact solution toL + under both PBC and OBC (with a total number of sites n = 2N and n = 2N − 1, respectively), using the eigenvectors of H S .

Majorana representation
As a first step, let us define the Liouville space K by a mapping from n spinless fermions to 2n Majorana particles: (2.9) For later convenience, we regroup c and d Majorana fermions into a whole set {w} under the vector notation: w = (w 1 , w 2 , . . . , w 2n ) T = (c 1 , . . . , c n , d 1 , . . . , d n ) T . Accordingly, in the Lindblad Eq. (2.1), the operators H = j,k w j H j,k w k and M ij = ν=g,l (l ν i,µ ) T (l ν µ,j ) * with L ν µ = j l ν µ,j w j take the following matrix forms: (2.10) H 0 and M 1,2 are n × n matrices holding entries, for instance, under OBC: γ i 's and η i 's stand for the sum and the imbalance of loss and gain dissipations: (2.12) It can be seen immediately that the Majorana representation in Eq. (2.9) is better adapted to diagonalize the Liouvillian in Eq. (2.3): in the adjoint fermion basis ϕ T = (ϕ T c , ϕ T d ), it incorporates the matrix blocks with X c = X d = −4iH 0 + 2M 1 . Notice that the damping matrix X is diagonal and depends only on the total strength of dissipations γ i . By contrast, the matrix Y is off-diagonal (thus couples ϕ c and ϕ d -fermions) and depends only on the imbalance of gain and loss η i . Taking into account the identical structure shared by X c and X d , one concludes that the rapidity spectrum determined by the full damping matrix is at least doubly degenerate, where β m = β c,m = β d,m . It is not difficult to find a unitary transformation U n×n = diag{1, i, 1, i, . . . , 1, i, 1} under which X c(d) is mapped to the generalized NH SSH Hamiltonian in Eq. (2.8): (2.15) The matrix form H S is given in Eq. (A8). We thus reveal one important relation aforementioned in Eq. (1.1) for the rapidity spectrum of our model: where E m represents the eigenvalues of H S . It can be checked directly that the equality in Eq. (2.16) holds true for the PBC spectrum as well.

Changing boundaries from PBC to OBC
We proceed to construct the complete set of eigenvectors of the damping matrix based on the mapping in Eq. (2.15). Let us write the generic eigenvalue equations of the NH SSH Hamiltonian: of which the exact solutions under different boundary conditions are derived in the Appendix. It renders that the pair of eigenvectors of X c(d) can be constructed as with corresponding eigenvalues β m in consistency with the relation (2.16): Under this construction, the biorthogonal normalization of the left and right eigenstates [18,34,75] are respected: When the boundary is switched from PBC to OBC as depicted in Fig. 1, we can extract the rapidity spectrum directly from E m solved in Eqs. (A3), (A9) and (A14): (2.23) In the OBC spectrum, the band index m ∈ {(±, q), 0} is assigned to n = 2N − 1 bands with discrete modes q = πm /N , m = 1, 2, . . . , N − 1. Given an odd number of sites, there emerge two right and left boundary modes at zero energy E 0 = 0 with exponential localization factors: As for the bulk spectrum, analogous to the simplified NH SSH model with only one asymmetric hopping term γ 1 [19,20,34], one finds up to a shift in q a general relation of rapidities between two boundary conditions: It is important to note that under OBC, in the region |r * L r R | < 1 or, equivalently, the complete set of eigenvectors exhibits a non-zero biorthogonal polarization [18] and holds a non-trivial non-Bloch topological invariant [19]. It further indicates at |r * L r R | = 1, the gap of E m closes [an alternative argument is given above Eq. (A13)]. As a result, when H S in Eq. (2.8) reduces to the Hatano-Nelson model at t 1 = t 2 , [76][77][78], a collapse of the exact bulk states is expected. One is nevertheless able to study the behavior of the system around these gap closing points via approximate variational states [34]. Meanwhile, in a dissipative quantum system composed of tight-binding bosons and non-linear (quadratic) Lindblad operators, a Liouvillian can be constructed in such a way that its diagonal subspace resembles the Hatano-Nelson model, thus giving rise to a similar Liouvillian skin effect [66].

Exactly solvable NMMs
We are now prepared to get an analytical set of NMMs for the Liouvillian in Eq. (2.3), in particular, under the open boundary condition. The essence is to remove in the upper triangular structure the off-diagonal block Y that entangles adjoint fermions ϕ c and ϕ d . More precisely, We find a solution for the transformation above where the covariance matrix C satisfies It is easy to discern that if η/γ = η 1 /γ 1 or η 2 /γ 2 , the covariance matrix holds a simple structure By definition in Eqs. (2.12), the solvable limit encompasses the following possibilities: A deeper understanding of the covariance matrix comes from the pairing function of Majorana fermions: C jk (t) = −Tr[w j w k ρ(t)] + δ jk . In its time evolution, by applying the Lindblad Eq. (2.1) in combination with the anti-commutation relations of Majorana fermions, one arrives at For a steady state, ∂ t C s = 0. Therefore, It can be interpreted that in the solvable limit of the Lindbladian in Eqs. (2.30), the covariance matrix encodes a stationary pairing pattern favoured by Majorana fermions at large times: A non-vanishing covariance matrix plays a role in counteracting the effect of the imbalanced gain and loss dissipations (η). Back to the original Hilbert space expanded by spinless fermions a, a more physical picture can be drawn as follows. Since n j,α = a † j,α a j,α = (1 − ic j,α d j,α )/2, the dissipative fermionic SSH chain has a tendency to evolve towards a uniformly distributed configuration with an occupation number: There is no correlation between a fermions. In general, one can express the steady state in the form of a product state, where θ j denotes an arbitrary phase difference.
In the end, with a covariance matrix fulfilling the transformation in Eq. (2.26), we can relate the NMMs in Eq. (2.14) to the left and right eigenvectors of the damping matrix in Eqs. In the same manner, under PBC, the covariance matrix takes a simple analytical form in the solvable limit in Eqs. (2.30), 37) and the set of NMMs associated with rapidities β PBC ν=± (q) is given by (2.38)

C. Spectrum and topology
In this section, we reveal the topology of the Liouvillian starting from the spectral winding of the NH damping matrix. Known for a NH system, the conventional bulkboundary correspondence is broken. Yet the prevalence of the exceptional points (EPs) in the OBC rapidity spectrum implies the Liouvillian skin effect, of which more physical consequences will be discussed in Section IV. We further analyze the formation of the Liouvillian gap that is found to be highly sensitive to the boundary conditions.

Spectral winding number and exceptional topology
Based on different classification schemes [73,79], open fermion matter falls into one of the ten NH Bernard-LeClair symmetry classes. Given a quadratic Liouvillian [73], the classification can be defined through the damping matrix X, or, in our case, X c(d) : The matrix Z resembling the Hamiltonian in the closed limit preserves the time-reversal, particle-hole, and pseudo-anti-Hermiticity (PAH, or generalized chiral) symmetries: Hence, each subspace of the damping matrix belongs to class BDI with a Z classification in 1D.
It should be noted that at the edge, two real Majorana fermions c and d recombine into one complex fermion a as indicated by Eq. (2.9). When the dissipations are turned on, this edge mode shares a finite lifetime with a contribution coming from the effective Liouvillian gap (non-vanishing as indicated by the purely imaginary total energy shift in Z) and another from the non-Hermiticity of the damping matrix [see Eq. (4.9)].
In the presence of the PAH symmetry, a Z classification is captured by the topological invariant, spectral winding number. By shifting the reference point to (0, −iγ) in the complex plane, it is equivalent to evaluate the winding of the Bloch Hamiltonian H S (q) [2,3]: In the framework of our model, the bond dissipators entail γ 1 ≥ 0, γ 2 ≥ 0. Figure 3 shows the dependence of the spectral winding number on one of the dissipation strengths with the signs of two symmetric hopping terms being either the same (t 1 = t 2 ) or the opposite (t 1 = −t 2 ). As soon as |ν| = 0, the left and right boundary modes exponentially pile up at different ends. This unique feature of H S can be applied to the design of the NH topological sensors exhibiting anomalous sensitivity that grows exponentially with the system size [38]. Inside the bulk spectrum, however, it is clear to see that the topological invariant ν obtained from the Bloch Hamiltonian fails to locate the boundary zero modes at . In a closed NH system, the concept of conventional bulk-boundary correspondence has thus been generalized to allow the reconstruction of topological quantities in the biorthogonal basis [18,19] such that the occurrence of the boundary modes are accurately predicted [see also Eq. (2.25)].
What happens to an open quantum system? Similarly, its relaxation dynamics have the Liouvillian skin effect once the NMMs of the Liouvillian pile up exponentially close to the boundary [60,66]. It is then crucial to look at the behavior around the EPs, arising naturally from the Lindblad master equation [80]: At EPs, the geometric multiplicity of the Liouvillian is smaller than the algebraic multiplicity which holds an order that scales with the system size [1,16,20]. We check that for t i = ±γ i , the rapidity spectrum β OBC m in Eqs. (2.23) indeed has one or three eigenvalues. Approaching one of the EPs, the set of NMMs merges into one or three linearly independent eigenstates. Moreover, after the mapping of the damping matrix to H S in Eq. (2.8), the adjoint fermions are only permitted to hop in one direction when t i = ±γ i , so all NMMs become exactly localized at that one end. Consequently, we envision the most drastic Liouvillian skin effect in close proximity to EPs, which, in terms of the momentum shift parameter r defined for linking two rapidities in Eq. (2.24), is manifested as As expected, with a sufficiently large system size, the Liouvillian skin effect is fully determined by the localization behavior of the bulk modes [see also Eqs. (A16)−(A18)] and the influence of the boundary modes is negligible. By varying the dissipation strength, alongside the spectral winding number, Fig. 3 compares the responses in different localization determinants: r j R (r j L ), r j (r −j ) for the piling of the right (left) boundary and bulk modes at unit cell j. To conclude, while not characterized by the Bloch topological invariant, the Liouvillian skin effect is embodied in the exceptional topology of the Liouvillian, or, more precisely, the damping matrix.

Liouvillian gap
We go on to study the development of the Liouvillian gap in relation to various hopping amplitudes and dissipation strengths together with its response to different boundary conditions. Apparently, the imaginary part of complex energy E in Eqs. (2.23) is bounded by ±(γ 1 +γ 2 ). Hence, ∆ = 2 min{Re[β m ]} ≥ 0. Figure 4 shows the real part of the rapidity spectrum as a function of t 1 under PBC and OBC.
Subjected to a periodic boundary, with finite bond dissipation inside unit cells γ 1 = 0, one observes drastically different gap closing behaviors with and without the second inter-unit-cell dissipation: (2.45) We find that after the introduction of the second bond dissipation, the gap closing line of the Liouvillian discovered in Ref. [60] becomes unstable and shrinks to a point in the phase diagram. It would lead to a collapse of the non-trivial (quasi-)NESS and a termination of the steady-state current in a wide parameter range (see more details in Section III B). Meanwhile, the gap-closing points vanish completely when the boundary opens up. In this circumstance, we resolve the Liouvillian gap from the exact spectrum in Eqs. (2.23): Here, θ(x) denotes the Heaviside step function which is defined as: Meeting the EPs at t i = ±γ i , a discontinuity appears in the derivatives of the gap with respect to the symmetric hopping term ∂ ti ∆ OBC (see Fig. 4) as well as the asymmetric hopping term ∂ γi ∆ OBC (see Fig. 7). Lastly, we mention briefly the effects of single-site loss and gain on the Liouvillian. Given a general set of loss and gain dissipators acting on individual sites, the only change to the Liouvillian becomes a modified constant term in the damping matrix that ultimately lifts the minimum of the Liouvillian gap, The gap now can no longer be removed by closing the boundary. Fortunately, the onsite dissipations will not add any NH term to the damping matrix, thus not altering the Liouvillian skin effect. Therefore, we implicitly assume γ l 0 = γ g 0 = 0.

D. Correlation function and comparison with the effective Hamiltonian
In this section, we derive a closed form of the singleparticle correlation function from the exact eigenmodes of the damping matrix. The rapidity spectrum determines the time-dependent part of the two-point correlator. Any observable consisting of even number fermionic operators can then be constructed by Wick's theorem. From the perspective of the correlation function, we compare the physics of the effective Hamiltonian that neglects Lindblad quantum jump operators with the picture of the full Liouvillian. Notably, when the bonds are subjected to purely loss dissipations, the two mechanisms become identical. (2.49) Starting from an arbitrary initial configuration that is not trivialC(0) = 0, we can integrate the above equation and implement the diagonalized damping matrix in the exponential: Taking into account the biorthogonality of the basis and the fact that the damping matrix is real, X * = X, we arrive at (2.51) At t = 0, without loss of generality, throughout the text we choose the system to be in a static configuration with each site completely filled: |Ψ 0 = ntot j=1 |1 j , which corresponds tõ Apparently,C(0) selects µ = µ . With no pairing between Majorana fermions of the same species in the initial state, c j c k t = d j d k t = 0. Next, we go back to the physical space and define the single-particle correlation function in the spinless fermion language: Q jk (t) = Tr[a † j a k ρ(t)]. After the mapping to Majorana fermions in Eq. (2.9), one rewrites it in terms of the pairing function, with n the total number of sites. The phase factor depends on whether the correlation resides on the same sublattice or not: (2.55)

Effective Hamiltonian in the absence of quantum jumps
For non-Gaussian Lindbladians, on the other hand, it is useful to study the short-time dynamics by ignoring quantum jumps in the Lindblad equation [55]: Here, we compare the effective Hamiltonian description with the full Lindblad master equation framework.
Without quantum jumps, the time evolution of the density can be described by where It turns out that the structure of the effective Hamiltonian becomes drastically different from the damping matrix [60]. We verify that rather than the total strength γ, the non-Hermiticity of H eff is related to the imbalance η between loss and gain dissipations: The purely imaginary energy shift scales with the size of the system: It is convenient to resolve the single-particle correlator directly from H eff according to the equation of motion in Eq. (2.56), with an effective rapidity spectrum: To work with the same basis as the damping matrix, we have applied the transformation in Eqs. (2.20). In the rapidity spectrum, s 0 prevents the Liouvillian gap from turning negative when η = i=1,2 (|γ l i | − |γ g i |)/2 < 0. By neglecting the quantum jumps in the Lindblad equation, we find that the dynamics of the dissipative system are not properly captured by the effective Hamiltonian at all times. For the short-time interval, regardless of the total dissipation strength γ, the imbalance η determines the Lindblad spectrum and makes the Liouvillian gap increase with the system size rather than stay a constant value as suggested by Eq. (2.16). In the long-time limit, the EPs of the exact solutions move to t i = ±η i , leaving the Liouvillian skin effect unpredictable compared with Eq. (2.44). Ultimately, the system always decays to an empty chain as there is no residual matrix Y in the equation of motion that can add up to a finite stationary occupation of fermions according to Eqs. (2.31)−(2.35).
In spite of all the discrepancies, however, once gain dissipators are suppressed on the bonds γ g 1 = γ g 2 = 0, we reach one special point where γ = η, Q jk,eff (t) = Q jk (t). (2.60) It infers that for the SSH chain with bond dissipations that only lead to losses, though truncated, the effective Hamiltonian encapsulates the full dynamics at arbitrary times.

III. NON-EQUILIBRIUM STEADY STATES
In Section II, we have already revealed the configuration of the trivial steady state in Eq. (2.35). It is a unique NESS under OBC and persists in PBC as long as the Liouvillian gap is not vanishing. With a gapless rapidity spectrum, NESS becomes degenerate and carries a stationary current that is independent of the dissipation strengths when the gain and loss contributions are in balance.

A. Open boundary: Uniqueness of NESS
First, we look at the steady state from the OBC rapidity spectrum. As shown in Fig. 5(a), every NMM has When |t 1 | ≥ γ 1 , all the bulk and edge modes stabilized by the open boundary share the same Liouvillian gap 2γ 1 . Whereas for |t 1 | < γ 1 , the Liouvillian gap decreases but stays positive as long as t 1 = 0. In this regime, the modes with the slowest decay rate appear at m = (ν, q) ∈ {(+, 0), (−, π)}.

B. Periodic boundary: Degeneracy, quasi-NESS, and stationary current
From Section II C, the Liouvillian gap can be closed by switching the boundary condition to PBC, thus lifting the degeneracy of NESS. For non-zero γ 1 and γ 2 , Fig. 4 shows the real part of the PBC rapidity spectrum holds a gapless point at t 1 = t 2 . It is easy to check that the rapidity of the bulk mode m * = (+, −π) vanishes completely: β m * = 0. Therefore, the right set of the steady states becomes three-fold degenerate: s R,0 = |NESS , s R,c = b c,m * |NESS , s R,d = b d,m * |NESS (accordingly, the left set is expanded by NESS |, NESS |b c,m * and NESS |b d,m * ). In contrast to OBC, after a long-time evolution, the final state of the periodic system now depends on the initial configuration and may appear as a superposition among different NESSs: ρ ss = ρ ss (ρ 0 ).
Compared to the trivial steady state in Eq. (2.35), the degenerate and quasi-NESSs can be viewed as plane waves of fermions with fixed momenta m * on top of a static uniform occupation. It is rather important to distinguish the degenerate and quasi-NESSs from the trivial one, especially on account of the former two being a direct consequence of the closing of the Liouvillian gap. We find the current flow [72,81] is such an ideal observable. Defined as j c (t) = ( i ntot ) j [ a † j a j+1 t − a † j+1 a j t ], the time-dependent current flow can be conveniently ob-tained from the single-particle correlator in Eq. (2.55): Our focal point is to study the behavior of the current in different NESSs. At larger times, only the Liouvillian gapless modes m * satisfying Re(β m * ) = 0 survive. It leads to Let us begin with the special limit γ 2 = 0 where all three types of NESSs coexist and assume t 2 > 0. For the trivial steady state m * = ∅, thus the current vanishes in the end j ss | γ2=0 = 0, |t 1 | > t 2 . (3.5) At phase transition points, the degenerate NESS generated by two bulk modes m * = (+, − π 2 [sgn(t 1 t 2 ) + 1]) supports a current flow with an amplitude In the intermediate gapless region, there emerges quasi-NESS from the modes m * = (±, ± arccos[−t 1 /t 2 ]). Analytically, the stationary current can be expanded in the orders of 1/N : with α = 1 − (t 1 /t 2 ) 2 . The time dependence in the steady steady current arises from the finite phase oscillation frequency, a unique feature possessed by quasi-NESS. Given sufficiently large system size, to the leading order O(N −1 ), we get For non-zero γ 1 and γ 2 , on the other hand, the steadystate current is carried by the gapless bulk mode m * = (+, −π) at t 1 = t 2 , while it vanishes elsewhere. Therefore, One immediately notices that when the gain and loss dissipations are in balance, namely η = 0, the prefactor (γ + η)/γ → 1. The steady state current is then independent of the dissipation strengths γ i . In Fig. 6, we confirm the analytical predictions on the steady state current for γ 2 = 0 in Eqs. (3.5)−(3.8) by a measurement of j c (t) at time γt = 10 5 . With balanced gain and loss bond dissipations, we verify the stationary current remains the same under different values of γ 1 . When γ 2 = 0, in the same manner, we still find a persistent current at the Liouvillian gap closing point regardless of the choices of γ 2 [see also Fig. 8(a)].
By contrast, in the description of the effective Hamiltonian, once γ = η the real part of the rapidity spectrum for a periodic chain in Eq. (2.59) is always gapped and a persistent current will not be observed in any allowed parameter regime.

IV. ANOMALOUS QUANTUM DYNAMICS
In this section, we search for dynamical signatures of the Liouvillian skin effect in dissipative quantum systems, originating from the piling up of the NMMs exponentially close to an open boundary. Compared with previous studies [60,63,67], we show the relaxation behaviors directly obtained from our exact solutions for an odd number of sites n = 2N − 1 and, at the same time, include the impact of the second bond dissipators. Apart from a diverging lifetime without gap closing [66], we find that other global observables such as a tail of a dynamical current flow and a chiral damping wavefront center can also serve as good probes of the Liouvillian skin effect, the nature of which will be related to the exceptional topology of the damping matrix.

A. Relaxation of current flow
As revealed in Section III B, for a general set of bond dissipators (γ 1 = 0, γ 2 = 0), a non-vanishing current in the steady state in Eq. (3.9) helps us to distinguish the degenerate NESS from the trivial one on a periodic chain. Figure 7 depicts the relaxation process of the current before reaching the steady state under different boundary conditions. For PBC, the closing of the Liouvillian gap at t 1 = t 2 shown by Fig. 7(a) sustains a stationary current in Fig. 7(c) which, after averaging over all sites, scales with the inverse of the system size ∼ 1/N . Deviating from the gap closing point, the current flow vanishes at relatively short times. Changing the boundary condition to OBC, the gapless mode in the rapidity spectrum disappears. All bulk and edge NMMs immediately pile up exponentially at one of the boundaries, thus terminating the current flow in Fig. 7(c). If we zoom in to look at the region where the Liouvillian gaps of two spectra are comparable [for instance, γ 2 = 0.5 in Fig. 7(b)], the behaviors of the current flow seen from Fig. 7(d) turn out to be less sensitive to the boundary conditions. It infers that from the perspective of the current, the Liouvillian skin effect is better captured when the gap of the PBC rapidity spectrum is closed. In Fig. 8, we thus fix t 1 = t 2 and study the dynamics of the current flow occurring at the Liouvillian gap-closing point in a wide range of dissipation strengths. Consistent with our earlier prediction in Eq. (3.9) for balanced gain and loss dissipators (η = 0), Fig. 8(a) shows that even in the presence of very weak dissipations (γ 1 = 0.02, γ 2 = 0.01), the current flow along a periodic chain saturates to a finite value identical to the limit of strong dissipations. However, once the boundary opens up, driven by weak dissi-pations, the current flow decays with oscillations and has a much shorter relaxation time as indicated in Fig. 8(b). It also represents the behaviors of those points far away from EPs (t i = ±γ i ), bringing about weak Liouvillian skin effect. For stronger dissipations or | ln(|r 2 |)| 0, the amplitude of the current is enlarged as the NMMs continue to pile up in the same direction. The current flow also shares a longer relaxation time before the system evolves to the trivial NESS hosting a static uniform distribution of fermions. The relaxation time or the tail of the current flow is determined by two factors, the effective Liouvillian gap ∆ eff and the correlation length ξ(r 2 ) of the system. It is equal to the lifetime of the particle at one boundary in favour of the Liouvillian skin effect, of which more detailed analysis can be found in Section IV C.

B. Density evolution under damping
Next, we focus on the chiral damping phenomena in the particle-number distribution of the dissipative chain [60]. We find the damping itself displays the Liouvillian skin effect beyond the gapless point t 1 = t 2 . By turning on the second bond dissipators γ 2 , we show the center of the chiral damping wavefront can be tuned via the localization parameter r 2 of the bulk NMMs. It enables us to associate the Liouvillian skin effect with the existence of When the open boundary builds up, the current vanishes immediately after reaching the peak. Its relaxation time, equal to the lifetime of the particle at the edge under the Liouvillian skin effect, depends on the system correlation length ξ(r 2 ) and the effective Liouvillian gap ∆ eff (4.9).
the EPs in the damping matrix at which |r 2 | → 0 or ∞. We further clarify the sensitivity of dynamical observables to the boundary conditions at longer times. Our results are consistent with one of the earlier studies [63], demonstrating that probes in the short-time domain are not sufficient to distinguish the relaxation processes of open quantum matter under different boundary conditions.

Link with exceptional topology
First, we associate the chiral damping of the particlenumber operator with the exceptional topology of the damping matrix. In the full Lindblad master equation framework, the evolution of the particle number to NESS, n j (t) = n j (t) − n j,ss can be built on our exact single particle correlator a † j a j t from Eq. (2.55): Under OBC, a large system size enables us to safely neglect the contributions from the boundary eigenmodes in Eq. (A9). Given exact solutions of the left and right bulk eigenmodes in Eqs. (2.21), (A17) and (A18), one derives an asymptotic scaling of particle number evolution for sites in the odd sublattice A: A similar expression can be obtained for the even sublattice B. Here, the damping factor r 2 comes from the localization of the left bulk eigenmodes. When |r 2 | → 0(∞), for instance, the left eigenmodes are localized on the right (left) end of the chain, pushing the chiral wavefront towards the same boundary while the right eigenmodes center on the opposite. In the time-dependent part of the particle number operator, we replace the original gap ∆ OBC with an effective Liouvillian gap ∆ eff . The comparison of two quantities is made in Fig. 13. Since we are interested in the full dynamics of the relaxation process, instead of taking the long-time limit and keeping track of the slowest decaying mode, we should include all bulk modes with various decaying rates. The summation then leads to an effective Liouvillian gap that is also useful in the evaluation of the particle lifetime in Eq. (4.9). The two extremities of the Liouvillian skin effect |r 2 | → 0 or ∞ occur at EPs: t i = ±γ i , when the adjoint fermions in the damping matrix are only allowed to hop in one direction after the mapping to H S in Eq. (2.8). It is noteworthy that with a large system size, the chiral damping condition can be further relaxed to irrespective of whether the PBC rapidity spectrum is gapless or not. So far, one discerns that within the framework of the full Lindbladian, the Liouvillian skin effect is closely related to the exceptional topology of the non-Bloch damping matrix, whereas it does not depend on the topology of the Bloch NH Hamiltonian H S (q) as illustrated in Fig. 3. Meanwhile, considering γ i 's and η i 's belong to two sets of independent parameters in Eqs. (2.12), the effective Hamiltonian in Eq. (2.57) does not play a role here neither. The EPs of H eff are shifted to t i = ±η i and the topological regime of H eff (q) for η 1 η 2 ≥ 0 lies in Eqs. (2.42) with a replacement: γ i → η i .

Modulation of chiral wavefront center
We embark on the characterization of the chiral damping behavior emerging in the dissipative SSH chain when subjected to an open boundary. Let us define the polarization of the damping process according to where, as depicted in Fig. 1, the length of chain is chosen to be odd n tot = 2N − 1. Without loss of generality, we restrict ourselves to the case when the gain and loss dissipations are in balance: η 1 = η 2 = 0. Starting from a completely filled chain, the NESS should then be half filling at each site in Eq. (2.34). Plugging the analytical expression for particle number in Eq. (4.1), the tendencies of polarization under different hopping amplitudes and dissipation strengths are shown in Fig. 9. Figure 10 further illustrates the motion of the chiral damping wavefront along the chain as time goes by. At the initial time, ∆P = 1/2, it refers to a uniform particle number distribution consistent with our initial condition. Around |r 2 | = 1, when the Liouvillian skin effect is quite weak, the polarization stays close to 1/2 and the chiral damping wavefront is absent in Figs. 10(e) and 10(b). As |r 2 | 1 ( 1), the NMMs of the Liouvillian start to pile up towards the right (left) end such that the particles closer to that boundary are granted a longer lifetime, thus pushing the polarization to 1 (0). From Figs. 10(a) and 10(f), the chiral wavefront also appears most distinguished in these two limits. In the intermediate parameter regime of |r 2 |, the polarization evolves to a finite value in between [0, 1/2] or [1/2, 1] in accompaniment with a damping wavefront growing obscure, as shown by Figs. 10(d) and 10(c). Therefore, by tuning the parameter r 2 , we are able to modulate the polarization or the center of the chiral damping wavefront regardless of the topology of the PBC rapidity spectrum before the boundary opens up. The left and right columns of Fig. 10 correspond to the gap closing (t 1 = t 2 ) and gap opening (t 1 = t 2 ) points in the PBC spectrum, respectively.
On a side note, in Fig. 9 even at relatively large times, ∆P will not converge to 1/2 for its measurement is targeted on the excess of the particle number over a uni-formly half-filled NESS, rather than the real occupation number.

Boundary sensitivity
From the chiral damping phenomena, the Liouvillian skin effect takes place upon changing the boundary from PBC to OBC. Since the system evolves from a completely filled initial state to a uniform steady state, the damping serves as a global effect. Nevertheless, we can also resolve the motion of a single particle and study the sensitivity of the open system to boundary conditions on broader grounds.
It is argued in Ref. [63] that in the thermodynamical limit, the single-particle Green's function of NH systems is independent of the boundary conditions. We show in the following within the Lindblad master equation framework, the argument is indeed true for the evolution of the density matrix at short times. However, as soon as the motion of a particle involves the edge, due to the gap closing and the Liouvillian skin effect, there emerge drastic differences in the relaxation process under PBC or OBC.
Let us start by putting one particle in the middle of the chain: n j (0) = δ j,N . The initial condition corresponds to a Majorana pairing configuration, where the n×n matrix D only holds one non-zero element D (N,N ) = 1. The time evolution of the particle number at the jth site can be described bỹ ψ Rm (l)ψ Rm (l) .
Also, we focus on the relaxation process at the same point t 1 = t 2 as Ref. [63], where the gap of the PBC rapidity spectrum closes for arbitrary dissipation strengths γ 1 and γ 2 . The inclusion of γ 2 = 0 offers more tunability on the Liouvillian skin effect parameter r 2 . As shown in Fig. 11, the motion of the particle does not depend on the boundary condition until its trajectory hits the edge. Once that long-time evolution is permitted, contrary to Ref. [63], the motion differs a lot under PBC or OBC. We observe a persistent current j ss = 1/N circulating along showing an exponential localization tendency near the boundary: |c| ∼ e O(L/ξ) . Here, ξ denotes the correlation length of the system. It leads to In our dissipative SSH model, we look at the Liouvillian skin effect region |r 2 | < 1 and take the dynamical observable as the occupation of a particle residing at the right boundary. More importantly, we make fewer assumptions than Ref. [66] by defining the lifetime of the non-equilibrium particle according to where l is taken as a positive integer. The precision of the lifetime can be improved by increasing the value of l.
The definition in Eq. (4.8) has the advantage of including the contributions from all decaying modes, staying closer to a real measurement. For simplicity, we start with a completely filled chain and balanced gain and loss dissipators η 1 = η 2 = 0. Then, |ñ 2N −1 (0)| = 1 − 1/2 = 1/2. From the asymptotic scaling of the particle number evolution in Eq. (4.2), we establish that The correlation length of the system in Eq. (4.7) thus satisfies ξ(r 2 ) 2 | ln(|r 2 |)| , (4.10) where the factor of 2 comes from the identification L = 2N − 1. Consistent with Ref. [66], we find the lifetime of the particle at the edge grows linearly with the system size without a closing of the effective Liouvillian gap. Next, based on our exact solutions, we verify numerically the above relations by varying the system size and the skin effect parameter r 2 . It is also important to explore two limits, the weak and strong dissipations where the formation of the effective Liouvillian gap differs.
At weak dissipations γ 1 < |t 1 | and γ 2 < |t 2 |, all bulk and boundary modes share the same Liouvillian gap [see Eqs. (2.23) and Fig. 4]. Hence, ∆ eff = ∆ OBC = 2(γ 1 +γ 2 ). Figure 12(a) plots the numerical scaling of the particle lifetime at the edge τ 2N −1 = τ under different system sizes at γ 2 = 0. The linear dependence on N in Eq. (4.9) becomes more visible as the length of the chain increases. With stronger Liouvillian skin effect |r 2 | 1, the linear scaling holds true for a relatively small system size N ∼ 8. The system correlation length can be extracted from the slope and Fig. 12(b) confirms our analytical prediction in terms of the skin effect parameter in Eq. (4.10). Approaching the EP, γ 1 → t 1 = 1, the small discrepancy results from the fact that the strong polarization of the bulk eigenmodes challenges the numerical precision of the damping matrix decomposition in Eqs. (2.50).
At strong dissipations, on the other hand, from Fig. 4 different bulk and boundary modes form distinct real rapidity spectra. As a consequence, the effective Liouvillian gap incorporating effects of all decaying modes differs from ∆ OBC in Eq. (2.46). By definition, (4.11) Still, we can estimate its value from the established linear relation for the particle lifetime in Eq. (4.9) and at the same time, fix the correlation length by the skin effect parameter r 2 in Eq. (4.10). In Fig. 13, we show the formation of the effective Liouvillian gap beyond weak dissipations: γ 1 > |t 1 |, γ 2 < |t 2 |. Indeed, ∆ OBC provides a lower bound for the effective Liouvillian gap.
On the contrary, when γ = η, the lifetime of a particle at the right boundary under the truncated H eff would scale as in the large N limit. Without the quantum jumps, the lifetime decreases with the system size. As expected, the effective Hamiltonian becomes problematic in describing the long-time dynamics.

V. DISCUSSION
We have studied the Liouvillian skin effect in a bonddissipative SSH model, of which the rapidity spectrum and NMMs are exactly solvable. This has illuminated the relation between the NH skin effect, two different effective Hamiltonians, H eff , H S , and the full quantum master equation description. A number of dynamical phenomena, such as diverging relaxation times, inherited in the quantum setup originate from the anomalous boundary sensitivity of the NH Hamiltonians and have been investigated in detail. This paper also provides a solid platform to resolve the entanglement spectrum and identify transitions between different NH topological phases [82] whose experimental realization is feasible, e.g., in ultracold atoms with a momentum lattice [59].
Dynamical probes of the quantum Fisher information [47] and the application as NH topological sensors [46] provide intriguing questions for future investigation. Meanwhile, the current model reveals the Liouvillian skin effect in the relaxation process while posing the question on the search for steady states that inherit the exceptional topology of the NH damping matrix. It is found that using monitored quantum circuits, periodic measurement allows access to biorthogonal steady state observables [29]. Future works could also be centered around the interaction effects on generic Lindbladians. For instance, the dynamical mean-field theory can be implemented to uncover the interplay of dissipation and FIG. 13. Effective Liouvillian gap (green dots) as a function of γ1 beyond weak dissipations: γ1 > |t1|, γ2 < |t2|. We take t1 = t2 = 1, γ2 = 0.2, and τ → |ñ2N−1(τ )| = e −5 |ñ2N−1(0)| = 0.0034. The fitting of the lifetime τ in Eq. (4.9) is performed with a finite system size N = 8 ∼ 14. The inset shows its linear dependence at γ1 = 2.5, |r 2 | = 0.286. The background corresponds to the real part of the rapidity spectrum generated by the bulk modes (blue curves) and boundary modes (red curve) in a system of N = 14 unit cells. ∆OBC (black curve) provides a lower bound for the effective Liouvillian gap. environmental fluctuations in the open quantum matter that consists of interacting fermions [79] or bosons [83]. A mixture of particle species is a good path to the discovery of further variations of the Liouvillian skin effect. From a broader theoretical point of view, effective field theories with bosonization prove powerful in studying dissipative Luttinger liquids [84]. In presence of weak symmetries, one may tackle even strongly interacting Liouvillians [85]. The inclusion of disorder and chaos may be possible provided the solvable limit in a disorder-averaging SYK Lindbladian [86]. γ 2 = 0. The spectrum recovers the known solution by setting γ 2 → 0 and retains the main features in the new limit.
On one hand, there exist two zero-energy boundary states that are fully suppressed on the sublattice B, with associated parameters r R = −(t 1 − γ 1 )/(t 2 + γ 2 ), r L = −(t 1 + γ 1 )/(t 2 − γ 2 ). The two eigenstates are biorthogonal to each other, which brings about a constraint on the normalization factors: It is easy to discern that, depending on the absolute values of r L and r R , the left and right boundary modes can be localized at different ends: sgn[ln(|r L |)] = sgn[ln(|r R |)]. We find the parameters satisfying this condition all stay in the topological regime in Eqs. (2.42) where the spectral winding number is non-trivial.
On the other hand, to derive the bulk spectrum, we present an intuitive approach by analogy to Refs. [19,34]. Our first step is the identification of the gap-closing points when the boundary is opened. One way is to find a transformation matrix R such that the NH SSH Hamiltonian is mapped to its Hermitian counterpartH S : A proper construction leads to R = R 1 R 2 , where R 1 = diag{1, r 1 , r 1 , r 2 1 , r 2 1 , . . . , r N −1 1 , r N −1 1 }, R 2 = diag{1, 1, r 2 , r 2 , r 2 2 , r 2 2 , . . . , r N −2 2 , r N −1 2 }, and r 1 = (t 1 − γ 1 )/(t 1 + γ 1 ), r 2 = (t 2 − γ 2 )/(t 2 + γ 2 ). The Hermitian SSH chain is embedded with the anisotropic hopping strengths with the gapless phase transition occurring at |t 1 | = |t 2 |. H S inherits from the transformation these gap closing points: An alternative method of reproducing the gapless points in the OBC spectrum is to study the biorthogonal polarization of the NH system [18,34], which changes its integer value at |r * L r R | = 1, in consistency with our result in Eq. (A13).
Under the transformation in Eq. (A11), one can also link the bulk eigenstates of H S to the ones ofH S byψ R = Rψ R andψ L =ψ L R −1 , which indicates the piling up of the right (left) states at one end with an exponential localization factor r j (r −j ), where r = r 1 r 2 = r R /r L . It is equivalent to a change of momentum in the Bloch phase factor e iqj from q to q − i ln(r), thus allowing us to build the OBC bulk spectrum from the PBC one [19,20]. More precisely, at a fixed momentum q, like the PBC spectrum the bulk energies always come in pairs (ν = ±): E m =0 = E OBC ν (q). The two spectra are related by E OBC ± (q) = E PBC ± (q − i ln(r)) = ± t 2 1 + t 2 2 − (γ 2 1 + γ 2 2 ) + 2 (t 2 1 − γ 2 1 )(t 2 2 − γ 2 2 ) cos(q).
Given an even number of sites n = 2N , the OBC spectrum obtained in this way applies only to the large N limit [19]. For an odd number of sites n = 2N − 1, with additional mirror symmetry [34,87], it can be shown that the bulk spectrum in Eq. (A14) becomes exact for any finite N . Let us verify by building the exact left and right eigenvectors. We choose q = πm /N , m = 1, . . . , N − 1, so that in total the index m ∈ {0, (±, q)} reproduces a complete set of 2N − 1 bands. A trial function for the right eigenstates can be written as where the component in the jth unit cell should be linked with the PBC eigenvector in Eq. (A4) by the same momentum shift: ψ Rν (q, j) ∼ r j e iqjψ PBC Rν (q − i ln(r)).