Dissipative Boundary State Preparation

We devise a generic and experimentally accessible recipe to prepare boundary states of topological or nontopological quantum systems through an interplay between coherent Hamiltonian dynamics and local dissipation. Intuitively, our recipe harnesses the spatial structure of boundary states which vanish on sublattices where losses are suitably engineered. This yields unique nontrivial steady states that populate the targeted boundary states with infinite lifetimes while all other states are exponentially damped in time. Remarkably, applying loss only at one boundary can yield a unique steady state localized at the very same boundary. We detail our construction and rigorously derive full Liouvillian spectra and dissipative gaps in the presence of a spectral mirror symmetry for a one-dimensional Su-Schrieffer-Heeger model and a two-dimensional Chern insulator. We outline how our recipe extends to generic noninteracting systems.

Here, we consider an alternative approach that crucially depends on the interplay of both coherent dynamics and dissipation.This approach alleviates both practical and fundamental challenges: with directly accessible ingredients, we can devise a combination of Hamiltonian dynamics and dissipation whose interplay uniquely prepares the system in what is the key hallmark of a topological phase, namely its boundary states.At long times all bulk modes vanish and only the boundary states prevail.A minimal example is shown in Fig. 1.Adding loss on a single site in a Su-Schrieffer-Heeger (SSH) chain evolves the system into a state localized at the same boundary.Our recipe is generic and can be applied to essentially any topological (or non-topological) noninteracting system.It also does not rely on fine tuning.Instead, the success of the construction is rooted in choosing lattices harboring boundary modes with the unique property of being the only eigenstates that vanish exactly on some sublattices [25][26][27][28]-where we engineer loss.In fact, in the SSH example of Fig. 1, the same boundary steady state is obtained whenever loss is applied to one or more of the B-sites.
We derive the aforementioned results analytically, and in the presence of a spectral mirror symmetry we obtain the full exact Liouvillian spectrum of the corresponding Lindblad master equation.We exemplify our construction for SSH chains with even and odd number of sites (Figs.1-3) and for a two-dimensional (2D) Chern insulator hosting exact chiral edge steady states (Fig. 4).We also consider the effect of gain yielding steady state currents and how our exact solutions provide key intuitions for non-solvable systems (cf.Fig. 3).

II. SETUP AND DISSIPATIVE SSH MODELS
We consider the Lindblad master equation [24,29] which captures quantum dynamics of Markovian dissipative systems [30][31][32][33][34], where ρ is the density matrix of the system and µ denotes the summation over all types of jump operators.As a minimal example, we study an SSH model of spinless fermions with an odd number of sites L = 2N − 1, with Hamiltonian H S = N −1 j=1 t 1 a † j,A a j,B + t 2 a † j+1,A a j,B + H.c. We add local loss on the first B site: Here, a † j,A (a j,B ) are creation (annihilation) operators on the sublattice A(B) in the jth unit cell satisfying anticommutation relations: The quadratic Hamiltonian and linear Lindblad dissipators yield a quadratic Lindbladian diagonalizable in the third quantization approach [35][36][37].The damping matrix X which governs the dynamics can be mapped to a NH tight-binding Hamiltonian encoding information of both the original Hamiltonian and Lindblad dissipators (see Appendix A), under a unitary transformation U = diag{1, i, 1, i . . ., 1}.For a single loss on the first B site, γ = γ B = |γ l 0,B |/2, I = diag{1, 1, 0, . . .} and Υ = diag{γ/2, −γ/2, 0, . . ., 0}.The damping matrix X = X c(d) denotes the equal contribution from two Majorana fermions species: a j,A = 1 2 (c j,A − id j,A ), a j,B = 1 2 (d j,B + ic j,B ).It can also be decomposed into its eigenmodes: X = m β m |ψ Rm ⟩⟨ψ Lm | with m the band index.The left and right eigenvectors satisfy the biorthogonal relations [12,38,39]: This decomposition enables us to derive dynamical observables by integrating the Lindblad master equation.The particle number ⟨n j (t)⟩ = Tr[a † j a j ρ(t)] reads (see Appendix A) where ⟨n j ⟩ s = ⟨n j (t = ∞)⟩ = 0 corresponds to the trivial non-equilibrium steady state (NESS) in the presence of loss, i.e. an empty chain, and the initial condition is chosen as a completely filled chain.The eigenvalues β m of the damping matrix coincide with the rapidity spectrum of the Liouvillian [8,35].Its real part encodes the decay rates of different modes towards the NESS and we define the dissipative gap ∆ = 2 min{Re[β m ]} ≥ 0.
For an odd number of sites, H S hosts a zero-energy boundary mode fully suppressed on the B sublattice [28]: with r = −t 1 /t 2 the localization factor and N , the boundary mode is exponentially localized at the left (right) end of the chain.Intuitively, the frustrated nature of the boundary mode indicates its robustness against local loss on any B site.Indeed, with a single loss on the first B site, we find a vanishing rapidity for this boundary mode: where β 0 = 0 and ψ R0 = ψ L0 = ψ 0 while all bulk modes have a finite dissipative gap.It implies that starting from all sites that are filled, through the dissipation the system always selects the boundary mode as the non-trivial steady state.For t → ∞, the particle number becomes To verify our observation, in Fig. 1(b), we calculate the time evolution of the particle number on individual sites by numerically diagonalizing the damping matrix.At sufficiently long times, the boundary mode (red line), with |r| = | − t 1 /t 2 | = 0.5 and site occupation ⟨n j ⟩ s predicted by Eq. ( 5), exponentially localizes (|r| < 1) even when the single B-site loss is placed close to the same left end of the chain.We also observe that with single weak loss (γ/2 ≤ ||t 1 | − |t 2 ||), the bulk dissipative gap is inversely proportional to the chain length: ∆ bulk ∝ γ/N .The single loss of our minimal model is particularly useful for small system size in modern experimental setups.
In the second example, while keeping the dissipationless feature of the boundary mode, we can make the dissipative gap of the bulk modes saturate in the large system size limit by adding loss on the entire B sublattice: Ll j,B = γ l 0,B a j,B , ∀j [red arrows in Fig. 2 (a)].The damping matrix X in Eq. ( 2) holds new entries with I = 1 L×L and Υ = (γ/2)×diag{1, −1, 1, . . ., −1, 1}.The rapidity for the boundary mode in Eq. ( 4) remains rigorously zero, β 0 = 0, and it is thus selected as the nontrivial steady state again.This steady boundary mode with a vanishing Liouvillian gap [denoted by the red dot in Fig. 2(b)] is not included in Ref. [40], while it is referred to as an 'edge dark state' in Ref. [41].Remarkably, we can obtain the bulk dissipative gap analytically by exploiting spectral mirror symmetry under open boundary condition (OBC) with odd sites.The rapidity spectrum of the (2N − 2) bulk modes comes in pairs β m = γ/2 ± iE m , where E m denote the eigenenergies of H NH from the mapping of Eq. ( 2) and we have adapted the notation for the band index m = (±, q) with q = πm ′ /N , m ′ = 1, 2, . . ., N − 1.The spectrum obeys spectral mirror symmetry β OBC ± (q) = β OBC ± (−q), establishing a direct link to the spectrum under periodic boundary con- dition (PBC) [28,39,42]: reflecting the absence of a NH skin effect.We thus obtain . Shown in Fig. 2(b), all bulk modes have a finite gap and the minimum ∆ bulk = 2 min{Re[β ± (q)]} determines the dissipative preparation time for the steady boundary mode: τ ∼ ∆ −1 bulk .In the large-N limit, we find the analytical solution to the bulk dissipative gap: 2 , otherwise.An exact complete set of right and left eigenvectors for the damping matrix X is obtained as well (see Appendix A), where the bulk modes can be viewed as a superposition of two Bloch waves with opposite momenta vanishing on the last B site (which is removed in the odd chain).We are able to analytically resolve the full time evolution of the particle number shown in Fig. 2(c).For weak dissipation with ∆ bulk = γ, the damping wavefront arises from the steady boundary mode, in contrast to previous studies with a major contribution from bulk skin modes [6][7][8][9] (here our bulk modes are delocalized).
Next, we study the effect of small gain Lg j,A = γ g 0,A a † j,A , ∀j on the A sublattice [green arrows in Fig. 2(a)].The damping matrix X in Eq. ( 2) is still exactly solvable with entries γ The boundary mode remains an eigenmode, yet it is associated with a non-zero rapidity β 0 = γ A leading to a finite gap: ∆ 0 = 2γ A .A small gain also elimi- nates the empty state as trivial NESS.In Appendix B, we identify the analytical structure of the non-trivial NESS in the presence of γ A as a mixed state continuously connected to the empty state and the boundary mode when γ A → 0. The system exhibits a steady state current J S = (i/L) j (⟨a † j a j+1 ⟩ s − ⟨a † j+1 a j ⟩ s ), which we can obtain in a closed form [see Eq. (B9)].We compare the analytical result with numerics in Fig. 2(d).Once 0 < γ A ≪ γ B , regardless of the initial conditions, the system eventually evolves to the non-trivial NESS with a localization structure approximating the boundary mode in Eq. ( 4).
We now address the non-solvable model by obtaining the eigenmodes of the damping matrix through exact diagonalization (ED).The simplest scenario is encountered when we consider the SSH chain with an even number of sites L = 2N .For L = 2N − 1, the model is symmetric under the exchange of t 1 and t 2 bonds [see Fig. 2(a protected under the B-sublattice loss [43,44].From the rapidity spectrum in Fig. 3(c), we observe this boundary mode only at |t 1 | < |t 2 | and find that its dissipative gap decays exponentially with the chain length [see Fig. 3(d)].The lifetime of the topological boundary state is enhanced exponentially by increasing the system size: τ ∼ ∆ −1 ∼ exp (αN ) with α > 0, a remarkable feature due to topological protection.

III. GENERALIZATION TO DISSIPATIVE LATTICE MODELS IN ANY DIMENSION
Our recipe for devising unique steady boundary states naturally generalizes to any dimension with localization on boundaries of any co-dimension including surfaces, corners, edges and hinges.This is achieved by inferring results on boundary state solutions that vanish exactly on certain sublattices [26,27].In presence of a spectral mirror symmetry the full Liouvillian spectrum may be obtained [28].Increasing both the bulk dimension and the codimension simultaneously is straightforward, yielding e.g.corner steady states on the breathing kagome lattice [26,45], in full analogy with the SSH chain.Increasing the bulk dimension is suitably done by dimensional extension which is slightly more technically involved.Here, we provide an explicit example of dimensional extension in the preparation of dissipationless chiral edge modes of a two-dimensional Chern insulator on the honeycomb lattice [28,46].
Figure 4(a) illustrates the corresponding cylinder geometry: we impose PBC on the x direction with an even number of sites L x and OBC along the y direction with an odd number of sites L y = 2N − 1.The real nearest-neighbor and next-nearest-neighbor hopping amplitudes are denoted as t and t ′ .The complex hoppings t ′ e iϕ with ϕ = π/2 are allowed only between unit cells along the x direction.For each k x , a mapping to the generic SSH-like model can be established x,0 compared with the actual value).With the B-sublattice loss γ = γ B , the rapidity spectrum consists of two copies (µ = ±1) for distinct Majorana fermions species (c, d).For the topological chiral mode, β µ 0 = −iµϵ A while for the bulk modes, the spectral mirror symmetry 4(c) and 4(d) show the real and imaginary parts of the rapidity spectrum of the chiral mode (red) and bulk modes (blue).The bulk modes have a vanishing dissipative gap in the large system limit at the momenta k * x which satisfies |r(k * x )| = 1 where the chiral mode becomes delocalized and switches sides.Nevertheless for momenta when the chiral mode is isolated within the real rapidity spectrum of Fig. 4(c), the system has a finite instantaneous ∆ bulk to separate the chiral and bulk modes for any system size, and thus inherits the localization structure to the boundary mode of the SSH model with |r(k This feature is illustrated by the steady state particle number of Eq. ( 5) in Fig. 4(b).

IV. DISCUSSION
In this work, we have explored an alternative to ground state cooling and Floquet engineering of topological phenomena [17,[47][48][49][50][51], as well as other finite-time methods to probe boundary states [52,53].Our generic approach utilizes both dissipation and coherent dynamics to prepare boundary states as steady states which can be operated at arbitrarily long times.It circumvents the problem of thermal excitations faced at finite temperature in existing platforms.Compared with other dissipative preparation schemes [1,15], in our recipe coherent dynamics from the Hamiltonian itself is introduced, while the form of dissipation is greatly simplified.
Our exact construction is not only simple, but also very general.By choosing appropriate lattice geometries it carries over directly to a plethora of topological and nontopological boundary states, including Fermi arcs and higher-order states at corners and hinges which all have the desired nodal boundary state structure [25][26][27][28].Go-ing beyond the idealized model, one can break spectral mirror symmetry in the SSH Hamiltonians by adding disorder on the nearest-neighbor hopping terms.We find that the unique localization structure of the boundary mode is a generic property of any odd length chain: it still vanishes on the entire B sublattice without its dissipative gap opening, and is robust against onsite-potential disorder on any B site (see Appendix D).In real experiments, weak disorder on the other sublattice and furtherneighour interactions directly interfere with the boundary mode and should be precisely controlled.The latter can be effectively suppressed on sufficiently deep optical lattices [54] that make the SSH models applicable [53,[55][56][57].
From the theoretical side, our analytical study offers promising outlooks for future exploration as well.One could test targeted boundary state distillation in open systems inherent with Liouvillian skin effects [6][7][8][9], critical phenomena [58,59] and phase transitions [60][61][62].In particular, the preparation of steady boundary modes with NH skin effects would bestow exponentially enhanced sensitivity on quantum sensing devices [63][64][65][66][67].It is equally exciting to extend the current framework from open fermions to bosonic and even hybrid systems via a similar third quantization approach [68][69][70].
Very recently, experiments on light in lossy optical waveguides showed that a similar boundary state preparation is possible at the level of NH effective Hamiltonians in classical settings [43].The present work shows that a similar NH phenomenology is also relevant in the quantum realm and opens up different avenues for quantum control.
We also note that related systems with staggered loss or weak measurements have been considered earlier in the context of quantum walks described by effective NH Hamiltonians [71][72][73].There, however, rather than the preparation of topological boundary states the focus was on other aspects such as defects, phase transitions and mean displacement.
Recently, we became aware of two independent works numerically observing the amplification of states at one boundary in Chern insulator models with gain and loss [41,74].In these studies, the boundary modes are classified and dynamically selected [43,44,75] as long-lived modes on even-length lattice models, but not as steady states once the system size becomes finite.By contrast, our predominantly analytical work focuses on oddlength lattice models and achieves this goal with generic boundary states hosting an exact zero dissipative gap.We also consider systems in which the entire family of boundary states is prepared, such as the Chern insulator model.This indicates a great flexibility and generality in harnessing structured dissipation for probing topological boundary state physics through its interplay with coherent Hamiltonian dynamics.Here we present the exact solutions to the Lindblad master equation for the SSH model with loss on the B sublattice at arbitrary odd number of sites using third quantization.The exactly solvable dissipationless boundary mode together with the damping bulk modes give us full access to the quantum dynamics of the relaxation process, as well as a non-trivial non-equilibrium steady state (NESS).

Derivation of the damping matrix in third quantization
We start by writing the Hamiltonian of an SSH chain of spinless fermions with an odd number of sites: . The operator a † j,A (a j,B ) creates (annihilates) a fermion on the sublattice A(B) in the j-th unit cell and they satisfy fermionic anticommutation relations: {a j,α , a † j ′ ,α ′ } = δ j,j ′ δ α,α ′ .The last unit cell is broken and the total number of sites becomes odd, L = 2N − 1. Adding particle loss on the B sublattice with the jump operator Ll j,B = γ l 0,B a j,B yields a quadratic Lindbladian from the Lindblad master equation in Eq. ( 1), of which the rapidity spectrum can be obtained via third quantization [35][36][37].We briefly review the approach and adopt the same Majorana representation as in Ref. [8].One spinless fermion can be mapped to two Majorana fermions per site: This particular choice of mapping decouples two sectors in the damping matrix belonging to different Majorana fermions species c and d.Let us group them into a whole set under a vector notation w = (w 1 , w 2 , . . ., w 2L ) T = (c 1 , . . ., c L , d 1 , . . ., d L ) T .Majorana fermions are their own anti-particles as evinced by w † j = w j , and they obey anticommutation relations {w j , w k } = 2δ j,k such that each Majorana fermion squares to one.The density matrix ρ in the original Hilbert space of dimension 2 L ×2 L is now embedded in the set of Majorana operators with the new form , where α j ∈ {0, 1}.In the Majorana representation, the Hamiltonian and Lindblad dissipators take the form: H = j,k w j H j,k w k and where where γ α 's and η α 's stand for the sum and the imbalance of the most generic loss and gain dissipators on two sublattices for Ll j,α = γ l 0,α a j,α , Lg j,α = γ g 0,α a † j,α , ∀j. (A4) Over the 2 2L -dimensional Liouville space K, we now define adjoint fermionic operators that annihilate and create Majorana fermions: φ j |P α ⟩ = δ αj ,1 |w j P α ⟩, φ † j |P α ⟩ = δ αj ,0 |w j P α ⟩, which obey {φ j , φ † k } = δ j,k .Since the Fermi parity P F = (−1) j φ † j φj is conserved, [ L, P F ] = 0 and (P F ) 2 = 1, we can represent the Liouvillian in the even parity sector P F = +1 as In the chosen Majorana fermions representation in Eq. (A1), the full damping matrix is decoupled in different Majorana fermions species.The diagonal and off-diagonal blocks read with the unitary transformation U = diag{1, i, 1, i, . . ., i, 1} and γ = |γ l 0,B |/2.Therefore, the rapidity spectrum satisfies the relation where E m denotes the eigenvalues of the matrix H NH = H S + iΥ.

Boundary mode and full rapidity spectrum with spectral mirror symmetry
Under open boundary conditions (OBC), for the SSH chain with an odd number of sites L = 2N − 1, we first identify the zero-rapidity boundary mode β m=0 = 0.
It comes from the zero-energy boundary state of the SSH Hamitonian: where the normalization factor reads which leads to zero rapidity from Eq. (A9), Physically, this is consistent with the observation that the wave function of the zero-energy boundary mode is fully suppressed on the B sublattice for an SSH chain with an odd number of sites.The boundary mode survives under the dissipation on the frustrated sublattice in the infinitely long-time limit.For a chain with an odd number of sites, the bulk spectrum fulfils a spectral mirror symmetry: The mirror symmetry allows one to establish a rigorous relation between the spectrum under OBC and the spectrum under periodic boundary conditions (PBC) [28,39,42], with q = πm ′ N , m ′ = 1, 2, . . ., N − 1.We plot the real part of the rapidity spectrum as a function of the hopping amplitude and the dissipation strength in Fig. 5.
From the exact solutions, we can also obtain the dissipative gap that separates the boundary and bulk modes, together with the relaxation time τ for the system to evolve to the boundary mode in the long chain limit where (A15)

Exact eigenmodes of the damping matrix and dynamical observables
The complete set of eigenvectors of the damping matrix can be used to obtain dynamical observables.Apart from the exact boundary mode in Eq. (A10), next we show how to construct the analytical solutions to the bulk eigenmodes.
To make the solutions more generic, we start from a Bloch form of the non-Hermitian tight-binding Hamiltonian H NH into which the damping matrix can be transformed through Eq. (A8), where h x , h y ∈ R and h z ∈ C. The non-Hermitian term enters into the sublattice potential.From the eigenvalue equations, H NH u R (q) = E(q)u R (q) and H † NH u L (q) = E * (q)u L (q), one obtains the right and left eigenmodes of the Bloch Hamiltonian, u R,± (q) = 1 2h(h∓h z ) Due to the spectral mirror symmetry, we establish the relation for energy eigenvalues switching from PBC to OBC, with q = πm ′ N , m ′ = 1, 2, . . ., N −1.Therefore, the bulk eigenmodes under OBC can be constructed as a superposition of the PBC eigenmodes at opposite momenta, Here, the minus sign is determined by the boundary condition, ensuring that the overall wavefunction vanishes on the B site of the last broken unit cell.
If we denote the band index as m ∈ {0, (±, q)}, it is straightforward to check that the boundary and bulk modes in Eqs.(A10) and (A19) satisfy the biorthogonal relations [12,38,39]: ψ * L,m • ψR,l = δ m,l .The eigenmodes of the damping matrix can be obtained with an additional unitary transformation according to Eq. (A8), and they inherit the biorthogonal relations: For the SSH model with loss on the B sublattice, we obtain the bulk modes by identifying It is then convenient to resolve the time evolution of the observables with the complete set of exactly solvable eigenstates of the damping matrix.Applying the anticommutation relations of Majorana fermions {w j , w k } = 2δ j,k to the Lindblad master equation in Eq. ( 1), we arrive at the equation of motion for the covariance matrix For the trivial steady state, Now, we define the expectation value of a local observable with respect to the trivial NESS: Starting from an arbitrary initial configuration that is not trivial C(0) ̸ = 0, we can integrate the equation of motion and implement the diagonalized damping matrix in the exponential, Due to the biorthogonality of the basis, one reaches a compact form, At t = 0, we choose the system to be a completely filled chain: ⟨n j (0)⟩ = 1, ∀j.It corresponds to a covariance matrix, that selects µ ̸ = µ ′ .Going back to the physical space consisting of spinless fermions, we define the single-particle correlator Q jk (t) = Tr[a † j a k ρ(t)].The mapping to Majorana fermions in Eq.( A1) yields where the phase factor depends on whether or not the correlation resides on the same sublattice, Combined with Eqs.(A26) and (A27), the single-particle correlator takes an explicit form in terms of the exact solutions of the damping matrix, The particle number operator then reads where ⟨n j ⟩ s = 0 denotes the trivial steady state.
Appendix B: Non-equilibrium steady state for generic loss and gain Below, we show the analytic structure of the non-trivial NESS with generic loss on the B sublattice and gain on the A sublattice.It also leads to a closed form for the steady state current.
To find the NESS that satisfies we begin by rewriting each matrix into the blocks, It simplifies to In the presence of generic loss and gain, an analytical solution to the steady state covariance matrix follows While the matrix 1 L×L in D s gives back the state of an empty chain due to loss, c 0 denotes the contribution from the boundary mode immune to loss and coefficients p m,m ′ denote the overlaps between the bulk modes arising from gain, We consider the scenario when γ B ≫ γ A and γ A → 0 + .The non-trivial NESS can be connected to a sum of the empty state and the boundary mode.Since the former FIG. 6. Illustration of a Chern insulator on a honeycomb lattice consisting of two sublattices A and B. We consider a system with PBC in the x direction and OBC in the y direction, yielding a cylinder geometry.We wrap the cylinder such that its are of zigzag type on the left and bearded type on the right.The vectors ⃗ a1 and ⃗ a2 indicate the Bravais lattice vectors of the honeycomb lattice, while ⃗ δ1 indicate the vectors in the three possible bond directions.The sites are coupled by nearest-neighbor hopping t (depicted in yellow) and next-nearest-neighbor hopping t ′ (depicted in green).We further require that the next-nearest-neighbor hoppings along the x direction acquire an additional phase factor e iϕ to open a Chern/Haldane insulator gap.In the text, we adopt ϕ = π/2.
gives zero particle number occupation, regardless of the initial conditions, the dissipative system always relaxes to a unique NESS with a localization structure approximating that of the boundary mode.
For the measurement of the steady-state current J S = (i/L) j (⟨a † j a j+1 ⟩ s − ⟨a † j+1 a j ⟩ s ), we reach a closed form by substituting Eqs.(A28) and (B7) to the above definition: p m,m ′ ψ Lm (j)ψ * Lm ′ (j + 1).(B9) Appendix C: Liouvillian spectrum of 2D Chern insulator Here, we present more details on the derivation of the exactly solvable rapidity spectrum of a 2D dissipative Chern insulator [28,46] through a mapping to the generic SSH model.
As shown in Fig. 6, we place the honeycomb lattice on a cylinder and adopt PBC along the x direction with M unit cells and OBC along the y direction with N unit cells.The last unit cell along the y direction is broken, leading to L y = 2N − 1 sites.We choose the zigzag edge on the left boundary and the bearded edge on the right boundary.The nearest-neighbor and next-nearestneighbor hopping strengths are denoted by real parameters t and t ′ .The complex hoppings t ′ e iϕ with ϕ = π/2 are allowed only within unit cells along the y direction.Setting the honeycomb lattice constant to 1, different unit cells are connected by the vectors (C1) First, we perform a Fourier transform on the x components of the fermionic annihilation and creation operators (α = A, B), ⃗ r = (l, j), a α (l, j) = 1 √ M kx e ikxl a α (k x , j). (C2) For each k x = 2πm ′ /M, m ′ = −M/2, . . ., M/2 − 1, the real-space Hamiltonian shares the form The matrix elements can be obtained as follows: where Next, let us absorb the complex phases appearing in the hopping terms by the transformation a A (k x , j) → e −i(k1−k2)/2 a A (k x , j). (C5) It allows us to map H(k x ) to a generic SSH model with real hopping terms as before, with For simplicity, we have adopted the renormalized k x,0 , compared with the original value k x,0 on the honeycomb lattice.
Adding the B sublattice loss Ll j,B = γ l 0,B a j,B , ∀j, we apply the Majorana fermions representation in Eq. (A1).The generic SSH Hamiltonian reads where And the damping matrix takes the form The eigenmodes of Xc(d) can be exactly solved by the method we have developed earlier in Eqs.(A8)-(A20), which requires the identification of their NH Bloch Hamiltonians in Eq. (A16).From the mapping of the damping matrix in Eq. (A8), one arrives at where µ = ±1 for c, d, respectively.For the chiral edge mode, while for the bulk modes (k y = πm ′ /N with m ′ = 1, . . ., N −1), the spectral mirror symmetry β ± (k x , k y ) = β ± (k x , −k y ) leads to = γ/2 ± i t 2 1 + t 2 2 + 2t 1 t 2 cos(k y ) + (iγ/2 − µϵ A ) 2 .Furthermore, we check that the trivial steady state remains the same as the 1D SSH model (an empty state of spinless fermions along the y direction), The time evolution of the particle number distribution can be resolved as well, × ψ L,µm ′ (j)ψ * R,µm ′ (l)ψ R,µm (l)ψ * L,µm (j).Starting from a completely filled lattice driven by the B-sublattice loss, the steady-state particle number distribution becomes that of the chiral edge mode, ⟨n j ⟩ s = r j−1 − r j+1 1 − r 2N (j odd), ⟨n j ⟩ s = 0 (j even).(C16) Appendix D: Robustness of the boundary state against sublattice and bond disorders In the following, we discuss the robustness of the boundary state in Eq. ( 4) hosted by the 1D SSH model with odd sites, and generalize its localization structure in the presence of bond disorders.
Further-neighbor interactions and weak onsitepotential disorder on the A sublattice would eliminate the boundary mode from being a steady state of the Liouvillian.However, with loss engineering, there are other types of disorders that the boundary mode is robust against, including strong random disorders in the B-sublattice potentials and in the t 1 and t 2 hopping terms.Let us include them explicitly in the model with B-sublattice loss,

2 FIG. 1 .
FIG. 1. Minimal example: boundary state from single site loss.An open SSH chain with local loss only on the first site of the B sublattice (a), exhibiting a dynamical damping towards the boundary state at the same boundary (b) as indicated by the particle density ⟨nj(t)⟩ for a chain of length L = 2N − 1 = 9 with t1 = 0.5, t2 = 1.0 and γB = 0.5.

2 FIG. 2 .
FIG. 2. Solvable SSH chain with odd number of sites and uniform gain and loss.(a) Illustration of an open SSH chain with loss rate γB on the B sublattice and gain rate γA on the A sublattice.(b) Rapidity spectrum for γB = 0.5 (green dots) and γB = 3.0 (blue dots).The red dot at zero corresponds to the steady boundary mode.(c) Dynamics of the particle density ⟨nj(t)⟩ for a chain of length L = 2N −1 = 27 with t1 = 2.0, t2 = 1.0, γA = 0.0, γB = 0.5.(d) Steadystate current JS as a function of γA.The other parameters are as in (c).

FIG. 3 .
FIG. 3. Topology and scaling in SSH chains with even number of sites and uniform loss.(a) Dynamics of the density ⟨nj(t)⟩ for a topologically trivial chain of length L = 2N = 28 with t1 = 2.0, t2 = 1.0, γB = 0.5.(b) Dynamics of the density ⟨nj(t)⟩ for a topological chain of length L = 28 with t1 = 0.5, t2 = 1.0, γB = 0.5.(c) Real part of the rapidity spectrum as a function of t1 for a system with L = 46 sites and t2 = 1.0, γB = 1.0.(d) Dependence of the Liouvillian gap ∆ = 2 min{Re[βm]} on the chain length L for an even number of sites.Blue dots are the numerical result, while the red line indicates an exponential fit f (N ) = C exp(−αN ) with C = 1.202 and α = 0.695.
FIG. 4. Preparation of chiral Chern insulator edge states.(a) Chern insulator on a cylinder geometry with (b) the distribution of the steady chiral mode prepared by Bsublattice loss as a function of kx along the y direction.We adopt Ly = 2N − 1 = 31 for OBC and Lx = 200 for PBC.The chiral edge state is localized at the left (right) boundary when |r(kx)| = |2 cos(kx/2)| < 1 (> 1).The localization position changes at k * x = 2π/3, 4π/3 when |r(k * x )| = 1.(c),(d) Complex rapidity spectrum of the Chern insulator for N = 16 and t ′ /t = √ 3/2, t = γ = 1 demonstrating bulk modes (blue lines) and exact chiral mode (red line).The real part (c) illustrates the dissipative gap and the imaginary part (d) is equivalent to the full energy spectrum.
Appendix A: Liouvillian spectrum of 1D SSH model at odd lengths ) with X c = X d = −4iH 0 + 2M 1 .The upper-triangular form of L+ indicates that the eigenvalues of the Liouvillian coincide with those of the damping matrix X [8, 35].After proper diagonalization, we are able to express the Liouvillian in terms of rapidities β and normal master modes (NMMs) b ′ , b: L+ = − L m=1 β m (b ′ c,m b c,m + b ′ d,m b d,m ).(A7) Here, β m = β c,m = β d,m with the band index m and NMMs satisfy the anticommutation relations {b ′ c,m , b c,l } = {b ′ d,m , b d,l } = δ m,l .Furthermore, we can map the damping matrix in the Liouvillian to the SSH Hamiltonian with additional non-Hermitian imbalanced chemical potential on the two sublattices.For the (uniform) loss on the B sublattice only, we arrive at

FIG. 5 .
FIG. 5. (a),(b) Real part of the rapidity spectrum as a function of t1 for t2 = 1 and N = 46.The gray lines show the structure of the periodic system while the blue and red lines correspond to the bulk and boundary modes under OBC.We set different loss dissipation strengths on the B sublattice: (a) γ = 0.5, (b) γ = 2, (c) γ = 5.(d) Real part of the rapidity spectrum as a function of γ for t1 = 2, t2 = 1, and N = 46.The green line denotes the half of the bulk dissipative gap estimated according to Eq. (A15) for N ≫ 1.
γ = γ B = |γ l 0,B |/2.It is convenient to introduce a new set of Pauli matrices ⃗ τ acting on different subspaces of Majorana fermions species c and d, in order to decouple them, X 0 = E 0 ψ 0 : E 0 = 0, ψ 0 = N (r 1 , 0, r 2 , 0, . . ., 0, r N ) T ,r 1 = 1, r m = factor N −2 = N m=1 r 2 m .It still shares a localization structure fully suppressed on the B sublattice, thus immune to both loss and on-site disorder on the same sublattice.Indeed, after a proper rotation in the subspaces of different Majorana fermions species c and d as has been performed in Eqs.(C6)-(C10), random on-site disorder acting as a real chemical potential enters the rotated damping matrix in the diagonal terms belonging to B sites and becomes purely imaginary (µ = +1, −1 for c, d), of the generalized boundary mode ensures an infinite lifetime as a steady state, Xc,d ψ 0 = β 0 ψ 0 , β 0 = 0.(D4) Casting the Hamiltonian into a matrix form, it is straightforward to check that due to local destructive interference arising from disordered nearest-neighbor hopping terms, the targeted zero-energy boundary state takes a new form, ,1 a † j,A a j,B + t j,2 a † j+1,A a j,B + H.c.,V = j ω j,B a † j,B a j,B , Ll j,B = γ l 0,B a j,B , ∀j. (D1)