Role of interactions in a dissipative many-body localized system

Recent experimental and theoretical efforts have focused on the effect of dissipation on quantum many-body systems in their many-body localized (MBL) phase. While in the presence of dephasing noise such systems reach a unique ergodic state, their dynamics is characterized by slow relaxation manifested in non-exponential decay of self-correlations. Here we shed light on a currently much debated issue, namely the role of interactions for this relaxation dynamics. We focus on the experimentally relevant situation of the evolution from an initial charge density wave in the presence of strong dephasing noise. We find a crossover from a regime dominated by disorder to a regime dominated by interactions, with a concomitant change of time correlators from stretched exponential to compressed exponential form. The strongly interacting regime can be explained in terms of nucleation and growth dynamics of relaxing regions - reminiscent of the kinetics of crystallization in soft matter systems - and should be observable experimentally. This interaction-driven crossover suggests that the competition between interactions and noise give rise to a much richer structure of the MBL phase than anticipated so far.


I. INTRODUCTION
Many-body quantum systems in the presence of quenched disorder undergo a transition between an ergodic phase and a many-body localized (MBL) phase [1][2][3][4][5] . While the transport properties of the MBL phase are still debated 6 it is generally accepted that it is characterized by a slow growth of entanglement entropy 7-10 , and ergodicity breaking which has been observed in numerical studies 5,11,12 and experiments [13][14][15][16] . Some nonergodic aspects of the MBL phase have also been argued to be present with translation invariance 17,18 .
While most literature has focused on closed quantum, the imperfect isolation of the cold atomic ensembles used in recent experimental observations of MBL calls for an understanding of the effect of dissipation on the MBL phase 19,20 . In Ref. 21 a chain of interacting fermions in contact with an infinite temperature dephasing bath was studied numerically. At conditions where the closed system would be in the MBL phase, a slow approach to the infinite temperature state was observed in the open system, characterized by a stretched exponential decay of self-correlations. Stretched exponential behavior was confirmed analytically in Ref. 22 in terms of a non-interacting (Anderson) system valid for large disorder. Similarly, in Ref. 23 the scaling properties of the same system were studied in the large disorder limit, finding independence of the dynamics from interactions.
A central question is therefore whether interactions play any role in the relaxation to the ergodic state due to dephasing in an otherwise MBL system. Here we address this question by studying the dissipative dynamics of a disordered XXZ chain in its MBL phase 21 . Our main result is that depending on the interaction strength the system explores two different regimes within the MBL phase. We show that in the dynamics of an initial spindensity wave, depending on interaction strength, there is a crossover from a disorder dominated regime to an interaction dominated regime, whose observable signature is a change of behaviour of self-correlators from a stretched exponential to a compressed exponential dependence with time. This latter behaviour is due to nucleation and growth of relaxing regions. A crossover of this sort is often a manifestation of non-equilibrium and aging behaviour in soft matter and glassy systems, see for example [24][25][26][27] .

II. MODEL
We consider a paradigmatic MBL system, the disordered XXZ chain in a spinless fermions description, where we denote withĉ † k the fermion creation operator, withn k =ĉ † kĉ k the number operator, and the random field h k ∈ [−h, h] is independently drawn for each site from a uniform distribution. This model exhibits an MBL transition for h c /J 7.2 5,10,28 .
Following 21 we couple the system to an infinite temperature Markovian dephasing bath. At weak coupling between system and bath, the dynamics can be described by a Lindblad quantum Master equation 29,30 where ρ is the system's density matrix and γ ≥ 0 sets the coupling to the bath. Eq. (2) has the advantage of being experimentally relevant, as it can be derived from microscopic principles for experiments on both cold fermionic 31 and bosonic 32 gases in the lowest band of an optical lattice. The decoherence is caused by off-resonant scattering of photons forming the lattice potential, and the dissipation rate γ is controlled by the detuning and intensity of the trapping laser. Eq. (2) conserves fermion number and in what follows we focus on the half-filling sector.

III. RATE EQUATION DESCRIPTION
For times t 1/γ, there are two situations described by the Master equation (2) which may be reduced to a classical rate equation. The first is the limit of large dephasing, γ J [33][34][35][36] . The second is that of large interactions and/or large longitudinal fields, V, h J 23 . This effective dynamics describes the evolution of the diagonal elements of the density matrix p α . If we express these in terms of the probability vector |p = α p α |α , where |α are the N !/(N/2)! 2 Fock states in the half-filling sector (see also 23 ) Eq. (2) reduces to where P k =n k +n k+1 − 2n knk+1 . Eq. (3) describes classical hopping of particles on the lattice, with a rescaled time τ = J 2 γt/h 2 , cf. 22,23 . The rate for hopping between site k and k + 1 is given by where ∆h k = h k+1 − h k 37 and n k = Tr(n k |p ) is the total probability of having an excitation on site k. In the following we will fix the energy scale to γ = 1. The rates Γ k are configuration-dependent, see Fig. 1(a). In the dynamics described by Eq. (3), the rates Γ k act as a kinetic constraint 33 , as often encountered in systems with complex relaxation like glasses 38 : the form of the rates Γ k does not determine the properties of the stationary state, but rather the relaxation pathways.
The rates Γ k are random through the field h k . Their distribution P (Γ) depends on the strength of the interactions and on the specific configuration under consideration. The analytical expression of P (Γ) is given in Appendix A and plotted in Fig. 1(b) for various values of V for the configurations of the left column of Fig. 1(a). When V < 2h, the distribution is bimodal, with a peak at Γ/h 2 = 1 from values of the field such that ∆h k = ±V , [cf. Fig. 1(a)], and another peak at values Γ/h 2 ∼ 4/(3h) 2 (for V 2h). The form of P (Γ) changes qualitatively when V > 2h when ∆h k = ±V is not accessible, and the distribution P (Γ) becomes unimodal retaining only the slower peak, which for V 2h is centered at Γ ∼ V −2 . The qualitative change of P (Γ) already hints at the different dynamical regimes depending on the strength of the interactions.

IV. DISTINCT DYNAMICAL REGIMES WITHIN THE MBL PHASE
To explore the relaxation dynamics we focus on the case in which the initial state is the charge density wave (CDW) state, where the corresponding probability is • denoting empty and occupied sites, respectively. This is relevant for recent experiments [13][14][15][16] , where ergodicity properties of were studied via the evolution of an initial CDW quantified by the imbalance which gives a direct readout of the self-correlations and thus accounts for the ergodicity properties of the system. In Fig. 2 we report our results on the imbalance, obtained by realising Eq. (3) via kinetic Monte Carlo, averaged over disorder,Ī. The decay ofĪ becomes slower for increasing interactions. We quantify this slowing down by defining the saturation time T such thatĪ(T ) = e −2 . As shown in Fig. 2(a) we observe two different regimes: For V < 2h the saturation time shows little dependence on V , while for V > 2h it increases with increasing interaction, signalling a slowdown of the dynamics. The inset shows that in the region V < 2h, while T is approximately independent of V , the shape of the relaxation function depends on the strength of the interaction.
Our data is well fitted by the functionĪ(τ ) ∼ exp − (τ /T ) β . This form is motivated by the analytical arguments below. The results on the exponent β and the time-scale T are reported in Fig. 2(b)-(c). We find that at V 2h the relaxation of the imbalance switches from a stretched exponential behavior (β < 1) to a compressed exponential behavior (β > 1) (see Fig. 2(b)). Despite the rapid increase in β at V 2h which suggests a sharp acceleration of the dynamics, the increase in the timescale T combines to give the slowing down observed in Fig. 2(a). The minimum in T at V ∼ 0.3h, and the large V behavior are compatible with the results in 22 . A finite size study for the exponent β is shown in the inset of Fig. 2(b). Although in the stretched exponential regime (V < 2h) finite size effects have a marginal impact, in the compressed exponential regime (V > 2h) they cause a saturation of the exponent to lower values. The origin of this behavior will become clear below.
The two regimes of stretched and compressed exponential decay, for V 2h and V 2h, respectively, can be argued as follows. When V 2h the dynamics is dominated by disorder, and we can set V = 0. In this case the long time dynamics is characterized by large portions of the chain in which the system has relaxed (giving null contributions to the imbalance), with isolated non-relaxed pairs of neighboring sites corresponding to the largest ∆h k . The approach of I(τ ) to equilibrium is then determined by those sites. Their dynamics can be studied by focusing on, say, sites k and k +1 with a single excitation between them with relaxed neighbors serving as a bath. That is, we set n k to the stationary average 1/2 for k > k + 1 or k < k. This setting is sketched in Fig. 3(a). We calculate the population as n k ≡ Tr(n k |p ) where the Tr operator is the trace, or sum, of the probabilities of the state i.e. Tr |p = α p α . We consider the two sites in the initial configuration (k, k + 1) = (1, 0), and since the dynamics conserves the number of excitations we can set p ↑↑ = p ↓↓ = 0. The evolution of n k is given asṅ such that, using Eq. (3), we obtaiṅ As anticipated we set n k+2 and n k−1 to their relaxed value n k−1 = n k+2 = 1/2, such that we can write the equations for the sites under consideration aṡ We are interested in the local imbalance I k = n k+1 − n k which can be obtained by integratinġ (9) In this case, the rates in Eq. (4) depend only on the difference of the random fields on the sites they are connecting. The rates associated to two contiguous links (e.g., Γ k and Γ k+1 ) are therefore not statistically independent, since they both depend on the field on the site they share, but those of links further apart are. When solving Eq. (9) we can treat Γ k−1 and Γ k+1 as independent. As an approximation we set them equal when averaging over the disorder Γ k−1 = Γ k+1 = Γ , leading tō where P (Γ k , Γ ) is the joint probability reported in Appendix A. A numerical integration of Eq. (10) gives a stretched exponential behavior. In Fig. 2(b) the results on β obtained by fitting are compared with the numerical data in the weak interaction regime, showing good agreement. In fact the local imbalance is a good approximation of the imbalance in the regime we are considering. This is because non-relaxed links are far-apart enough to be considered independent, and they give the same average contribution.
In the opposite limit V 2h, in contrast, the first step costs V + ∆h k when starting from the CDW state. This sets the time-scale τ ∼ V 2 to observe transitions of the kind ...
. We refer to these as nucleation events, happening at homogeneous rate Γ n 2(h/V ) 2 . After a nucleation event has occurred relaxation can proceed via transitions like .., whose rates are independent of V [see e.g. Fig. 1(a)]. This dynamics is conveniently described by the following coarse-grained approximation. Since the imbalance is a quantity with a period of two sites it is natural to divide the chain in the following way: (1, 2)|(3, 4)|...| (N − 1, N )|. We focus then on a set of new degrees of freedom labeled by the contribution that the pairs of original sites bring to the imbalance: The CDW configuration corresponds to the |1, 1, 1, 1, ..., 1, 1 state, and a nucleation event creates either two 0 sites or a -1 site as shown in Fig 3(b). In what follows we will focus on the case in which two 0 sites are created. Once a nucleation has happened the events |0, 1 ↔ |−1, 0 , and |1, 0 ↔ |0, −1 are possible with the non-interacting rate Γ k = h 2 /(1 + ∆h 2 k ). This pair of reversible processes implies that the 0 sites can be treated as random walkers, which when moving away from each other create a growing region of −1 sites. The site dependent differences between these rates are small for the time-scales we are considering (V 2h), and we will assume a constant rate Γ e . When averaged over random realizations the extension of the region between the 0 sites expands following the lawḠ e (τ ) ∼ √ Γ e τ , contributing with net zero imbalance since the sites falling in this region are now equally likely to be a 1 or −1. This growth dynamics together with the initial nucleation events -reminiscent for example of a crystallization process -is well described by the so-called Avrami law 39-43 . Here we give a sketch of the derivation in our case. The average number of nucleation events up to a given timeν(τ ) can be found by integratingν = N Γ n /2. Not accounting for overlap of the expanded regions, the total number of transformed sites can be expressed as This dynamics is sketched in Fig. 3(c) for a single expanding region of transformed sites. Overlaps can be excluded by assuming the increment in transformed sites dN tr is proportional to dN multiplied by the probability of not having an already transformed site (1 − 2N tr /N ), giving Initializing our dynamics in the untransformed state the imbalance at a given time is given byĪ(τ ) = 1 − 2N tr (τ )/N , leading to the compressed exponential behavior with exponent β = 3/2 observed in Fig. 2(b). Equation (12) also yields the functional dependence of the time-scale T ∼ (V /h) 4/3 for large V /h, which is confirmed by our numerical results in Fig. 2(c). This picture breaks when the distance between nucleation events becomes comparable to the system length. In this case we can consider the expansion of a nucleated region as instantaneous and the imbalance as fully relaxed after a single nucleation event. In a single realization we can then model the imbalance as I(τ |τ ) = 1 − θ(τ − τ ), where τ is the time at which the first nucleation event happens. The probability of nucleation at this time is given as π(τ ) = N exp −N τ /V 2 /V 2 , such that the imbalance averaged over realizations is This is the origin of the strong size dependence of the dynamics for large V , such as the saturation of the exponent β to 1 in the inset of Fig. 2(b) for e.g. a system of length N = 10.

V. PARTICLE LOSS
We finish by considering another situation that can occur experimentally, that of particle loss. This corre-sponds to the processes |• → |• , which can be modelled in our effective description by adding to the r.h.s. of Eq.(3) the operator where κ is the loss rate. Loss acts to relax the local imbalance to 0 in a non-collective manner. For κ γJ 2 h −2 the imbalance decays as I ≈ e −κt , and none of the above features survive. In contrast, for κ γJ 2 h −2 only the nucleation-expansion dynamics is significantly modified: decay can act as a nucleation event, and is dominant when κ γJ 2 /V 2 . This affects marginally the value of the compressed β, but results in a saturation of the time scales T , T for large enough V , see insets of Fig. 2(b)-(c).

VI. CONCLUSION
We have considered the effect of interactions on the dynamics of a MBL system subject to dephasing noise. We found two relaxation regimes, one dominated by disorder, and one dominated by interactions. The physical manifestation is a crossover in the decay of time correlators, from stretched to compressed exponential in time. While the stretched exponential regime was expected for weak interactions 21,22 , the crossover to compressed exponential is novel. Our effective classical approximation suggests that this regime would also exist with particle loss. This latter behaviour is due to nucleation and growth dynamics dominating relaxation, a regime which should display strong finite size effects. The above described dynamics should be observable in current experiments on MBL under controlled noise.
Here we obtain analytically the distributions of rates P (Γ) depicted in Fig. 1(b). We consider the rate associated with a fermion hopping in an "interacting" config- uration, namely the ones depicted in the left column of Fig. 1(a). In particular we will focus on the bottom case in the left column of Fig. 1(a), where the interaction V comes with a plus sign in the rate, the other case being trivially deducible by this case. The "non-interacting" cases can be extracted easily from the results below by setting the interactions V = 0. The quantity of interest is the rate for a hopping event involving sites k and k + 1 (which we normalize by h 2 for simplicity), corresponding to where ∆h k = h k+1 − h k . In Eq. (A1) we made the inverse function explicit since it will be used below. Since we focus on a single rate we will drop the site dependency for all the quantities at hand. The rates Γ are random variables, since they depend on the difference of the random field on two contiguous sites. The distribution of the difference ∆h can be easily extracted form the random fields' one since both h k and h k+1 are identically distributed with the same probability between −h and h.
The index s in Eq. (A3) is summed over the inverse functions in the region under consideration. Considering the case where V < 2h, as depicted in Fig. 4 the inverse function is multivalued in the region ∆h ∈ [−2h, 2h − 2V ], so that in Eq. (A3) we have to sum on both the solutions (namely g −1 V (Γ k ) ± ). In the region ∆h ∈ (2h − 2V, 2h] on the other hand the inverse function is single-valued having as single contribution the branch g −1 V (Γ k ) + . The probability density for the rates Eq. (A3) is defined then as In the case V > 2h the function g −1 V (Γ k ) is never multivalued in the region of parameters considered. This results in a probability density function which is defined by the first line of Eq. (A5).
The joint probability P (Γ k , Γ ) used in Eq. (10) is more involved and cumbersome to write in a closed form. We can though express it in a form that allows for direct numerical integration. Calling h k , h k+1 , h the random fields on respectively the k, k + 1, and relaxed sites (see main text for an explanation) we can define This expression can be readily plugged into Eq. (10) giv-ing the results presented in Fig. 1(b) in the main text.