The bosonic skin effect: boundary condensation in asymmetric transport

We study the incoherent transport of bosonic particles through a one dimensional lattice with different left and right hopping rates, as modelled by the asymmetric simple inclusion process (ASIP). Specifically, we show that as the current passing through this system increases, a transition occurs, which is signified by the appearance of a characteristic zigzag pattern in the stationary density profile near the boundary. In this highly unusual transport phase, the local particle distribution alternates on every site between a thermal distribution and a Bose-condensed state with broken U(1)-symmetry. Furthermore, we show that the onset of this phase is closely related to the so-called non-Hermitian skin effect and coincides with an exceptional point in the spectrum of density fluctuations. Therefore, this effect establishes a direct connection between quantum transport, non-equilibrium condensation phenomena and non-Hermitian topology, which can be probed in cold-atom experiments or in systems with long-lived photonic, polaritonic and plasmonic excitations.


Introduction
Transport phenomena are of relevance for almost all areas of physics and technology with transport of electric currents and heat conduction in solids being two prototypical examples.While electric currents are carried by electrons, i.e., massive fermionic particles, heat transfer can be understood as the emission and reabsorption of quantized lattice vibrations, i.e, non-conserved bosonic excitations.However, despite relying on very different microscopic mechanisms, both transport scenarios share many similarities.For example, depending on the mean free path, transport can either be ballistic or diffusive, where in the latter case Ohm's law and Fourier's law describe a similar linear relation between the current and the applied voltage or temperature gradient.Therefore, a general question of interest is under which conditions 'anomalous transport' with a qualitatively very different phenomenology can be observed.
In this paper, we consider the setup shown in Fig. 1 (a) as an elementary model to study dissipative transport of bosons.Here, bosons injected from a hot reservoir on the right can incoherently hop between neighboring sites of a one dimensional lattice, before being dumped into a second reservoir on the other end.This process has two key features: First, in the  Bosons injected from a thermal particle reservoir with mean occupation number nr on the right can incoherently hop along the lattice with asymmetric rates Γ l and Γ r , before being emitted into a second reservoir with occupation number nl on the left.A directional hopping can be imposed, for example, by applying a potential gradient with an energy offset U between neighboring sites.(b) Under stationary conditions, this hopping asymmetry combined with the bosonic particle statistics results in the bosonic skin effect, i.e., the formation of a finite boundary region with a staggered density profile.The two insets show sketches of the Wigner distribution for individual lattice sites, indicating that within this boundary region, the odd sites are in a condensed state with broken U(1) symmetry, while all other lattice sites exhibit a thermal distribution.See text for more details.
presence of a bias, the hopping rates to the left and right, Γ l and Γ r , are in general different, in which case the transport is asymmetric, i.e., directional.Second, the hopping rates toward sites that are already occupied are enhanced by the bosonic particle statistics.Therefore, this process can be seen as the bosonic counterpart of the celebrated asymmetric simple exclusion process (ASEP) [1][2][3][4]-a common model for directed transport of fermions or classical hardcore particles-and one speaks of an asymmetric simple inclusion process (ASIP) instead. 1ompared to fermions as the carriers for electric currents, the dissipative transport of bosonic particles has attracted considerably less attention so far.This can be attributed to a lack of conventional solid-state systems where this physics could be observed.However, this situation has changed recently and a variety of experimental platforms have now become available where non-equilibrium processes with bosonic particles can be probed.This includes, for example, cold atoms in optical lattice potentials, where different techniques to study transport have already been demonstrated [9][10][11][12][13].Furthermore, it has been shown in various experiments that long-lived photonic [14][15][16], polaritonic [17][18][19][20][21][22][23] or plasmonic [24] excitations can behave as massive bosonic particles and equilibrate with the surrounding material, before they eventually decay.The ongoing experimental advances in these platforms naturally raise the question of how transport in such settings is affected by the bosonic particle statistics of the carriers.
In the following analysis, we investigate the properties of the ASIP in a thermal transport scenario, where we focus primarily on the stationary current and the density profile along the lattice.In the absence of asymmetry, we recover the usual diffusive transport in this model as well, characterized by a linear population gradient and a Fourier law for the current.However, as soon as a finite degree of asymmetry is introduced, the transport becomes ballistic and particles accumulate in a finite boundary region near the drain.Moreover, as the total current through the system increases, we observe a transition from a smooth pile-up to a zigzag structure, as depicted in Fig. 1 (b), with odd (even) sites being highly (weakly) populated.This phase represents a rather unusual non-equilibrium configuration, where the particle distribution alternates on every lattice site between a thermal distribution and that of a coherent state with broken U(1)-symmetry.This emergence of coherences in a purely dissipative and thermal transport scenario is very surprising and related to non-equilibrium condensation phenomena [14,19,[25][26][27][28] that have no counterpart in fermionic transport.Therefore, we identify this boundary condensation as a unique feature of the ASIP model and call it the bosonic skin effect.
The observed accumulation of particles in a dissipative transport scenario is indeed very reminiscent of the so-called non-Hermitian skin effect (NHSE) [29][30][31][32][33][34][35][36][37][38][39][40][41][42][43].This effect refers to the boundary localization of the eigenfunctions of certain non-Hermitian lattice Hamiltonians and is thus frequently discussed in connection with their topological classification [32,35,42,44,45].However, such non-Hermitian models do not conserve the norm of the wavefunction nor the particle number.Therefore, beyond their mathematical interest, the relevance of the NHSE and other spectral features of non-Hermitian systems for actual quantum transport processes is not immediately clear and still a subject of ongoing investigations [31,[36][37][38][39][40][41]43].Here, by mapping the dynamics of density fluctuations in our system onto the paradigmatic Hatano-Nelson model (HNM) [29,35,42,45], we establish a direct correspondence between the eigenvalue structure of this non-Hermitian Hamiltonian and the stationary states of the ASIP transport problem.This correspondence relies on a subtle difference between Dirichlet and Neumann boundary conditions for the HNM and provides important additional insights into the nature of the predicted boundary transition.In particular, we find that the onset of condensation coincides with the appearance of a higher-order exceptional point in the HNM and occurs without a closing of the dissipative gap.This distinguishes the bosonic skin effect from other dissipative quantum phase transitions [46,47] and, in summary, reveals an unexpectedly rich interplay between transport, non-equilibrium condensation effects and non-Hermitian physics.
The remainder of the paper is structured as follows.In Sec. 2, we introduce the ASIP model and the main transport equations that we use to describe it.In Sec. 3, we present the bosonic skin effect and discuss the onset of the zigzag phase within mean-field theory, before investigating the full particle distribution and condensation effects in Sec. 4. Finally, in Sec. 5, we discuss the connection between the ASIP and the HNM, before summarizing our main findings in Sec. 6.Additional details about the analytic derivations and numerical methods are presented in the appendices.

Model
We consider the transport of bosons in a 1D lattice, as depicted in Fig. 1 (a).Here, the bosons are injected from a thermal reservoir on the right and propagate along a chain of L lattice sites through incoherent hopping processes, before being emitted into a second reservoir on the left.In the following we are primarily interested in asymmetric transport, Γ l > Γ r , where Γ l and Γ r denote the hopping rates to the left and to the right, respectively.

The ASIP master equation
We model the dynamics of this system by the Lindblad master equation where ρ is the system density operator.Here, the first term describes the incoherent hopping of bosons along the lattice.This process is described by the Liouville superoperator [40,[48][49][50]] where âp (â † p ) are the bosonic annihilation (creation) operators for lattice site p and we have introduced the short notation In Eq. ( 2), the jump operator â † p+1 âp (â † p−1 âp ) destroys a boson at site p and creates a boson at site p+1 (p−1) instead.This process conserves the total particle number and it is thus different from particle loss or gain.As a direct consequence of this particle number conservation, each jump operator is quadratic in â and â † , and therefore the hopping process is nonlinear.
The second and the third term in Eq. ( 1) represent the coupling to the thermal particle reservoirs to the left and to the right, which we model by Here κ l and κ r denote the coupling rates to the two reservoirs and nl and nr are the corresponding thermal occupation numbers.Note that while we will only consider thermal baths in this work, other pumping mechanisms, such as incoherent gain, would result in a behavior that is qualitatively very similar to what is discussed below.

Asymmetric hopping
Before we proceed, let us briefly comment on the physical motivation behind this asymmetric transport model.A very generic scenario is depicted in Fig. 1 (a), where bosons are confined to a lattice with an energy gradient, for example, an optical lattice for cold atoms [10,12], a nanophotonic lattice for exciton polaritons or plasmons [21,24], etc.In this case, due to a large energy offset U > 0 between neighboring sites, coherent tunneling is suppressed, but in the presence of a phononic bath, the bosons may still transition between neighboring sites by emitting or absorbing vibrational excitations.Such a process can be modelled by a phononassisted tunneling term of the form where the bosonic operators bp represent local bath excitations.Roughly speaking, for a particle to jump to the left, it must lose the energy ∼ U by emitting it into the environment.
Conversely, to jump to the right, it must absorb the same amount of energy.Therefore, a bath at low temperature, where emission processes are more likely than absorption, favors hopping to the left.More precisely, under the assumption that the bath is sufficiently Markovian, its dynamics can be eliminated to derive an equation of motion for the reduced system density operator ρ only.While some details may depend on the specific implementation (see Appendix A for a more detailed derivation), this master equation will be, quite generically, of the form given in Eq. ( 1), with hopping rates satisfying where T phon is the temperature of the phononic bath, which determines the asymmetry in this setting.Apart from such naturally occurring dissipative hopping mechanisms, there are also many systems where this asymmetric hopping processes can be engineered.For example, in optical lattices, directed dissipative hopping can be implemented via Raman processes [40,[51][52][53], which involve atomic or cavity decay as a source of dissipation and directionality.Ideas for realizing number-conserving dissipation processes for photons have also been discussed for optomechanical systems [54,55] and circuit QED [56], and can be readily adapted for the implementation of directed hopping processes as well.In the following we do not consider any of these possible implementations specifically, but rather address the general properties of the transport model given in Eq. (1).

Transport
In this work we focus primarily on the stationary transport of particles between two thermal reservoirs.In the absence of asymmetry, transport would be solely driven by the temperature gradient between the reservoirs, i.e., by the difference between nr and nl .For asymmetric rates, Γ l ̸ = Γ r , a directed particle flow develops even without any external temperature bias.
To characterize transport in different parameter regimes, we consider the average stationary current J as well as the stationary density profile n p = 〈n p 〉 = 〈â † p âp 〉 along the chain.Throughout this paper we adopt the convention that symbols with hats represent quantum operators, while symbols without hats denote their averages.Starting from the master equation in Eq. (1), the mean occupation number n p of any of the sites changes in time as This equation has the form of a conservation law, where, for any p ∈ [1, N − 1], is the average particle current between sites p and p + 1.Note that we have adopted the convention that a positive J p,p+1 implies a current flowing from right to left, i.e., from site p + 1 into site p.From Eq. ( 7) we already see that the current depends non-linearly on the density, due to bosonic bunching: indeed, the probability for a particle on site p + 1 to jump to site p is enhanced by a factor 1 + np , which depends on the population of the target site.On the boundaries, the currents represent the flow of particles into the left bath and from the right bath, respectively.In the steady state, the particle current is conserved along the chain and we obtain Note, however, that this uniformity of the current does not imply a uniform profile for the density n p .

Mean-field dynamics
Although Eq. ( 1) contains only dissipative terms and no additional coherent interactions between the bosons, these incoherent processes are nonlinear and therefore do not permit a closed set of equations for the mean occupation numbers.In addition, since the number of possible bosonic configurations scales exponentially with the number of lattice sites L, brute-force numerical solutions of the master equation are also inaccessible for the parameter regimes of interest.Therefore, to proceed we resort to a mean-field decoupling of the equations of motions by factorizing expectation values as 〈n p np+1 〉 ≈ 〈n p 〉〈n p+1 〉.Under this approximation, the average current reads The system is then described by a set of L nonlinear differential equations, which can be solved efficiently numerically and also permit exact analytical solutions in the steady state.
To benchmark the validity of the mean-field approximation, we compare these predictions with exact Monte-Carlo simulations for small systems sizes and low occupation numbers n p ≲ 1 and with phase-space simulations based on the Truncated Wigner Approximation (TWA) [57] for larger occupation numbers.Within their respective regimes of validity, we find almost perfect agreement between the numerical results and the stationary distributions obtained from mean-field theory.Further details about these numerical methods and some of the benchmarks can be found in Appendix B.

Hydrodynamic limit
Additional insights about the transport dynamics in our system can be obtained by considering the continuum (or hydrodynamic) limit.To do so, we rewrite the mean-field equations of motion as where Γ A = Γ l − Γ r and Γ S = (Γ l + Γ r )/2.Then, under the assumption that the n p vary slowly between neighboring sites, we can replace them by a continuous field n(x, t), where x is the dimensionless position along the lattice.Away from the edges, this field obeys the partial differential equation This is, in essence, the well-known Burgers' equation [58][59][60], a simplified version of Navier-Stokes equation in hydrodynamics.The parameters Γ A and Γ S can thus be interpretated as non-linear advection and diffusion rates, respectively.With the left side of the lattice being initially empty, the possible solutions of Eq. ( 12) include propagating shock fronts of the form [60] where nsw is the height, c sw = Γ A (n sw + 1) the speed and w sw = 2Γ S /(n sw Γ A ) the width of the wavefront.These solutions clearly illustrate how the bosonic enhancement factor affects transport.First, the velocity of the density wave scales with the typical density nsw .Second, the bosons in the high density region propagate faster than the bosons at the front, which leads to a compression of the wave and w sw going to 0 for very large nsw .
While the Burgers' equation provides valuable intuition about the transport dynamics in our system, it is based on a continuum approximation.As such, it is only expected to hold in a 'laminar' regime, i.e., when the effective Reynolds number associated with a typical occupation number nsw , is small. 2In the opposite limit, the characteristic length scale, w sw ∼ O(1), becomes of the order of the lattice spacing and new features can arise from the discreteness of the lattice and the presence of boundaries.

Relation to the ASEP
By replacing the bosonic operators in Eq. ( 1) by operators âp that obey fermionic anticommutation relations, i.e., {â p , â † p } = 1, we obtain the master equation describing the ASEP.In this case, the site occupation numbers n p obey the same equation as in Eq. ( 6), but with a fermionic current Here, rather than being enhanced, the hopping to neighboring sites is prohibited by the Pauli exclusion principle, if the site is already occupied.The properties of the ASEP have been extensively studied in the literature [1][2][3][4].This includes, most notably, the scaling of current fluctuations [61,62] in infinite lattices, which falls into the Kardar-Parisi-Zhang (KPZ) universality class [4,63,64].The ASEP is thus closely connected to surface growth and related non-equilibrium phenomena.It is therefore interesting to understand how the change from an exclusion to an inclusion process affects these properties.These aspects, however, will be discussed in more details elsewhere [65].Instead, here we focus on novel effects that are unique to the ASIP and reveal themselves already at the mean-field level.

The bosonic skin effect
In the following section, we explore in more details the stationary state of the transport master equation in Eq. ( 1), which we describe in terms of the mean occupation numbers n p and the current J.

Transport regimes
In a first step, we show in Fig. 2 examples of the stationary density profile n p for a lattice of L = 15 sites, together with the scaling of the current J as a function of L. From these plots we identify three qualitatively different transport regimes.

Diffusive transport
In Fig. 2 (a) we first consider the symmetric case Γ l = Γ r , where the stationary density profile along the chain is simply a linear interpolation between nl and nr .This is also expected from 2 In the literature, the Burgers' equation is often defined as x n, with Re the Reynolds number, and n ∼ O(1).In our case, Eq. ( 12) includes an additional linear term ∼ Γ A ∂ x n and the field assumes values between 0 and n ∞ .However, we can rescale the field as n → n/n ∞ and time as t → 2Γ A n ∞ t and remove the linear term by a Galilean transformation.After these transformations we recover the standard form of the Burgers' equation with Re given by Eq. ( 14).
For all plots nr = 10 and two different values of nl = 0 (blue lines) and nl = 20 (yellow lines) have been considered.The insets show the current J versus the lattice size L, in log-log scale and for three different values of nl = 0, 5, 9.For Γ A = 0 we recover a linear population gradient and the Fourier law for the current, as expected for diffusive transport.For any Γ A > 0 and large L, the current becomes independent of both L and nl , indicating ballistic transport.In this regime, we observe the formation of a finite boundary region of size ξ, as indicated by the shaded area.As the asymmetry increases, the width ξ shrinks and vanishes for Γ A /Γ l ≃ 0.17.Beyond this point, a finite boundary region, but with an oscillating density profile, reappears.For all plots, we have set κ r = κ l = Γ l .
Burgers' equation in the continuum limit, Eq. ( 12), which for symmetric hopping describes pure diffusion.In this regime, the current obeys the Fourier law and decreases with system size, i.e., Interestingly, this diffusive transport is independent of the particle statistics and it is the same for bosons, fermions and noninteracting classical particles.

'Laminar' asymmetric transport
For a sufficiently large lattice, L ≫ 1, the diffusive transport turns into directional transport for any finite hopping imbalance, Γ A ̸ = 0.In this case, the stationary density profile is flat and assumes a constant value of n p ≃ n ∞ across most parts of the lattice.The exception is a region of size ξ close to the left reservoir, where the density gradually adjusts to a boundary value, which depends on the occupation number of the left bath, nl .Most importantly, for a lattice size L ≫ ξ, the stationary current J > 0 is completely independent of both nl and the length of the chain [see the inset of Fig. 2 (b)].This is true even though Γ r is still finite.This is in contrast to ballistic transport in coherent systems [66][67][68][69], where the stationary current depends on the properties of both reservoirs.
While the quantitative details in this regime are already affected by the bosonicallyenhanced hopping rates, the population profile is still qualitatively similar to what one would obtain for asymmetric hopping of independent classical particles.Moreover, since the effective Reynolds number introduced in Eq. ( 14) is still small, this behavior is well described by the continuous Burgers' equation in Eq. ( 12) and we can draw a close analogy with the regime of laminar flow in fluid dynamics.

'Turbulent' asymmetric transport
When either the asymmetry or the right bath occupation nr are further increased, the size of the boundary region, ξ, decreases and reaches ξ = 0 at a critical value Γ c A ≡ Γ c A (n r ).At this specific point, the density profile is completely flat, with the exception of site p = 1, which is coupled to the left reservoir.As shown in the inset of Fig. 2 (c), since the relevant length scale vanishes, the current at this critical value is independent of the system size and adopts the value Remarkably and somewhat unexpectedly, this situation occurs already for finite Γ r , i.e., under conditions where particle flow in both directions is still possible.
As the directed particle flow is further increased, a boundary region of finite size ξ reappears.In this regime, however, the occupation numbers vary strongly between neighboring sites and we observe a zigzag configuration with a decaying envelop.Counter-intuitively, as we keep increasing Γ A , we find that the extent of this zigzag configuration increases in the direction opposite to the propagation.The transport in this regime is ballistic as well, i.e., for sufficiently large L the current J ≈ κ r nr > 0 , (18) is independent of both nl and the system size.However, in contrast to the smooth pile-up observed above, this rapidly oscillating density profile is no longer captured by the Burgers' equation.This behavior is found for high effective Reynolds numbers and in analogy with turbulent flow in fluid dynamics, we observe a build-up of excitations at small length scales.
In our discrete lattice setting, this leads to a breakdown of the continuum approximation. 3his staggered accumulation of particles in alternating lattice sites, rather than being distributed smoothly across the lattice, does not appear in analogous models for directed transport of fermions or classical particles.Since it arises from a purely dissipative process, this pattern must also be distinguished from the formation of standing waves in coherent channels [69].It is thus a unique consequence of bosonic bunching.

Stationary density profile
Let us now proceed with a more in-depth analysis of the stationary density profile.In the steady state, the current J is uniform across the lattice and we can use Eq.(10) to relate the occupation numbers between neighboring sites by for all p.For a large enough lattice, L ≫ 1, and p large, the occupation numbers near the right reservoir approach a constant value n p ∼ n p+1 = n ∞ , which is determined by the fixed point of this equation.This leads to the following general relation, between the stationary current and the asymptotic particle density.The boundary condition for the reservoir on the right also gives us J = κ r (n r − n ∞ ), which allows us to compute explicitly the asymptotic density, and from it the stationary current J.Note that both quantities are smooth functions of all the system parameters and don't exhibit any sharp features.For large nr we obtain n ∞ ∼ nr and a current J ≈ κ r nr , which is limited by the influx of particles from the right reservoir.
The left boundary condition imposes J = κ r (n 1 − nr ), meaning that n p ̸ = n ∞ for small site numbers p.In Appendix C we show in more details how the relation in Eq. ( 19) can be used to determine the full density profile n p in the limit L → ∞, which for Γ l > Γ r can be written in the form Here, µ is a constant that depends on the properties of the left reservoir, but its precise dependence is not important for the following discussion.In Eq. ( 22) we have also introduced the parameter which is the bosonically-enhanced speed of propagation.Indeed, c is closely related to the speed of the shockwaves discussed in connection with the Burgers' equation ( 12), but determined by the self-adjusted, stationary density n ∞ .By looking at the first term on the right side of Eq. ( 22), we see an exponential decay of the excess population, which can be re-expressed as Therefore, in both regimes, we can define the characteristic decay length As we increase Γ A or nr , ξ decreases, and goes to zero for c = Γ S .This allows us to identify the critical value of the hopping imbalance, at which point ξ = 0 and the system changes between the smooth and the zigzag boundary configuration observed above.Beyond this point, we acquire an extra phase π, which explains the alternating occupation numbers for values of The full dependence of ξ on Γ A and nr is plotted in Fig. 3, which clearly shows a sharp drop to zero along the transition line Γ A = Γ c A .For all points in this plot a value of κ r = Γ l has been assumed, and the results are independent of both nl and κ l .

Nonlinear transport and the Fibonacci sequence
While the first term in Eq. ( 22) defines the characteristic size of the boundary region, it is important to keep in mind that the full density profile is not described by a simple exponential decay.This deviation, represented by the second term in Eq. (22), is due to the nonlinear nature of transport arising from the bosonic particle statistics.To obtain additional insights about this profile, we show in Appendix C that the stationary occupation numbers can be rewritten in the form n p = a y p−1 where the new quantities y p obey the recursion relation This reformulation shows that rather than being described by an exponential decay, the mathematical structure of n p is given by the ratio of successive coefficients of a generalized Fibonacci sequence defined by Eq. ( 28), also known as a Lucas sequence.For example, in the special case of Γ r = 0 and a current J = Γ l , we obtain a = b = 1 and d = 0 and the populations n p then oscillate toward n ∞ = (1 + 5)/ 2 in the same way that the ratio of successive coefficients of the Fibonacci sequence oscillates towards the golden ratio.This observation is not just a purely mathematical curiosity, but a very generic feature of nonlinear transport.Indeed, any transport model with a next-neighbor nonlinear recursion relation of the type αn p n p+1 + β n p + γn p+1 = δ will lead to a density profile of the form given in Eq. (27).By contrast, recursion relations of the type β n p + γn p+1 = δ, as encountered in linear transport models, give rise to a simple exponential population profile.

Boundary condensation
The strong bunching of the bosons in certain lattice sites, as observed for Γ A > Γ c A , is somewhat similar to the formation of a Bose-Einstein condensate, where at low temperatures bosons tend to accumulate in a single momentum mode.However, in our setting this effect is observed under conditions where a large thermal current passes through the system, and locally one would expect a thermal distribution of particles instead.To resolve these two conflicting physical pictures, we must go beyond mean-field theory and take a closer look at the full particle number distributions and the coherence properties of our system.

Density fluctuations
To study effects beyond mean-field theory, we use numerical simulations based on the TWA.Within the TWA, the Wigner distribution is sampled by complex phase-space variables α p that follow stochastic trajectories.Symmetrically-ordered expectation values of the form 〈â †n p âm q 〉 sym are then approximated by the corresponding stochastic averages 〈α * n p α m q 〉.We refer to Appendix B for more details about this method.In Fig. 4 (a) we use the TWA to evaluate the equal-time two-particle correlation function g (2)  p (0) = for each of the lattice sites, and once the system has reached a steady state.The phase space plots below this curve show the corresponding distributions of the α p , as obtained from the individual trajectories in the numerical simulation.These sample the Wigner distribution of that site.We see that, near the right reservoir, the value of this correlation function is g (2) (0) ≃ 2, as expected for a thermal state [71].The corresponding Wigner distributions are very close to a Gaussian distribution centered around α = 0. Near the left boundary, however, g (2) (0) decreases for all odd sites and approaches a value of g (2) (0) ≈ 1, which indicates a coherent state.In this case the corresponding phase-space distribution has the shape of a symmetric ring with a maximum at a finite value of |α p | ≈ n p .In contrast, on all even sites the distribution remains Gaussian-like and centered around α p = 0, although values of g (2) (0) > 2 indicate small deviations from an exact thermal distribution.This overall behavior is further confirmed by the probability distributions P(|α p | 2 ) plotted in Fig. 4 (b).

U(1) symmetry breaking and phase coherence
The ASIP describes a purely incoherent hopping process.This means that the full master equation given in Eq. ( 1) is diagonal in the number basis and it is invariant under the local U(1) symmetry transformations âp → âp e iφ p .This symmetry is also clearly visible in the phasespace plots in Fig. 4 (a), which are fully symmetric under rotation.However, these results only describe an ensemble average in the steady-state, while within a given experimental realization, or for finite time, the U(1) symmetry can still be spontaneously broken.
To analyze potential symmetry-breaking effects in our system, we are interested in how long information about the phase in a given site is preserved.This is quantified by the coherence function g (1)  p (τ) = lim In Fig. 4 (c), we show the evolution of g (1)  p (τ) as a function of the delay time τ.We see that for odd sites near the left reservoir, this correlation function decays over a timescale  (1) (τ) (30) for odd (left) and even (right) sites near the boundary (all the even sites have extremely similar behavior, here we show only p = 2 for simplicity).For all plots we have set κ = Γ l , Γ r = 0, nr = 30 and nl = 0.For the plots in (c), we have used a reference time of t = 10Γ −1 l , which is sufficient to reach the steady state.
which is multiple times longer than the typical relaxation timescales in this system.In contrast, for even sites, no such extended phase correlation can be observed and the coherence vanishes on timescales much faster than Γ −1 l .Note that we do not observe any significant cross-correlations between any of the lattices sites, either.
To understand this emergence of coherence in more details, we consider the totally asymmetric case, Γ r = 0, and also assume nl = 0 for simplicity.Under these assumptions, the phase-space variable α 1 of the first lattice site obeys the stochastic equation (see Appendix B.2) where dW is a Wiener process.For the current discussion we have also adopted the convention n 2 ≡ |α 2 | 2 − 1/2 to be consistent with symmetrized expectation values, 〈â † p âp 〉 sym = n p + 1/2 = 〈|α p | 2 〉, even on the level of a single trajectory.Close to the steady state, the occupation number of the second site can be expressed in terms of the recursion relation in Eq. ( 19) and approximated by After reinserting this results into Eq.( 31), we obtain a closed diffusion equation for the variable α 1 , which is of the form From the deterministic part of this equation, we see that the nonlinear hopping process acts like an effective saturable gain.For sufficiently large J, this leads to a growth of the initial amplitude, which then saturates at a value n 1 = |α 1 | 2 − 1/2 ∼ J/κ l , consistent with the steady state result obtained from the mean-field analysis in this regime.
For J/κ l ≫ 1 and once the amplitude α 1 has been amplified to a large value, it can be approximately written as α 1 (t) ≃ n 1 e iφ 1 (t) , with a fixed n 1 and a phase φ 1 (t) that obeys This phase diffusion equation predicts a decay of the ensemble-averaged amplitude according to with a coherence time of τ coh = 4J/κ 2 l ≃ 4n r /κ l .In Fig. 5 we consider a scenario in which the lattice is initialized in a symmetry-broken state, i.e., a coherent state with a very small but finite amplitude on each site.The displacement direction can change from site to site, but remains the same from one trajectory to the next.For this initial configuration, the plots in Fig. 5 (a) show the successive evolution of the Wigner distribution of site p = 1.We clearly see that the small initial displacement is quickly amplified to its steady-state value, after which the phase diffuses on a much longer timescale.Eventually, we recover the ring-shaped profile shown in Fig. 4 (a).In Fig. 5 (b) we plot the evolution of |〈α 1 〉| for different values of nr .The long-time decay of this quantity agrees very well with the analytic prediction in Eq. (35).Note also that, initially, each trajectory breaks the symmetry in the same direction, set by the initial perturbation.Hence, for intermediate times, the symmetry breaking is present even at the level of the density matrix.

Summary
In summary, the results presented in this section show that the zigzag structure observed at the mean-field level is consistent with the picture of an alternating lattice of condensed and thermal-like bosonic states.Consistently with other non-equilibrium condensation phenomena or closely related lasing effects [14,19,[25][26][27][28], the Bose-condensed sites in our system are characterized by a spontaneously broken U(1)-symmetry with a phase coherence time that is long compared to the typical relaxation timescales in this system.The most surprising finding in our setting is that this effect occurs only in every other site near the boundary, while neighboring sites and other parts of the lattice remain close to a thermal state.This configuration is specific to the current transport scenario, where the stationary populations are determined by the nonlinear recursion relations discussed in Sec.3.3, rather than by energetic considerations or an external gain mechanism.
Note that condensation effects have also been discussed for zero-range [5] and other attractive transport processes [5,6], where even on a periodic lattice all particles eventually Figure 5: (a) Evolution of the Wigner distribution of site p = 1, when the system is initially prepared in a symmetry-broken state with |〈α p 〉| = 3 and a random phase.The three plots show the resulting phase-space distributions obtained in a TWA simulation for times Γ l t = 0, 2, 40 and for nr = 80.On a short timescale, the initial displacement is amplified, while phase diffusion is observed over much longer times.(b) Logarithmic plot of the ensemble-averaged amplitude of the first site for the same initial conditions, but assuming different thermal occupation numbers of the right reservoir.After a short amplification, we observe an exponential decay of the average amplitude due to phase diffusion.The dashed lines represent the analytic prediction for this decay, as given in Eq. (35).For all plots, Γ l = κ r , Γ r = 0, L = 10 and nl = 0. accumulate in a single site.This is not the case for the ASIP considered here, where for periodic boundary conditions the system would simply evolve into an infinite-temperature state with all particle configurations being equally likely.Therefore, the presence of a boundary is essential to observe this type of condensation, which would not follow from an analysis of bulk properties only.

Asymmetric bosonic transport and the Hatano-Nelson model
As already pointed out in the introduction, the accumulation of particles near one end of the lattice in a dissipative transport model shares many similarities with the NHSE.This effect refers to the fact that the eigenfunctions of certain non-Hermitian lattice Hamiltonians, which are extended over the whole lattice for periodic boundary conditions, become exponentially localized when open boundary conditions are introduced.A prominent example where this effect occurs is the HNM [29], which, indeed, has originally been introduced to describe directional transport of bosons.However, in contrast to the full ASIP master equation considered here, the HNM is formulated in terms of a tight-binding Hamiltonian with asymmetric tunnelling amplitudes.Such a Hamiltonian is necessarily non-Hermitian, meaning that it does not preserve probabilities, particle numbers or operator commutation relations.Therefore, despite a considerable interest in the spectral properties of the HNM and related non-Hermitian Hamiltonians, their relevance for actual quantum transport problems often remains unclear.
While consistent embeddings of the HN Hamiltonian into a proper master equation have been discussed [31, 37-39, 41, 43], these works have considered linear jump operators, which describe particles being exchanged with the environment.In this case, the evolution does not obey a conservation relation with a well-defined current, and the connection to the original dissipative hopping problem is lost.In the following we show instead how an explicit connection between the HNM and the ASIP transport problem can be established at the level of density fluctuations.This discussion complements the single-particle analysis of Ref. [40], and provides a new interpretation of the HNM at the level of a many-body transport problem.It also reveals a surprising relation between the dynamics of fluctuations and the stationary state of this system.

The non-Hermitian skin effect
The HNM is the simplest model to study boundary localization in non-Hermitian systems.It is described by the lattice Hamiltonian where the âp represent non-interacting bosons or fermions, whose dynamics is then fully described by the tunneling matrix Here, x = 1 for periodic boundary conditions and x = 0 for an open chain.For J r = J l we recover the usual tight-binding Hamiltonian with real-valued single particle eigenenergies, E k = 2J r sin(k).The corresponding momentum eigenstates are extended over the whole lattice, both for open and periodic boundary conditions.For J r ̸ = J l , by contrast, the tunneling to the left and to the right is no longer the same, and Ĥ † HN ̸ = ĤHN .Still, when assuming periodic boundary conditions, the eigenfunctions of h HN remain plane waves, ψ k (p) ∼ e −ikp , where k ∈ [−π, π), but with a complex spectrum which describes an ellipse in the complex plane.In contrast, for open boundary conditions, all eigenmodes are exponentially localized near one end of the chain [72], and are no longer orthogonal to each other.The corresponding spectrum is given by Thus, the spectrum changes from a closed loop to a line in the complex plane (see Fig. 6 and the discussion below).This transition from an extended to a localized set of wavefunctions when changing from periodic to open boundary conditions occurs in many other related lattice models, and has been dubbed NHSE [29][30][31][32][33][34][35][36][37][38][39][40][41][42][43].
From Eq. ( 40) we see that when J r and J l have the same sign, i.e., J r J l > 0, the singleparticle energies are real and therefore describe solutions that oscillate in time.However, when J r J l < 0, the spectrum is purely imaginary, i.e., it describes decaying or amplified solutions.These two regimes are separated by a so-called exceptional point (EP) at J r = 0, where the Hamiltonian of Eq. ( 36) becomes defective and cannot be diagonalized anymore.Instead, the tunneling matrix adopts a Jordan normal form which has only a single eigenmode with energy E EP = 0 and a wavefunction ψ EP (p) = δ p1 , which is fully localized on the first site.The other basis elements are so-called generalized eigenvectors, i.e., they are transformed into ψ EP through the action of h HN .The NHSE and the presence of exceptional points have recently attracted a lot of attention, in particular in connection with the classification of topological properties of non-Hermitian lattice systems [32,35,42,44,45].

Linearized boson transport
Let us now return to our mean-field model in Eq. ( 6) and consider a situation where at some initial time t = 0 the whole lattice is prepared in a state with a flat density distribution n p (0) = n ∞ .For the successive evolution we make the Ansatz and assume that the fluctuations ε p remain small compared to n ∞ .This is justified for short times and, more generally, under the condition nr + nl ≈ 2n ∞ .We can then linearize the mean-field equations of motion and obtain with m = (n l + nr − 2n ∞ ) and δ i j the Kronecker delta, and we have set κ l = κ r = κ for simplicity.
To connect this result to the HNM discussed above, we introduce the vectors with a non-Hermitian Hamiltonian Therefore, ignoring the coupling to the reservoirs for now, i.e. κ → 0, we see that the density fluctuations ε p obey an effective Schrödinger equation with a non-Hermitian Hamiltonian h, which, by identifying J r ↔ c − Γ S and J l ↔ c + Γ S , is very similar but not identical to h HN .In particular, the diagonal elements of h are shifted by a constant imaginary part −2iΓ S and, compared to h HN , there is an additional term c + Γ S in the first entry of h.The first change merely shifts all the eigenenergies in the complex plane towards negative imaginary values, enforcing stable dynamics.The second change, as we will see, arises from the boundary conditions.

Neumann boundary conditions and the steady state
To understand the differences between h and h HN , we emphasize that the dynamics of the fluctuations ε p in Eq. ( 43) can still be written as a continuity equation, with currents This set of currents, ⃗ j = ( j 1,2 , j 2,3 , . . . ) T , then obeys the equation of motion We see that, up to a global shift, it is the dynamics of current fluctuations that is governed by the non-Hermitian lattice Hamiltonian h HN with Dirichlet boundary conditions In other words, the linearized dynamics in our system is indeed governed by the HNM, but imposing Neumann boundary conditions for the density fluctuations ε p .This is physically consistent with the assumption κ = 0 made in this analysis.This subtle change in the boundary conditions has an important consequence for the spectrum of h, namely the existence of a steady state.More precisely, in Appendix D we show that Spec{h} where Spec{A} L is the spectrum of matrix A in L dimensions.This means, first of all, that the spectrum of density fluctuations in the ASIP model shares all the spectral features of the HNM, which we discussed in Sec.5.1 above.In addition, there exists a unique steady state with E ss = 0 and a wavefunction Up to nonlinear corrections, which have been omitted in the current analysis, this wavefunction agrees with the stationary density profile derived in Eq. (22).Note that the existence and the shape of this steady state does not change when the coupling to the reservoirs is no longer neglected, since the term ∼ κ(ε 1 − m) in Eq. ( 43) merely fixes the magnitude of the fluctuation at the first site and ε L ∼ ψ ss (L) → 0.

Discussion
In Fig. 6, we plot the eigenvalues of h HN −2iΓ S , for Dirichlet and periodic boundary conditions, and compare them with the spectrum of h.These plots confirm that the eigenvalue structure of h mimics that of the shifted HNM, except for the existence of a steady state with E ss = 0.For open lattices, the non-zero eigenvalues coalesce near c = Γ S , which corresponds to the (L − 1)-th order exceptional point EP for J r = 0 in the HNM.By contrast, the steady-state mode remains well isolated and pinned at the origin.The explicit form of the steady-state solution in Eq. ( 51) confirms that this exceptional point coincides with the transition point into the zigzag phase in the full ASIP master equation.Hence, we have found a situation in which an EP for the higher-energy modes is directly connected with an observable configuration change in the steady-state.The spectrum of h coincides with the one of hHN , plus the isolated steady-state at the origin.The two insets at the top depict the shape of the steady-state eigenmode ψ ss before and after crossing the EP at a value of c = Γ S .These results show that the transition into the zigzag structure of the steady state of h coincides with the EP for its higher-energy modes.
In Ref. [40], this exponentially localized steady-state was also obtained, but only in the "smooth" phase, by considering the full spectrum of the hopping Liouvillian L hop restricted to the single-particle subspace.For a single boson, there are no nonlinearities, which corresponds to the limit n ∞ → 0 and c = Γ A /2.In this case, even in the fully asymmetric limit, Γ r → 0, the EP can be approached, but not crossed.A closely related analysis has also been performed for generalizations of the HNM with purely linear jump operators [39].Here, the many-body stationary states for bosons and fermions exhibit an accumulation of particles near one boundary, but these features are rather broad and the direct connection to the NHSE has not been found there.We conclude that while many of these models show formally similar excitation spectra, the crossing of the EP and the transition into the zigzag phase is related to the nonlinearity of the underlying transport equations, which allows us to fulfill the condition c > Γ S through a bosonically enhanced propagation speed.At the same time, while being a many-body effect, for Γ r → 0 the zigzag pattern can already be observed deep in the quantum regime, i.e., for an average density of n ∞ < 1 (see Fig. 7).
The correspondence between the EP in the fluctuation dynamics and the transition in the stationary density profile is actually quite surprising.Naively one would expect that the EP, which occurs at an imaginary offset of −2iΓ S , mainly influences the transient dynamics of decaying fluctuation modes.Instead, it signifies a sharp transition in the stationary fluctuation mode, which remains spectrally well isolated from the EP.This is in stark contrast to what is usually assumed for non-equilibrium phase transitions in dissipative systems, where the phase transition point coincides with a closure of the dissipative gap [46,47].The origin of this paradoxical situation can be traced back to the conservation of fluctuations, which, when decaying in site p, reappear in site p − 1. Fluctuations thus propagate across the chain, and fully decay only when they reach the edge.Therefore, while near the EP there is only a single eigenvalue that sets the timescale of the dynamics, it can still take a (diverging) time τ relax ∝ L for the system to fully relax.As pointed out in Refs.[36,40,43], it corresponds to the time required for the excitations to propagate along the chain.This distinguishes the analysis of such transport transitions from other non-equilibrium phase transitions in unbiased systems.

Summary and conclusions
In summary, we have studied the dissipative thermal transport of bosons through a lattice with asymmetric hopping rates, as described by the ASIP.Compared to analogous models for fermions or distinguishable particles, dissipative transport of bosons is characterized by hopping events that are accelerated by the presence of other particles.Our analysis showed that despite the simplicity of this process and without including any additional coherent interactions, this bosonic enhancement already gives rise to a highly non-standard transport phenomenology including ballistic currents, the formation of a boundary region with coexisting thermal and Bose-condensed sites, as well as the spontaneous development of coherence in a purely dissipative system.
In contrast to other condensation mechanisms that have been investigated for various (classical) inclusion processes [5,6], the predicted transition for bosonic transport relies on the presence of a boundary and would be absent in an infinite or periodic lattice.This creates a natural connection to the HNM and related non-Hermitian lattice models, which we discussed in full details in Sec. 5.This analysis establishes a direct correspondence between the EP in the complex excitation spectrum of the HNM and the transition point in the stationary density profile of the ASIP.It also shows that, while closely related to the NHSE, the formation of the zigzag phase is a genuine many-body effect that does not appear for single-particle or linear transport models.
In conclusion, this bosonic skin effect creates an interesting connection between transport physics, non-equilibrium phase transitions and non-Hermitian physics.For highly asymmetric hopping rates, which can be engineered, for example, for cold atoms in optical lattices, the predicted zigzag phase is already observable at the level of a few atoms.Instead, in nanophotonic lattices for optical photons or exciton polaritons, where one might only achieve a small, temperature-induced bias, the necessary condition c ∼ Γ S can still be reached by assuming pumped reservoirs with a considerably higher density.Therefore, since the main features associated with this transition are rather robust with respect to the details of the model, they should be observable in a variety of bosonic lattice systems, whenever hopping is predominantely incoherent.
Due to the exponentially growing configuration space, the exact dynamics of P({n p }, t) can only be calculated for very small lattices and low occupation numbers.Instead, for larger lattices we sample the probability distribution via a Monte-Carlo simulation.To do so, the boson numbers n p (t) for each site are treated as stochastic variables, which during an infinitesimal time step d t evolve according to Here, the d N l,r p = 0, 1 are independent random variables and indicate that a boson has hopped to the left (right) when d N l p = 1 (d N r p = 1).The probabilities for these events are . By starting from a given initial configuration, {n p (t = 0)}, and evolving a total number of N t stochastic trajectories in time, we can approximate the expectation value of any function of operators np by an ensemble average.For example, This method becomes exact in the limit N t → ∞, and therefore also accounts for cross-site correlations, C pq (t) = 〈n p nq 〉(t) − 〈n p 〉(t)〈n q 〉(t), which are neglected in mean-field theory.It cannot, however, be used to predict quantities such as cross-site coherences of the form 〈â † p âp+1 〉, because those involve off-diagonal elements of the density matrix.Furthermore, this method is limited to low average occupation numbers, since otherwise the rate of jumps, and therefore also the total simulation time, increases significantly.

B.2 High density regime: Truncated Wigner approximation
The TWA is a technique for simulating the dynamics of bosons in phase space, which is spanned by complex amplitudes α p and α * p defined on each site p.The state of the full lattice is then fully described by a multi-mode Wigner distribution W ({α p }, t) on this space and expectation values of symmetrically-ordered operator products can be obtained from the moments of this function.For example, To obtain the equation of motion for W ({α p }) we use the substitutions [57] â etc., to convert the master equation ( 1) for the density operator into a partial differential equation for W .To illustrate this approach, let us consider only a single term, where we have used the short-hand notation ∂ i = ∂ α i , and . Note that the same equation was derived in [73], where the two bosonic modes represented Schwinger bosons describing a d-level system.
The TWA consists in neglecting in this equation all third-order derivatives.This approximation is expected to be accurate when the number of bosons in the chain is high (see [57,74] for a more detailed discussion).Hence, this method provides a complementary treatment to the one presented in the previous section.After we performed the TWA, we obtain a Fokker-Planck equation, governed by a drift vector ⃗ A and a diffusion matrix D, where we have used Einstein's sum convention and the 2L greek indices run over all α p and α * p .For the example given in Eq. (B.6) above, the corresponding diffusion matrix is given by where we have ordered the four independent variables as (α 1 ,α * 2 ,α * 1 ,α 2 ).We have also omitted the constant terms ±1/4, which cancel when adding the contributions from all lattice sites, expect at the boundaries.For any other site, the diffusion matrix is positive semi-definite and can therefore be written as D = BB † with Therefore, it is possible to unravel the Fokker-Planck equation in terms of stochastic trajectories in phase space, which follow the (Ito) equations , where the dW i are complex-valued Wiener processes satisfying 〈dW i dW * i 〉 = 1 and 〈dW i dW i 〉 = 〈dW i 〉 = 0.By defining d V = (dW 1 − dW * 2 )/ 2, we can write the stochastic equations as This derivation can be generalized in a straightforward manner to all lattice sites and including the hopping to the right and the coupling to the reservoirs.Altogether we end up with the following set of stochastic differential equations where all the d V i are independent complex Wiener processes.Note that in the equation for α 1 , the diffusion rate in the last term can become negative, when the coupling to the left reservoirs is too weak.This problem does not occur in any of the presented results, where we assume κ l = Γ l .In this case, the noise processes ∼ d V 1 and ∼ d V l can be combined in a single stochastic process, and for Γ r = nl = 0 we obtain Eq. ( 31).

B.3 Benchmarking the mean-field approximation
In Fig. 7, we compare the results obtained with these two numerical methods with the predictions from mean-field theory in the limits of low and high occupation numbers.These plots show that all the features in the stationary density profile discussed in the main text are accurately reproduced by both methods, within their respective range of applicability.In particular, the exact results from the Monte-Carlo simulations demonstrate that the predicted density patterns are already visible in parameter regimes where there is on average less than one boson per site.We also find that the mean-field prediction for the transition point Γ c A is well reproduced by both methods (not shown here).

C Derivation of the stationary density profile
In this section we provide additional details about the derivation of the steady-state occupation numbers n p within the mean-field approximation.The starting point for this derivation is Eq.(19), which for L → ∞ already determines the relation between the current J and the asymptotic occupation number n ∞ , as given in Eq. (20) where µ = β/α.By rewriting the above result in terms of the ratio (n p − n ∞ )/(n 1 − n ∞ ) we obtain Eq. ( 22), from which the decay of the zigzag structure becomes more obvious.At this point, the parameters n ∞ and µ are still unknown and must be determined by the boundary conditions.Since in the steady-state the current is constant we obtain In the limit of a large lattice, we can set n L = n ∞ , which gives us a quadratic equation for n ∞ , Γ A n 2 ∞ + (Γ A + κ r )n ∞ = κ r nr , with a solution displayed in Eq. ( 21).Finally, from the current into the left reservoir and the result for n p=1 in Eq. (C.5) we can determine the value of µ.For nl = 0 and Γ r = 0, its explicit expression is In the most general case, its precise functional dependence on all the system parameters is complicated and of limited interest.

D Eigenstates of HNM with Neumann boundary conditions
In this appendix we derive the relation between the spectra of h HN and h, which correspond to the HNM with Dirichlet and Neumann boundary conditions, respectively.An alternative derivation, and further results on these kinds of matrices, can be found in [72].We introduce here the L-dimensional current vector ⃗ j = ( j 0,1 , j 1,2 , j 2,3 , ...) T , which includes the component j 0,1 = 0.Then, according to Eq. ( 47), we obtain the linear relation This means that for each of the L − 1 eigenfunctions ⃗ ψ k of h HN we obtain an eigenvector of h, with a modefunction ∇ ⃗ Φ k , and eigenenergy E k = E HN k − 2iΓ S .Moreover, since det(V ) = 0, there is one additional eigenstate ψ ss with energy E ss = 0, which satisfies V ψ ss = hψ ss = 0.It is straightfoward to check that this eigenstate is of the form given in Eq. (51).Putting these two sets of eigenstates together, we finally obtain the result of Eq. ( 50).

Figure 1 :
Figure 1: Asymmetric bosonic transport.(a) Sketch of the ASIP setup studied in this work.Bosons injected from a thermal particle reservoir with mean occupation number nr on the right can incoherently hop along the lattice with asymmetric rates Γ l and Γ r , before being emitted into a second reservoir with occupation number nl on the left.A directional hopping can be imposed, for example, by applying a potential gradient with an energy offset U between neighboring sites.(b) Under stationary conditions, this hopping asymmetry combined with the bosonic particle statistics results in the bosonic skin effect, i.e., the formation of a finite boundary region with a staggered density profile.The two insets show sketches of the Wigner distribution for individual lattice sites, indicating that within this boundary region, the odd sites are in a condensed state with broken U(1) symmetry, while all other lattice sites exhibit a thermal distribution.See text for more details.

Figure 2 :
Figure 2: Plots of the steady-state occupation numbers n p for a lattice of L = 15 sites and different degrees of asymmetry:Γ A /Γ l = 0 (a), Γ A /Γ l = 0.05 (b), Γ A /Γ l = 0.17 (c), and Γ A /Γ l = 1 (d).For all plots nr = 10 and two different values of nl = 0 (blue lines) and nl = 20 (yellow lines) have been considered.The insets show the current J versus the lattice size L, in log-log scale and for three different values of nl = 0, 5, 9.For Γ A = 0 we recover a linear population gradient and the Fourier law for the current, as expected for diffusive transport.For any Γ A > 0 and large L, the current becomes independent of both L and nl , indicating ballistic transport.In this regime, we observe the formation of a finite boundary region of size ξ, as indicated by the shaded area.As the asymmetry increases, the width ξ shrinks and vanishes for Γ A /Γ l ≃ 0.17.Beyond this point, a finite boundary region, but with an oscillating density profile, reappears.For all plots, we have set κ r = κ l = Γ l .

Figure 3 :
Figure 3: Dependence of the skin length ξ as defined in Eq. (25) on the hopping asymmetry Γ A and on the thermal population of the right reservoir, nr .When Γ A is exactly zero (thick dark line at the bottom of the diagram), we recover the usual diffusive behavior.The dashed line corresponds to Γ A = Γ c A , at which point ξ = 0. Below (above) this line, the steady-state population exhibits a smooth (zigzag) profile near the left boundary.The inset shows ξ along the horizontal green line at Γ A = 0.5Γ l .For all points in this plot a value of κ r = Γ l has been assumed, and the results are independent of both nl and κ l .

Figure 4 :
Figure4: (a) Plot of the second-order correlation function g(2) (0) for a lattice of L = 10 sites, as obtained from a TWA simulation with 5000 trajectories.The phasespace distributions below each point indicate the distributions of the amplitudes α p , in the complex plane, at the final time of the simulation.(b) Distributions of the values of |α p | 2 and (c) plots of the coherence function g(1) (τ)(30) for odd (left) and even (right) sites near the boundary (all the even sites have extremely similar behavior, here we show only p = 2 for simplicity).For all plots we have set κ = Γ l , Γ r = 0, nr = 30 and nl = 0.For the plots in (c), we have used a reference time of t = 10Γ −1 l , which is sufficient to reach the steady state.

Dirichlet periodicFigure 6 :
Figure 6: Complex spectrum of shifted HNM, hHN = h HN −2iΓ S 1, for different boundary conditions.The green squares ('periodic') and the blue dots ('Dirichlet') represent the eigenvalues of hHN on a lattice of L = 19 sites with periodic and open boundary conditions, respectively.The red triangles ('Neumann') are the eigenvalues of h as given in Eq. (45) for a lattice of L = 20 sites.The four lower panels show the complex spectra of these Hamiltonians for different values of c, increasing counter-clockwise.The spectrum of h coincides with the one of hHN , plus the isolated steady-state at the origin.The two insets at the top depict the shape of the steady-state eigenmode ψ ss before and after crossing the EP at a value of c = Γ S .These results show that the transition into the zigzag structure of the steady state of h coincides with the EP for its higher-energy modes.

Figure 7 :
Figure7: Top row: Comparison between the steady-state populations as predicted by mean-field (MF) theory and by the TWA, for nr = 10 and for Γ r /Γ l = 0.95 (left) and Γ r = 0 (right).Bottom row: Comparison between the steady-state populations as predicted by mean-field theory and by exact Monte-Carlo (MC) simulations, for nr = 1 and Γ r /Γ l = 0.5 (left) and Γ r = 0 (right).For these results we simulated 5 × 10 3 trajectories for the TWA and 5 × 10 5 trajectories for the Monte-Carlo method.For all plots, κ r = κ l = Γ l , nl = 0, and L = 10.