Thermalization of strongly interacting bosons after spontaneous emissions in optical lattices

We study the out-of-equilibrium dynamics of bosonic atoms in a 1D optical lattice, after the ground-state is excited by a single spontaneous emission event, i.e. after an absorption and re-emission of a lattice photon. This is an important fundamental source of decoherence for current experiments, and understanding the resulting dynamics and changes in the many-body state is important for controlling heating in quantum simulators. Previously it was found that in the superfluid regime, simple observables relax to values that can be described by a thermal distribution on experimental time-scales, and that this breaks down for strong interactions (in the Mott insulator regime). Here we expand on this result, investigating the relaxation of the momentum distribution as a function of time, and discussing the relationship to eigenstate thermalization. For the strongly interacting limit, we provide an analytical analysis for the behavior of the system, based on an effective low-energy Hamiltonian in which the dynamics can be understood based on correlated doublon-holon pairs.


I. INTRODUCTION
An important fundamental question in many-body quantum mechanics is to what extent and under which conditions an isolated system perturbed away from equilibrium will undergo thermalization, in the sense that at long times the system will reach a steady state where simple observables equal the values for a thermal distribution [1][2][3][4][5]. Recently it has been possible to observe integrable dynamics for strongly interacting cold gases confined to move in one dimension [6], and thus reach regimes where systems either do not thermalize in this standard sense [7] or otherwise undergo generalized thermalization [8,9].
These fundamental questions also have an important impact on the application of cold atoms as quantum simulators [10,11]. Development in experiments with these systems has reached a stage where it is possible to tailor interesting many-body Hamiltonians, and study the many-body physics of corresponding ground states. However, some of the most interesting physical regimes require the realization of states with small energy gaps, requiring exceptionally low temperatures and entropies, which provide a key challenge for current experiments [12]. This has been particularly true in attempts to observe quantum magnetism in optical lattices within strongly interacting regimes |U | J, where U is the on-site interaction energy, and J/ indicates the tunneling amplitude for particles moving between neighboring sites. There, magnetic order in multi-species mixtures is driven by terms of the order of J 2 /U , which is typically small for current experiments [12][13][14][15][16][17][18][19].
In this context, it is very important to understand and control the many-body effects of competing heating processes in experiments. Previously, it has often been as-sumed that all of the energy added to the system will be thermalized, causing an effective increase in temperature. However, in regimes in which the system does not thermalize excitations entirely (or for some types of excitations, at all), the behavior is more complex. A simple example where excitations cannot thermalize on typical experimental timescales is given by spontaneous emission events (incoherent light scattering). In such processes, the main contribution to the increase in the average energy of the system as a function of time comes from atoms being excited to higher bands of an optical lattice [20][21][22][23]. Because the bandgap is usually much larger than the typical energy scales U and J of dynamics in the lowest band, it is not possible for the system to thermalize most of the energy added to the system as a function of time.
At the same time, understanding this process is made more complex by the fact that processes exciting particles to higher bands are rare compared with processes leaving particles in the lowest band in typical experiments [20,21]. As a result, the dynamics can actually be dominated by heating in the lowest band of the lattice, the basic effects of which have been studied in several recent articles for bosons [20,21,[24][25][26][27][28] and for fermions [29][30][31]. It is then natural to ask, in particular, whether the system can thermalize after spontaneous emission events that leave atoms in the lowest band.
In this article, we study the interplay between spontaneous emissions and thermalization in detail, asking under what conditions key simple quantities such as the quasimomentum distribution and the kinetic energy relax to values given by a thermal distribution. In each case, we ask this question in a practical context by fixing a thermal distribution such that we take a canonical ensemble with the temperature T chosen so that the ex-pectation value of the energy in the thermal distribution matches the expectation value of the energy after a spontaneous emission event. In the case of strong interactions, we might expect that the system reaches integrable limits where the system does not thermalize. On the other hand, one might also expect that as in general, the Bose-Hubbard Hamiltonian is non-integrable, that outside the special cases of strong interactions (U/J → ∞) and noninteracting systems (U = 0), thermalization would essentially always occur.
Instead, we find that this generalization is not universally applicable. It is important to note that spontaneous emissions correspond to local quenches for the many-body state, leading to population of low-lying excited states. As we showed previously [24], after a spontaneous emission, depending on the parameter regime and the corresponding different characteristics of the lowenergy spectrum: either (i) the system relaxes over short times to thermal values of the quasimomentum distribution and kinetic energy, or (ii) on short timescales, the system relaxes to states that are clearly non-thermal, even if all atoms remain in the lowest Bloch band. The latter regime was found to occur if the interactions are strong enough so that the ground-state of the system corresponds to a Mott insulator state. Below we elaborate on these results especially for the quasimomentum distribution, and discuss their interpretation in terms of eigenstate thermalization.
In the limit where the system is close to an ideal Mottinsulator (a product state with a fixed number of particles per site) and in the case of unit filing, one can describe the dynamics of the system in an effective model, where one restricts the local Hilbert space to having a maximum of two particles per site. The low-lying excitations can then be described as correlated pairs of holes ("holons") and doubly occupied sites ("doublons") [26,32]. Here we use this approximation to show that no thermalization is present for spontaneous emission in this limit. We compare our result to exact numerical calculations either using exact diagonalization in small systems, or by making use of t-DMRG techniques [33][34][35][36].
This article is organized as follows. We begin by introducing spontaneous emissions and thermalization in these systems in Section II. Then, in Section III we review the explanation of the basic behavior we find in terms of the eigenstate thermalization. In Section IV we discuss how the behavior can be understood in the strongly interacting limit, making predictions that could be observed for the propagation of excitations when the system does not thermalize in that limit. In Section V we then provide a summary and outlook.

II. THERMALIZATION OF SPONTANEOUS EMISSIONS
Spontaneous emissions are a fundamental source of heating in optical dipole potentials [22,23], and one of the key contributing heating sources in current experiments with cold atoms in optical lattices [20,21]. As shown in Ref. [21], a multi-band many-body master equation can be derived to describe this process in the case of far detuned lattice light. An adiabatic elimination of the electronic excited states leads to an equation for atoms in their electronic ground states. Taking this, and assuming the typical experimental case in which the lattice spacing a is comparable to the wavelength of scattered photons, a = λ/2, the dynamics of the many-body density operator ρ follows a master equation of the form Here H MB is a multi-band Bose-Hubbard Hamiltonian and Lρ includes dissipative Lindblad terms with on-site jump operators ∝ (b is the bosonic annihilation operator at site i, for a particle in band m. For relatively deep optical lattices, transition rates for inter-band processes coupling neighboring Bloch bands are strongly suppressed by the square of the Lamb-Dicke parameter, η = 2πa T /λ. For a typical experiment with a lattice depth V 0 ≈ 8E R [E R = 4π 2 2 /(2mλ 2 ), where m is the mass of the atom], η 2 ∼ 0.1. In the usual case of a red-detuned optical lattice the dominant processes are thus intra-band processes, which return the atoms to their initial Bloch band. Provided we initially consider the system to be in the lowest band, processes accessing higher bands are suppressed by a factor of the order η 4 or greater, and the master equation for the 1D model simplifies to:ρ where n i = b † i b i is the particle number operator for the lowest band, and γ is the scattering rate. The Hamiltonian is the usual single band Bose-Hubbard Hamiltonian In the following we will focus on a system with M sites and N = M particles, so that the filling isn = N/M = 1.
In this case the system undergoes a ground-state quantum phase transition between a superfluid (SF) and a Mott insulating (MI) state at the critical value, which was estimated to be U c ≈ 3.25(5)J [37,38]. The underlying physics of the lattice-photon scattering process can be summarized as the environment obtaining information about the position of an atom. Thereby the atoms are localized on the length-scale of a single lattice-site size. Numerically we can fully simulate the master equation (1) also for large 1D systems by combining t-DMRG techniques with a Monte-Carlo sampling over quantum trajectories [39][40][41][42]. There, the density matrix is obtained from a statistical average over many trajectories with randomly applied jump operators [39]. To give an example, in Fig. 1 we plot the dynamics for the evolution under equation (1). We start in the ground-state for a system with U = 2J (SF regime, M = 48). Then we turn on light scattering for a short time up to tJ = 1 and for different values of γ (illustrated in Fig. 1a). This leads to an increase in energy (Fig. 1b). Once the heating is switched off, the system then relaxes. This relaxation is for example seen when looking at the evolution of the quasi-momentum distribution During the heating period the localization of particles leads to a "lifting" of tails of this distribution, as seen in Fig. 1c. Once the heating is switched off, the high quasi-momentum components become redistributed and the system is relaxing to a broadened distribution, depicted in Fig. 1d. The question whether the system thermalizes or not depends now on the shape of this distribution. In the case that the steady state distribution coincides with a Boltzmann distribution for a temperature corresponding to the same mean energy, we call the system thermalized, otherwise we don't. Obviously, the question of thermalization is connected to the observable under consideration. In Ref. [24] it was found that for example the tails of the quasi-momentum and the kinetic energy thermalize on experimentally relevant time-scales in the SF regime and that this thermalization breaks down in the MI.
Here, we demonstrate a specific example of this in Figs. 1e,f. There we show relaxation of a high-q component of the quasimomentum distribution, q = 40π/48, in a system with M = 48 lattice sites. We can clearly see that after a short period of time, this component relaxes to its value in a corresponding thermal distribution provided that we initially start in a superfluid state with U/J = 2. In the MI limit, we see that the behavior is similar, in so far as there is a rapid relaxation of this value of quasimomentum on short times. However, this leads to a value that is markedly different to the corresponding thermal value. It is naturally possible that in this regime, there are processes on very long timescales that will eventually lead to thermalization of this quantity. But it is clear that on the MI side of the phase transition, the system exhibits a qualitatively different behavior in terms of how close it is to thermal distributions on the scale of a few tunnelling times. This will significantly impact the ability of the system to thermalize energy added by heating on the timescales of the heating process itself.
It is important to note that in the regime where the system is relaxing, the thermalization and all of the corresponding dynamics is determined by unitary closed system dynamics alone. In the remainder of this paper we will simplify the question of thermalization to the coherent dynamics after a single spontaneous emission event, i.e. we ask whether expectation values of observables relax to thermal values in a unitary time-dynamics after we apply a single quantum jump. Thus, we consider a unitary time-evolution of an initial state |ψ l (t = 0 + ) = n l |ψ g ||n l |ψ g || , where |ψ g denotes the ground-state of the system, and then take a stochastic average over the sites l, which is weighted proportional to the square of the density on site l, ψ g | n 2 l |ψ g . The mechanism of this thermalization of closed quantum systems can be explained via eigenstate thermalizaton, a mechanism which we will review in the next section.

III. EIGENSTATE THERMALIZATION
The question whether unitary time-dynamics leads to thermalization or not can be re-phrased in terms of the "eigenstate thermalization" [1][2][3], which we will briefly review here. Consider the situation where the initial state is not an eigenstate of the Hamiltonian (such as it is the case after a spontaneous emission). The state can be expanded into energy eigenstates defined via H |α = E α |α , |ψ = α c α |α , and the time-evolution is thus given by Then the time-dependent expectation value of any ob-servableÔ is If the system relaxes into a steady state, this steady state must be identical to the long-time average Thus, the value of the steady state expectation (if such a steady state is developed) only depends on two variables: i) The expansion coefficients of the initial state, c α , and ii) the expectation values of the energy-eigenstates with the particular observable, α|Ô |α . The eigenstate thermalization hypothesis states, that if the the eigenenergy expectation values of an observable vary smoothly in the energy window defined by the c α (and the off-diagonal α|Ô |β , α = β are small), the system relaxes into a state for which the observable can be described with a micro-canonical ensemble [1,2]. In our situation, spontaneous emission only adds a small amount of energy to the system. Thus, the question whether this energy thermalizes or not will be determined by the low-energy spectrum. Thus, we can use an exact numerical diagonalization to obtain the lowest energy eigenstates in systems with N = M = 10 sites and particles. As an example, in panels (a)-(c) of Fig. 2 we show the matrix elements α|Ô|β for the kinetic energy in the system. As observable we use the kinetic energyÔ = H kin . The eigenvalues are ordered ascending according to their energy. We find that only in the SF case of U = 2J, the diagonal elements vary smootly as function, while the off-diagonal elements are zero. This verifies the findings of Ref. [24], and is consistent with the breakdown of the low-energy thermalization when entering the MI regime.
In addition we plot a direct comparison of the two expectation values in Fig. 2d. We can calculate the expectation value on the one hand from the diagonal ensemble and on the other hand from a canonical ensemble with the same mean-energy. As initial state we choose the state after a single spontaneous emission on site i = M/2. Although finite size effects play an important role in this small system, it is already obvious that a discrepancy between the two ensembles arises when entering the MI regime. As was found by t-DMRG calculations in [24], this also holds for larger systems. In the next section we will calculate the low-lying energy eigenstates exactly analytically in the limit of large interactions in the thermodynamic limit.

IV. DYNAMICS IN THE STRONGLY INTERACTING LIMIT
In this section we will discuss how we can understand the dynamics analytically for the case of strong interactions. In the extreme case of infinite U , at unit filling the particles will be in an idealized Mott insulator state, i.e., |ψ U →∞ = n |1 n . This state is an eigenstate of the jump operator with eigenvalue one, so that light scattering leaving particles in the lowest band does not change the state or introduce any extra energy to the system.

A. Hard-core bosons
Going beyond this trivial limit, one can use an approximation in which we assume that the bosons are hard-core bosons (HCBs) [43]. The idea is to restrict the Hilbert space to local states with a maximum of one particle per site. Mathematically, this can be achieved by using the standard Bose-Hubbard model (2) while additionally imposing fermionic anti-commutation relations to the creation/annihilation operators for bosons on the same site, {b i , b † i } = 1. These particles are not real fermions, since they still commute on different sites, [b i , b † j ] = 1. However this can be fixed by a Jordan Wigner transformation [44] b where the new quasi-particles now are real fermions, In this picture the Hamiltonian becomes a Hamiltonian for non-interacting fermions: and the ground-state takes the product form |ψ G = N n N i P i,n |i (where |i denotes the state of a particle at site i). The jump-operators -as is true for all site-local operators -are unaffected by the Jordan-Wigner transformation and so in this transformed nota- If the system has fewer particles than lattice sites, spontaneous emissions will give rise to heating in the lowest band, as the states are not eigenstates of c † i c i , which does not commute with the Hamiltonian. Simple results in this case are discussed in Ref. [39]. However, in the special case of unit filling, which we would like to treat here, the only state described by this simple HCB form is again the trivial state with a single particle on each site, which is an eigenstate of all operators c † i c i and c † i c i+1 . It is then clear that we need to go beyond this treatment in order to understand heating in a Mott Insulator state with finite U .

B. Doublon/Holon calculation
To study thermalization of spontaneous emission we have to go to the next order in the approximation and allow for states with 0, 1, and 2 particles per site. It turns out that also in this case approximate analytical calculations can be derived [26]. Assuming that the state of the system remains close to the Mott insulator state withn = 1, to first order excitations will be given by doubly occupied sites on the one hand, and holes on the other. It is now possible to introduce creation/annihilation operators with a vacuum given by the ideal MI state |vac i = |1 i as Since there can only be one hole or doublon on each site these quasi-particles must obey the hard-core constraint, which can again be expressed with on-site anti- Note that in principle one also has to add a constraint in order not to have a doublon and holon on the same site. However, based on the assumption that the total number of doublyoccupied sites and holes remains low we will neglect these terms, as was originally done in Ref. [26]. As in the case of the HCB model one can now turn the quasi-particles into proper fermions via the Jordan-Wigner transformationd In the picture of these fermionic quasi-particles, the Hamiltonian reads A discrete Fourier transformation into quasi-momentum space (x j = 1 √ M q e ijqa x q , for the operator x j , q ∈ {0, 2π/M, 2 × 2π/M . . . 2π}) gives the quadratic Hamiltonian (up to a constant) The corresponding dispersion relations are given by E d (q) = −4J cos(q) + U , E h (q) = −2J cos(q), and E dh (q) = 2J √ 2 sin(q). This Hamiltonian can be diagonalized straightforwardly by a Bogoliubov transformation where the new fermionic quasi-particles c d/h are linear combinations of the d q and h q , i.e. are correlated pairs of doublons and holons. The dispersion relations of the particles are given by (assuming U > 6J) The ground state in this approximation is the vacuum |vac , which is defined by c d/q |vac = 0 for all q.
C. Spontaneous emission in the doublon/holon picture Using the inverse transformation, (c † d,q , c h,−q ) U † B (q) = (d † q , h −q ), we can express the the local particle number operator for site m, n m in the Bogoliubov frame. Applying this operator to the ground-state, we find n m |GS = |vac where we introduced the matrix elements of the transformation as u q = U B (q) [1,1], and v q = U B (q) [2,1]. The time-evolved state after the jump is therefore √ N |ψ(t) = |vac where N ≡ GS|n 2 m |GS . The initial probablitydistribution of the wavepacket w(q, q ) ≡ u q v q − u q v q is proportional to and examples for U = 7J and U = 10J are depicted in Fig. 3. The wavepackets are peaked at opposite momenta and therefore, the jump creates holes and double occupations which move in the opposite direction. For example, for U/J = 10, the wavepackets are peaked at qa = −q a ≈ 0.30π and we can calculate the groupvelocity from (20) as d dqa d (q)| 0.3π ≈ 3.75a/J. We can test this analytical result with an exact t-DMRG time-evolution calculation after a single jump. We compare the dynamics after a jump in the center of a system with N = M = 48 site in Fig. 4a. There we plot the time-evolution of the number of doubly occupied sites as a function of time. The white dashed line in Fig. 4a indicates the light-cone according to the analytically expected group-velocity and is in excellent agreement with the numerical calculations.

D. Behaviour of the kinetic energy and similar variables
We now consider a specific type of observable, such as the kinetic energy, which in quasi-momentum space and in the Bogoliubov frame can be written as with matrix elements O ij for each q. Calculating the expectation value of this observable with the time-evolved state gives The important points to note are that all cross-terms corresponding to the O 12 and O 21 part of the operator disappear. This is because terms such as vac|c h,−p c d,p c † d,q c † h,−q |vac = δ p,q δ p,q become zero due to the symmetric form of the specific wave-packet that is created in the jump, w(q, q ) ≡ u q v q −u q v q , w(p, p) = 0. Similarly, Kronecker deltas in the terms of the form vac|c h,−r c d,r c † d,p c d,p c † d,q c † h,−q |vac = δ p,q δ r ,p δ r,q lead to a cancellation of time-dependent terms. Thus, we find that any observable of this type is time-independent after the spontaneous emission. This is a remarkable result, given the fact that the state (22) is not an eigenstate of the Hamiltonian. The effect of a jump operator here is that it puts the system into a state whose density matrix is diagonal in the space of single doublon-holon pair excitations, for which case these specific observables don't evolve in time. In particular, the expectation value of observables of the form (24) can be effectively expressed by a density matrix -O = tr(ρO) -of the form ρ = λ vac |vac vac| + q,q λ q,q |q, q q, q |, where |q, q ≡ c † d,q c † h,q |vac . Thus, we conclude that after a single spontaneous emission, observables of the form (24) can immediately be described by the diagonal ensemble, i.e. the effective density matrix for those observables takes a diagonal form in the basis low-lying energy states.
As an example, let us now consider the kinetic energy operator, which can be written as In the thermodynamic limit, we can replace the sum over quasi-momenta with the integral, (1/M ) q → d(qa)/(2π), which yields We compare this analytical result to our time-dependent t-DMRG calculations This is in excellent agreement with our numerical results even for relatively small U/J as shown in Fig. 4b.

E. Eigenstate expectation values
We can check the eigenstate expectation values of the kinetic energy analytically. In Ref. [24], it was found that, once one enters the MI regime, these do not coincide with the expectation values for Boltzmann distributions with the corresponding energies anymore (cf. Section III). Now we can check whether these agree with our analytical calculation. Therefore, using exact diagonalization we compute the 500 lowest eigenstates in a system of N = M = 10 sites. These are shown as dots in Fig. 5a,b for the case of U = 10J and U = 20J, respectively. In addition we add the values we obtain analytically for the ground state (grey dot) and the branch of single doublonholon pair excitations (grey lines). We find good agreement for large interactions. Since the first branch converges towards the analytical calculation, we conclude that this branch indeed due to the single doublon-holon pair excitations, which is further verified by the fact that it contains M 2 − M points. The second branch must thus correspond to two pair excitations. Already soon after entering the MI U ∼ 5 the gap and branches start to appear [24]. This indicates that the doublon-holon calculations are indeed a good approximation even for regimes that are not particularly deep in the MI phase.

V. SUMMARY AND OUTLOOK
In summary, we have investigated the non-equilibrium dynamics of bosons in an optical lattice in the presence of spontaneous emission events. We show clearly parameter regimes in which the system relaxes to thermal distributions for simple quanitites (especially the quasimomentum distribution), and others where it relaxes on short timescales to non-thermal values. We can understand this behavior by applying eigenstate thermalization considerations to small systems for which we can perform exact diagonalization calculations. We also show that we can understand the dynamics in the strongly inter-acting limit well in terms of propagating doublon-hole pairs, which are analytically tractable. All of these results are expected to be directly observable in ongoing experiments. In addition to the possibility to measure quasimomentum distributions, the propagation of doublon-holon pairs as investigated here is a key signature for the effects of spontaneous emissions that could be measured in quantum gas microscope experiments with single-site resolution in ways that were previously implemented for quenches within the MI regime [32]. Measurements of this type could be used as a diagnostic tool for heating of many-body states in optical lattices, which could in turn be used to improve the robustness of quantum simulators. In the future, these studies can be continued towards thermalization for fermions in optical lattices undergoing spontaneous emissions [29][30][31], and the incorporation of partial thermalization from excitations to higher Bloch bands.