Dynamics of Negativity of a Wannier-Stark Many-Body Localized System Coupled to a Bath

An interacting system subjected to a strong linear potential can host a many-body localized (MBL) phase when being slightly perturbed. This so-called Wannier-Stark or `tilted-field' MBL phase inherits many properties from the well-investigated disordered MBL phase, and provides an alternative route to experimentally engineer interacting localized systems without quenched disorder. In this work, we investigate the dynamics of entanglement in a Wannier-Stark MBL system coupled to a dephasing environment. As an accessible entanglement proxy, we use the third R\'{e}nyi negativity $R_3$, which reduces to the third R\'{e}nyi entropy in case the system is isolated from the environment. This measure captures the characteristic logarithmic growth of interacting localized phases in the intermediate-time regime, where the effects of the coupling to the environment are not yet dominating the dynamics. Thus it forms a tool to distinguish Wannier-Stark MBL from non-interacting Wannier-Stark localization up to intermediate time-scales, and to quantify quantum correlations in mixed-state dynamics.


I. INTRODUCTION
Over the past decade there has been a huge research interest in disordered interacting many-body systems. It was realized that these systems could host a manybody localized (MBL) phase provided that the disorder is strong enough [1][2][3][4]. This MBL phase is robustly non-thermalizing, and should be contrasted to a single-particle Anderson localized phase occurring in non-interacting systems [5]. While both phases are characterized by the absence of transport, there are also notable differences, most prominently in the entanglement dynamics under a quench. In the Anderson localized system, quantum correlations cannot propagate through the system. Hence the entanglement saturates after a fast ballistic initial growth resulting from local rearrangements of particles. In the interacting MBL system, on the contrary, quantum correlations can propagate in the system but that happens only logarithmically in time [6][7][8]. This slow growth can be understood in terms of effective exponentially decaying interactions between socalled 'localized integrals of motion' (LIOMs), that form a phenomenological picture to understand MBL [9,10].
The existence of MBL has been experimentally confirmed [11][12][13][14], and the logarithmic spread of quantum correlations has been observed as well [15][16][17][18]. Recently, there have been many proposals to establish disorderfree types of localization [19][20][21][22][23][24][25][26][27][28][29], including for instance lattice gauge theories [26], or mixtures of two types of particles where the light ones are localized on the heavy ones [21-23, 30, 31]. In particular, it has been realized recently that many features of MBL are also inherited by interacting systems subjected to a strong linear potential [28,29,[32][33][34][35], yet also differences have been identified and understood in terms on Hilbert-space fragmentation [36,37]. In the non-interacting case this phenomenon is referred to as Wannier-Stark localization [38,39]; in the interacting case it is referenced to as Wannier-Stark MBL (or shortly Stark MBL in the literature). This type of localization has the advantage that it can potentially be induced solely by the tuning of an external electric field, without the need of engineering internal properties in the system. Experimental signatures of non-ergodic dynamics in systems subjected to a tilted field have been observed in Refs. [40][41][42], while Ref. [43] was the first experiment that investigated the effect of a tilted potential on the approach to equilibrium.
Typical experiments are never fully isolated from the surrounding environment. This implies that it is hard to distinguish interacting types of localization from noninteracting types of localization. The most prominent difference between the two types is the logarithmic spreading of quantum correlations in the interacting (and isolated) case after a quench. In this work we therefore focus on the entanglement dynamics of interacting and non-interacting Wannier-Stark systems that are coupled to a dephasing Markovian environment. This extends our previous work [44] where we considered a similar setup for disorder-induced MBL. Related setups have also been considered to mostly investigate transport properties in Refs. [44][45][46][47][48][49][50] in the context of disorder-induced MBL, and in Ref. [51] in the context of Wannier-Stark MBL coupled to a dephasing environment.
Measuring quantum correlations in open systems is challenging both theoretically as well as experimentally. From the theoretical side, it is hard to find generically computable measures that do not rely on a full diagonalization of the (partially transposed) density matrix. We circumvent this problem by considering moments of the partially transposed density matrix and calculate the third Rényi negativity R 3 [44,[52][53][54][55][56], which is not an exact entanglement monotone but which captures the relevant dynamics as it behaves quantitatively similar as the negativity [44,54]. From the experimental side, it is challenging to measure non-local correlations as full-state tomography is exponentially expensive in the system size, and as joint measurements on multiple copies of the state are also hard to engineer [54,[57][58][59][60]. Recently there has been a lot of progress in measuring Rényi entropies in synthetic quantum matter by random unitary measurements [15,[61][62][63][64][65], and this toolbox naturally includes the measurement of mixed state entanglement via Rényi negativities [56].

II. MODEL AND SETUP
We consider the XXZ Hamiltonian with an on-site potential where S x,y,z i are the spin-1 2 operators. As on-site potential we take a linear potential with small quadratic corrections of the form to induce Wannier-Stark localization. Here γ gives the slope of the linear potential, and α describes small quadratic corrections to it. The parameter ∆ describes the strength of the nearestneighbor interactions which can be seen by writing the Hamiltonian (1) in terms of spinless fermions by means of the Jordan-Wigner transformation. When there are no interactions present, ∆ = 0, the system is 'single-particle' Wannier-Stark localized. In the presence of interactions, it is Wannier-Stark many-body localized. We will focus on the weakly interacting regime where ∆J < γ as to avoid resonance effects where e.g. a particle could move up the linear potential without energy cost.
Localization also occurs in the same model (1) for randomly disordered on-site potentials. In the interacting case, there is an MBL phase if the disorder is strong enough, i.e. higher than a critical disorder strength under which the system is thermalizing.
The Wannier-Stark model exhibits many similar properties as the model with disorder, e.g. most importantly a slow logarithmic growth of entanglement under a quench, if the linear field γ is sufficiently strong, and if there is some non-uniformity to this linearity [28,29,33,35]. When only a linear field is applied, the (non-interacting) system contains a lot of degeneracies, and therefore properties like the level spacing statistics, can deviate from the ones of the disorder-induced localized systems. By adding a quadratic gradient most of these degeneracies are resolved, and the level-spacing statistics of the two cases are very similar [28]. We couple the system to a simple, yet realistic [13], Markovian dephasing environment which is modeled by the jump operators L i = √ ΓS z i , such that the timeevolution of the state is described by the following Lindblad master equation [66] Note that this type of dynamics also conserves the U (1) symmetry of the Hamiltonian leading to a fixed total magnetization M = i S z i . In the limit of a purely dephasive coupling, the off-diagonal matrix elements of the density matrix decay exponentially with decay rate Γ. Hence it is clear that this type of environment coupling destroys all quantum correlations as it drives the system into a trivial infinite temperature state, that is the unique steady state of Equation (3). A sketch of our setup is shown in Figure 1.
For simulating the time evolution according to the Lindblad master Equation (3), we use the time-evolving block decimation (TEBD) algorithm on density matrices [67,68] where the (non-unitary) time-evolution operator is given by U (t) = exp(Lt) with L the Lindblad superoperator This time-evolution operator acts on a (vectorized version) of the density matrix |ρ as |ρ(t + dt = exp(Ldt) |ρ(t) . The matrix-product decomposition of |ρ reads with χ the maximal bond dimension we allow for in the simulation. Hence the so-called operator entanglement entropy of |ρ is limited to log χ during the evolution towards the steady state. Because of the localized nature of the isolated system and the addition of dephasing noise, the operator entanglement entropy stays moderate at all times in the considered setups. For more details about the numerical implementation, we refer to our previous work Ref. [44].
We study a global quench setup by starting from an initial (pure) product state ρ 0 = |ψ 0 ψ 0 | where |ψ 0 is a product state in the zero magnetization sector. Our goal is to measure the dynamics quantum of quantum correlations in such a setup. However, in a mixed state manybody context it is very challenging to separate quantum correlations from classical correlations in a generic way. There exist experimentally accessible quantities like the quantum Fisher information [12,[69][70][71], however these so-called entanglement witnesses rely on system specific properties and are hence not as generic.  3 . We bipartition the system into two equal subsystems A and B, and take the transpose of the degrees of freedom of partition B. After which we take the trace of the system (green lines).
A famous generic criterion is the so-called positive partial transpose (PPT) criterion of Peres [72]. It states that if a bipartite density matrix ρ AB is separable, then ρ T B AB has only positive eigenvalues, where T B denotes the operation of partial transposition. The negativity is then a generic entanglement monotone that quantifies the violation of this criterion [73] however, because of the trace norm, this quantity relies on the full diagonalization of ρ T B AB which is exponentially hard in the system size. Therefore the negativity is not computable for large systems.
However, as detailed in Ref. [44,54,55] we can use the third Rényi negativity to estimate the amount of bipartite entanglement in the system Note that the equivalent quantities using the first and second power of the density matrix are trivially zero as Tr ρ T B AB = Tr ρ AB = 1 and Tr ρ T B AB 2 = Tr ρ 2 AB . The quantity R 3 is however not an entanglement monotone as the negativity (for the definition of an entanglement monotone see e.g. Ref. [74]) as it is always zero for simple product states, but not necessarily for separable ('classically' correlated) states. However, in case of disordered MBL we have shown in Ref. [44] that it has the same dynamical properties as the negativity. This is manifested for instance by very similar stretching exponents. The contraction of the network to compute Tr (ρ T B AB ) 3 is sketched in Figure 2.
For a pure state, the density matrix has the form ρ AB = |ψ AB ψ AB |, hence we have that and thus reduces to the third Rényi entropy [53]. Here, ρ B = Tr A ρ AB denotes the reduced density matrix of subsystem B.
In what follows we will always assume that we have two equal partitions in the system and shortly denote R 3 (ρ AB ) ≡ R 3 .

III. RESULTS
In this section, we discuss our results. In the closed case, we first look at the entanglement dynamics by making a quench from a Néel state, and secondly at the dynamics by making a quench from a random product state in the zero magnetization sector and average over the results. Afterwards, we discuss the results in the open case Γ > 0.

A. Isolated system
We start by investigating the dynamics of R 3 in the isolated system Γ = 0. In this case we have that R 3 is directly related to the third Rényi entropy as discussed in the previous section.

Quench from a Néel state
When we start initially from a Néel state |↓↑↓↑ . . . , we can track how fast the Wannier-Stark MBL system 'loses' information about this initial state pattern. This can be realized by considering a quantity like the imbalance where N z e /N z o is summing over the occupation numbers N z i = S z i +1/2 for even/odd sites. For localized systems I decays to a highly non-thermal (i.e. non-zero) value, and this forms a simple and accessible experimental probe for localization [11]. However, a priori the imbalance decay does not directly allow one to distinguish interacting types of localization from non-interacting types of localization. This is shown in Figure 3, where we show both the imbalance and entanglement dynamics in the interacting (∆ = 0.5) and non-interacting (∆ = 0) cases for various slopes of the linear potential. From this it is clear that only the entanglement dynamics provides a striking difference between the two cases.

Quench from a random initial state
We now consider a random product state in the zero magnetization sector and average the results over the different realizations. This can be seen in Figure 4(a), where we show the logarithmic growth of R 3 for different interaction strengths at field parameters α = γ = 2J. Without interactions ∆ = 0 the system is single-particle Wannier-Stark localized and entanglement cannot propagate through the system. Fluctuations are stronger in the non-interacting limit, as the level spacing is only inversely proportional to the system size. The time-scale at which the effect of the interactions becomes dominant is set by t int ∼ (∆J) −1 . In the data obtained for finite interactions, a cross-over time-scale t cross becomes apparent, which is absent in the case of disorder-induced MBL, where there is a faster logarithmic growth, after the initial ballistic growth up to times tJ ∼ 1. Beyond this cross-over regime, there is then a slower, logarithmic growth. We attribute the existence of this cross-over regime to the quadratic contribution of the potential of strength α in (2). When increasing the quadratic deviations, the cross-over regime becomes less pronounced, as we show in Figure 4(b).
In Figure 4(c) we show the logarithmic growth at one particular interaction strength ∆ = 0.5 for various system sizes. Also, t cross is doubled when the system size is doubled which confirms that the parameter α of our local potential (2) is indeed governing t cross . When increasing system size, α thus has to be rescaled accordingly when considering the in the literature commonly used form of the local potential [28].

B. Open System
We now turn to the investigation of the dynamics of R 3 in the open system Γ > 0. When the system is coupled to a dephasing enviroment, all entanglement structure will be eventually lost as the system heats up to the infinite temperature state. The time scale on which the depasing starts to dominate the dynamics is set by the coupling strength t deph ∼ 1/Γ. In order to allow that the interactions can still dominate the dynamics at intermediate times, we need to have that t int t t deph which implies that we must have that Γ/J ∆. Hence, the dephasing strength must be sufficiently weak compared to the interaction strength to be able to still observe signatures of the logarithmic entanglement growth. In Figure 5 we show the entanglement dynamics for various coupling strengths. One of the main advantages of looking at the Rényi negativity R 3 is that we can capture the logarithmic growth of quantum correlations even if the system becomes slightly mixed. This onset of logarithmic growth can indeed be still observed for sufficiently weak dephasing strengths in Figure 5(a). In principle, for characterizing interacting dynamics it is sufficient that the entanglement reaches a value that is higher than the maximal value of the oscillations in the non-interacting case. The decay of R 3 in the non-interacting case is shown in Figure 5(b) for comparison.
From disordered MBL it is known that the tails of the decay of the imbalance [46,47] and the negativity or R 3 [44] are stretched exponentials ∼ e −(Γt/a) b with b < 1. These stretched exponentials are in that case understood as a superposition of many local exponential decays [46], originating from very broad distributions of exponentially decaying couplings in the phenomenological LIOM picture of MBL [9,10,75]. However, the Wannier-Stark case is disorder free, which could therefore lead to a different behavior of the tails of the decay. In the non-interacting case, such a difference has been reported in the decay of the imbalance in Ref. [51].
In Figures 6 and 7, we show the decay of respectively R 3 and the imbalance under a quench from the Néel state with Γ = 0.1. These figures have a double logarithmic scale on the y-axis and single logarithmic scale on the xaxis, such that a stretched exponential ∼ e −(Γt/a) b would appear as a straight line with b the slope, and a related to the offset. We have fitted stretched exponentials to the data in the intermediate-time regime. However, in the late-time regime, the functional form of the decay changes. For small system sizes, exponential decay can be observed at late times as we illustrate in the Appendix A. In this Appendix we also demonstrate that averaging over initial states does not lead to a quantitatively different behavior in the decay of R 3 .
In Ref. [51] they report that the decay of the imbalance happens according to an exponential in the noninteracting case, while in the interacting case it happens according to a stretched exponential. They however consider much stronger tilts and interactions. As we focus here on a more moderate regime of weaker tilts and weak interactions as relevant for experiments [33], we can not distinguish qualitative differences between both cases in the late-time dynamics. We averaged over about 100 initial states in the M = 0 sector for Γ > 0, and over about 300 states for Γ = 0 (exact diagonalization was used in this case).

IV. CONCLUSION
In this work we have investigated the dynamics of entanglement in a non-interacting Wannier-Stark localized and in an interacting Wannier-Stark many-body localized system. In the closed case, we have observed a cross-over regime that is absent in the case of disordered MBL and that is related to the quadratic corrections from linearity of the field. In the open system with dephasing noise, it is possible to still observe parts of the characteristic logarithmic growth for sufficiently weak dephasing strengths.
Our results confirm that the Wannier-Stark MBL system indeed inherits many properties from the disordered MBL system in particular when considering the entanglement dynamics coupled to an environment. However, we do not observe a robust stretched-exponential functional form of the entanglement decay in the Wannier-Stark case.
For future work, it would be interesting to investigate the entanglement dynamics under different types of dissipation. In particular for non-hermitean Lindblad operators, it could be interesting to investigate whether a physically relevant entanglement structure could still exist in the non-trivial steady state. In this Appendix, we present some exact diagonalization data, obtained by numerically integrating the Lindblad Equation (3) for a small system of L = 8 sites. In Figure 8 we show the collapse of the tails for various coupling strengths Γ for a quench starting from the Néel state. We have fitted stretched exponentials ∼ e −(Γt/a) b in the intermediate-time regime, and exponentials ∼ αe −Γt/β in the late-time regime. In Figure 10 we show the same data in a different scale that makes the exponential form of the late-time tails clearly visible. In Figures 9 and 11, we show the data when we average over different initial states. From this, it is clear that the averaging does not alter the behavior of the tails quantitatively, e.g. the best-fit stretching exponents in the intermediate-time regime stay very similar as can be seen by comparing Figures 8 and 9.     11. The same data as in Figure 9 but in a different scale to make the exponential form of the tails visible.