Probing local relaxation of cold atoms in optical superlattices

# Probing local relaxation of cold atoms in optical superlattices

A. Flesch, M. Cramer, I. P. McCulloch, U. Schollwöck, and J. Eisert 1 Institut für Theoretische Physik C, RWTH Aachen University, 52056 Aachen, Germany
2 Institute for Mathematical Sciences, Imperial College London, SW7 2PE London, UK
3 School of Physical Sciences, The University of Queensland, Brisbane, QLD 4072, Australia
4 Institute for Physics and Astronomy, University of Potsdam, 14476 Potsdam, Germany
###### Abstract

In the study of relaxation processes in coherent non-equilibrium dynamics of quenched quantum systems, ultracold atoms in optical superlattices with periodicity provide a very fruitful test ground. In this work, we consider the dynamics of a particular, experimentally accessible initial state prepared in a superlattice structure evolving under a Bose-Hubbard Hamiltonian in the entire range of interaction strengths, further investigating the issues raised in Ref. [Phys. Rev. Lett. 101, 063001 (2008)]. We investigate the relaxation dynamics analytically in the non interacting and hard core bosonic limits, deriving explicit expressions for the dynamics of certain correlation functions, and numerically for finite interaction strengths using the time-dependent density-matrix renormalization (t-DMRG) approach. We can identify signatures of local relaxation that can be accessed experimentally with present technology. While the global system preserves the information about the initial condition, locally the system relaxes to the state having maximum entropy respecting the constraints of the initial condition. For finite interaction strengths and finite times, the relaxation dynamics contains signatures of the relaxation dynamics of both the non-interacting and hard core bosonic limits.

July 5, 2019

## I Introduction

Both in classical and quantum physics, a complete framework for the description of equilibrium properties of arbitrary physical many-body systems exists, although the explicit calculation of many-body equilibrium properties is an often unsolved challenge. The situation is much less satisfactory when it comes to the study of the non-equilibrium properties of many-body systems, where such a general framework is missing and may not even exist. Research in this field has therefore focussed on relatively specific issues and types of non-equilibrium. One of the issues taking center stage is whether quantum many-body systems in non-equilibrium evolving coherently under a local Hamiltonian equilibrate or not. If so, one may ask whether the equilibrium states can in some sense be described by a thermal state. This old and fundamental question of the equilibration of quantum many-body systems has enjoyed quite a renaissance recently [1–36].

A specific setting of coherent non-equilibrium quantum dynamics is provided by quantum quenches, where one starts from an eigenstate of some initial Hamiltonian and pushes the system out of equilibrium by a sudden change (or “quench”) of system parameters, leading to a new Hamiltonian. One then considers the evolution of the system under the new Hamiltonian. A further restriction is provided by the assumption that the quantum system under consideration is closed, i.e., has no coupling to a bath of degrees of freedom that might assist the relaxation process. Time evolution (and hence the potential relaxation to some equilibrium) will obviously be constrained by the constants of motion, i.e., Hermitian operators commuting with the new Hamiltonian whose expectation values are fixed by the initial state; in that sense, any relaxed state will to a certain degree show some memory of the initial state.

It has been conjectured that—in some sense— the quantum system should relax to the maximum entropy state consistent with the expectation values of the constants of motion fixed by the initial state (see, e.g., Refs. Rigol (); Olshanii (); Cazalilla ()), also referred to as a generalized Gibbs ensemble Jaynes (). This may be attractive because it reminds of Jaynes’ derivation of equilibrium statistical mechanics.

This observation is in conflict with the fact that obviously, if the system can be meaningfully treated as a closed quantum system, one cannot expect the whole system to relax: Initially pure states will remain so in time under a unitary time evolution [2–6]. After all, the entire information about the initial condition is still stored in the system, albeit in a dilute fashion. Yet, this observation is by no means in contradiction with the possibility that in any local observation, the system may appear perfectly relaxed, even without invoking a time average Relax (); Tegmark (); Barthel (); Kollath (). The key point is that locally, one may well expect the relaxation to be true Momentum (): For any subset of sites in a sufficiently large lattice system, the reduced state may well converge to the reduced state of the maximum entropy state, given the conserved quantities of motion, and stay relaxed for an arbitrary long time. Indeed, for such a subset it is, under suitable assumptions about its interactions with the rest of the world, easy to make contact to Jaynes’ formulation of statistical mechanics.

Ref. Relax () considers a variant of the question in which this local relaxation of subsystems, referred to as local relaxation conjecture, can in fact be rigorously proven to hold exactly: This is the one where one evolves a state deep in a Mott phase according to a Bose-Hubbard Hamiltonian corresponding to the deep superfluid phase—treated as a non-interacting system. In this setting, the local relaxation is in fact true: The reduced state of a block of consecutive sites indeed converges to a maximum entropy state

 ^ϱI(t)→^ϱmax (1)

(in trace norm) for large systems and large times, having maximum entropy consistent with the constants of motion Relax (). Note that there is no time average and the initial state is not a Gaussian state.

Also, for free bosonic and fermionic systems, and for Gaussian initial states, it has been shown rigorously in Refs. Tegmark (); Barthel () under which conditions the local relaxation conjecture is true. It turns out that in these cases, the presence and speed of local relaxation depends crucially on the single-particle spectrum and the dimensionality of the systems.

The in some sense inverse case of a quench from the superfluid phase to the Mott insulating phase, but generally at finite interaction strengths, has been studied in Ref. Kollath (). In this non-integrable system, numerically two distinct non-equilibrium regimes have been found where equilibrated states resembled thermal states or states with memory.

The physical intuition why local relaxation happens is the following: If one switches to a new Hamiltonian, the system is no longer in equilibrium. Hence, one has local excitations at each point [2, 3, 12–16]. They cannot travel arbitrarily fast, however, as there is a finite speed of information transfer in lattice systems LR (); LRC (). At each site, in time more and more “waves” of farther and farther away sites can possibly have a significant influence on this site. This is related to a finite speed of sound or of information transfer in a lattice system LRC (). Then, local relaxation may be a consequence of the incommensurate influence of these excitations generating mixing and thermalization. The “thermalization” time scale is hence governed by the speed of information transfer Note ().

This also links to kinematical approaches to the problem [37, 47–51], arguing that most states anyway look locally very much relaxed in that they have large entropy. A random pure state (as taken from the Haar measure) will locally have a large entropy. Specifically, in interacting systems one should expect such a local relaxation to be true, too, an aspect that will be studied in this work.

So far the discussion has not taken into account the recurrence happening in any closed quantum or classical system Relax (). For finite, but large system sizes, recurrence times will however become so long that recurrence effects become negligible, in that the quantum many-body system can be effectively locally equilibrated on a much shorter time scale than the recurrence time.

We will see that while settings exhibiting local relaxation may be generated in various fashions, it is a greater challenge to actually probe signatures of local relaxation. This apparent dilemma—that to demonstrate local relaxation appears to necessitate local addressing—will be resolved in this work, by making use of a period- setting, further developing the idea of Ref. Letter (). The taken path opens up a way to experimentally explore local relaxation effects using atoms in optical superlattices. We systematically investigate various aspects of local relaxation in such a setting, both using analytical as well as numerical methods, based on a time-dependent density-matrix renormalization-group (t-DMRG) approach.

## Ii Experimental setup: Ultracold atoms in optical superlattices

The current surge of interest in the relaxation of quantum systems after a quench is mainly motivated by the advent of ultracold atoms in optical lattices [52–56]. These systems are highly attractive as they allow for sudden controlled manipulations of system parameters, hence quenches, are strongly interacting, hence non-trivial, are on experimental time scales essentially closed quantum systems, and therefore show coherent quantum dynamics. Systems are also sufficiently large to show non-trivial many-body behavior.

In the present context, however, the major drawback of ultracold atoms is that despite the unprecedented possibilities of manipulation the study of local relaxation provides an experimental challenge. This is due to the fact that local (i.e., site-resolved) measurements on ultracold atoms in optical lattices are still not satisfactory, albeit rapid progress is being made. In this work, we propose to study instances of local relaxation even in a setting where strictly local quantities can not easily be studied by using optical superlattices [57–60].

Following the setup very recently realized by Bloch and coworkers Wells (); TimeResolved (), we consider bosonic Rubidium-87 atoms in a period- optical superlattice geometry: Two standing-wave laser fields at wavelengths nm and nm are superimposed with fixed relative location to provide a superlattice geometry of lattice constants nm and nm. It is experimentally possible, among other things, to (i) change the relative strength of the two optical lattices and (ii) shift their relative position, by altering phase and detuning of the optical superlattice. Assuming the strength of the lattice with lattice constant to be fixed, (i) allows to couple and decouple this lattice into an array of double well potentials, combining one odd (o) and one even (e) site of the original lattice. In such double well potentials, (ii) allows to introduce a bias between the chemical potentials of the odd and even sites, .

Isolating double wells and biasing odd vs. even sites, in turn, allows for the preparation of patterns of atoms and for extracting local quantities to the degree that odd and even sites can be distinguished from each other:

• Preparation of periodic patterns is achieved by loading the superlattice while introducing a bias between odd and even sites such that due to a shifting of particles all particles are on either odd or even sites after loading. Using further experimental techniques Wells (), multiple occupancies on the occupied sites can be eliminated, leaving a sequence of empty and single-occupied sites.

• Period- local measurements can be obtained by mapping odd and even sites to different Brillouin zones: Assuming completely decoupled double wells, each part of the double well has multiple bands separated by well-defined energies. Biasing, say, the odd sites relative to the even sites by an energy in excess of the separation energy of the band-separation energy, odd-site particles are reloaded into the higher band of the even sites, whereas the even-site particles stay in the lower band. A standard time-of-flight mapping then shows the even-site particles in the first Brillouin zone, the odd-site particles in the higher Brillouin zones.

• More sophisticated measurements with period- can also be performed in principle. When letting the system evolve for some defined hold time, and then freezing the time evolution by ramping up the barrier, one can also measure nearest-neighbor correlations in the lattice Wells (); TimeResolved (); Superlattice (); Superlattice2 (). There are several ways to proceed here: On the one hand, one can isolate double wells, let them dephase. Then upon free expansion, the average correlations can be measured as in a double slit experiment. One the other hand, correlations can be mapped to densities. To the extent that the barrier between the wells is sufficiently high such that the time evolution can be described by a collection of double well systems (and higher Wannier bands do not have to be taken into account), we can investigate the time evolution in each individual double well. Up to the presence of a confining potential, the total time evolution then reflects the time evolution in each of the double wells. Let us label the two modes in any of the wells by and , one may apply an appropriate free Hamiltonian to map correlations onto on-site densities. Specifically, for we find that

 2im⟨^a†1^a2⟩=⟨^a†2(t)^a2(t)⟩−⟨^a†1(t)^a1(t)⟩), (2)

when letting evolve under the free Hamiltonian

 ^h=−(^a†1^a2+^a†2^a1), (3)

until . One can hence measure the imaginary part of the appropriate correlators, and hence the period- correlators in the full lattice. In settings where in the Heisenberg picture the phase map for each right well is feasible, one can also measure the real part. If one has small interactions (because of the Wannier band problem—a description in terms of Wannier functions still has to be valid, accompanied by a non-vanishing ), then this leads to a small dephasing in this mapping. Direct numerical simulation can take this properly into account, allowing for a correct interpretation of the experimental observation. Such a technique allows to measure correlators from a time-resolved observation of on-site densities. Similar time-resolved measurements have recently been performed in optical superlattices TimeResolved ().

As shown in Fig. 1, we propose to start from a two-periodic initial state prepared by the superlattice setup as described above, where all odd sites are occupied by exactly 1 boson and all even sites are empty. The initial state vector of the entire lattice is hence

 |ψ(0)⟩=|1,0,1,0,…1,0⟩. (4)

Tools how to experimentally achieve that sites are to a good approximation occupied by a single atom, and not two or more atoms, are described in Refs. [57–59]. If the -lattice is suddenly switched off, i.e., the system quenched, the state vector will evolve in time according to the conventional Bose-Hubbard Hamiltonian

 ^H=−JL∑i=1(^b†i+1^bi+^b†i^bi+1)+U2L∑i=1^ni(^ni−1), (5)

where and are the standard interaction and hopping parameters of the Bose-Hubbard model that can be calculated microscopically from the lattice parameters. The system size is given by , whereas the boson number is from this setup. will always be even in the following, in line with the proposed setup; various boundary conditions will be imposed. We work in units where and .

We will not consider the effect of the occupation of higher order bands, but will stay within the limit of applicability of the Bose-Hubbard model. We will for the purposes of the present paper also neither consider an additional harmonic confining potential—except when explicitly stated otherwise—nor additional dephasing effects due to statistical fluctuations of local fields. Instead, we will systematically flesh out what behavior is expected in this slightly idealized setting, to see how local relaxation manifests itself here, to comment on the impact of imperfections later.

As is not an eigenstate of , we are in a non-equilibrium situation. The initial state is uncorrelated and shows inhomogeneous density of periodicity . Coherent quantum dynamics is now expected to homogenize densities locally and to build up non-local correlations.

Local relaxation is now monitored by the global measurement of the total occupation of even, , and odd, , sites. In a translationally invariant setting, this gives access to local observables as

 ⟨^Ne(t)⟩=L/2∑i=1⟨^n2i(t)⟩=L2⟨^n2i(t)⟩, (6)

and similarly for odd sites.

If one has experimental access to the variances , density-density correlations (see also Ref. NoiseNoise ()) between all even and all odd sites may be obtained through

 σ2tot−σ2e(t)−σ2o(t)2=⟨^Ne(t)^No(t)⟩−⟨^Ne(t)⟩⟨^No(t)⟩ (7) = L/2∑i,j=1(⟨^n2i(t)^n2j−1(t)⟩−⟨^n2i(t)⟩⟨^n2j−1(t)⟩).

In fact, the value of is fixed to for the proposed pattern, so in repeated experiments of same length, this quantity will have variance zero, .

Moreover, the quasi-momentum distribution defined as

 S(q,t)=1LL∑i,j=1eiq(i−j)⟨^b†i(t)^bj(t)⟩,q∈[0,2π], (8)

is also accessible as a global quantity, from time-of-flight measurements.

We now proceed as follows. In order to obtain analytical results we set and respectively, which are both essentially non-interacting limiting cases, can be solved exactly and show very similar, but not generally identical results (compare also Ref. LightCone ()). These results will be extended to the interacting finite- case which will be studied using t-DMRG method. For the time-scales considered it is an effectively quasi-exact method.

## Iii Exact solutions for limiting cases

Both limiting cases and are or can be mapped to free models. This case has been considered for Gaussian initial conditions in Refs. Tegmark (); Barthel (), as well as for non-Gaussian product initial conditions in Ref. Relax (). In these cases, the reduced state of a subsystem

 ^ϱI(t)=trL∖I[^ϱ(t)]→^ϱmax. (9)

will converge for large systems and long times to , being the maximum entropy density operator as constrained by the initial conditions (or, more precisely, if recurrences are considered, this will be true in 1-norm to an arbitrarily small approximation error for an arbitrarily long time). Also, fermionic Gaussian initial states have been considered in Refs. Relax (); Barthel (). Local relaxation will therefore happen in these cases.

The physical mechanisms behind the relaxation processes in the cases of and is sightly different, however: In the case , it is due to dephasing in the sense that freely propagating bosons lead to reduced density operator contributions of quickly oscillating phases that average out. In Ref. Relax (), this intuition leads even for non-Gaussian initial states to maximum entropy states, by invoking a quantum version of a central limit theorem, and exploiting the finite speed of information transfer LRC ().

In the case , real scattering processes happen, albeit of a very specific form that allows for a formal mapping to a non-interacting fermionic problem (and to the XX model, compare also Ref. LightCone ()): The scattering is simply relegated to the internal sign structure of the wave function. We will also see that the two relaxation processes lead to different time-evolutions of most physical observables.

For all cases of nonzero finite , the situation should be somehow intermediate. For the specific setup chosen here, interacting particles will learn of each others existence only after some initial time they need to come into contact. We would therefore expect that for very short time scales, observables should evolve as in the limit, with a crossover in behavior for longer time scales when they start interacting. The question is whether for which interaction strengths the limiting pictures remain essentially valid and whether there is a genuinely different intermediate interaction regime with different relaxation behavior. Compared to Refs. Relax (); Tegmark (); Barthel (), in this work we focus to a lesser extent on rigorous mathematical convergence statements—like invoking quantum versions of central limit theorems—but instead put more emphasis on the phenomenology of the physical relaxation process as such in the -periodic setting.

### iii.1 Non-interacting bosons: U=0

In the translationally invariant case of non-interacting bosons, the Hamiltonian can be diagonalized by Fourier transforming to new bosonic operators to yield

 ^HU=0=L∑k=1λk^a†k^ak, (10)

where

 λk=−2cos(2πk/L) (11)

are the eigenvalues of the circulant Hamilton matrix with entries . While the real experiment will not have periodic boundary conditions, our calculations show that for realistic system sizes and time scales the difference between open and periodic boundary conditions is negligible.

In the Heisenberg picture, the initial state vector remains time-independent, whereas the operators evolve in time as

 ^bi(t)=L∑j=1Vi,j(t)^bj(0),V(t)=e−itH,Vi,j(t)=1LL∑k=1e−itλke2πik(i−j)/L. (12)

Straightforward algebra then yields the exact time-evolution of two-point correlations, see Appendix,

 fi,j=⟨^b†i(t)^bj(t)⟩=12δi,j−(−1)i2Vj,i(2t). (13)

In the thermodynamic limit this turns into

 fi,j=12δi,j−(−1)iij−i2Jj−i(4Jt), (14)

where is a Bessel function of the first kind.

From the above expression, one can derive the quasi-momentum distribution for finite , see Appendix. For it is constant, , while one finds

 (15)

for all other .

With slightly more effort, density-density correlations as in Eq. (7) emerge as (see Appendix)

 ⟨^n2i(t)^n2j−1(t)⟩−⟨^n2i(t)⟩⟨^n2j−1(t)⟩ (16) =f2i,2j−1f2j−1,2i−2L/2∑k=1|V2i,2k−1(t)V2j−1,2k−1(t)|2,

which, as they are local quantities, relax for large systems and long times. For large one finds for the global density-density correlator (see Appendix)

 ⟨^Ne(t)^No(t)⟩−⟨^Ne(t)⟩⟨^No(t)⟩→−L16(3+J0(8Jt)−4[J0(4Jt)]2). (17)

We also show in detail the exact local relaxation to a maximum entropy state for the case of the initial condition , largely following Ref. Relax (): For every block of consecutive sites , every small error and any desired, arbitrarily long recurrence time there exists a system size and a local relaxation time such that

 ∥^ϱI(t)−^ϱmax∥1<ε (18)

for all times , see Appendix. Here, denotes the trace norm. Hence, after the quench for , for -periodic initial conditions, the local state becomes in time a maximum entropy state, to an arbitrarily good approximation.

### iii.2 Hardcore bosons: U→∞

In the limit , the interaction manifests itself exclusively in that bosonic occupation numbers are upper bounded by . This leads to a well-known mapping in case of quantum ground states: The hard core limit is equivalent to the XX spin model and a model of spinless free fermions. But even in time evolution, in the limit of large , the population of sites with particle number larger than unity is dynamically suppressed to an arbitrary extent: The expectation value of the new Hamiltonian is obviously preserved under time evolution, , , where

 ^H0=L∑i=1(^b†i+1^bi+^b†i^bi+1),^H1=12L∑i=1^ni(^ni−1). (19)

For initial states that are supported on span one has . Furthermore , where

 ^H2=L∑i=1∞∑k=2|k⟩⟨k|, (20)

where the latter projector defined on each site . This, in turn, means that

 tr[^H0(^ϱ(t)−^ϱ(0))]U≥% tr[^H2^ϱ(t)]≥0. (21)

Hence, one can ensure to arbitrary accuracy for large that is only supported on . One therefore arrives in a perfectly meaningful way at the familiar hard core limit of the Bose-Hubbard model, which is equivalent to a free fermionic spinless system.

This mapping to non-interacting spinless fermions is done through the familiar Jordan-Wigner transformation

 ^bn=e−iπ∑m

Under this transformation the initial state vectors turn into

 |ψ⟩=^c†1(0)^c†3(0)^c†5(0)…|0⟩ (23)

and the Hamiltonian reads upon a Fourier transformation to new fermionic operators

 ^HU=∞=L∑k=1λk^d†k^dk (24)

with as in Eq. (11) and the time-evolution of the operators is given by

 ^ci(t)=L∑j=1Vi,j(t)^cj(0) (25)

with as in Eq. (12). This is formally identical to the case. However, there are two differences: Periodic boundary conditions in the original bosonic model map to (anti-)periodic boundary conditions . Hence, boundary conditions stay periodic if the particle number is odd, to which we restrict ourselves in the following.

The second, more important difference is that differences in results will show up due to the Jordan-Wigner transformation of operators and the difference between bosonic and fermionic (anti)commutators. It is easily shown that local densities and correlations between nearest neighbors translate directly (see XIII.4), such that the respective results for are identical to those for . For longer-ranged two-point correlations, including structure functions, as well as density-density correlations results differ. Density-density correlations are given by (see XIII.4)

 ⟨^n2i(t)^n2j−1(t)⟩−⟨^n2i(t)⟩⟨^n2j−1(t)⟩=−f2i,2j−1f2j−1,2i, (26)

where is as in Eq. (13). For large the global density-density correlator is then given by (see Appendix)

 ⟨^Ne(t)^No(t)⟩−⟨^Ne(t)⟩⟨^No(t)⟩→−L16(1−J0(8Jt)). (27)

### iii.3 Interacting softcore bosons: Finite U

In this case, we are no longer facing an integrable model. In order to study the relaxation dynamics, we will therefore turn to the time-dependent variant VidalTEBD (); Daley (); WhiteFeiguin () of the DMRG (density-matrix renormalization group) method White (); RMP (). This method allows to follow the coherent time-evolution of strongly interacting quantum systems very precisely. Its reach in time is, however, limited by the growth of entanglement in quantum systems: Linear entanglement growth, for example, leads to an exponential growth in numerical resources needed. As we will see, interestingly, the system under study is characterized by very strong linear entanglement growth. In the free instances this linear increase is indeed provably correct Relax (); Schuch (), see Sec. IX. Hence, the reachable times are quite short () with up to states kept in the simulations. For most quantities, this turns out to be sufficiently long to make contact to the non-interacting results and to read off long-time behavior. In particular, the results allow for a quantitative comparison to experiments.

A subtlety arises from the fact that DMRG prefers open boundary conditions, and is hence closer to experiment. However, we would like also to compare to analytical results, where periodic boundary conditions are preferable. As it will turn out, for system sizes and times considered, the difference is negligible.

## Iv Time evolution of densities

For the time evolution of local densities both exactly solvable cases and give the same result

 ⟨^ni(t)⟩=12−(−1)i2LL∑k=1e4itJcos(2πk/L)→12−(−1)i2J0(4Jt)(L→∞). (28)

Odd- and even-site densities relax symmetrically about the axis to , with an asymptotic decay as .

All -DMRG results for finite are perfectly compatible with relaxation of local densities to . On very short time scales () particles have typically not collided yet and are not yet sensitive to the values of , hence, finite results are identical to the limiting cases . The relaxation behavior for intermediate times, however, deviates quite strongly. Interaction effects become visible right after particles make contact. This can be seen in Figs. 2-6 (all calculated for ), where we compare the time evolution of local densities for various finite values of to the special cases . For small (exemplified by ) and large (exemplified by ), the comparison indicates that the relaxation of local densities is governed by the behavior of the limiting non-interacting cases. For intermediate values of (exemplified by ), scattering seems to be most effective and lead to much faster damping and relaxation, much beyond the above square root time dependence. This is plausible, as close to a quantum critical model there is a limiting point in the spectrum of the new Hamiltonian, leading to specifically effective relaxation. Deviations from the limiting behavior are sufficiently strong that they should be visible experimentally.

Under the assumption that for the time scales achievable by DMRG an asymptotic power-law can already be read off, one finds the exponents given in Fig. 7. While for very small and all one finds a slope similar to , as for the limiting cases, there is an intermediate regime where relaxation is much faster. It must be stated, however, that in this regime the precise slope is hard to extract, even an exponentially fast decay cannot be excluded completely, but we are not aware of a physical reason why there should be a qualitative change of decay behavior from power-law to exponential. Qualitatively, stronger interactions should lead to a more distinct dispersion, leading in turn to a quicker decay.

Concerning experimental implementations, a number of further issues require a discussion: The experimentally accessible systems will be of finite size, leading to recurrence effects. Moreover, there will be a parabolic confining potential, which will modify results. Furthermore, one will encounter initial states at non-zero temperature, as well as possible additional fluctuations due to inhomogeneities in applied fields, different from the apparent local relaxation.

So far, we have confronted analytical results in the limit with DMRG results for . In order to check whether on the time scales reachable by DMRG recurrence effects can be seen for this system size, we have rerun selected calculations for , observing no relative change in results above 1 percent. In particular, the shapes of oscillatory behavior remained completely unchanged. In the limit, finite system DMRG results and infinite as well as finite system analytic results agree completely on the finite times reachable by DMRG, even despite the different boundary conditions (open boundary conditions (OBC) for DMRG, periodic boundary conditions (PBC) for analytical statements). Recurrences, where the difference between finite and infinite system becomes obvious, happen much later, see Fig. 8.

The effect of the trap, as shown in Fig. 10, is to generate effective reflections from the edges of the system, leading to much earlier recurrences. However, in experimental setups, this effect would for realistic traps set in late enough to allow for sufficiently long observation.

## V Time evolution of local correlators

Let us now consider the nearest-neighbor correlator . This quantity is specifically interesting, and can in principle be measured by means of exploiting the tuning of the double well potential of the superlattice, and appropriate timing (see above, and compare also Refs. Wells (); TimeResolved (); Superlattice ()). It is of interest as it goes beyond local densities: The build-up of correlations in time starting from the uncorrelated initial state becomes visible.

In the limiting free cases, identical results are found:

 ⟨^b†i+1(t)^bi(t)⟩=(−1)i2LL∑k=1e4Jticos(2πk/L)e−2πik/L→−(−1)i2iJ1(4Jt)(L→∞). (29)

The real part of the correlator is strictly zero for all times, whereas the imaginary part relaxes to with an asymptotics of after a quick growth to a maximal value of about at time . This quick growth reflects the buildup of correlations due to particle motion with speed linear in .

For finite , the scenario has marked similarities and differences. Considering Figs. 11 and 12, one sees that on short time scales the buildup of the imaginary part of correlations is identical for arbitrary . This simply reflects the fact that due to the distance between particles at , no collisions have yet happened on these time scales. Only when the interaction strength becomes visible, the relaxation to 0 follows different paths. As for local densities there is a clear trend that relaxation is fastest around , reflecting the particularly efficient scattering there. In the real part, convergence to a finite value occurs for all finite , but not such a clear picture of relaxation speeds occurs.

If one plots the asymptotic values of the real part of the nearest neighbor correlators (Fig. 13), one sees that there is a maximum around , reflecting the efficient scattering there. Interestingly, the large- behavior can be very well approximated by a curve for all and larger.

Indeed, this dependence is exactly what one would expect in the thermal or Gibbs state of the Bose-Hubbard Hamiltonian : We hence look at local correlations in the Gibbs state

 ^ϱβ=e−β^H/Z (30)

with . We will take

 ^H0=U2L∑i=1^ni(^ni−1). (31)

as our unperturbed Hamiltonian including interactions and look at leading orders in the hopping term

 ^V=−JL∑i=1(^b†i+1^bi+^b†i^bi+1). (32)

Then we find, using standard thermal perturbation theory, up to first order in

 (33)

This means that we have ( if and are nearest neighbors, zero otherwise)

 ⟨^b†i^bj⟩=δi,j⟨^b†i^bi⟩ (34) +δ⟨i,j⟩JU∑n,me−β(En+Em)n(m+1)eβU(n−m−1)−1z2(n−m−1)

where and are the local energies of the unperturbed Hamiltonian. So, indeed, within the validity of perturbation theory, we do find the anticipated linear dependence on , as seen also in DMRG simulations in the time-dependent scenario. By definition, the correlators merely probe local quantities. This corroborates the intuition that locally, the system is indistinguishable from the situation as if globally the system was in a state maximizing the entropy, under the constraints of motion Relax (); Tegmark (); Barthel ().

To complete the picture, let us finally remind ourselves of how the maximum entropy states of the global system locally look like. The constants of motion of interest here are the occupation numbers in bosonic or fermionic momentum space, respectively. The global maximum entropy state—under the constraints of motion—for is then found to be

 ^ϱ=L⨂i=123e−ln(3)^b†i^bi. (35)

This result is completely consistent with our result for the relaxation of the characteristic function (see Appendix XIII.2). One can now calculate expectation values in equilibrium for the local observables considered above (by construction in agreement with our earlier analytical findings) as

 ⟨^b†i^bj⟩=12δi,j,⟨^ni^nj⟩−⟨^ni⟩⟨^nj⟩=34δi,j. (36)

In the case of one finds

 ^ϱ=12L\mathbbm1, (37)

which means that

 ⟨^b†i^bj⟩=12δi,j,⟨^ni^nj⟩−⟨^ni⟩⟨^nj⟩=14δi,j. (38)

We have also numerically investigated the relaxation behavior of longer-ranged correlators . As an example, we depict in Fig. 14 the relaxation of the real part of , the imaginary part being zero. Relaxation here is of monotonically increasing effectiveness with . The relaxation dynamics for and is different, revealing the fundamentally different character of the non-interacting limiting cases.

## Vi Local relaxation

Let us summarize the properties of local relaxation that can be extracted in this experimental setup. In all quantities, there are three time regimes:

1. The first time regime is associated with the building up of correlations: This is governed directly by the coupling strength to the nearest neighbor, and on the shortest times is identical for all , before collision processes become important. The time scale is , as this is when the first collisions happen and establish correlations.

2. The second time regime is associated with local relaxation. There is fast oscillating dynamics between neighboring sites, happening at a time scale dictated by the hopping , which also governs the speed of sound. Local relaxation as such is slow, in the exact analysis a slow polynomial decay: The Bessel function fulfills as a strict bound, with an asymptotic envelope

 J0(x)=x−1/2+o(x−1/2). (39)

one finds relaxation (not due to decoherence but due the dilution of information over the lattice) to the true local equilibrium. For finite , numerics is consistent with polynomial decay, which for intermediate seems to be much faster, but still polynomial.

The intuition is that this local relaxation is due to influences of excitations travelling with the speed of information propagation from farther and farther lattice sites, broadened by dispersion. This means that this time scale, even in the interacting case of , should be defined by the speed of information propagation in the system. Note that to date, there is no rigorous bound known to the speed of information propagation in the fully interacting Bose-Hubbard model (this being a consequence of the lattice sites having an unbounded number of particles, see also comments in Refs. Nachtergaele (); HarmonicLiebRobinson (); Buer ()). On physical grounds, yet, it should be expected to be similar to , i.e. ”ballistic” transport at the high energies provided by the quench, slightly modified by (possibly the resulting bounds can only be formulated for specific initial states having small local particle number). This is also what the numerics shows. Note that the low-energy speed of sound need not be relevant here.

3. The third (very large) time is the recurrence time, which seems to be shortened substantially in the presence of a realistic trap. As in the presence of a trap, the excitations no longer travel with a constant speed (the speed of sound), but are slowly reflected, one should expect a quicker relaxation, and this is also the behavior the analytics exhibits. Generally, for moderate to large system sizes it is already beyond the reach of our simulations. It would be interesting to see, possibly in exact diagonalization, whether the recurrences are weakened compared to the free solutions due to interactions in the system, see also Ref. LightCone ().

While these findings are essentially independent of the interaction strength , we also find three local relaxation regimes for different :

• For very small values of (up to ) relaxation dynamics is quite close to the non-interacting bosonic limit.

• For larger values of the system seems to show a behavior very similar to the hard core bosonic or free fermionic limit with similar relaxation exponents of order . Local correlations are consistent with locally equilibrated subsystems perturbatively coupled in .

• The case of intermediate appears to be a special case for various observables, marking the “boundary” between the “free bosonic case” and the “hard core boson case” , so the fermionic one. Collisions lead to the most efficient relaxation in this case.

To summarize, the dynamics of local quantities shown is consistent with the limiting and cases, but shows a richer phenomenology in particular for intermediate .

## Vii Quasi-momentum distribution

In this subsection, we briefly compare what is seen in the above local relaxation with the situation in momentum space. The quasi- momentum distribution,

 S(q,t)=1LL∑i,j=1eiq(i−j)⟨^b†i(t)^bj(t)⟩, (40)

obtainable from time-of-flight experiments, is no longer a local quantity, but a global one probing the state of the entire lattice, as well as the boundary conditions and possibly a confining additional potential. In fact, in the free model for , one finds that the quasi-momentum distribution will be very little time-dependent, and if so, in a chaotic fashion showing little structure, see Fig. 17. More specifically, we find (see Appendix) for , , and

 (41)

for all other . In case of , the quasi-momentum distribution hence indeed does not show any signatures of the local dynamics, see Fig. 17. This is, again, not inconsistent at all with the notion of local relaxation: This quantity is a global one (and local only insofar as correlators of far away lattice sites are suppressed via a quickly rotating exponential function). For any subblock, the dynamics drives the system towards the values of equilibrium. In momentum space, this is not necessarily the case; at least not on these time scales. In the free situation of , this leads to an absence of relaxation of this quantity.

It would still be interesting to measure this quantity for close to free cases for short times, to demonstrate the very point that the information about the initial condition is fully contained in the system, and merely diluted. Hence, obviously, one cannot expect true global relaxation to happen, unless a further mechanism of decoherence due to additional external degrees of freedom is present.

In Fig. 18 we show the quasi-momentum distribution for the situation including a trapping potential. The effect of the additional harmonic potential is very significant. Yet, again, note the absence of signatures of the time scales of local dynamics as compared to Figs. 9, 10.

How does this compare to the case of having a finite on-site interaction ? Again, the behavior reminds of the situation in the free case of . The quantity is too non-local to grasp the effect of local relaxation. Still, for longer times, one may also expect relaxation in some translationally invariant properties in momentum space, and the quasi-momentum distribution may well relax. Such a situation is observed in a flow equations approach in case of an interacting model in Ref. Moeckel (). For the DMRG calculations we use a slightly different definition of the quasi-momentum distribution

 S(q,t)=2L+1L∑i,j=1sin[πqiL+1]sin[πqjL+1]⟨^b†i(t)^bj(t)⟩. (42)

This ensures that is a constant of motion for in the free system in case of OBC. Fig. 15 and Fig. 16 show this quasi-momentum distribution for for different values of . Indeed, signatures of such a behavior are also observed in Fig. 15, where the momentum distribution is depicted for different values of time . For larger times, the quasi-momentum distribution appears to relax in the interacting case of . This effect is even more distinct in the case of stronger interactions, as depicted in Fig. 16.

To reiterate, the very absence of relaxation in close to free settings on short times gives rise to an interesting situation: One could experimentally see signatures of local relaxation. One could, however, also see that the information becomes more and more dilute in position space, but is still preserved in the system as such, as seen in the quasi-momentum distribution.

## Viii Global density-density correlator

One potentially experimentally accessible global correlator is the connected density-density correlator

 ⟨^Ne(t)^No(t)⟩−⟨^Ne(t)⟩⟨^No(t)⟩ (43)

which we show for different in Fig. 19. Numerics indicates that this quantity looks effectively relaxed after some time as long-range correlators contributing will be small. For the free case of , one analytically finds for large (see Appendix)

 ⟨^Ne(t)^No(t)⟩−⟨^Ne(t)⟩⟨^No(t)⟩→−L16(3+J0(8Jt)−4[J0(4Jt)]2), (44)

which relaxes to for large times. For hardcore bosons the global density-density correlator for large is given by

 ⟨^Ne(t)^No(t)⟩−⟨^Ne(t)⟩⟨^No(t)⟩→−L16(1−J0(8Jt)), (45)

which relaxes to for large times.

## Ix Entanglement and entropy growth

The creation of excitations at all points of the lattice will in general create a significant degree of entanglement in the time evolving state. Since at each site, an excitation starts to travel through the lattice, a linear increase of the entanglement entropy of subblocks in time is to be expected Entanglement0 (); Entanglement (); Entanglement2 (). More precisely, it has been shown in Refs. Entanglement (); Entanglement2 () that for any time evolution under a local Hamiltonian with finite-dimensional constituents, starting from a product state, the entanglement entropy of a subblock grows in one-dimensional systems in at most as

 E(^ϱI(t))≤c0t+c1, (46)

for suitable constants . This means that for any constant time, the entanglement entropy satisfies what is called an area law. For larger and larger block sizes, the entanglement entropy will eventually saturate in (and, e.g., not logarithmically diverge).

In turn, this upper bound is saturated, in the following sense: In Ref. Relax () and explicitly in the Appendix -norm convergence is shown to products of Gaussian states

 ∥^ϱI(t)−^ϱmax∥1=∥^ϱI(t)−^ϱ⊗sG∥1→0, (47)

when considering initial conditions and time evolution under , where is a Gaussian product state. This observation immediately implies a bound to the entanglement entropy that is linear in time, if the block size can only be chosen appropriately. This has also been made explicit in Ref. Schuch (), in that for the fermionic instance of this problem (or, equivalently, a spin model), there exist constants such that

 E(^ϱI(t))≥c2t+c3, (48)

for and and , again for . In other words, for a larger and larger block size, one encounters a linear increase of the entanglement entropy of that block. In both Ref. Relax () and Schuch () the local states converge to maximal entropy states under the constraints of motion, in the latter case starting from a fermionic Gaussian state.

This means that eventually, one will have to use matrix-product states in the DMRG approach that make use of exponentially many parameters in time. This linear increase, provably correct in the above cases, is consistent with a wealth of numerical findings in quenched systems, as well as with arguments using conformal field theory Daley (); Entanglement3 (); Entanglement4 (); Calabrese ().

Knowing that we have to asymptotically deal with a linear increase, we have plotted the entanglement entropy as obtained from a t-DMRG approach. We depict here the entanglement entropy as a function of time in a slightly different geometrical setting, which still has implications on the approximatability of a state with a matrix-product state in DMRG. This is the setting of the symmetrically bisected half-chain, for and . The above linear lower bound is also rigorously true for this geometry (at least for finite-dimensional constituents), whereas the linear upper bound is certainly expected to be valid.

We numerically find an initial sublinear regime in time, followed by a linear regime, see Fig. 20. The linear regime is plausible when considering the linear propagation of the excitations in the lattice. This sequence of a sublinear regime followed by a linear one is plausible when considering the observation that eventually, excitations travel through the lattice at a finite speed. This is also consistent with the above linear upper and lower bounds. This increase also eventually limits the time to which the time evolution can be traced using t-DMRG. It is interesting to observe that the entanglement growth depends only very weakly on interaction strength , whereas the (low-energy) speed of sound depends quite strongly on KollathDensity (). This reflects the fact that after the quench, we are dealing with a very high energy state of the new Hamiltonian, where excitations are expected to be essentially ballistic with a propagation velocity , such that the conventional low-energy speed of sound is not relevant.

## X Relationship to kinematical approaches

We only very briefly touch this issue here. In this work, similarly to Refs. Relax (); Tegmark (); Barthel (), we have considered local relaxation to an apparent equilibrium state. The overall system undergoes time evolution under a local Hamiltonian , whereas subsystems appear relaxed. This observation is in an interesting relationship to kinematical approaches Jaynes (); Page (); Page2 (); Kine1 (). Indeed, if one draws a random pure state of a large system, this will be—in the limit of a large environment—almost always maximally mixed. Random here means, drawn from the unitarily invariant measure, so the Haar measure Page (); Page2 (); Kine1 (); Kine2 () (or a measure on the energy shell for Gaussian states Kine3 ()). So one could argue that “most states are locally relaxed”. This might render a relaxation to high entropy states plausible. Yet, ironically, the image of the positive time axis under the time evolution corresponds, of course, to a one-dimensional manifold in state space: So in order to show that relaxation follows from such a kinematical argument, one has to demonstrate that for a given local Hamiltonian, time evolution ensures that the state of the global system stays within the typical subspace. This is an interesting programme and the proof constitutes an interesting challenge in case of interacting systems, which has not yet been completed. In case of free systems as in Ref. Relax (), the findings may indeed be interpreted in this way.

## Xi Summary

In this work, we have introduced a setting in which local relaxation in quantum many-body systems can be probed, without the need of actually addressing single sites: This is the setting of a Bose-Hubbard model, in which preparation and read-out can be done with period- translational invariance, as can be achieved by exploiting optical superlattices. We have approached this idea by deriving analytical expressions for the relevant quantities in case of the free limits of and in the Bose-Hubbard model, as well as using a systematic t-DMRG approach in the time-dependent interacting case. For the case, we presented in detail a true convergence of subsystems to the maximum entropy state in -norm. In several ways, the interacting setting reminded of the non-interacting case, certainly concerning the mechanism of relaxation. The time-scales are different, however, showing a faster relaxation compared to the inverse square-root dependence. Also, the dependence on the interaction strength reflects the same dependence in the corresponding Gibbs state of the Hamiltonian the system is quenched to.

In this way, signatures of local relaxation can be measured with present technology: Local densities are found to relax to those of a quasi-thermal state on well-defined time scales, often with a similar to inverse square-root time dependence. More sophisticated measurements would reveal correlators and density-density correlations, again showing local relaxation in a characteristic fashion. Hence, by exploiting period- translational invariance, one has a tool at hand that effectively can be viewed as if one looked at a small local subsystem, showing all signatures of the local relaxation. In turn, the quasi-momentum distribution is a global quantity, one that shows in the free limit of no relaxation at all, and merely detects the boundary conditions quickly.

Hence, this is a quite exciting situation that both the presence of the memory of the initial condition could be experimentally probed—in the absence of relaxation in “too global quantities”—as well as local relaxation: Locally, the system “appears relaxed”, for very long times, until recurrences become relevant. The technology offered by recent advances in the use of optical superlattices should hence open up the way to experimentally access such fundamental and old questions as the mechanism of local relaxation in quantum many-body systems.

## Xii Acknowledgements

We would like to thank T. Barthel, P. Calabrese, R. Fazio, M.B. Hastings, C. Kollath, T.J. Osborne, B. Paredes, and M.M. Wolf for discussions related to the subject of local relaxation in quenched quantum many-body systems. We specifically thank I. Bloch for several fruitful discussions on non-equilibrium dynamics and optical superlattices, and for the kind hospitality in Mainz. This work has been supported by the EU (QAP, COMPAS), the DFG (FOR 801), the EPSRC, the EURYI, and Microsoft Research.

## Xiii Appendix: Proofs

We largely follow Ref. Relax () to show local relaxation. For the special situation at hand—one spatial dimension and — the proofs simplify significantly and we state them here for completeness.

### xiii.1 Preliminaries

The bosonic operators evolve in time according to (this can be shown by solving Heisenberg’s equation of motion or applying the Baker-Hausdorff formula)

 ^bi(t)=eit^H^bie−it^H=L∑j=1Vi,j^bj,V(L,t)=e−itH, (49)

where the entries of the Hamilton matrix are given by . For the translationally invariant setting one has

 Vi,j(L,t)=1LL∑k=1e−itλke2πik(i−j)/L, (50)

where

 λk=−2cos(2πk/L) (51)

are the eigenvalues of .

### xiii.2 Proof of local relaxation

In this subsection, we rigorously show that for large times and large system size, the state of a subblock becomes exactly a maximum entropy state. The key ingredients to this proof are Lieb-Robinson ideas and a central limit-type theorem. The proof is not too technically involved, but also far from being trivial, and we present it for completeness. It is, after all, quite astonishing that one does arrive at maximum entropy states without a time average. This proof largely follows Ref. Relax (), adapted to our situation of an alternating sequence of bosons and no bosons per site as initial condition.

We can now no longer merely think in terms of moments, as we want to prove convergence of the state itself. Instead of studying the quantum state, we investigate its characteristic function in phase space. Pointwise convergence of the characteristic function implies convergence in trace norm for the state. For subsets , the local state on is given by a partial trace

 (52)

and the corresponding characteristic function to represent the state is defined as

 χ(α;t)=tr[^ϱI∏i∈I^Di(αi)],^Di=eαi^b†i−α∗i^bi. (53)

Here the are the complex phase space coordinates and

 |ψ(t)⟩=e−it^H|1,0…,1,0⟩ (54)

is the time evolved state vector. We aim at showing that the state of any subblock of consecutive sites locally relaxes. For any such block, we find from Eq. (49) that

 χ=⟨1,0,…,1,0|∏i∈L^Di(βi)|1,0,…,1,0⟩, (55)

where we defined

 βi(t)=s∑j=1αjV∗j,i(t). (56)

We thus find the following explicit form of the characteristic function of :

 χ=L/2∏i=1⟨0|^D2i(β2i)|0⟩L/2∏i=1⟨1|^D2i−1(β2i−1)|1⟩=e−β†β/2L/2∏i=1(1−|β2i−1|2), (57)

where we made use of the explicit matrix elements of the displacement operator in the Fock basis. Now, follows from unitarity of and we proceed by proving that the above product converges pointwise in to .

#### xiii.2.1 The causal cone

The key ingredient for most of what follows is the fact that the entries become arbitrarily small for sufficiently large and . This can be seen as follows: The may be thought of as Riemann-sum approximation to the integral

 12π∫2π0dϕe2iJtcos(ϕ)eilϕ=ilJl(2Jt), (58)

where is a Bessel function of the first kind, for which one has for all and all . The error involved in this approximation is

 ∣∣Vl−ilJl(2Jt)∣∣≤π(l+2Jt)L, (59)

i.e., we have the bound

 |Vl|≤π(l+2Jt)L+(2Jt)−1/3, (60)

which converges to zero if we fix , let and then . However, we need a bound on the entries of for all . To this end we complement the above bound with a Lieb-Robinson type bound on the influence of sites with large : As for , we have for , i.e.,

 Vi,j=∑n≥d(it)nn!(Hn)i,j. (61)

Now, for any matrix one has , where indicates the operator norm. Furthermore, . We thus find

 ∣∣Vi,j∣∣≤∑n≥d(6Jt)nnn≤∑n≥d(6Jt)ndn=(6Jt/d)d1−6Jt/d, (62)

independent of . Hence, matrix entries with are exponentially suppressed in , defining a ”causal cone” as the influence of sites with is exponentially small in . For given we have if , where is given by the solution to . A crude bound on may be obtained from noting that

 (6Jt/d)d1−6Jt/d≤6Jt/d1−6Jt/d. (63)

Combining this with the bound for inside the cone from above, we find for given and all that

 |Vi,j|<ε for all t:4ε2

The entries of are thus arbitrarily small for sufficiently large and . In particular .

#### xiii.2.2 Central limit-type argument

We have from the previous section that becomes arbitrarily small for sufficiently and . Thus, for given the become arbitrarily small, in particular we may assume . Now, for one has , i.e.,

 L/2∑i=1∣∣log(1−|β2i−1|2)+|β2i−1|2∣∣≤L/2∑i=1|β2i−1|4≤supj|βj|2L/2∑