Bistabilities and domain walls in weakly open quantum systems

Weakly pumped systems with approximate conservation laws can be efficiently described by a generalized Gibbs ensemble if the steady state of the system is unique. However, such a description can fail if there are multiple steady state solutions, for example, a bistability. In this case domains and domain walls may form. In one-dimensional (1D) systems any type of noise (thermal or non-thermal) will in general lead to a proliferation of such domains. We study this physics in a 1D spin chain with two approximate conservation laws, energy and the $z$-component of the total magnetization. A bistability in the magnetization is induced by the coupling to suitably chosen Lindblad operators. We analyze the theory for a weak coupling strength $\epsilon$ to the non-equilibrium bath. In this limit, we argue that one can use hydrodynamic approximations which describe the system locally in terms of space- and time-dependent Lagrange parameters. Here noise terms enforce the creation of domains, where the typical width of a domain wall goes as $\sim 1/\sqrt{\epsilon}$ while the density of domain walls is exponentially small in $1/\sqrt{\epsilon}$. This is shown by numerical simulations of a simplified hydrodynamic equation in the presence of noise.


Introduction
Generalized Gibbs ensembles are powerful concepts to (approximately) describe complex interacting systems with a set of exact or approximate conserved quantities Q i , i = 1, ..N . The density matrix ρ obtains the simple form ρ ∼ e − i λ i Q i , where the Lagrange parameters λ i are in one-to-one relation to the expectation values of the Q i [1][2][3]. This approach can even be used if the conservation laws are only approximately valid and if the system is weakly driven out of equilibrium as long as scattering processes which conserve the Q i dominate the dynamics. For example, to describe the Bose-Einstein condensation of exciton-polaritons or photons [4][5][6][7][8], it is useful to introduce a chemical potential for these particles despite the fact that particle number is not exactly conserved in the systems. The value of the chemical potential is then determined by balancing loss and pumping rates. Similarly, in solid state materials driven out of equilibrium, e.g., by a short laser pulse, one can use the weak coupling of phonons to electrons to introduce two different temperatures for the two subsystems. Here the relevant approximately conserved quantities Q 1 and Q 2 are the energies of the phonon and electron system, respectively. The corresponding λ i are identified with their inverse temperatures. Simple rate equations then describe the time-evolution within such two-temperature models [9]. Recently, we have generalized this notion also to approximately integrable systems with an infinite number of -approximate -conservation laws, where we could show that one can create giant heat and spin currents in driven spin chains [10]. Similar concepts can also be used to describe many-body localized systems coupled to phonons and an external drive [11].
V m m Figure 1: Left: Single particle in a 0D system whose deterministic dynamics is governed by a symmetric double-well potential. In the absence of noise, the symmetry of the potential is broken by the initial conditions. However, any finite noise triggers consecutive transitions between the two minima of the potential and eventually restores the symmetry in the longtime limit. Right: A 1D generalization of the aforementioned case for a field theory. Similarly to the 0D case there is no symmetry breaking at finite noise strength. Instead, depending on the noise strength a non-zero density of domain walls forms.
In this work, we want to study in a controlled way a weakly driven system with approximate conservation laws where the concept of a (generalized) Gibbs ensemble breaks down. Starting from a 1D system with just two exact conservation laws (energy and magnetization), we add weak perturbations of strength which break the corresponding symmetries and drive the system out of equilibrium. We choose the perturbations in such a way that they induce a bistability in the magnetization and argue that noise terms naturally generate domains and domain walls in such systems. The emergence of bistabilites in non-equilibrium systems in the limit of strong drive and dissipation has recently gained increased attention due to experimental observations, in e.g. driven Rydberg ensembles [12,13], nonlinear photon lattices [14][15][16], semiconductor microcavities [17], and QED with cold atoms [18,19]. In low-dimensions a symmetry cannot be spontaneously broken due to thermal fluctuations. A well known zero-dimensional example is the supercritical pitchfork bifurcation [20] For the deterministic part of the dynamics (α = 0) one obtains for µ > 0 two stable solutions at ± √ µ, see Fig. 1. In the absence of noise, the symmetry of the underlying symmetric doublewell potential is broken by the initial conditions. However, at any finite noise strength α = 0, the symmetry is restored in the long-time limit and the corresponding Focker-Planck equation yields P (x) ∼ exp(−V (x)/α) as a stationary probability distribution. Ref. [21] discusses, for example, such a zero-dimensional case by investigating a Dicke model with non-linear noise. Similarly, in 1D systems with short-ranged interactions arbitrarily weak noise will generically induce domain walls thus rendering any description in terms of noiseless (generalized) Gibbs ensembles invalid. Here the (finite) energy cost of a domain wall plays a similar role as the potential barrier of the zero-dimensional example, see Fig. 1. In dimensions larger than one, in contrast, an Ising symmetry can be spontaneously broken even in the presence of (sufficiently weak) noise [22,23]. In the following, we will investigate a simple 1D model which allows one to study the role of approximate conservation laws, the validity of Gibbs ensembles and the relations of bistabilities and noise in a controlled way. We discuss how an effective description in terms of (noisy) hydrodynamics can be obtained and solve a simplified version of these equations numerically.

Model
We consider a antiferromagnetic (J > 0) one-dimensional XXZ spin chain The next-nearest neighbor interaction J renders the model non-integrable. The unperturbed Hamiltonian H 0 therefore has only two conservation laws: the magnetization in z-direction and the energy, We assume that the system is driven out of equilibrium by the weak coupling to a Markovian bath. The dynamics of the density matrix ρ is thereby governed by the Liouville equatioṅ where, importantly, is assumed to be small. We aim to construct the coupling such that the dynamics exhibits a (local) bistability in the presence of noise. To achieve this goal we consider two competing perturbationsD i (i=1,2) whose relative strength is controlled by the parameter γ ∈ [0, 1]. As Lindblad operators we choose While the first Lindblad operator flips spins which leads to noise and heating, the second Lindblad operator aligns neighboring spins by transforming ↑↓↑ to ↑↑↑ and ↓↑↓ to ↓↓↓. Therefore it naturally induces a bistability in the total magnetization of the system. For γ = 1, i.e., in the absence of the σ x term, the fully polarized states | ⇑ = | ↑ . . . ↑ and | ⇓ = | ↓ . . . ↓ are the two unique dark states of the system and the steady-state density matrix is simply . Importantly, this steady state is completely noiseless. In the following we will only consider the situation where the noise is finite, γ < 1.

Hydrodynamic Approximations
For = 0, in the absence of any coupling to an environment, the steady-state density matrix in the thermodynamic limit is simply given by Here the parameters λ 1 and λ 2 are not fixed by the dynamics but only by initial conditions. Scattering processes of the non-integrable interacting system are essential to establish this steady state. For a finite, but tiny value of it is clear that the system will remain locally close to such states (for a quantitative discussion of corrections we refer to Ref. [24]). If such stationary states are not unique (e.g, due to a bistability), we can, however, not expect that locally the same values of λ i are obtained as we will show in detail below. Instead, we should parametrize the system with space-dependent Lagrange parameters λ i (r). This leads to the following ansatz for the density matrix Here we integrate (in the functional integral sense) over smoothly varying space-dependent Lagrange parameters λ i (r). Theq i (r) are the (coarse-grained) local charge density operators with Q i = drq i (r), Z[λ i ] is the partition sum for a fixed configuration of λ i (r), and δρ is a correction to the density matrix arising from gradients of λ i (r), briefly discussed below. The (yet unknown) functional P t [λ i (r)] describes the (classical) probability for a given configuration of Lagrange parameters defined by λ 1 (r) and λ 2 (r). In general, P t [λ i (r)] depends on time. It describes the dynamics on a time scale of order 1/ , which is assumed to be much larger than all internal equilibration times [25]. Instead of developing directly a theory for the probability distribution P t [λ i (r)] in the spirit of a Fokker-Planck equation, we will use a description in terms of a (generalized) Langevin equation for the fields λ i (r) or, equivalently, the corresponding local expectation values of the charge densitiesq i (r). This approach has the advantage of being much more intuitive. In the following we use q i (r, t) to denote the expectation values of the coarse-grained local densities for one realization of the underlying Langevin process. Technically, we will perform a gradient expansion around the homogeneous solutions q i (r) = const. To zeroth order in the gradient expansion, we can locally approximate the density matrix close to the position r 0 by We use this density matrix to compute the change of the conserved charge densities linear in Here L is the size of the system and f i is the averaged, deterministic force which depends on the local Lagrange parameters λ i (r 0 ), or equivalently on the local densities q i (r 0 ) (i = 1, 2) evaluated at r 0 . Within our Langevin approach we expect that the coupling to the bath also leads to a noise term ξ i (r), with ξ i = 0 and ξ i (r, t)ξ j (r , t ) ≈ N ij δ(t − t )δ(r − r ). As we describe in Appendix A.2, the noise matrix N ij can be computed from the time evolution of Q i Q j [26]. Importantly, both the forces f i and the noise matrix N ij are linear in as they arise both from the coupling to the Lindblad operators. Both are also functions of the local charges q i (r, t). We compute f i and N ij using exact diagonalization of small systems. For the parameters investigated by us the effective temperatures are rather high and thus finite size effects turn out to be small. Fig. 2 displays the forces for γ = 0.9. As expected from our construction, we obtain two stable fixed points at a magnetization of approximately m ≈ ±0.7. We denote the values of conservation laws at the fixed points as q F P 1 i and q F P 2 i with q F P 1 1 = −q F P 2 1 and q F P 1 2 = q F P 2 2 by symmetry. In the absence of the noise term, these two solutions would lead to spontaneous symmetry breaking. We also find an unstable fixed point at m = 0 and e ≈ 0.55. It is important to note that in the presence of approximate conservation laws even a tiny coupling to a non-equilibrium bath can strongly modify the system. In our example a state with a large magnetization and high energy is approached in the long-time limit even for very small perturbations .
As a next step, we have to compute contributions to ∂ t q i arising from terms proportional to gradients of the local charges ∂ r q i . Due to the space-reflection symmetries all linear gradients vanish on average. The other gradient terms can be calculated for = 0 as they are finite in this limit. Their form is well known from standard hydrodynamics and we obtain Here, D ij is the matrix of diffusion constants of the unperturbed model H 0 defined by j i = −D ij ∂ r q j where j i is the current of the conserved densities q i . Technically, they arise from the correction δρ in Eq. (7) which induces gradients of the Lagrange parameters. The matrix of diffusion constants depends on q 1 and q 2 and can at = 0 be computed using Kubo formulas evaluated in the corresponding thermal Gibbs state. The first two terms on the right-hand side have been copied from Eq. (10). The last term, again computed for = 0, is the usual thermal noise with where are the susceptibilities of the Q i . Note that the thermal noise ∂ r ξ th i obeys the conservation laws as it is proportional to a derivative while the nonequilibrium noise ξ i does not. The equations (11) describe the hydrodynamics of our driven system and we expect that they are exact in the limit of small as they have been derived in a systematic expansion in and gradients, keeping always the leading corrections. To understand their properties in the limit of small , it is useful to rewrite the equations using rescaled variables. Employing that the forces are linear in , we introduce rescaled variables, τ = t , x = r √ ,f = f / ,ξ = ξ/ ,ξ th = ξ th / √ . In these variables, we obtain equations which have exactly the same form as Eqs. (11), The only difference is that nowf is independent of and the only dependence arises from the two noise terms which both turn out to be proportional to √ , This immediately shows that both noise terms are of equal importance for our hydrodynamic theory. Furthermore, the analysis justifies a posteriori the gradient expansion underlying the derivation of our equation: higher order terms are suppressed by powers of √ . All parameters of our hydrodynamic equations can in principle be calculated from correlation functions of the unperturbed system H 0 only. By far the most difficult part of the calculation is the numerical determination of the diffusion constants D ij of the unperturbed system as function of the q i . While there has been an enormous recent progress in the numerical calculation of transport coefficients in 1D systems [27], this is still a challenging problem suffering from huge finite size effects. As all of our qualitative results do not depend on the numerical values and functional dependence of the transport coefficients, we are not trying to calculate those. Instead, we will use in the following mainly the scaling arguments given above in combination with a numerical investigation of a strongly simplified version of Eqs. (11). First, instead of tracking the dynamics in the two-dimensional space q 1 and q 2 , we concentrate on the magnetization m as this is the variable which shows a bifurcation. Second, we approximate the q i dependent matrix of diffusion coefficients by a single constant D. Finally, we adjust the forces of the right-hand side accordingly and obtain a strongly simplified toy model for the fluctuation induced domain-wall formation As we are only interested in the dependence of our result, we approximate the force by ) (a is the lattice constant of the microscopic model), and ξ th ξ th = 2DχT δ(r − r )δ(t − t ), where we simply set χT = a/2 for our toy model. The functional form used for f and the non-thermal noise is motivated by the infinite temperature limit where one can easily calculate all terms analytically, see App. A.1 and App. A.2. Within our toy model, a bistability is obtained for γ > 4/5 in the noiseless case. To analyze the properties of Eqs. (11) (and its simplified version Eq. (15)), we first consider the noiseless limit by neglecting ξ and ξ th . In this case, two trivial solutions are given by the fixed points, q i = q F P 1 i and q i = q F P 2 i . More importantly, there is also a 'domain wall' solution obtained from the boundary condition lim r→−∞ q i = q F P 1 i and lim r→∞ q i = q F P 2 i . As it is obvious from our scaling analysis, the width of the domain wall is proportional to 1/ √ . For our toy model, one can calculate the shape of the domain wall also analytically,  Fig. 3 shows such a domain wall for the toy model where it is compared to our numerical results. In our numerical simulations we use rescaled variable where length and time are measured in units of a and a 2 /D, respectively. Equivalently, one can set D = J = a = 1 and replace by˜ = (Ja 2 /D). We discretize space in steps of size 0.25 and time in steps of 0.001, using Heun's method for integration [28]. As discussed in the introduction, we expect that for any finite strength of fluctuations, a finite density of such domain walls occurs in the steady state. This is confirmed by simulations of our simplified model shown in Fig. 4 for different values of . The figure shows m(r, t) after some initial waiting time in which the system obtains its (fluctuating) steady state. For = 0.1, domains are huge but their size drops rapidly when˜ is increased. The time scale which governs a reversal of the local magnetization depends also strongly on . In Fig. 5 we show the density of domain walls, or equivalently, the inverse distance of domain walls. To estimate the density of domain walls analytically, we can use well-known results obtained for equilibrium systems. Here it is important to note that our effective theories Eqs. (11) and also Eq. (15) are not equivalent to an equilibrium theory (they will, for example, not fulfill the second law of thermodynamics) as the two noise terms do not encode thermal noise of a single temperature. If we, however, switch off the noise contribution ξ th in Eq. more precisely The effective temperature T eff = J(1 − γ) is set by the strength of fluctuations of ξ and therefore linear in . Hence, we expect that the density of domain walls is proportional to e −E DW /T eff or n DW ∼ e −c/ √ for 1 (18) where c = √ 2D(5γ−4) 3 3a √ Jγ(1−γ) if we only include fluctuations from ξ, ignoring corrections from ξ th . Our scaling analysis, Eqs. (14), strongly suggests that these results also hold when the second noise term ξ th is switched on as it has the same scaling properties. Only the prefactor c should become smaller when an extra source of noise induces more domain walls. This is confirmed by our numerical results. The solid line in Fig. 5 is a fit to Eq. (18). The fit works very well for 0.08 ˜ 0.4. Deviation for very small values of˜ arise from finite size effects when the distance of domains 1/n DW becomes of the order of the system size (L = 10.000 in our simulation). Within our numerics we obtain when including ξ th a value of c =c ≈ 0.77 (D = J = a = 1) that is indeed smaller than our analytical prediction c ≈ 0.99 obtained for the model without thermal noise. We have also compared the analytical result to the numerics when considering only fluctuations due to ξ, thus validating our numerical approach.

Discussion
Weakly driven classical and quantum systems can exhibit properties with no equilibrium analogy. Our example shows, that even a very weak driving term can induce ferromagnetism in an antiferromagnetic system. In contrast, very large Hamiltonian perturbations are needed to transform an antiferromagnet to a ferromagnet. Nevertheless, phase transitions in the driven system share many similarities with finite-temperature phase transitions, at least in cases where the stationary points of the Lindblad evolution are not noiseless absorbing dark states [29][30][31][32]. We have shown that for a weakly-driven system with approximate conservation laws one can describe the physics at large length scales by noisy hydrodynamic equations which are similar (but not identical) to corresponding equations for thermal systems. An important consequence of the noise is that phase transitions only occur in dimensions larger than one.
In the one-dimensional example analyzed by us one finds instead that at each finite noise strength a finite density of domain walls with a density proportional to e −c/ √ arises. The width of the domain walls are not determined by energetic arguments but instead by the interplay of diffusion and the drive with strength . Therefore the width scales with 1/ √ . The origin of this peculiar behavior is the (approximate) conservation of magnetization in the system, which ultimately allows one to drive the system far from equilibrium by only weak perturbations.
Our analysis can also be seen as an example of a weakly driven system which can not be described simply by a (generalized) Gibbs ensemble as used by us, e.g., in Ref. [25]. Due to the existence of several fixed points and due to strong fluctuations effects in low-dimensional systems, it is necessary to consider instead ensembles of (generalized) Gibbs ensembles with fluctuating Lagrange parameters. In more simple situations, where only a single attractive Gibbs state exists, one can instead expect that large fluctuations are sufficiently rare to allow for systematic expansions around (generalized) Gibbs states [24]. We expect that the notion of fluctuating hydrodynamics will also be very useful to explore the physics of driven approximately integrable systems with an infinite number of conservation laws.

A.1 Generalized Forces
The generalized force f = (f 1 , f 2 ) can to leading order in be calculated with the formula r 0 (i = 1, 2) which yields In the infinite temperature limit the force simplifies to f 1 (m) = J γ 4 (1 − m 2 ) − (1 − γ) m.

A.2 Noise
Here, we derive the generalized Einstein relation Next we use the approximation O β (t) − O β (t − ∆t) = t t−∆t dt Ȯ β (t ) to write This finally gives We can use this relation to determine the noise-noise correlation matrix. As an example we calculate ξ 1 ξ 1 which is used in the numerical simulation of our toy model. While the unitary part of the dynamics and the second Lindblad term do not yield a contribution, the first Lindblad term yields where Γ = J (1 − γ).