Probing local relaxation of cold atoms in optical superlattices
Abstract
In the study of relaxation processes in coherent nonequilibrium 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 BoseHubbard 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 timedependent densitymatrix renormalization (tDMRG) 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 noninteracting 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 manybody systems exists, although the explicit calculation of manybody equilibrium properties is an often unsolved challenge. The situation is much less satisfactory when it comes to the study of the nonequilibrium properties of manybody 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 nonequilibrium. One of the issues taking center stage is whether quantum manybody systems in nonequilibrium 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 manybody systems has enjoyed quite a renaissance recently [1–36].
A specific setting of coherent nonequilibrium 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 BoseHubbard Hamiltonian corresponding to the deep superfluid phase—treated as a noninteracting 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
(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 singleparticle 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 nonintegrable system, numerically two distinct nonequilibrium 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 manybody 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 timedependent densitymatrix renormalizationgroup (tDMRG) 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 nontrivial, are on experimental time scales essentially closed quantum systems, and therefore show coherent quantum dynamics. Systems are also sufficiently large to show nontrivial manybody 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., siteresolved) 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 Rubidium87 atoms in a period optical superlattice geometry: Two standingwave 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 singleoccupied 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 welldefined energies. Biasing, say, the odd sites relative to the even sites by an energy in excess of the separation energy of the bandseparation energy, oddsite particles are reloaded into the higher band of the even sites, whereas the evensite particles stay in the lower band. A standard timeofflight mapping then shows the evensite particles in the first Brillouin zone, the oddsite 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 nearestneighbor 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 onsite densities. Specifically, for we find that
(2) when letting evolve under the free Hamiltonian
(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 nonvanishing ), 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 timeresolved observation of onsite densities. Similar timeresolved measurements have recently been performed in optical superlattices TimeResolved ().
As shown in Fig. 1, we propose to start from a twoperiodic 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
(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 BoseHubbard Hamiltonian
(5) 
where and are the standard interaction and hopping parameters of the BoseHubbard 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 BoseHubbard 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 nonequilibrium 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 nonlocal 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
(6) 
and similarly for odd sites.
If one has experimental access to the variances , densitydensity correlations (see also Ref. NoiseNoise ()) between all even and all odd sites may be obtained through
(7)  
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 quasimomentum distribution defined as
(8) 
is also accessible as a global quantity, from timeofflight measurements.
We now proceed as follows. In order to obtain analytical results we set and respectively, which are both essentially noninteracting 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 tDMRG method. For the timescales considered it is an effectively quasiexact 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 nonGaussian product initial conditions in Ref. Relax (). In these cases, the reduced state of a subsystem
(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 1norm 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 nonGaussian 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 noninteracting 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 timeevolutions 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 Noninteracting bosons:
In the translationally invariant case of noninteracting bosons, the Hamiltonian can be diagonalized by Fourier transforming to new bosonic operators to yield
(10) 
where
(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 timeindependent, whereas the operators evolve in time as
(12) 
Straightforward algebra then yields the exact timeevolution of twopoint correlations, see Appendix,
(13) 
In the thermodynamic limit this turns into
(14) 
where is a Bessel function of the first kind.
From the above expression, one can derive the quasimomentum distribution for finite , see Appendix. For it is constant, , while one finds
(15) 
for all other .
With slightly more effort, densitydensity correlations as in Eq. (7) emerge as (see Appendix)
(16)  
which, as they are local quantities, relax for large systems and long times. For large one finds for the global densitydensity correlator (see Appendix)
(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
(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:
In the limit , the interaction manifests itself exclusively in that bosonic occupation numbers are upper bounded by . This leads to a wellknown 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
(19) 
For initial states that are supported on span one has . Furthermore , where
(20) 
where the latter projector defined on each site . This, in turn, means that
(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 BoseHubbard model, which is equivalent to a free fermionic spinless system.
This mapping to noninteracting spinless fermions is done through the familiar JordanWigner transformation
(22) 
Under this transformation the initial state vectors turn into
(23) 
and the Hamiltonian reads upon a Fourier transformation to new fermionic operators
(24) 
with as in Eq. (11) and the timeevolution of the operators is given by
(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 JordanWigner 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 longerranged twopoint correlations, including structure functions, as well as densitydensity correlations results differ. Densitydensity correlations are given by (see XIII.4)
(26) 
where is as in Eq. (13). For large the global densitydensity correlator is then given by (see Appendix)
(27) 
iii.3 Interacting softcore bosons: Finite
In this case, we are no longer facing an integrable model. In order to study the relaxation dynamics, we will therefore turn to the timedependent variant VidalTEBD (); Daley (); WhiteFeiguin () of the DMRG (densitymatrix renormalization group) method White (); RMP (). This method allows to follow the coherent timeevolution 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 noninteracting results and to read off longtime 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
(28) 
Odd and evensite 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. 26 (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 noninteracting 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 powerlaw 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 powerlaw 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 nonzero 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 nearestneighbor 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 buildup of correlations in time starting from the uncorrelated initial state becomes visible.
In the limiting free cases, identical results are found:
(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 BoseHubbard Hamiltonian : We hence look at local correlations in the Gibbs state
(30) 
with . We will take
(31) 
as our unperturbed Hamiltonian including interactions and look at leading orders in the hopping term
(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)
(34)  
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 timedependent 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
(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
(36) 
In the case of one finds
(37) 
which means that
(38) 
We have also numerically investigated the relaxation behavior of longerranged 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 noninteracting 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:

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.

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
(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 BoseHubbard 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 lowenergy speed of sound need not be relevant here.

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 noninteracting 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 Quasimomentum 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,
(40) 
obtainable from timeofflight 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 quasimomentum distribution will be very little timedependent, 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 quasimomentum 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 quasimomentum 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 onsite interaction ? Again, the behavior reminds of the situation in the free case of . The quantity is too nonlocal 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 quasimomentum 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 quasimomentum distribution
(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 quasimomentum 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 quasimomentum 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 quasimomentum distribution.
Viii Global densitydensity correlator
One potentially experimentally accessible global correlator is the connected densitydensity correlator
(43) 
which we show for different in Fig. 19. Numerics indicates that this quantity looks effectively relaxed after some time as longrange correlators contributing will be small. For the free case of , one analytically finds for large (see Appendix)
(44) 
which relaxes to for large times. For hardcore bosons the global densitydensity correlator for large is given by
(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 finitedimensional constituents, starting from a product state, the entanglement entropy of a subblock grows in onedimensional systems in at most as
(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
(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
(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 matrixproduct 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 tDMRG 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 matrixproduct state in DMRG. This is the setting of the symmetrically bisected halfchain, for and . The above linear lower bound is also rigorously true for this geometry (at least for finitedimensional 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 tDMRG. It is interesting to observe that the entanglement growth depends only very weakly on interaction strength , whereas the (lowenergy) 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 lowenergy 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 onedimensional 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 manybody systems can be probed, without the need of actually addressing single sites: This is the setting of a BoseHubbard model, in which preparation and readout 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 BoseHubbard model, as well as using a systematic tDMRG approach in the timedependent 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 noninteracting case, certainly concerning the mechanism of relaxation. The timescales are different, however, showing a faster relaxation compared to the inverse squareroot 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 quasithermal state on welldefined time scales, often with a similar to inverse squareroot time dependence. More sophisticated measurements would reveal correlators and densitydensity 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 quasimomentum 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 manybody 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 manybody systems. We specifically thank I. Bloch for several fruitful discussions on nonequilibrium 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 BakerHausdorff formula)
(49) 
where the entries of the Hamilton matrix are given by . For the translationally invariant setting one has
(50) 
where
(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 LiebRobinson ideas and a central limittype 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
(53) 
Here the are the complex phase space coordinates and
(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
(55) 
where we defined
(56) 
We thus find the following explicit form of the characteristic function of :
(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 Riemannsum approximation to the integral
(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
(59) 
i.e., we have the bound
(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 LiebRobinson type bound on the influence of sites with large : As for , we have for , i.e.,
(61) 
Now, for any matrix one has , where indicates the operator norm. Furthermore, . We thus find
(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
(63) 
Combining this with the bound for inside the cone from above, we find for given and all that
(64) 
The entries of are thus arbitrarily small for sufficiently large and . In particular .
xiii.2.2 Central limittype 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.,