Dynamics of many-body localization in the presence of particle loss

TL;DR

This paper explores how particle loss affects many-body localization in cold-atom systems, revealing that loss combined with interactions suppresses entanglement and alters the system's dynamics, with implications for experimental studies.

Contribution

It introduces a symmetry-based method to analyze particle loss effects on many-body localization, extending understanding beyond dephasing to include particle loss in open quantum systems.

Findings

Particle loss suppresses operator space entanglement entropy.

Interactions combined with loss lead to dissipation and localization suppression.

The method enables long-time analysis of pure and mixed states under loss.

Abstract

At long times residual couplings to the environment become relevant even in the most isolated experiments, creating a crucial difficulty for the study of fundamental aspects of many-body dynamics. A particular example is many-body localization in a cold-atom setting, where incoherent photon scattering introduces both dephasing and particle loss. Whereas dephasing has been studied in detail and is known to destroy localization already on the level of non-interacting particles, the effect of particle loss is less well understood. A difficulty arises due to the `non-local' nature of the loss process, complicating standard numerical tools using matrix product decomposition. Utilizing symmetries of the Lindbladian dynamics, we investigate the particle loss on both the dynamics of observables, as well as the structure of the density matrix and the individual states. We find that particle loss in the presence of interactions leads to dissipation and a strong suppression of the (operator space) entanglement entropy. Our approach allows for the study of the interplay of dephasing and loss for pure and mixed initial states to long times, which is important for future experiments using controlled coupling of the environment.

Figures

Equations

$$\frac{\partial}{\partial t}{\rho}=\mathcal{L}[\rho]=-i[H,\rho]+\sum_{i}\gamma_{i}\Big{(}L^{\phantom{{\dagger}}}_{i}\rho L_{i}^{\dagger}-\frac{1}{2}\{L_{i}^{\dagger}L_{i}^{\phantom{{\dagger}}},\rho\}\Big{)},$$

$$H=-J\sum _i(c_i^{†}{c}_{i+1}^{}+h.c.)+\sum _i{V}_i{n}_i+U\sum _i{n}_i{n}_{i+1}.$$

$$|\psi ⟩=\sum _{i_1,i_2,\dots ,i_L}\text{Tr}(A^{i_1}\cdots A^{i_L})|i_1\cdots i_L⟩,$$

$$\mathcal{O}=\mathbb{I}\times \mathbb{I}\times \cdots {\mathcal{O}}_{i,i+l}\cdots \times \mathbb{I},$$

$$\rho =\sum _{n_1,n_2,\cdots }\sum _{n_1^{\prime },n_2^{\prime },\cdots }{\rho }_{\{n_i^{\prime }\}}^{\{n_i\}}|n_1\cdots ⟩⟨\cdots n_1^{\prime }|.$$

$$c_i|\cdot n_i\cdot ⟩⟨\cdot n_i^{\prime }\cdot |c_i^{†}=(-1{)}^{N_{<i}+N_{<i}^{\prime }}|\cdot {\hat{n}}_i\cdot ⟩⟨\cdot {\hat{n}}_i^{\prime }\cdot |,$$

$$\mathcal{L}[\rho ]={\mathcal{L}}^{(2)}[\rho ]+{\mathcal{L}}^{(1)}[\rho ],$$

$$\mathcal{L}^{(1)}_{i,\pm}[\rho]=\gamma_{i}\Big{(}\pm c^{\phantom{{\dagger}}}_{i}\rho c_{i}^{\dagger}-\frac{1}{2}\{c_{i}^{\dagger}c_{i}^{\phantom{{\dagger}}},\rho\}\Big{)},$$

$$\rho =\sum _{j_1,j_2,\dots }\text{Tr}[{\lambda }^0{\Gamma }^{j_1}{\lambda }^1\cdots {\Gamma }^{j_L}{\lambda }^L]{\sigma }^{j_1}\otimes \cdots \otimes {\sigma }^{j_L}.$$

$$e^{\mathcal{L}t}\approx\Big{(}\Pi e^{\frac{\Delta t}{2N}\mathcal{L}^{(1)}}\Pi e^{\frac{\Delta t}{2N}\mathcal{L}^{(e)}}\Pi e^{\frac{\Delta t}{N}\mathcal{L}^{(o)}}\Pi e^{\frac{\Delta t}{2N}\mathcal{L}^{(e)}}\Pi e^{\frac{\Delta t}{2N}\mathcal{L}^{(1)}}\Big{)}^{N}.$$

$$\frac{d\mathcal{I}}{dt}\approx -{\gamma }_{\mathrm{eff}}(t)\mathcal{I}(t)\approx -\gamma n(t)\mathcal{I}(t)=-\tilde{\gamma }{e}^{-{\gamma }_{l}t}\mathcal{I}(t),$$

$$\mathcal{I}(t)\propto\exp\Big{\{}-\Big{[}\frac{\tilde{\gamma}}{\gamma_{l}}(1-e^{-\gamma_{l}t})\Big{]}^{\beta}\Big{\}}.$$

Full text

Keywords: Many-body localization, Time evolution, Open quantum systems, Matrix product operators, Quantum Simulation, Cold atoms in optical lattices

1 Introduction

Advances in synthetic quantum systems allow for the simulation of paradigmatic model Hamiltonians in highly controlled settings, closely corresponding to ideal closed systems. As a residual coupling to the environment persists even in these systems, it is crucial to address the effects of such non-zero coupling, in particular when studying properties connected to slow intrinsic timescales. Such simulations are well established in the context of quantum-optics [1], where dissipative effects of the environment can be captured in well-controlled microscopic models [2, 3]. Recently, experiments with cold atoms in optical lattices have demonstrated controlled coupling to an environment, through incoherent light scattering and particle loss [4]. Thus, it is important to have a means to benchmark and better understand these processes through numerical simulations. A renewed interest for this problem for the case of many-body systems has focused mostly on bosonic or spin systems, but the efficient simulation of dissipative dynamics for fermionic many-body systems has only started to attract attention [5]. A particular problem arises for fermionic systems due to their statistics, introducing signs whenever particles are exchanged with the environment. When working in an occupation-number basis, we can formally write these signs as the appearance of string operators keeping track of the operator commutations to reach a given site, which renders the formally short-ranged problem long ranged.

Many-body localization presents an example of a model system that has attracted considerable attention in recent years [6, 7]. Isolated many-body-localized systems fail to thermalize even at infinite time scales, providing a novel paradigm contrasting ergodic behavior. Recently, such systems have been studied in trapped-ion settings [8] and with cold atoms in optical lattices [9, 10]. However, in the latter both particle loss and dephasing due to incoherent light scattering are naturally present [4] and destroy localization at long times. Several theoretical works have studied the influence of a dephasing bath on the relaxation dynamics in a many-body-localized system in the Markov approximation [11, 12, 13], which is well suited for typical experimental parameters. The effect of particle loss, however, has only been addressed on a phenomenological level providing insight for special initial states [12].

The time evolution of a system coupled to a Markovian bath is governed by the Lindblad equation for the density matrix

[TABLE]

where $\mathcal{L}$ is referred to as the Lindbladian operator and in the last equality the first term describes the unitary time evolution with Hamiltonian $H$ and the second term the coupling to the environment. The system operators $L_{i}$ are called jump operators and directly couple to the bath with rate $\gamma_{i}$. We will focus on the case provided by spinless fermions on a one-dimensional lattice with random on-site energies and nearest-neighbor interactions,

[TABLE]

Here, $n_{i}=c_{i}^{{\dagger}}c^{\phantom{{\dagger}}}_{i}$, where $c^{\dagger}_{i}$ creates a fermion on site $i$, and $V_{i}\in[-h,h]$ are on-site energies, which are independent random variables. This Hamiltonian has been extensively studied as a model (closed) system and (for $U=2J$) is expected to be localized for $h\gtrsim 7.2J$ [14, 15, 16]. We are then interested in processes that occur naturally in experiments, namely particle loss, i.e., $L_{i}=c_{i}$, as well as dephasing or local density measurement with $L_{i}=n_{i}$ [17].

A powerful approach for the simulation of both closed and open one-dimensional systems is the matrix product decomposition of states and operators [18]. In general, a state can be written as

[TABLE]

where each $A^{i_{n}}$ is a matrix of maximal dimension $\chi\times\chi$ per physical index $i_{n}$. For spinless fermions, $i_{n}$ corresponds to the local occupation, i.e., the Hilbert space is spanned by all $|n_{1}n_{2}\dots\rangle=\Pi_{\{i\}}c^{\dagger}_{i}|\Omega\rangle$, where the product runs over the set $\{i\}$ of occupied sites and $|\Omega\rangle$ represents the vacuum state with no particles. Equivalently, a general operator $\mathcal{O}$ can be represented in matrix product form [19]. For even combinations of single-particle operators with a finite support between sites $i$ and $i+l$, such as given in short-ranged Hamiltonians, these operators can be written as

[TABLE]

where $\mathbb{I}$ denotes a local identity operator and the operator $\mathcal{O}_{i,i+l}$ only acts on the physical indices of sites $i$ to $i+l$. With only local updates required, this allows for the efficient computation of time evolution through time-evolving block decimation (TEBD) [20].

This locality is, however, lost, if we want to annihilate (create) a single particle from a site $i$, since the fermionic statistics requires the introduction of a non-local string of phases (c.f. the Jordan-Wigner transformation). In the following, we show how using an appropriate symmetry, the ‘bond parity’, a matrix-product-operator (MPO) approach for the density matrix allows for the simulation of the relaxation dynamics of an MBL system coupled to a Markovian bath. Simulating the full density matrix simplifies the implementation of mixed initial states and further gives us access to longer times in the presence of a bath. We investigate this latter point by calculating the operator-space entanglement entropy (OSEE), which measures the factorizability of the density matrix [21] and is a measure for the efficiency of the matrix product formalism. While this entropy shows a logarithmic growth for both closed systems and dephasing bath, thereby resembling the entanglement entropy of the pure-state evolution, particle loss leads to a rapid decrease. Finally, we employ a quantum trajectory approach as an alternative view on the entanglement generated during the dynamics and to compare computational efficiency.

2 Methods

2.1 Matrix-Product-Operator Formalism

In general, we can write the density matrix as

[TABLE]

In this basis, the term $L_{i}^{\phantom{{\dagger}}}\rho L_{i}^{\dagger}$ with jump operators $L_{i}=c_{i}$, which describe the loss of particles, will acquire a sign depending on the number of particles to the left of site $i$ both in the bra and the ket element,

[TABLE]

where $N_{<i}=\sum_{j<i}n_{j}$ and the state on the right side is equal to the state on the left side apart from particles at site $i$ missing. The sign of the (local) term thus depends on the total sum of particles up to site $i$, which we will refer to as bond parity $\mathcal{P}_{i}=(-1)^{N_{<i}+N_{<i}^{\prime}}$. This dependence on the particle number renders the operators non-local, and standard TEBD approaches can not be used anymore. Note that for all the other terms in (1), there is no additional sign, as the operators always come in pairs, thus canceling out the strings.

In order to deal with the additional string or sign associated with the loss operators, we first separate the full Lindbladian into bond and single-site terms

[TABLE]

where $\mathcal{L}^{(2)}$ contains the unitary evolution, and $\mathcal{L}^{(1)}$ contains only single-site terms. Only these latter terms can suffer from a relative sign, such that for particle loss, i.e., $L_{i}=c_{i}$, we split it further into

[TABLE]

acting only locally. Importantly, while the loss term changes the number of particles, both single-site and bond terms preserve the bond parity on the left and right bond from where they act and hence, all updates can be performed in the subspaces of $\mathcal{P}_{i}=\pm 1$.

We have implemented the time evolution within the superoperator generalization of TEBD for density matrices [22], where the density matrix (5) is written as

[TABLE]

Here, $\sigma^{j}$ is a basis for the $2\times 2$ matrices, e.g., the Pauli matrices together with the identity $\sigma^{0}$. The physical indices then come with respective (bond parity) quantum numbers $+1$ for the diagonal $\sigma^{0}$ and $\sigma^{3}$ and $-1$ for $\sigma^{1}$ and $\sigma^{2}$. In (9), we follow the notation of Vidal [23], where the singular values are stored in diagonal matrices $\lambda^{i}$. Finally, we separate the bond operators into even and odd bond operators $\mathcal{L}^{(e)}$ and $\mathcal{L}^{(o)}$ and use a Trotter decomposition for the time-evolving block decimation,

[TABLE]

Similarly to other implementations for conserved quantities, we keep track of the parity by enriching each bond index $\alpha$ with a map $\mathcal{P}:\alpha\rightarrow\pm 1$. The bond operators are applied as shown in figure 1(a) for updating the MPO matrices represented by $\Gamma$ and $\lambda$. However, unlike standard TEBD implementations respecting quantum numbers [24], the correct operator has to be chosen for the single-site term in order to preserve the parity subspaces, see figure 1(b). Note that while one could further simplify the calculation by working with occupation quantum numbers for unitary evolution and dephasing, this is no longer true when considering particle loss.

3 Results

We have employed the above formalism to calculate the effect of a general bath including both particle loss and dephasing (with respective coupling rates $\gamma_{l}$ and $\gamma_{d}$) on a system described by (1) and (2). Following previous literature, we calculate the time evolution of an initial occupation imbalance, $\mathcal{I}(t)=[n_{e}(t)-n_{o}(t)]/[n_{e}(t)+n_{o}(t)]$ with $n_{e/o}(t)$ the density of particles on even/odd sites to access the system’s relaxation dynamics. Unless noted otherwise, we start from a pure state given by the perfect density wave, i.e., $|\psi_{0}\rangle=|0101\dots\rangle$ and use numerical parameters $\chi=100$ (bond dimension) and $\Delta t=0.1\;[1/J]$.

3.1 Pure loss

The case describing a bath that induces no dephasing but only loss can be understood in both the weak- and strong-interaction limit from a phenomenological perspective based on integrals of motion [12], and thus provides a good benchmark to study first. In the absence of interactions loss does not change the normalized imbalance since the occupancy on even and odd sites follows the same exponential decay. However, upon introducing a nearest-neighbor interaction $U\neq 0$ the imbalance is expected to decrease. This is due to the change of the local potential a particle experiences as a neighboring particle is removed, which may enable hopping. Only after the particle density has been reduced sufficiently does the imbalance stabilize as the system becomes effectively non-interacting.

This behavior can be qualitatively described in a mean-field picture, where the time evolution of the imbalance can be approximated by [12, 13]

[TABLE]

where $n(t)=n_{0}\exp(-\gamma_{l}t)$ is the particle density and $\tilde{\gamma}=\gamma n_{0}$ is an effective decay rate. Integrating (11) yields the time dependence for the imbalance

[TABLE]

Here, we have introduced an exponent $\beta$ phenomenologically to account for the fact that the mean-field picture is a simplification [12] and to yield better agreement with the numerical results.

Figure 2(a) shows the imbalance relaxation due to loss for a system of length $N=20$ sites, $h=10J$, and various interactions. As expected, the imbalance relaxes initially with a rate increasing with the interaction strength before saturating once the density has dropped sufficiently (see inset). For small interactions $U\ll h$, the initial decay can be calculated perturbatively and is expected to scale as $\gamma_{l}(JU/h^{2})^{2}$ [12]. For large interactions, on the other hand, the initial decay becomes independent of the interaction strength. This dependence of the effective decay rate $\tilde{\gamma}$ on the interaction strength is shown in figure 2(b). Equation (12) yields a robust fit to the numerical time traces with $\beta\approx 0.9$ for all interactions $U\geq 10J$ and allows for a direct extraction of the decay rate. For comparison, the error bars denote the difference when fitting to a fixed $\beta=1$. Unfortunately, the fitting does not work as robustly for smaller interactions and instead, we confirm the initial decay rate $\tilde{\gamma}\sim U^{2}$ via a data collapse of the time traces shown in the inset of figure 2(b).

3.2 Dephasing and loss: weak interactions

Having established the dynamics of the system when only particle loss is present, we next address the interplay of both dephasing and loss. As can be seen in figure 3 the inclusion of dephasing processes for weak to intermediate interactions has a detrimental effect on the localization. In this regime, the dephasing processes directly lead to dissipation, and additional particle loss provides only a subleading effect. Regardless, the relative effect of loss is largest close to the critical disorder strength and decreases for stronger disorder. This is in agreement with recent experiments with spinful fermions [4]. The inset of figure 3 emphasizes this trend by comparing the effective initial decay rate for dephasing only with dephasing and loss combined, extracted using a stretched exponential decay for both time traces for various disorder strengths.

3.3 Dephasing and loss: strong interactions

As the interaction strength increases, the initial density-wave state at half filling has more and more overlap with an eigenstate of the system. In such a situation, the dephasing-induced hopping of particles is suppressed as $1/U^{2}$ and dephasing has vanishing effect. Hence, the particle-loss term is expected to dominate, which is clearly shown in figure 4 for $U=50J\gg h$ comparing the dynamics of loss only, dephasing only, and a combination of the two. Only after sufficient particles have been lost does the dephasing become relevant, leading to additional diffusion.

For a more general investigation of the dynamics at large interactions, it is thus necessary to not only focus on the perfect density wave with $n=0.5$ as the initial state, but also on density-wave states with lower overall density. The matrix product operator approach is well suited for implementing these states, and figure 4 additionally shows the imbalance decay for pure dephasing, pure loss, and the combination for additional densities $n=0.3,0.4$. For these densities, the initial states have less overlap with a single eigenstate and the imbalance drops to a lower value at short times correspondingly. In addition, the influence of dephasing increases since for the initial state not every second site is occupied and thus the effect of the nearest-neighbor interaction is decreased. In the limit of vanishing initial density, we would thus recover the non-interacting result, where loss has no effect on the imbalance dynamics.

3.4 Evolution of entanglement entropy

To investigate the efficiency of our approach we consider the density matrix’s entanglement with respect to the middle bond of the state with $n=0.5$. The OSEE provided by the MPO representation of the density matrix (9), i.e., $S_{\sharp}=-2\sum_{\alpha}(\lambda^{i}_{\alpha})^{2}\log\lambda^{i}_{\alpha}$ for bond $i=N/2$ is a measure of the factorizability of the density matrix and thus its unbound growth signals the breakdown of the MPO approach.

Figure 5 compares the OSEE evolution for various baths starting from the perfect density wave, a product state. For $U=50J$, the OSEE stays close to zero until particle loss removes the system from being close to an eigenstate, confirming the observations of the previous section for $n=0.5$. For $U=2J$, the closed system OSEE shows a $\log t$ growth, in accordance with the (pure state) entanglement entropy, which limits the accessible times to $t\approx 100/J$. Pure dephasing leads to a decrease of the OSEE at intermediate times, but a characteristic $\log t$ growth sets in at long times [11]. With the inclusion of loss, however, the OSEE quickly drops to zero after times $t\sim 1/\gamma_{l}$ irrespective of interaction strength, since at this timescale the loss of particles becomes significant and effectively prevents long-range correlations. Note, however, that there is still a significant number of particles left in the system at this time.

Finally, we compare the evolution of the OSEE with the entanglement entropy obtained via a quantum trajectory method [2, 25] as shown with dashed lines in figure 5. For comparison, we calculated the entanglement entropy using a logarithm with base equal to the local Hilbert space dimension for both the quantum trajectory method ($d=2$) and the density matrix evolution ($d=4$). While the overall behavior of the two entropies is similar, the entropy from the trajectory method shows both a slower growth and a longer decay. This difference can be understood as follows: First, for a trajectory the local density measurement (representing the dephasing) does not affect the entanglement by resetting it to a product state. This behaviour is demonstrated by the blue dashed line in figure 5. Second, the probabilistic nature of the quantum trajectory approach includes trajectories where jumps, i.e., coupling to the bath, occur very early, or not until late in the time evolution. This resembles more closely the actual experimental situation, where only averaging over many runs results in the statistical behavior described by the density matrix. However, for disordered systems, where the initial entropy growth is slow, our MPO approach provides a more efficient method for the calculation of the dynamics induced by a Markovian bath.

4 Conclusion

We have implemented the Lindbladian time evolution for a fermionic system coupled to a bath through dephasing and loss employing a matrix-product-operator formalism. For this purpose, we have identified the ‘bond parity’ as a symmetry of the Lindbladian evolution and implemented this symmetry in a quantum-number-conserving TEBD calculation. Our approach allows for the efficient calculation of the time evolution of any (mixed) initial state in the presence of a bath to long times and thus for the comparison to experimental setups that have recently become available.

As a specific example of interest, we have used our approach to investigate the effects of loss on the dynamics of an otherwise many-body-localized system. We have confirmed the behavior of an initial charge imbalance expected from a phenomenological model of pure loss in the limits of weak and strong interaction. Importantly, we have shown how dephasing and loss have different significance for inducing diffusion at weak compared to strong coupling, as well as close to the localization transition compared to deep in the localized phase.

The interplay between interactions, loss, and dephasing that we study here could be realised in ongoing experiments. Recent advances in cooling and observing atoms that have a large intrinsic dipole-dipole interaction such as Erbium allow for the realisation of lattice models with significant off-site interactions [26]. Beyond these specific systems, our matrix-product-operator implementation that retains parity information provides a new general tool for simulating particle loss and other dissipative processes where we have to account for such signs, significantly expanding the available possibilities for efficient simulation of dynamics induced by a Markovian bath.

EvN gratefully acknowledges funding by the Swiss National Science Foundation through grant P2EZP2 172185. Work at Strathclyde is supported in part by the US Air Force Office of Scientific Research grant number FA2386-14-1-5003.