# Out-of-Equilibrium Dynamics of the Bose-Hubbard Model

Authors: Malcolm P. Kennett
Journal: ISRN Condensed Matter Physics
 (2013)

---

## Abstract

The Bose-Hubbard model is the simplest model of interacting bosons on a lattice. It has recently been the focus of much attention due to the realization of this model with cold atoms in an optical lattice. The ability to tune parameters in the Hamiltonian as a function of time in cold atom systems has opened up the possibility of studying out-of-equilibrium dynamics, including crossing the quantum critical region of the model in a controlled way. In this paper, I give a brief introduction to the Bose Hubbard model, and its experimental realization and then give an account of theoretical and experimental efforts to understand out-of-equilibrium dynamics in this model, focusing on quantum quenches, both instantaneous and of finite duration. I discuss slow dynamics that have been observed theoretically and experimentally for some quenches from the superfluid phase to  the Mott insulating phase and the picture of two timescales, one for fast local equilibration and another for slow global equilibration, that appears to characterize this situation. I also discuss the theoretical and experimental observation of the Lieb-Robinson bounds for a variety of quenches and the Kibble-Zurek mechanism in quenches from the Mott insulator to superfluid. I conclude with a discussion of open questions and future directions.


---

## Body

## 1. Introduction

The Bose-Hubbard model (BHM) is a minimal model of interacting bosons on a lattice. The original focus of work on the BHM [1] was in the context of experiments on superconductor-insulator transitions in granular superconductors [2] and Josephson junction arrays [3] and for 4He in porous media [4]. The proposal by Jaksch et al. [5] that the BHM could be realized by cold atoms in an optical lattice and the subsequent experimental demonstration of a superfluid to Mott insulator transition in this system by Greiner et al. [6] has lead the focus of work on this model to shift to cold atoms. The tunability of parameters in cold atom systems [7, 8], particularly as a function of time, naturally leads to interest in the out-of-equilibrium dynamics of the BHM especially in the vicinity of a quantum critical point [6].The out-of-equilibrium dynamics of interacting quantum systems is an area of very active research, due to the challenge of understanding such a nontrivial problem. Unlike equilibrium statistical mechanics, where there are clear prescriptions for determining the state of a system, in out-of-equilibrium systems, not only do the current parameters of the Hamiltonian play a role, but also its history. Finding general principles and widely applicable calculational approaches to obtain greater insight into these problems are important goals of  experimental and theoretical work in this field. In this context, the BHM offers up an example of a system in which there is strong interplay between theory and experiment and hope for advancing these goals, both for the BHM itself, and potentially in a broader context.In addition to the interest in understanding the out-of-equilibrium dynamics of interacting bosons, there has been much excitement about the prospects for using cold atoms in optical lattices for quantum simulations [9]. Any quantum simulator should be characterized and reliable on well-understood problems before being used to simulate poorly understood problems. The temporal evolution of quantum systems might well be a future task for a quantum simulator, and hence it is worthwhile to try to characterize the out-of-equilibrium dynamics for one of the simplest models of interacting bosons.There have been a number of extensive reviews on cold atoms in optical lattices, which include discussions of the BHM [10–12] and a recent review on the out-of-equilibrium dynamics of closed interacting quantum systems [13]; hence, one might ask whether a review of out-of-equilibrium dynamics in the BHM is required. The large number of publications on this topic in the past ten years and the current activity in the field suggest (at least to this author) that it may be valuable to provide a narrower focus to take stock of what has been learned and what directions may be fruitful to explore in the future.On the theoretical side, there are a number of challenges to studying the out-of-equilibrium dynamics of the BHM. In particular, even though there are a wide array of theoretical approaches that have been used to study the BHM, all have their drawbacks. Of analytic approaches, most are accurate only in some portion of parameter space, either for weak interactions or strong interactions, or have limitations that there are a large average number of bosons per site. Numerical approaches can be essentially exact but are limited in the size of system that can be simulated due to finite computational resources and for the same reason may be restricted in the dimensions for which they are useful. As a consequence, especially in situations such as a quantum quench, which may require a description of both weakly and strongly interacting limits of the BHM within the same calculation, it may be challenging to find methods that are fully trustworthy for all parameters in the calculation. Nevertheless, this should not discourage such efforts but does imply that the results that should generally be most trusted are those that can be derived from a variety of different methods, with differing strengths and weaknesses.On the experimental side, as mentioned above, the main focus of work on the BHM is in the context of cold atoms in optical lattices. The flexibility in engineering the dimensionality of the lattice and parameters in the BHM, as well as their temporal dependence through appropriate choices of laser configurations and parameters, makes these systems well suited to studying out-of-equilibrium dynamics. Additionally, to a very good approximation, cold atoms in optical lattices can be treated as isolated quantum systems, so the interaction of the system with an environment is not a primary consideration and the observed dynamics is intrinsic to the system. However, there is the complicating feature that cold atoms require a harmonic confining potential as well as the optical lattice potential, as opposed to the homogeneous systems often preferred by theorists. The trap introduces an inhomogeneous density profile, so that there can be multiple phases present in the same trap, which allows a broad range of parameter space to be probed in a single experiment but places constraints on the ability to observe quantum critical behaviour. Out-of-equilibrium effects in imaging methods and thermometry also complicate the interpretation of some experiments.In addition to the quantum quench protocol mentioned above, several other protocols have been considered for studying out-of-equilibrium physics in the BHM, both theoretically and experimentally. These include periodic driving in time [14–20], setting up states with finite superfluid current [21–23], adding dissipation [24, 25], and setting up states with a linear potential (artificial electric field) [26–28]. I have chosen to focus primarily on out-of-equilibrium dynamics associated with quantum quenches, which have attracted the bulk of the attention of work on out-of-equilibrium dynamics in the BHM, but wish to highlight these other efforts in the interest of completeness.This paper is structured as follows: in Section2, I review the Bose Hubbard model and discuss its phases and phase transitions, along with the effects of a harmonic trapping potential. In Section 3, I survey theoretical techniques that have been used to study equilibrium and out-of-equilibrium dynamics in the BHM, discussing some of the findings within various methods and also strengths and weaknesses of individual methods. In Section 4, I review experiments that have probed the out-of-equilibrium dynamics of the BHM, particularly those that have involved a crossing of the quantum critical region associated with the Mott insulator-superfluid phase transition. Finally, in Section 5, I conclude and discuss future directions that may be of interest.In writing this paper, I have tried to represent the wide variety of approaches and results that have been obtained on the BHM. In so doing, I have tried to include all relevant work. In an area that has seen as many publications as this one, it is easy to overlook a piece of work; I apologise to any author whose work I may have missed.


## 2. The Bose Hubbard Model

The Hamiltonian for the Bose Hubbard model takes the form(1)H^BH=-∑〈ij〉Jij(a^i†a^j+a^j†a^i)+U2∑in^i(n^i-1)-μ∑in^i,=H^J+H^0,
where a^i and a^i† are annihilation and creation operators for bosons on site i, respectively, with the usual commutation relations [a^i,a^j†]=δij, and where n^i=a^i†a^i is the number operator on site i. The Jij are the hopping parameters (the choice that will generally be considered here is that Jij=J for nearest neighbour sites and zero otherwise on a cubic lattice), U is the strength of interactions (assumed to be positive, corresponding to repulsive interactions), and μ is the chemical potential. The Hamiltonian can be written as a sum of a piece H^0 that contains only single site terms, with
(2)H^0=H^U-μN^=U2∑in^i(n^i-1)-μ∑in^i,
and H^J corresponding to the kinetic energy. The term H^U in H^0 contains on-site interactions between bosons (longer range interactions can lead to additional orderings, such as supersolidity [29], but I will consider on-site interactions exclusively). The chemical potential determines the average filling per site for a given U and J.A phase diagram which is qualitatively correct in dimensions higher than one can be determined from simple considerations and mean field theory. IfU≪J, then the model describes weakly interacting bosons hopping on a lattice, which will condense and lead to a superfluid at low temperatures. In the opposite limit of J=0, the ground state can be determined by considering the single site Hamiltonian H^0. Working in the grand canonical ensemble (results in the canonical ensemble can be obtained by a Legendre transformation), the ground state energy is minimized for n0 bosons on each site when the chemical potential satisfies n0≤μ/U≤n0+1. This phase with an integer number of bosons per site is the incompressible Mott insulator. For finite J/U, there is generically a transition between the Mott insulator and superfluid phases. It is relatively straightforward to calculate the mean field phase diagram as a function of J/U and μ/U [1, 30–32] as described in Section 2.1 which leads to the well-known Mott insulator lobes illustrated in Figure 1. Considerable effort has been expended to obtain more accurate determinations of the phase boundaries, using quantum Monte Carlo [33–40], series expansions [41–43], and the density matrix renormalization group [44, 45]. The mean field phase diagram is qualitatively similar to more exact calculations in dimensions 2 and higher but differs greatly from accurate determinations of the phase boundary in 1 dimension, as illustrated in Figure 2.Figure 1
The mean field phase diagram of the BHM illustratingμ± (9) for the first five Mott lobes. A vertical line indicates a possible range of the local chemical potential μ(r)=μ-V(r) in a trap.(a) Phase diagram of the one dimensional BHM as a function of hopping(t/U) and chemical potential (μ/U) as calculated using series expansions [41], Quantum Monte Carlo simulations [35, 36], and density matrix renormalization group (DMRG) calculations [45]. (b) The same data as in (a) but with an expanded view of the tip of the lobe where the Kosterlitz-Thouless transition takes place [44]. Adapted figure with permission from T. D. Kühner, S. R. White, and H. Monien, Phys. Rev. B 61, 12474 (2000). Copyright 2000 by the American Physical Society.
(a)
(b)In general it is not simple to write down the ground state when bothU and J are finite, but in the limits J=0 and U=0, the ground state takes the form of a particularly simple product state. When U=0, one may write the ground state as a product of identical single particle states [46]:
(3)ΨU=0(x1,…,xN)=∏i=1N‍ϕ(xi),
where ϕ(x) satisfies the Gross-Pitaevski equation, and when J=0, the ground state is a product of Fock states:
(4)|Ψ〉J=0=∏i‍|n0〉i,
where n0 is the number of bosons per site in the Mott insulator.


### 2.1. Mean Field Phase Diagram

In the simplest mean field theory [1, 31, 47], one introduces the superfluid order parameter as ψ=〈a^i†〉=〈a^i〉, then the hopping terms in the Hamiltonian simplify to
(5)a^i†a^j≃ψ(a^i†+a^j)-ψ2,
and the Hamiltonian may be approximated as a sum of single site terms:
(6)H=zJψ2+∑i[12Uni(ni-1)-μni-zJψ(a^i†+a^i)].
Take the hopping term to be a perturbation to the Mott insulator phase, which for J=0 has the ground state given in (4) when there are n0 bosons per site. One may calculate the appropriate Landau expansion of the ground state energy as a series in powers of ψ to obtain [47]
(7)Egs=a0(n0,U-,μ-)+a2(n0,U-,μ-)ψ2+a4ψ4,
where U-=U/zJ, and μ-=μ/zJ. An estimate for the phase boundary of the Mott insulator may be found from requiring that
(8)a2(n0,U-,μ-)=n0U-(n0-1)-μ-+n0+1μ--U-n0+1=0,
which yields the expressions [47]
(9)μ-±=12[U-(2n0-1)-1]±12U-2-2U-(2n0+1)+1,
with the ± referring to the upper and lower branches of the phase boundaries. The zero temperature phase diagram is illustrated in Figure 1.


### 2.2. Phase Transitions

The form of the mean field phase diagram illustrated in Figure1 is qualitatively correct for dimensions d≥2 but is progressively less accurate in lower dimensions. For the transition between the Mott insulator and superfluid at the tip of the Mott lobe, where 〈n〉=1 on both sides of the transition, mean field theory predicts (zJ/U)c≃0.17 [1, 30, 42, 43, 48], where z is the coordination of the lattice (z=2d for a cubic lattice). Quantum Monte calculations and series expansions have found the following positions for the transitions in dimensions up to 3: (J/U)c≃0.29 in d=1 [44], (J/U)c≃0.05974 in d=2 [41], and (J/U)c≃0.034 in d=3 [40].The nature of the transition from the Mott insulator to superfluid depends on dimensionality and density. At fixed density, with the average number of bosons per site an integer, then the transition in ad-dimensional BHM is in the universality class of the d+1 dimensional XY model [1, 31, 32, 49] for which the lower and upper critical dimensions are dl=1 and du=3, respectively [1]. Hence, for d=1, one expects there to be a Kosterlitz-Thouless (KT) transition, and it was observed that in the Mott insulator the gap, Δ, scales as [41, 44]
(10)Δ~e[(J/U)c-J/U]-1/2.
As illustrated in Figure 2, the first Mott lobe in the 1d BHM has a nontrivial structure in the vicinity of the KT transition point, and in fact there can be a reentrant superfluid if one crosses the transition at fixed chemical potential [44].Ind=2 and 3, one finds a scaling form for the gap in the Mott insulator of
(11)Δ~[(JU)c-JU]zv,
where for d=2, z=1, ν=0.67 [50] and for d=3, z=1 and ν=1/2 [1, 31]. For the generic transition, where the density is not fixed, for example, at fixed chemical potential, then zν=1, in d=1,2, and 3. The compressibility and the superfluid density are predicted to manifest quantum critical behaviour. Quantum Monte Carlo simulations of the 1-dimensional BHM in [33] found z≃2.04 and ν≃0.48, in good agreement with the scaling predictions of Fisher et al. [1] of z=2 and ν=1/2.


### 2.3. Cold Atoms in Optical Lattices

Advances in laser cooling of atoms that led to the achievement of Bose Einstein Condensation (BEC) were needed in order to obtain quantum degenerate cold atoms in optical lattices [11]. The localization of cold atoms in an optical lattice takes advantage of the polarizability of the atoms. In a standing wave laser field from counterpropagating beams, the atoms are trapped either at the nodes or antinodes of the laser intensity, depending on the sign of their polarizability. The form of the subsequent lattice potential, for example, for a 3-dimensional cubic lattice, is
(12)V(x,y,z)=V0[sin2(kx)+sin2(ky)+sin2(kz)].
The strength of the optical lattice V0 is usually expressed in units of the recoil energy Er=ℏ2k2/2m, where k is the wavevector of the laser. The lattice constant of the optical lattice is a=λ/2, where λ is the laser wavelength.It was realized by Jaksch et al. [5] that the combination of an optical lattice and cold bosons would allow for the realization of the Bose Hubbard model. In particular, if one starts from the Hamiltonian for interacting bosons in the presence of a lattice potential Vlattice(x) and a trap potential Vtrap(x), then the Hamiltonian for a 3-dimensional system is
(13)H=∫‍d3xψ^†(x)[-ℏ22m∇2+Vlattice(x)+Vtrap(x)]ψ^(x)+124πasℏ2m∫‍d3xψ^†(x)ψ^†(x)ψ^(x)ψ^(x),
where ψ^(x) is a boson field operator for atoms in a given internal atomic state and as is the s-wave scattering length. When bosons are cooled to lie in the lowest Bloch band of the periodic potential, the Hamiltonian can be expressed in the tight binding form of (1) by expanding the field operators in the Wannier functions ψ^(x)=∑i‍a^iw(x-xi) (the role of other Bloch bands, especially in dynamics, has been investigated by several groups [51, 52]). The parameters in the tight-binding model may be related to integrals over the Wannier functions, for which analytic expressions may be given for a deep lattice (V0≫Er) [46]
(14)J=∫‍d3xw*(x-xi)[-ℏ22m∇2+Vlattice(x)]w(x-xj)≃4πEr(V0Er)3/4e-2(V0/Er)1/2,U=4πasℏ2m∫‍d3x|w(x)|4≃8πkasEr(V0Er)3/4.
Note that the hopping parameter J is much more sensitive than the interaction parameter U to the depth of the well. This is to be expected, since the hopping is related to wavefunction overlap between different sites and will be exponentially suppressed in a deep lattice, whereas the interactions are at a single site and should be much less sensitive to the depth of the lattice. The values of J and U will generally be assumed to be spatially uniform throughout the trap, with the effects of the trapping potential being incorporated into the tight binding term by adding the term
(15)Htrap=∑iVini
to the BHM for a homogeneous system, where Vi is the strength of the trap at site i. The solutions found for the homogeneous BHM can be used to understand the behaviour of the BHM in a spatially inhomogeneous trapping potential by assigning a local chemical potential at radius r in the trap (assuming rotational symmetry)
(16)μ(r)=μ-V(r).
Within the local density approximation (LDA), the phase of the system at radius r may be determined to be that of the homogeneous system found at μ=μ(r). This implies that for a sufficiently deep trap and not too large a value of J/U, there can be concentric shells of the Mott insulator and superfluid phases, as follows from the schematic line μ(r) shown in Figure 1. Quantum Monte Carlo (QMC) simulations [53] of bosons in a quadratic trap (where there is no LDA assumption) also show concentric domains of superfluid and Mott insulator, which, due to the incompressible nature of the Mott phase, leads to steps in the density as a function of radius, giving rise to the so-called wedding cake structure [53, 54]. The presence of these density and local chemical potential variations within a trap complicates efforts to relate the properties of the homogeneous and trapped BHM, respectively [55].QMC simulations have proven to be a very useful tool to obtain results that are in principle exact for both the homogeneous BHM and bosons on an optical lattice in a trap. Unlike many other theoretical methods, finite temperature effects both in the homogeneous system and in traps have also been taken into account. Mahmud et al. [39] used QMC simulations to note that the presence of a trap affects the correlations in the superfluid at finite temperature. At zero temperature, the correlations for a two-dimensional system in a trap were observed to be intermediate between those expected for a uniform one-dimensional superfluid and those of a uniform two-dimensional superfluid. At finite temperatures, the correlations were observed to be similar to those expected in one dimension, but at a lower temperature. Direct comparisons between the results of QMC simulations and experiments [39, 56–58] have indicated good agreement for equilibrium systems. In a particularly detailed recent study, Trotzky et al. [56] compared measurements of the critical temperature for the finite temperature phase transition between a normal fluid and the superfluid phase as a function of U/J for a 3-dimensional system of bosons in an optical lattice with QMC simulations of the same system. Good agreement between the measured and calculated values was found for U/J≲20 with less satisfactory agreement for deeper lattices (stronger interactions). This and previous work confirm that QMC is an important tool for calculating the equilibrium properties of the BHM and of bosons in optical lattices, giving a baseline to which out-of-equilibrium dynamics can be compared.


#### 2.3.1. Quantum Criticality in a Trap

As discussed by Zhang et al. [59], there are two ways that quantum criticality can manifest itself that are accessible in a cold atom experiment. One possibility is that it can be seen in the scaling of thermodynamic quantities [31, 60, 61]. The other is that it may be seen via the dynamic passage through a quantum critical region. A signature of quantum criticality is the diverging of the correlation length at the critical point. For cold atoms in a trap, the size of the atom cloud in the trap is the largest lengthscale in the problem, leading to a cutoff in the correlation length. The inhomogeneity in the local chemical potential also limits the spatial extent of different phases. At higher temperatures the coherence length ξ shrinks to be less than the cloud size so the role of the trap is lessened [60]. Recently, Hazzard and Mueller [60] and Zhou and Ho [61] proposed schemes to observe quantum critical scaling in a trap. Hazzard and Mueller suggested comparing the density and compressibility by writing the density in the form n=n0+nu, and the compressibility as κ=κ0+κu, where n0 and κ0 are the zero temperature values of density and compressibility at given μ, T, and U for J=0. Then a plot of κuT1-d/z against nuT-d/z should collapse onto a single curve. On the other hand, Zhou and Ho suggested comparing the density and compressibility to the chemical potential, also separating out regular and singular contributions to n and κ. Unlike Hazzard and Mueller’s approach, which uses only local observables, Zhou and Ho’s approach requires a global quantity, the equilibrium chemical potential. A drawback for both methods is that at higher temperatures there can be overlapping critical regimes that affect the scaling of thermodynamic quantities [60]. For example, as illustrated in Figure 3, at temperatures above the zero temperature Mott insulator, along the line where there is particle-hole symmetry, the particle dilute Bose gas (p-DBG) and hole DBG (h-DBG) critical regions overlap and the resulting critical behaviour is that of the finite density O(2) rotor model. Such overlapping critical regions may complicate attempts to extract critical exponents from experiment. An additional complication that may affect future experiments was pointed out by Pollet et al. [62], that a finite correlation length can lead to violations of the LDA as different radii in the trap cannot be treated as independent. They estimated this effect to be not currently observable within experimental error bars. However, there may be ways to circumvent this issue by using finite size scaling [61].(a) An illustration of universality classes of the BHM [60]. Transitions via paths (1) and (2) are described by the dilute Bose gas (DBG) universality class, and path (3) is described by the O(2) universality class. The shaded regions indicate where each universality class applies: blue indicates DBG physics, orange indicates O(2) physics, and the green region is governed by the finite density O(2) universality class. (b) The schematic finite temperature phase diagram along path (4) [60]. Adapted Figure with permission from K. R. A. Hazzard and E. J. Mueller, Phys. Rev. A 84, 013604 (2011). Copyright 2011 by the American Physical Society.
(a)
(b)A manifestation of quantum criticality more in line with the theme of this paper is in the dynamic passage through a quantum critical point [59]. This should leave signatures in density and entropy currents and may also give access to the Kibble-Zurek mechanism (KZM) [63, 64]. The KZM predicts defect formation after a system crosses a second-order thermodynamic phase transition and has been generalized to quantum phase transitions [65–67]. If one focuses on a time-dependent coupling g(t) which is close to the critical coupling gc, then close to the critical point the gap Δ scales as Δ∝|g(t)-gc|zν. When g(t) gets sufficiently close to gc, then Δ becomes too small to maintain adiabaticity at some time t=t*. At this timescale, the energy scale for excitations will be Δ*~|g˙(t*)|zν/(zν+1). The correlation length at time t* will be ξ*~Δ*-1/z, which implies a density of excitations of nex~ξ*-d~|g˙(t*)|dz/(zν+1) [66]. This derivation of the number density of excitations after a quench assumes a picture in which there is adiabatic evolution up to some distance |g(t*)-gc| from the transition at which point the relaxation time diverges and the defects are frozen until the system has passed through the critical region. As recently pointed out [68, 69] for a sufficiently slow quench, there can still be evolution of the defect configuration as the system passes through the critical region, which can modify the Kibble-Zurek predictions and may be relevant to slow quenches in the BHM. These considerations are all for a uniform system and do not consider the possible presence of a trap.


### 2.4. Other Proposals to Realize the BHM

In addition to cold atoms, there have also been proposals to realize the BHM and the Mott insulator to superfluid transition in photonic [70, 71] and polaritonic systems [72–75].


## 2.1. Mean Field Phase Diagram

In the simplest mean field theory [1, 31, 47], one introduces the superfluid order parameter as ψ=〈a^i†〉=〈a^i〉, then the hopping terms in the Hamiltonian simplify to
(5)a^i†a^j≃ψ(a^i†+a^j)-ψ2,
and the Hamiltonian may be approximated as a sum of single site terms:
(6)H=zJψ2+∑i[12Uni(ni-1)-μni-zJψ(a^i†+a^i)].
Take the hopping term to be a perturbation to the Mott insulator phase, which for J=0 has the ground state given in (4) when there are n0 bosons per site. One may calculate the appropriate Landau expansion of the ground state energy as a series in powers of ψ to obtain [47]
(7)Egs=a0(n0,U-,μ-)+a2(n0,U-,μ-)ψ2+a4ψ4,
where U-=U/zJ, and μ-=μ/zJ. An estimate for the phase boundary of the Mott insulator may be found from requiring that
(8)a2(n0,U-,μ-)=n0U-(n0-1)-μ-+n0+1μ--U-n0+1=0,
which yields the expressions [47]
(9)μ-±=12[U-(2n0-1)-1]±12U-2-2U-(2n0+1)+1,
with the ± referring to the upper and lower branches of the phase boundaries. The zero temperature phase diagram is illustrated in Figure 1.


## 2.2. Phase Transitions

The form of the mean field phase diagram illustrated in Figure1 is qualitatively correct for dimensions d≥2 but is progressively less accurate in lower dimensions. For the transition between the Mott insulator and superfluid at the tip of the Mott lobe, where 〈n〉=1 on both sides of the transition, mean field theory predicts (zJ/U)c≃0.17 [1, 30, 42, 43, 48], where z is the coordination of the lattice (z=2d for a cubic lattice). Quantum Monte calculations and series expansions have found the following positions for the transitions in dimensions up to 3: (J/U)c≃0.29 in d=1 [44], (J/U)c≃0.05974 in d=2 [41], and (J/U)c≃0.034 in d=3 [40].The nature of the transition from the Mott insulator to superfluid depends on dimensionality and density. At fixed density, with the average number of bosons per site an integer, then the transition in ad-dimensional BHM is in the universality class of the d+1 dimensional XY model [1, 31, 32, 49] for which the lower and upper critical dimensions are dl=1 and du=3, respectively [1]. Hence, for d=1, one expects there to be a Kosterlitz-Thouless (KT) transition, and it was observed that in the Mott insulator the gap, Δ, scales as [41, 44]
(10)Δ~e[(J/U)c-J/U]-1/2.
As illustrated in Figure 2, the first Mott lobe in the 1d BHM has a nontrivial structure in the vicinity of the KT transition point, and in fact there can be a reentrant superfluid if one crosses the transition at fixed chemical potential [44].Ind=2 and 3, one finds a scaling form for the gap in the Mott insulator of
(11)Δ~[(JU)c-JU]zv,
where for d=2, z=1, ν=0.67 [50] and for d=3, z=1 and ν=1/2 [1, 31]. For the generic transition, where the density is not fixed, for example, at fixed chemical potential, then zν=1, in d=1,2, and 3. The compressibility and the superfluid density are predicted to manifest quantum critical behaviour. Quantum Monte Carlo simulations of the 1-dimensional BHM in [33] found z≃2.04 and ν≃0.48, in good agreement with the scaling predictions of Fisher et al. [1] of z=2 and ν=1/2.


## 2.3. Cold Atoms in Optical Lattices

Advances in laser cooling of atoms that led to the achievement of Bose Einstein Condensation (BEC) were needed in order to obtain quantum degenerate cold atoms in optical lattices [11]. The localization of cold atoms in an optical lattice takes advantage of the polarizability of the atoms. In a standing wave laser field from counterpropagating beams, the atoms are trapped either at the nodes or antinodes of the laser intensity, depending on the sign of their polarizability. The form of the subsequent lattice potential, for example, for a 3-dimensional cubic lattice, is
(12)V(x,y,z)=V0[sin2(kx)+sin2(ky)+sin2(kz)].
The strength of the optical lattice V0 is usually expressed in units of the recoil energy Er=ℏ2k2/2m, where k is the wavevector of the laser. The lattice constant of the optical lattice is a=λ/2, where λ is the laser wavelength.It was realized by Jaksch et al. [5] that the combination of an optical lattice and cold bosons would allow for the realization of the Bose Hubbard model. In particular, if one starts from the Hamiltonian for interacting bosons in the presence of a lattice potential Vlattice(x) and a trap potential Vtrap(x), then the Hamiltonian for a 3-dimensional system is
(13)H=∫‍d3xψ^†(x)[-ℏ22m∇2+Vlattice(x)+Vtrap(x)]ψ^(x)+124πasℏ2m∫‍d3xψ^†(x)ψ^†(x)ψ^(x)ψ^(x),
where ψ^(x) is a boson field operator for atoms in a given internal atomic state and as is the s-wave scattering length. When bosons are cooled to lie in the lowest Bloch band of the periodic potential, the Hamiltonian can be expressed in the tight binding form of (1) by expanding the field operators in the Wannier functions ψ^(x)=∑i‍a^iw(x-xi) (the role of other Bloch bands, especially in dynamics, has been investigated by several groups [51, 52]). The parameters in the tight-binding model may be related to integrals over the Wannier functions, for which analytic expressions may be given for a deep lattice (V0≫Er) [46]
(14)J=∫‍d3xw*(x-xi)[-ℏ22m∇2+Vlattice(x)]w(x-xj)≃4πEr(V0Er)3/4e-2(V0/Er)1/2,U=4πasℏ2m∫‍d3x|w(x)|4≃8πkasEr(V0Er)3/4.
Note that the hopping parameter J is much more sensitive than the interaction parameter U to the depth of the well. This is to be expected, since the hopping is related to wavefunction overlap between different sites and will be exponentially suppressed in a deep lattice, whereas the interactions are at a single site and should be much less sensitive to the depth of the lattice. The values of J and U will generally be assumed to be spatially uniform throughout the trap, with the effects of the trapping potential being incorporated into the tight binding term by adding the term
(15)Htrap=∑iVini
to the BHM for a homogeneous system, where Vi is the strength of the trap at site i. The solutions found for the homogeneous BHM can be used to understand the behaviour of the BHM in a spatially inhomogeneous trapping potential by assigning a local chemical potential at radius r in the trap (assuming rotational symmetry)
(16)μ(r)=μ-V(r).
Within the local density approximation (LDA), the phase of the system at radius r may be determined to be that of the homogeneous system found at μ=μ(r). This implies that for a sufficiently deep trap and not too large a value of J/U, there can be concentric shells of the Mott insulator and superfluid phases, as follows from the schematic line μ(r) shown in Figure 1. Quantum Monte Carlo (QMC) simulations [53] of bosons in a quadratic trap (where there is no LDA assumption) also show concentric domains of superfluid and Mott insulator, which, due to the incompressible nature of the Mott phase, leads to steps in the density as a function of radius, giving rise to the so-called wedding cake structure [53, 54]. The presence of these density and local chemical potential variations within a trap complicates efforts to relate the properties of the homogeneous and trapped BHM, respectively [55].QMC simulations have proven to be a very useful tool to obtain results that are in principle exact for both the homogeneous BHM and bosons on an optical lattice in a trap. Unlike many other theoretical methods, finite temperature effects both in the homogeneous system and in traps have also been taken into account. Mahmud et al. [39] used QMC simulations to note that the presence of a trap affects the correlations in the superfluid at finite temperature. At zero temperature, the correlations for a two-dimensional system in a trap were observed to be intermediate between those expected for a uniform one-dimensional superfluid and those of a uniform two-dimensional superfluid. At finite temperatures, the correlations were observed to be similar to those expected in one dimension, but at a lower temperature. Direct comparisons between the results of QMC simulations and experiments [39, 56–58] have indicated good agreement for equilibrium systems. In a particularly detailed recent study, Trotzky et al. [56] compared measurements of the critical temperature for the finite temperature phase transition between a normal fluid and the superfluid phase as a function of U/J for a 3-dimensional system of bosons in an optical lattice with QMC simulations of the same system. Good agreement between the measured and calculated values was found for U/J≲20 with less satisfactory agreement for deeper lattices (stronger interactions). This and previous work confirm that QMC is an important tool for calculating the equilibrium properties of the BHM and of bosons in optical lattices, giving a baseline to which out-of-equilibrium dynamics can be compared.


### 2.3.1. Quantum Criticality in a Trap

As discussed by Zhang et al. [59], there are two ways that quantum criticality can manifest itself that are accessible in a cold atom experiment. One possibility is that it can be seen in the scaling of thermodynamic quantities [31, 60, 61]. The other is that it may be seen via the dynamic passage through a quantum critical region. A signature of quantum criticality is the diverging of the correlation length at the critical point. For cold atoms in a trap, the size of the atom cloud in the trap is the largest lengthscale in the problem, leading to a cutoff in the correlation length. The inhomogeneity in the local chemical potential also limits the spatial extent of different phases. At higher temperatures the coherence length ξ shrinks to be less than the cloud size so the role of the trap is lessened [60]. Recently, Hazzard and Mueller [60] and Zhou and Ho [61] proposed schemes to observe quantum critical scaling in a trap. Hazzard and Mueller suggested comparing the density and compressibility by writing the density in the form n=n0+nu, and the compressibility as κ=κ0+κu, where n0 and κ0 are the zero temperature values of density and compressibility at given μ, T, and U for J=0. Then a plot of κuT1-d/z against nuT-d/z should collapse onto a single curve. On the other hand, Zhou and Ho suggested comparing the density and compressibility to the chemical potential, also separating out regular and singular contributions to n and κ. Unlike Hazzard and Mueller’s approach, which uses only local observables, Zhou and Ho’s approach requires a global quantity, the equilibrium chemical potential. A drawback for both methods is that at higher temperatures there can be overlapping critical regimes that affect the scaling of thermodynamic quantities [60]. For example, as illustrated in Figure 3, at temperatures above the zero temperature Mott insulator, along the line where there is particle-hole symmetry, the particle dilute Bose gas (p-DBG) and hole DBG (h-DBG) critical regions overlap and the resulting critical behaviour is that of the finite density O(2) rotor model. Such overlapping critical regions may complicate attempts to extract critical exponents from experiment. An additional complication that may affect future experiments was pointed out by Pollet et al. [62], that a finite correlation length can lead to violations of the LDA as different radii in the trap cannot be treated as independent. They estimated this effect to be not currently observable within experimental error bars. However, there may be ways to circumvent this issue by using finite size scaling [61].(a) An illustration of universality classes of the BHM [60]. Transitions via paths (1) and (2) are described by the dilute Bose gas (DBG) universality class, and path (3) is described by the O(2) universality class. The shaded regions indicate where each universality class applies: blue indicates DBG physics, orange indicates O(2) physics, and the green region is governed by the finite density O(2) universality class. (b) The schematic finite temperature phase diagram along path (4) [60]. Adapted Figure with permission from K. R. A. Hazzard and E. J. Mueller, Phys. Rev. A 84, 013604 (2011). Copyright 2011 by the American Physical Society.
(a)
(b)A manifestation of quantum criticality more in line with the theme of this paper is in the dynamic passage through a quantum critical point [59]. This should leave signatures in density and entropy currents and may also give access to the Kibble-Zurek mechanism (KZM) [63, 64]. The KZM predicts defect formation after a system crosses a second-order thermodynamic phase transition and has been generalized to quantum phase transitions [65–67]. If one focuses on a time-dependent coupling g(t) which is close to the critical coupling gc, then close to the critical point the gap Δ scales as Δ∝|g(t)-gc|zν. When g(t) gets sufficiently close to gc, then Δ becomes too small to maintain adiabaticity at some time t=t*. At this timescale, the energy scale for excitations will be Δ*~|g˙(t*)|zν/(zν+1). The correlation length at time t* will be ξ*~Δ*-1/z, which implies a density of excitations of nex~ξ*-d~|g˙(t*)|dz/(zν+1) [66]. This derivation of the number density of excitations after a quench assumes a picture in which there is adiabatic evolution up to some distance |g(t*)-gc| from the transition at which point the relaxation time diverges and the defects are frozen until the system has passed through the critical region. As recently pointed out [68, 69] for a sufficiently slow quench, there can still be evolution of the defect configuration as the system passes through the critical region, which can modify the Kibble-Zurek predictions and may be relevant to slow quenches in the BHM. These considerations are all for a uniform system and do not consider the possible presence of a trap.


## 2.3.1. Quantum Criticality in a Trap

As discussed by Zhang et al. [59], there are two ways that quantum criticality can manifest itself that are accessible in a cold atom experiment. One possibility is that it can be seen in the scaling of thermodynamic quantities [31, 60, 61]. The other is that it may be seen via the dynamic passage through a quantum critical region. A signature of quantum criticality is the diverging of the correlation length at the critical point. For cold atoms in a trap, the size of the atom cloud in the trap is the largest lengthscale in the problem, leading to a cutoff in the correlation length. The inhomogeneity in the local chemical potential also limits the spatial extent of different phases. At higher temperatures the coherence length ξ shrinks to be less than the cloud size so the role of the trap is lessened [60]. Recently, Hazzard and Mueller [60] and Zhou and Ho [61] proposed schemes to observe quantum critical scaling in a trap. Hazzard and Mueller suggested comparing the density and compressibility by writing the density in the form n=n0+nu, and the compressibility as κ=κ0+κu, where n0 and κ0 are the zero temperature values of density and compressibility at given μ, T, and U for J=0. Then a plot of κuT1-d/z against nuT-d/z should collapse onto a single curve. On the other hand, Zhou and Ho suggested comparing the density and compressibility to the chemical potential, also separating out regular and singular contributions to n and κ. Unlike Hazzard and Mueller’s approach, which uses only local observables, Zhou and Ho’s approach requires a global quantity, the equilibrium chemical potential. A drawback for both methods is that at higher temperatures there can be overlapping critical regimes that affect the scaling of thermodynamic quantities [60]. For example, as illustrated in Figure 3, at temperatures above the zero temperature Mott insulator, along the line where there is particle-hole symmetry, the particle dilute Bose gas (p-DBG) and hole DBG (h-DBG) critical regions overlap and the resulting critical behaviour is that of the finite density O(2) rotor model. Such overlapping critical regions may complicate attempts to extract critical exponents from experiment. An additional complication that may affect future experiments was pointed out by Pollet et al. [62], that a finite correlation length can lead to violations of the LDA as different radii in the trap cannot be treated as independent. They estimated this effect to be not currently observable within experimental error bars. However, there may be ways to circumvent this issue by using finite size scaling [61].(a) An illustration of universality classes of the BHM [60]. Transitions via paths (1) and (2) are described by the dilute Bose gas (DBG) universality class, and path (3) is described by the O(2) universality class. The shaded regions indicate where each universality class applies: blue indicates DBG physics, orange indicates O(2) physics, and the green region is governed by the finite density O(2) universality class. (b) The schematic finite temperature phase diagram along path (4) [60]. Adapted Figure with permission from K. R. A. Hazzard and E. J. Mueller, Phys. Rev. A 84, 013604 (2011). Copyright 2011 by the American Physical Society.
(a)
(b)A manifestation of quantum criticality more in line with the theme of this paper is in the dynamic passage through a quantum critical point [59]. This should leave signatures in density and entropy currents and may also give access to the Kibble-Zurek mechanism (KZM) [63, 64]. The KZM predicts defect formation after a system crosses a second-order thermodynamic phase transition and has been generalized to quantum phase transitions [65–67]. If one focuses on a time-dependent coupling g(t) which is close to the critical coupling gc, then close to the critical point the gap Δ scales as Δ∝|g(t)-gc|zν. When g(t) gets sufficiently close to gc, then Δ becomes too small to maintain adiabaticity at some time t=t*. At this timescale, the energy scale for excitations will be Δ*~|g˙(t*)|zν/(zν+1). The correlation length at time t* will be ξ*~Δ*-1/z, which implies a density of excitations of nex~ξ*-d~|g˙(t*)|dz/(zν+1) [66]. This derivation of the number density of excitations after a quench assumes a picture in which there is adiabatic evolution up to some distance |g(t*)-gc| from the transition at which point the relaxation time diverges and the defects are frozen until the system has passed through the critical region. As recently pointed out [68, 69] for a sufficiently slow quench, there can still be evolution of the defect configuration as the system passes through the critical region, which can modify the Kibble-Zurek predictions and may be relevant to slow quenches in the BHM. These considerations are all for a uniform system and do not consider the possible presence of a trap.


## 2.4. Other Proposals to Realize the BHM

In addition to cold atoms, there have also been proposals to realize the BHM and the Mott insulator to superfluid transition in photonic [70, 71] and polaritonic systems [72–75].


## 3. Theory

In Section2.1, the mean field behaviour of the BHM was reviewed. In this section, I summarize a number of approaches that have been used to study the out-of-equilibrium behaviour of the BHM and results that have been obtained within each of these approaches.There is only one way for a system to be in equilibrium, but there are many ways to prepare out-of-equilibrium states; indeed, one of the characteristics of out-of-equilibrium states is that their history matters. The class of out of equilibrium states that I focus on primarily here is states that arise as a result of a quantum quench corresponding to the dynamic traversal of a quantum phase transition. In the context of the BHM, this refers to either or both ofJ and U being time-dependent parameters, usually such that the ground states for the initial and final values of U/J are in different phases. In order for the system to be out of equilibrium, the change of U/J in time must be sufficiently quick that the system is not able to evolve adiabatically from one phase to the other. If the path crosses the critical region, then the presence of diverging timescales at the quantum critical point ensures that this will be the case. Both instantaneous quenches and a continuous change in parameters will be considered. The quantum-quench protocol has received considerable recent interest beyond the BHM [13, 66, 67, 76–85] as the resulting systems give examples of out-of-equilibrium dynamics in interacting quantum systems, a class of problem that is still under active investigation. In the context of the Bose Hubbard model, quenches have generally taken one of two forms. Either a quench from the Mott insulating phase to the superfluid phase or a quench from the superfluid phase to the Mott insulating phase. Extensive theoretical effort has been expended on the effects of time-dependent J/U in the BHM in the form of quantum quenches [66, 86–117].There are several features that are desirable for any theory of equilibrium or out-of-equilibrium dynamics of the Bose-Hubbard model. First, it should allow for both superfluid and Mott insulating ground states (the Bogoliubov theory fails in this instance). Second, a variety of methods concur on the nature of the collective excitation spectrum in the two phases. In the superfluid phase, there is a gapless sound mode that arises from phase fluctuations of the superfluid order parameter. Near the fixed density quantum critical point (QCP) there is also a gapped (Higgs) mode corresponding to amplitude fluctuations of the superfluid order parameter that has been predicted in a variety of calculations [118–122] and observed experimentally for bosons in two-and-three dimensional optical lattices [123, 124]. In the Mott phase, the collective modes are particle-hole excitations. Third, the approach should be capable of describing nonlocal correlations in space and time.The methods discussed here mostly have the features outlined above and are: the weak interactions Bogoliubov approach (Section3.1); an infinite dimensional mean field approach (Section 3.2); Gutzwiller mean field theory (Section 3.3); a method for 1/Z corrections to mean field theory (Section 3.4); a slave particle approach (Section 3.5); numerical methods, including time-dependent density matrix renormalization group (t-DMRG), time evolving block decimation (TEBD), and Exact Diagonalization (Section 3.6); fermionization approaches (Section 3.7); closed time path approaches (Section 3.8); and exact results (Section 3.9). In addition to the methods discussed here, a number of other methods, such as bosonic dynamical mean field theory (DMFT) [125–129], the variational cluster approximation (VCA) [130, 131] and Quantum Monte Carlo (QMC) have been used to study the equilibrium properties of the BHM, but to the best of my knowledge these have not been used to address out-of-equilibrium dynamics in the BHM.


### 3.1. Weak Interactions: The Bogoliubov Approach

The Bogoliubov approach to the BHM is a mean field approach based on the weak interaction limit [47], similar to the Bogoliubov approach to a weakly interacting superfluid in the absence of a lattice. This suggests working in momentum space, reflecting the extended nature of the state and the expectation of a condensate with k=0. Writing the boson annihilation operator in momentum space:
(17)a^i=1Ns∑ka^ke-ik·ri,
with Ns the number of sites, and similarly for the creation operator, one may rewrite the BHM Hamiltonian as [47]
(18)H=∑k(ϵk-μ)a^k†a^k+U2Ns∑k∑k′a^k†a^-k†a^k′a^-k′,
with ϵk=-2J∑j=1dcos(kja) on a d-dimensional cubic lattice. The approximation in the Bogoliubov approach comes from noting that the number of condensate atoms N0≫1, in which case N0=〈a^0†a^0〉≃〈a^0a^0†〉, and hence one can write N0=〈a^0†〉〈a^0〉. Thus, 〈a^0†〉=〈a^0〉=N0, where the phase is chosen so that the expectation values are real. This observation is used to replace the annihilation and creation operators with their average value plus a fluctuation:
(19)a^0†→N0+a^0†;a^0→N0+a^0.
Keeping only terms to quadratic order in the boson operators in the Hamiltonian after this substitution leads to a Hamiltonian that can be diagonalized with a Bogoliubov transformation to
(20)Heff=-12Un0N0+12∑k[ℏωk-(ϵk+Un0)]+∑kℏωkb^k†b^k,
where n0=N0/Ns, and quasiparticle excitations b^k† are related to the a^k† by
(21)(b^kb^-k†)=(ukvkvk*uk*)(a^ka^-k†),
with |uk|2-|vk|2=1 and their energies are given by
(22)ℏωk=ϵk2+2Un0ϵk.
Limitations of the Bogoliubov theory are that it is only valid when the number of atoms in the condensate N0≫1 and hence does not indicate the presence of a Mott insulator phase when the condensate vanishes [47]. This suggests that it will not be useful for giving a quantitative description of out-of-equilibrium dynamics which involve crossing the quantum phase transition from superfluid to Mott insulator. However, as pointed out by Fischer et al. [86], it is still possible to infer some features of these out-of-equilibrium dynamics using a Bogoliubov approach. Fischer et al. [86] took the hopping term in the BHM to have the form
(23)Hhop=J∑ijMija^i†a^j,
where Mij takes into account both the connectivity of the lattice and the possibility of spatially nonuniform hopping. They studied the Heisenberg equation of motion for the operator a^i (setting ℏ=1):
(24)i∂ta^i=J∑jMija^j+Un^ia^i,
focusing on large fillings n≫1, in which case 1/n may be used as a small expansion parameter. They wrote the operator a^ as
(25)a^i=(ψ0+χ^i+ζ^i)A^N^,
where A^=a^Σ(a^Σ†a^Σ)-1/2N^1/2, with a^Σ=∑ia^i, and N^=A^†A^=∑in^i. The leading term in the expansion is ψ0, which describes the condensate and is determined from the solution of the Gross-Pitaevskii equation. The quantum corrections to this mean field are split into the linear corrections, which correspond to quasiparticle excitations of the superfluid, χ^i and nonlinear corrections ζ^i. In the superfluid phase well away from the Mott insulator phase, the non-linear corrections ζ^i may be ignored, provided Utn≪1. Fischer et al. [86] decomposed the operator
(26)a^i=eiϕ^in^i,
into phase ϕ^ and number n^ operators and then related the fluctuations χ^ to phase δϕ^ and number δn^ fluctuations, respectively. The dynamical equations for these fluctuations when expanded in the basis of eigenvectors of the hopping matrix (with eigenvalues λκ where κ are generalized momenta) are
(27)∂tδn^κ=2n0Jλκδϕ^κ,∂tδϕ^κ=-(Jλκ2n0+U)δn^κ,
where U and J may both be time dependent. If U and J are both independent of time, these equations lead to the usual Bogoliubov spectrum
(28)ωκ2=J2λκ2+2n0UJλκ,
found above in (22). When U is taken to be constant and J is allowed to vary as a function of time, Fischer et al. defined a new time variable via dτκ=λκJdt, leading to an equation of motion for number fluctuations
(29)(∂2∂τκ2+1+UJ2n0λκ)δn^κ=0.
The solutions of this equation depend sensitively on the way in which J is taken to zero as a function of time. If J(t) is of the form J(t)=J0(γt)-α, then for α<2 the number fluctuations oscillate in time with decreasing amplitude and frequency; for α=2 the fluctuations decay to zero, and for faster decay of J(t), either as a power law or exponential form J(t)=J0e-γt, the fluctuations end up being frozen in at a finite value. This can be understood as α=2 providing the limit for horizon formation, in analogy to cosmological horizons. For the exponential sweep, the equation of motion for number fluctuations becomes scale invariant (i.e., independent of κ) and the dimensionless parameter governing the physics is ν=Un0/γ, the ratio of an internal energy scale Un0, to an external one, the sweep rate, γ. This allows a notion of slow (ν≫1) and fast (ν≪1) sweeps. The leading order behaviour of the number fluctuations and phase fluctuations was determined to be [86]
(30)〈δn^κ2〉≃n0(1-e-2πν2πν),〈δϕ^κ2〉≃ν(1-e-2πν2πn0)γ2t2.
A weakness of these results is that they come from extrapolation of the mean field theory for the weak coupling superfluid phase, where ψ0 is dominant over linear and non-linear fluctuations, into the strong coupling phase—they also apply only in the limit n0≫1. The approach taken by Fischer et al. to address this issue is to note that if the quench takes place sufficiently quickly such that the strong coupling regime is reached for times t≪1/nU and t≪1/γ, then the state reached via Bogoliubov dynamics can be viewed as an initial state for evolution after the quench with the Hamiltonian H^0. From these considerations, they found that the correlation function usually considered to determine the presence of off-diagonal long range order (ODLRO) takes the form
(31)〈a^i†a^j〉≃n0e-U2t2Δ2(n0),Δ2(n0)=n0(1-e-2πν2πν).
This prediction for the decay of the ODLRO with time could be tested in the momentum distribution function in time of flight experiments either in the peak at k=0, or perhaps in the visibility as a function of time.As mentioned above, Fischer et al. [86] noted the analogy between the dynamics of number fluctuations under a quantum quench and dynamics in the early universe. In the long wavelength limit, where the lattice structure is ignored, that is, for wavelengths λ≫a, the lattice spacing, the equation of motion for number fluctuations may be written as
(32)(∂∂t1J(t)∂∂t-Un0a2∇2)δn^(x,t)=0,
where the speed of sound cs2(t)=Un0a2J(t) decreases as a function of time. If this decrease is sufficiently quick, then there is a horizon beyond which excitations of a given wavelength cannot propagate, given by
(33)Δn(t)=∫t∞dt′cs(t′),
where for J(t)=J0e-γt, the integral converges and one has Δn(t)=(2J0Un0/γ)e-γt/2. This horizon shrinks as t→∞, so at a given time t, modes with wavelengths shorter than Δ will oscillate, but longer wavelength modes with λ≫Δ will be frozen out. Identifying J˙/J with the Hubble constant, the frozen fluctuations of the scalar field are analagous to the quantum fluctuations of the inflaton field [132]. The presence of such a horizon would imply that after a quench to the Mott insulator, there will be frozen in number fluctuations.


### 3.2. Infinite Dimensional Limit

The infinite dimensional limit of the BHM is somewhat artificial, but it is a situation where an exact solution can be obtained, which may have relevance to the finite dimensional case. Sciolla and Biroli [87, 88] took this approach and studied quantum quenches in the infinite dimensional BHM. They took advantage of the fact that the ground state of the BHM is site permutation symmetric for any value of U/J to note that after a quench the state is also site permutation symmetric. This allows one to parametrize the states with the parameters xi, where xi is the fraction of sites with i bosons, i=0,1,2,…. Considering the particular example where the maximum number of bosons per site is 2 (a constraint that is easily relaxed), one has the conditions,
(34)x0+x1+x2=1;N=V(x1+2x2),
where N is the total number of bosons and V is the volume of the system. From these relations, one may specify the state in terms of x=x1. Sciolla and Biroli considered quantum quenches in which the interaction strength is taken from Ui to Uf at time t=0 and studied the subsequent time evolution of x and the superfluid order parameter |Ψ0|2. They found the interesting result that there is a dynamical transition at a value of Ufd=(Ui+Uc)/2, where Uc=3+22 (in units of J), at which there is exponential relaxation to the Mott state. This is illustrated in Figure 4 where it can be seen that in addition to the dynamical transition noted above, the value of 〈|Ψ0|2〉 tends to a finite value at long times, as opposed to the value of 0 expected if the system is in equilibrium. Sciolla and Biroli noted that this behaviour is an artifact of the infinite dimensional limit but suggest that the mean field approximation may capture aspects of the short-time dynamics in finite dimensional systems, so for large quenches one may expect to see freezing in superfluid correlations, and a failure to relax to the Mott state. A study by Snoek [89] using Gutzwiller mean field dynamics also found a dynamical transition of the sort identified by Sciolla and Biroli.Figure 4
Evolution of the microcanonical average of the superfluid order parameter〈|Ψ0|2〉 as a function of Uf for a quench with Ui=0. The mean field value of |Ψ0|2 is shown for comparison [87]. Adapted figure with permission from B. Sciolla and G. Biroli, Phys. Rev. Lett. 105, 220401 (2010). Copyright 2010 by the American Physical Society.


### 3.3. Gutzwiller Mean Field

The Gutzwiller mean field approach [133, 134] is a widely used way to look at the equilibrium phase diagram and out-of-equilibrium dynamics in the BHM in part because it is relatively straightforward to implement numerically (dynamics in systems with as many as 653 sites have been investigated [90]). The approach involves approximating the many-body state of the system as a product of single-site states in the form
(35)|Ψ〉=∏i‍∑n=0nmaxfn(i)|n〉(i),
with the constraint that ∑n‍|fn(i)|2=1 [12, 134]. The state |Ψ〉 is to be viewed as a variational ansatz, and so the fn(i) are to be found by minimizing
(36)〈Ψ|H^|Ψ〉.
For an average of one boson per site, the transition from superfluid to the Mott insulator is found for zJ/U=1/5.8, as in regular mean field theory.The Gutzwiller technique generalizes to time-dependent problems fairly straightforwardly [90, 135] by requiring that
(37)〈Ψ|iddt-H^|Ψ〉=0,
which leads to the equations [12]
(38)iddtfn(i)=U2n(n-1)fn(i)-JΦi*n+1fn+1(i)-JΦinfn-1(i),
where
(39)Φi=∑〈i,j〉〈aj〉=∑〈i,j〉∑nnfn-1(j)*fn(j).As with any mean field theory, the estimate of the critical coupling for the phase transition is not particularly accurate within the Gutzwiller approach. A more important failing is that due to the assumption of a variational state which has a product form, two site correlations factorize into single-site quantities and so the method does not capture correlations involving different sites. This is most problematic for intermediate values ofU/J where the ground state is not close to product form, unlike the two limits U=0 and J=0. Improvements to the Gutzwiller mean field have been made using a variety of approaches [91–93, 136] and can allow for perturbative corrections to short-range correlations [136]. Nevertheless, because of its simplicity, it is a very useful method for gaining an understanding of physics for regimes where exact numerical results are not easily obtained. For example, for out-of-equilibrium dynamics in dimensions higher than 1, for example, [94] (Section 4.2.3), or in the presence of complicated space and time-dependent potentials such as the “Gaussian spoon” consisting of a usual harmonic trap V0 with a Gaussian with velocity v:
(40)V(r,t)=V0(r)+V1e-((x-vt)2+y2)/W2,
suggested by Lundh [95] as a means to excite Mott insulating states.


#### 3.3.1. Projection Operator Approach

One suggestion to improve upon Gutzwiller mean field theory that allows the effects of quantum fluctuations to be included both in equilibrium and in dynamics is to use a projection operator approach [91, 92]. The particular projection operator introduced by Trefzger and Sengupta [91] was
(41)Pl=|n-〉〈n-|r⊗|n-〉〈n-|r′,
which lives on the link l between r and r′ and projects on to the manifold of states for which the occupation on sites r and r′ is n-. With this projection operator, the hopping term in the BHM may be written as (introducing the hopping operator Tl and replacing the sum over sites with a sum over links)
(42)-J∑〈rr′〉br†br′=∑lTl=∑l[PlTl+TlPl+Pl⊥TlPl⊥],
where Pl⊥=1-Pl. The terms PlTl+TlPl take the system out of the low energy manifold in Mott state where there are n- bosons on each site. The strategy employed in [91] to eliminate these terms from the Hamiltonian was to use a canonical transformation
(43)S=i∑l[Pl,Tl]U,
which leads to a transformed Hamiltonain
(44)H*=eiSHe-iS,
in which terms are kept up to order z2J2/U:
(45)H*=H0+∑lPl⊥TlPl⊥-1U∑l[PlTl2+Tl2Pl-PlTl2Pl-TlPlTl]-1U∑〈ll′〉[12PlTlTl′-TlPlTl′+12(TlPlPl′Tl′-PlTlTl′Pl′)+H.c.].
A Gutzwiller ansatz was then applied for the canonically transformed state |ψ′〉=eiS|ψ〉,
(46)|ψ′〉=∏r‍∑nfn(r)|n〉.
The phase diagram obtained using this state in [91] is significantly improved from mean field theory—in 2 dimensions the Mott insulator superfluid phase boundary is closer to QMC results than the mean field result is, and in 3 dimensions, the projection operator result is in very close agreement.As with Gutzwiller mean field theory, it is straightforward to generalize the method to allow time dependence. In particular, ifJ(t) is time dependent then the canonical transformation becomes time dependent and one has the following evolution equation:
(47)(iℏ∂t+∂S∂t)|ψ′〉=H*[J(t)]|ψ′〉,
provided J(t)/U≪1. This restricts the initial and final values of J to satisfy Ji/U,Jf/U≪1. On the other hand, it does not place restrictions on the ramp rate τ-1, which may be either fast or slow.Trefzger and Sengupta considered both quenches from the superfluid phase to the Mott insulator and from the Mott phase to the superfluid phase. In the quench from superfluid to the Mott insulator, they did not observe the expected universal scaling of either the residual energy(48)Q=〈ψf|H[Jf]|ψf〉-EG[ψf],
or the defect formation probability (which formed a plateau at large ramp times)
(49)P=1-|〈ψfΨG|ψf〉|2,
where |ψf〉 is the state at time tf after the ramp and |ΨG〉 is the ground state for the final value of J/U. EG[ψf] is the ground state energy determined by minimizing the energy functional for J=Jf. They observed that the timescale on which one might expect to see universality in P and Q can be estimated as the time for a boson to cross the system of size L, which may well exceed the system lifetime, so that the observed behaviour of P and Q is governed by local physics.On quenching from the Mott insulator to the superfluid phase, they found oscillations in the amplitude of the order parameter, with a period that scales likeT~(δJ)-α, where (δJ)=|Jf-Ji| and α=0.35±0.05 for quenches with a final value of J close to the critical J (higher values of J led to multiple frequency components being present). They also found P~(δJ)0.89 and Q~(δJ)1.90, both in disagreement with the exponents expected for a QCP with z=1 [13].The projection operator and canonical transformation used in this method have an advantage over Gutzwiller mean field theory in the Mott state in that they retain low energy excitations while removing high energy excitations via canonical transformation. In the superfluid phase close to the Mott insulator where number fluctuations are not too large, it is likely that this method may still be an improvement over regular Gutzwiller mean field theory, but deeper into the superfluid phase, where number fluctuations are much larger than in the Mott phase, it seems unlikely that this method will provide any improvement.


### 3.4.1/Z Corrections to Mean Field Theory

Mean field theory should be exact in the limit of infinite co-ordinationZ, but for finite dimensions, it is known to be inaccurate. Navez and Schützhold [93] developed a systematic expansion in 1/Z about mean field theory and used this expansion to explore a quench from a Mott insulator to superfluid. They started with reduced density matrices for 1 site (with site index μ):
(50)ρ^μ=Trμ{ρ^},
and 2 sites (site indices μ and ν)
(51)ρ^μν=Trμν{ρ^},
and similarly for higher numbers of sites. The correlated parts of these reduced density matrices, ρ^μνc, may be obtained from
(52)ρ^μν=ρ^μνc+ρ^μρ^ν.
Navez and Schützhold’s approach is based on the observation that the correlations in the Mott insulator have the scaling hierarchy
(53)ρ^𝒮c=O(Z1-|𝒮|),
where |𝒮| is the number of lattice sites in the set of lattice sets 𝒮. They showed that if this hierarchy is satisfied by the reduced density matrices ρ^𝒮c, then it is also satisfied by ∂tρ^𝒮c, so the hierarchy will remain valid for some finite time (of O(lnZ)) [93]. In the Mott phase with one particle per site ρ^μ=|1〉μ〈1|+O(1/Z). By making use of the hierarchy in (53) to neglect higher order reduced density matrices, they obtained an equation of motion for ρ^μνc, which leads to the following equations (with ℏ set to unity):
(54)(i∂t-U+3JTk)fk12=-2JTk(fk11+fk22+1),i∂tfk11=i∂tfk22=2JTk(fk12-fk21).
In (54), the fkij are the spatial Fourier transforms of the particle and hole correlation functions fμν11=Tr{ρ^h^μ†h^ν}, fμν12=Tr{ρ^h^μ†p^ν}, fμν21=Tr{ρ^p^μ†h^ν}, and fμν22=Tr{ρ^p^μ†p^ν}, where the local particle and hole operators are p^μ=|1〉μ〈2| and h^μ=|0〉μ〈1|, respectively, and Tk is (1/z)∑μ‍Tμνeik·(rμ-rν)≃1-k2/2m*, where Tμν is the connectivity matrix of the lattice.Navez and Schützhold considered a quench fromJ/U≪1 to a final much larger J, and hence assumed initial conditions fμν11=δμν, and all other fμνij=0. Writing the boson operators of the BHM as
(55)a^μ=h^μ+2p^μ+r^μ,
where  r^μ includes terms involving occupation numbers larger than 2, which decouple due to the scaling hierarchy assumption, they calculated 〈a^μ†(t)a^ν(t)〉 for an instantaneous quench from J≃0 to a final value of J=Jc(1+ϵ), where 0<ϵ≪1 and Jc is the critical value of J. In this limit, they estimated
(56)〈a^μ†(t)a^ν(t)〉≃𝒩(t)eγt2-(rμ-rν)2/c2,
where 𝒩(t) is weakly time dependent compared to the exponential, γ~J-Jc, and c2=3J(U-J)/m* with c a velocity scale. This form gives a constant propagation speed of the correlations, as found in the Lieb-Robinson bound [137]. The exponent in (56) also shows scaling behaviour in the vicinity of the critical point: as γ→γ′=λγ, the form is invariant under t→t′=t/λ and r→r′=r/λ. In addition to the growth of off-diagonal long range order, Navez and Schützhold also considered the growth of phase coherence by looking at the condensate fraction within a region 𝒮 with |𝒮|≪1. Defining the coarse-grained operator A^𝒮=∑μ∈𝒮a^μ/|𝒮|, which is the homogeneous mode which will give the largest eigenvalue of 〈a^μ†a^ν〉, this mode will have a macroscopic occupation N𝒮=〈A^𝒮†A^𝒮〉≫1 corresponding to the condensate fraction in the region 𝒮 after time t. For c2t2≪|𝒮|,  N𝒮/|𝒮|~1/|𝒮| and for |𝒮|≪c2t2, N𝒮/|𝒮|~eγt. This suggests that for distinct regions 𝒮 and 𝒮′ with 1≪|𝒮|,|𝒮′|≪c2t2, the relative phase between the regions correlates as
(57)〈ei(φ^𝒮-φ^𝒮′)〉≃e-γr2/2c2t,
where r is the distance between the regions and φ^𝒮 is defined via A^𝒮=N^𝒮eiφ^𝒮, and since 〈N^𝒮〉≫1, one can safely write N^𝒮≃N𝒮. This gives a growth of phase correlations which is diffusive, as opposed to ballistic growth of ODLRO found in (56).


### 3.5. Slave Particle Approach

Another approach in which the Hilbert space of the BHM is truncated to three states per site was put forward by Altman and Auerbach [96] who considered the BHM in the form
(58)H=U2∑i(n^i-n-)2-J∑〈ij〉(a^i†a^j+a^j†a^i)-μ∑i(n^i-n-),
and restricted the Hilbert space to states with n-,  n--1 or n-+1 bosons per site, allowing for a pseudospin 1 representation of the problem. This approach is most appropriate in and near the Mott insulating phase where number fluctuations are small. Altman and Auerbach introduced creation operators t^αi†, where |n-+α〉i=t^αi†|0〉i, with α=1,0,-1, subject to the constraint ∑αt^αi†t^αi=1. The original bosons of the BHM are represented as a^i†=n-t^0i†t^-1i+n-+1t^1i†t^0i. They also introduced modified coherent states of the form |Ω(θ,η,ϕ,χ)〉, with
(59)|Ω〉={cos(θ2)t^0†+eiηsin(θ2)×[eiϕsin(χ2)t^1†+e-iϕcos(χ2)t^-1†]cos(θ2)}|0〉,
and then chose a variational wave function of the form
(60)|ΦMF〉=∏i‍|Ω〉i.
The Mott phase is described by θ=0. When θ>0 there is a superfluid order parameter Ψ=n-sinθe-iϕ. Minimizing the variational energy leads to the usual mean field boundaries for the Mott lobes. To study excitations above the ground state, Altman and Auerbach performed the canonical transformation
(61)b^0i=cos(θ2)t^0+12sin(θ2)(t^1i+t^-1i),b^αi=sin(θ2)t^0-12cos(θ2)(t^1i+t^-1i),b^φi=12(t^1i-t^-1i),
subject to the constraint ∑mb^mi†b^mi=1. The mean field variational state may be written as |ΦMF〉=∏i‍b0i†|0〉, with bαi† and bφi† creating excitations about this state. They then used the constraint to eliminate b0 from the Hamiltonian and truncated the resulting Hamiltonian at quadratic order in the operators bα† and bφ†. After a Bogoliuibov transformation one obtains
(62)Hfluc=∑mkωmkβmk†βmk.
The modes associated with the two types of fluctuations in the superfluid phase are a massive amplitude (Higgs) mode and a phase (sound) mode [96]
(63)ωα(k)≃c2k2+Δ2,ωφ(k)≃c|k|,
where c and Δ may be expressed in terms of J,  U,  z, and n-. In the Mott phase, the two modes are degenerate and gapped and represent particle and hole excitations.Altman and Auerbach also constructed a path integral using the modified coherent states, and writing the action in terms ofΨ=n-sinθe-iφ, they obtained the action
(64)S=18zJn-2∫0t‍dt′∫‍ddr{|Ψ˙|2-(2Jn-)2z|∇Ψ|2-(2Jn-z)2(u-1)|Ψ|2-(zJ)2n-u|Ψ|4|Ψ˙|2}.
For large, n- the critical point is at u=U/4zJn-=1. From this action they obtained an equation of motion for the amplitude of the superfluid order parameter (after rescaling)
(65)Ψ¨=c2∇2Ψ+12Δ2Ψ(1-|Ψ|2).
The solution of this equation when Ψ is assumed to be spatially uniform shows oscillations as a function of time for a system prepared in the Mott state, as illustrated in Figure 5 (although these oscillations were not observed in a recent quench experiment [138])(a) Mean field energy in the superfluid state as a function ofΨ with the equiilbrium ground state indicated at (2); (b) solution of (65) as a function of time, compared to the equilibrium value of Ψ (dashed line) [96]. Adapted figure with permission from A. Altman and A. Auerbach, Phys. Rev. Lett. 89, 250404 (2002). Copyright 2002 by the American Physical Society.
(a)
(b)Altman and Auerbach noted some caveats to their results: the assumption of a uniform system may be spoiled by topological defects that get trapped by the Kibble-Zurek mechanism [63, 64]. For initial states near the transition, the correlation length should be large enough that such defects are widely spaced. However, for a quench from deep in the Mott phase, there will be a much higher initial density of defects. Altman and Auerbach argued that the k=0 mode will be the fastest growing, so that as J is increased to near the transition, defects will generally be widely separated by ξ=c/Δ. The question of the effects of the remaining vortices on the evolution of Ψ is an open question. The oscillations shown in Figure 5 are undamped—in general, one would expect that such oscillations would damp as a function of time. Altman and Auerbach argued that such damping could come from the coupling of the amplitude mode to low-energy phasons, which would lead to phason pair emission.The theory discussed in [96] is primarily for n-≫1, motivated by the experiments of Orzel et al. [139], for which n-~50. Further work by Huber et al. [118] extended the slave particle approach to allow for n- of order unity. Using a similar procedure in which they constructed the low energy excitations above the mean field variational ground state, they found qualitatively similar features to Altman and Auerbach, namely, a massive Higgs mode and a gapless sound mode in the superfluid and gapped particle hole excitations in the Mott insulator. The theory has similar caveats to the large n- theory that J/U not be too large and the density of excitations not be too large (which will certainly be satisfied in the Mott phase).


### 3.6. Numerical Methods

Numerical simulations give the opportunity to obtain quantitative results for both equilibrium and out-of-equilibrium dynamics both for the homogeneous BHM and for cold bosons in an optical lattice. Methods that have been used extensively for out-of-equilibrium dynamics include exact diagonalization and the two closely related methods: time-dependent density matrix renormalization group (t-DMRG), and time-evolving block decimation (TEBD). The t-DMRG [140, 141] and TEBD [142] methods are primarily limited to use in one dimension, due to the exponential increase in computational resources required in two dimensions [143], whereas exact diagonalization may be used in any dimension, but the Hilbert space grows exponentially with the number of bosons, limiting the number of sites that may be simulated. In both t-DMRG and TEBD, the Hilbert space is truncated, but low lying energy states are retained, allowing for the simulation of larger systems in 1 dimension than with exact diagonalization. Reviews by Schollwöck [143, 144] give a comprehensive discussion of details of these methods—my focus here is on the physics that has been obtained from numerical calculations for quenches from (i) the superfluid phase to the Mott insulating phase and (ii) the Mott insulator phase to the superfluid phase.


#### 3.6.1. Superfluid to Mott Insulator Quench

Kollath et al. [97] used a combination of exact diagonalization and t-DMRG to study the behaviour of the homogeneous BHM with unit filling after a sudden quantum quench from the superfluid phase with an initial value of the interaction strength Ui to a final value Uf≫Ui. Their simulations encompassed up to 18 sites in two dimensions and up to 64 sites in 1 dimension. They calculated the time-and space-dependent correlation functions 〈a^x†(t)a^0(t)〉 (note that Kollath et al. denote the boson creation operators of the BHM by b† rather than a† as used here) to characterize the dynamics after such a quench. This correlation function is seen not to decay to zero, but after an initial relaxation on a timescale of order 1/J is seen to have small amplitude oscillations, with period 2π/Uf about a finite value after a quench, as illustrated in Figure 6.Time dependence of the nearest neighbour correlations forUi=2J and Uf=40J obtained using (a) exact diagonalization and t-DMRG methods in 1 dimension and (b) exact diagonalization in 2 dimensions [97]. Adapted figure with permission from C. Kollath, A. M. Lauchli, and E. Altman, Phys. Rev. Lett. 98, 180601 (2007). Copyright 2007 by the American Physical Society.
(a)
1D chain,Ui=2J,  Uf=40J
(b)
2D square,Ui=2J, Uf=40JKollath et al. investigated the decay of the time-dependent correlation functions with distance after the decay of the initial transient and found that there were two different behaviours, depending on the final value ofUf. For Uf/J not too large, the correlations fit reasonably well to thermal correlations determined with QMC and suggest that the final state can be viewed as thermalized. However, for large Uf, there was a strong memory of the initial state, and the decay of correlations with distance was much slower than found either for a thermal state or in the ground state at U=Uf, as illustrated in Figure 7, which they characterized as a nonthermal steady state.Figure 7
Average spatial correlations as a function of site separation after a quench fromUi=2J to Uf=40J in (a) 1 dimension and (b) 2 dimensions. Squares indicate the results from exact diagonalization and triangles indicate results from t-DMRG. The diamonds show the result of equilibrium QMC calculations, and the circles indicate the zero temperature correlations in the ground state for Ui and Uf [97]. Adapted figure with permission from C. Kollath, A. M. Läuchli, and E. Altman, Phys. Rev. Lett. 98, 180601 (2007). Copyright 2007 by the American Physical Society.Kollath et al. [97] rationalized the slow relaxation they found for large Uf by focusing on the particle and hole excitations of the Mott insulator. If quasiparticle interactions are ignored, then the effective Hamiltonian in the Mott regime may be written in the form H0=∑k,αωkβαk†βαk, with βαk† a creation operator for a quasiparticle in the Mott insulator. Thermalization then proceeds through the relaxation of the quasiparticle distribution given by the initial conditions. This requires processes in which quasiparticle number is not conserved, which are strongly suppressed if the gap Δ is larger than W/2 where W~4zJ is the quasiparticle bandwidth. Higher order processes could still lead to thermalization, but on much longer timescales than accessible numerically. This picture relies on having a dilute quasiparticle population, that is, an initial state close to the transition, but at least in their numerics, Ui=2J, appear to have been sufficiently close to the transition (Uc=3.37J in 1 dimension [44] and Uc=16.7J in 2 dimensions [41]) for this physics to be observable.The single particle correlations studied in [97] and also density-density correlations 〈n^jn^j+r〉 were calculated as a function of time and position with similar Ui and Uf=40J for average fillings 〈n〉=1 and 0.5 by two of the authors of [97] in [98]. They observed that the single particle and density-density correlations showed a “light-cone” effect in that there was a propagating front in the correlation with velocity vs. This behaviour is similar to that demonstrated by Lieb and Robinson [137] for one dimensional spin systems with finite range interactions, that at long distances there is a maximal velocity at which correlations can spread. The values of vs for 〈n〉=1, for which Uf corresponds to the Mott phase, were on the order of 5-6J, whereas the values of vs for 〈n〉=0.53 (for which the system remains superfluid) after the change in U were of the order of 4J. Additionally, Lauchli and Kollath [98] calculated the change in von Neumann entropy for blocks of length l as a function of time: S(l,t)-S(l,0), and found that this showed linear growth for short times and saturated at a time t*(l)≃l/2ve for periodic boundary conditions and t*(l)≃l/ve for open boundary conditions with ve≃vs in accord with predictions of Calabrese and Cardy [145].The dynamics of local observables both in homogeneous systems and in the presence of a parabolic trapping potential were studied by Bernier et al. [99] using a combination of exact diagonalization calculations and t-DMRG to investigate quenches of the form U(t)=Ui+(t/τ)(Uf-Ui) for the 1d BHM. By focusing on local observables, such as the local compressibility κj=〈nj2〉-〈nj〉2, and the probability of occupation of a site by n bosons, Pn, they observed that the dynamics split into two regimes: short ramp times and long ramp times as illustrated in Figure 8. For short ramp times (τ≲ℏ/J), the behaviour near the centre of a homogeneous system and a trap is essentially identical, and indicative of local dynamics, since the ramp timescale is less than the hopping timescale. At longer times, the inhomogeneity of the trap becomes apparent as density redistribution takes place, as seen in different values for the local quantities in the presence and absence of the trap. This behaviour was seen both for a quench from superfluid to the Mott insulator (Ui=2J and Uf=4J) and for a quench within the Mott insulator (Ui=4J and Uf=6J).Figure 8
Dependence of local observables on ramp timeτ for both homogeneous systems and at the centre of a parabolic trap, n: average density per site; κ: compressibility; and Pn: the probability that a site has an occupation of n particles. (a) For a protocol with Ui=2J and Uf=4J, (b) for the same protocol with Ui=4J and Uf=6J. Adapted figure with permission from J. S. Bernier, G. Roux, and C. Kollath, Phys. Rev. Lett. Lett. 106, 200601 (2011). Copyright 2011 by the American Physical Society.In subsequent work investigating local observables and local currents in a trap, Bernier et al. [100] observed that regions with density of order 1 boson per site (even if not in a Mott insulator state) can act as “Mott barriers” that impede particle transport, slowing down dynamics associated with density redistribution in the trap. These ideas are discussed further in the context of quench experiments by Hung et al. [146] in Sections 4.2.2 and 4.2.3.


#### 3.6.2. Mott Insulator to Superfluid Quench

Clark and Jaksch [101] used TEBD for the one dimensional BHM to study the effect of varying the lattice depth V0(t) with time on the coherence in ramping from the Mott insulator to the superfluid phase. They considered slow ramps in which the lattice depth was evolved from the superfluid state with U/2J=2 to the Mott insulator with U/2J=40. Comparing the exact numerical solution with TEBD for a system of M=7 bosons, they established that their slow ramp had an infidelity of 1-F<10-4 so that their slow ramp could be viewed as adiabatic. They then considered larger systems with M=49 prepared with a slow ramp and studied how coherence was reestablished as the lattice depth was reduced. To do this, they studied a time-dependent correlation length ξc(t). They found that for a ramp in which V0 was reduced slowly in a similar manner to how it was increased, that the restoration of coherence was on a similar timescale either for a system that had been prepared with a slow ramp or for a system that had been prepared in the ground state for U/2J=20, with the timescale for single atom hopping governing the establishment of coherence. In contrast, for a system prepared with a slow ramp and then subjected to a fast ramp, Clark and Jaksch found that ξc increased quickly for short ramp times and reached the value expected in the superfluid on much shorter timescales than that of single-particle hopping, leading them to suggest that higher-order correlations must play a role in the establishment of coherence during a fast quench.


#### 3.6.3. Entanglement Growth

Numerical simulation of the time dependence of quantum systems using approaches based on matrix product states, such as t-DMRG and TEBD, is limited in evolution times due to the growth of entanglement [145, 147]. Motivated by an experiment in which there was a quantum quench for which relaxation dynamics persisted to timescales not accessible with DMRG [148], Daley et al. [149] focused on the growth of entanglement during a quench from a Mott insulator to a superfluid and proposed a protocol to measure this entanglement growth by obtaining the arbitrary order Renyi entropies
(66)Sn(ρ)=11-nlogTr{ρn}.
They tested the scheme with t-DMRG to show its robustness against measurement imperfections. Implementation of this scheme experimentally could lead to quantitative insights into the performance of quantum simulators [149].


### 3.7. Fermionization

In the limit thatU/J→∞, there can be no double occupancy of bosons on a site, which leads to the notion of hard-core bosons, which obey bosonic commutation relations if they are not on the same site, that is, [a^i,a^j†]=[a^i,a^j]=[a^i†,a^j†]=0 if i≠j, but have fermionic commutation relations {a^i,a^i†}=1, a^i2=[a^i†]2=0 on the same site. In one dimension, the hard-core boson model can be mapped onto fermions by the use of the Jordan-Wigner transformation [150]
(67)a^i†=f^i†∏j=1i-1‍e-iπn^jf,a^i=∏j=1i-1‍eiπn^jff^i,
where the f^i† and f^i are fermionic creation and annihilation operators on site i, respectively, and n^if=f^i†f^i is the fermionic number operator on site i. This mapping removes the hard-core constraint and allows for the efficient simulation of out-of-equilibrium dynamics [150, 151] of hard-core bosons.An alternative fermionization scheme which allows for a treatment of out-of-equilibrium dynamics of the BHM whenU is finite was recently introduced by Barmettler et al. [152] for the 1-dimensional BHM in the strongly interacting regime, where there are small fluctuations around integer filling. They truncated the Hilbert space to three states on each site, similarly to the slave particle approach, but introduce different auxiliary bosons b^j±†, such that
(68)a^j†=n-+1b^j+†+n-b^j-,
and for the three states on site j: |n-+m〉j with m=±1,0, the operators satisfy
(69)b^j+†|n-〉j=|n-+1〉,b^j-†|n-〉j=|n--1〉j,bj±|n-〉j=0,
with additional hard core constraints to ensure no double occupancy: (b^j±†)2=b^j±2=0. In addition, there is not allowed to be double occupancy of different species of bosons, which implies the constraint
(70)b^j,+†b^j,+b^j-†b^j-=n^j+n^j-=0.
Barmettler et al. used a Jordan-Wigner transformation similar to (67), to introduce two species of fermions (reflecting two species of bosons):
(71)b^j,σ=Zj,σc^j,σ,
where
(72)Zj,+=eiπ∑σ,j′<jn^j′,σ,Zj,-=Zj,+eiπn^j,+,
are nonlocal string operators which satisfy Zj,σ†=Zj,σ and Zj,σ2=1. This form removes the hard-core constraints on b^j+ and b^j- individually, but not on the product n^j+n^j-. This condition may be satisfied with the projector 𝒫=∏j‍Pj, where Pj=1-n^j+n^j-, in which case the BHM Hamiltonian may be represented by
(73)H=∑j‍𝒫{U2-J(n-+1)c^j+†c^j+1,+-Jn-c^j+1,-†c^j--Jn-(n-+1)(c^j+†c^j+1,-†-c^j-c^j+1,+)+H.c.+U2(n^j++n^j-)}𝒫.
In order to circumvent the complications of calculating with the projector in the Hamiltonian, Barmettler et al. studied the problem with the approximation of unconstrained fermions (UF), that is, 𝒫=1. With this assumption, the Hamiltonian may be diagonalized with a Bogoliubov transformation to the form
(74)HUF=∑k,σϵσ(k)γ^k,σ†γ^k,σ,
where ϵσ(k)=-σJcos(k)+ℏω(k) and
(75)ℏω(k)=12[E+(k)+E-(k)]2+4|Δ(k)|2,
with
(76)E+(k)=-2J(n-+1)cos(k)+U2,E-(k)=-2Jn-cos(k)+U2,Δ(k)=2iJn-(n-+1)sin(k).
Barmettler et al. found that for U/J<4(n-+1), the UF approximation breaks down.Quench in the Mott Phase. In order to study the out-of-equilibrium dynamics, Barmettler et al. chose an initial Fock state with filling n- and considered an instantaneous quench from this J/U=0 state to a final J/U which remains within the Mott phase. They calculated the time-dependent density-density correlation function
(77)Cd(t)=〈n^j(t)n^j+d(t)〉-〈n^j(t)〉〈nj+d(t)〉,
for sites separated by d and found that for d>1,
(78)Cd(t)≃-(n-(n-+1)JdU)2(Jd(J~t)J~t)2,
where J~=2J(2n-+1)/ℏ and Jd is a Bessel function.Focusing in particular on the case ofn-=1, they studied density correlations Cd(t) with both DMRG and the UF approximation for U/J=18 and U/J=9. These two methods both showed close agreement and a propagating front of correlations, corresponding to doublon-holon pairs, with the same constant velocity in both methods. By approximating the expression (78) in the d≫1 limit as
(79)Cd(t)~-(2d2/321/3ℏ3Ut)2Ai2[-(2d)1/3(6Jtℏ-d)],
Barmettler et al. made use of the fact that the Airy function Ai(-21/3z) has a peak for z0≃1.02 to infer that there is a propagation front with correlations exponentially suppressed for t<tpeak where
(80)Jtpeakℏ≃16[d+z0(d2)1/3],
from which it can be inferred that v∞=6J/ℏ is the velocity at which correlations spread as d→∞ for large U/J. Barmettler et al. provided a similar analysis for noninteracting bosons to deduce v∞=4J/ℏ when U/J=0. In the regime of validity of the UF approximation, they considered the relative velocity of quasiparticle pairs analytically, which takes a maximum value
(81)vmax=maxk|v(k)|=maxk|2ddkω(k)|≃2J(2n-+1)ℏ(1-8n-(n-+1)J2(2n-+1)2U2)+O(J4U3),
which corresponds to vmax=6J/ℏ for n-=1 and decreases monotonically with increasing J/U from the J/U=0 limit.The fermionization procedure introduced by Barmettler et al. gives very nice results that go beyond previous analytical work on the one-dimensional BHM and, as discussed in Section4.4, also give good agreement with experiment. The derivation of a velocity which governs the spreading of correlations at large distances is an important result, as it gives a concrete value for a Lieb-Robinson-like bound on the speed at which correlations spread. The Lieb-Robinson bound was originally formulated for interacting spins on a lattice [137, 153–155], but these results indicate that it also appears to hold in the Mott insulating phase of the one-dimensional BHM.


### 3.8. Closed Time Path Approaches

The closed time path (CTP) or Schwinger-Keldysh technique [156–161] is an approach that allows for the calculation of both equilibrium and out-of-equilibrium quantum dynamics within the same formalism. In CTP methods the problem is formulated in real time (rather than imaginary time, as for the Matsubara formalism), allowing out-of-equilibrium problems to be tackled without analytic continuation. The price to pay is that the number of fields in the theory doubles, with a second copy of each field propagating backwards in time. Consequently in calculations of Green’s functions, the notion of time ordering needs to be replaced by that of contour ordering [159] on a contour such as illustrated in Figure 9. There are several reviews that provide a detailed introduction to CTP techniques [158–161], so few details will be given here.Figure 9
Contour for the Schwinger-Keldysh technique.Several authors have used these approaches to study the BHM both from the weakly interacting limit (largeJ/U) [162–166] and the strongly interacting limit (small J/U) [112, 120, 167, 168]. The Schwinger-Keldysh approach has also been used to study out-of-equilibrium dynamics for the BHM with periodic driving [20], including allowing for the possibility of ohmic dissipation [20, 169]. I first discuss the two-particle irreducible formalism which has a weak interaction starting point and then discuss the strongly interacting approach.


#### 3.8.1. Two Particle Irreducible Formalism

The two particle irreducible (2PI) closed-time-path (CTP) formalism for the BHM has been developed in detail by Rey et al. [162] and Temme and Gasenzer [165] who studied the 1/𝒩 expansion.The 2PI formalism can be obtained by starting from the action for the complex fieldsΦ and Φ* in the Bose Hubbard model in the form
(82)S[Φ]=∫‍dt∑iΦia(t)τabyℏ∂tΦib(t)+∫‍dt{U4𝒩(JΦi+1a(t)τabxΦ1b(t))-U4𝒩[τabxΦia(t)Φib(t)]2},
where a compact notation has been introduced, with Φi1=Φi and Φi2=Φi*, and the matrices τx and τy have the same form as the corresponding Pauli matrices but are labelled so as not to be confused with spin. The quantities studied are ϕia(t)=〈Φia(t)〉 and Gijab(t,t′)=〈TCΦia(t)Φib(t′)〉-ϕia(t)ϕib(t′), where TC  indicates contour ordering. The theory is specified by the effective action Γ[ϕ,G], which takes the form
(83)Γ[ϕ,G]=S[ϕ]+i2Tr[lnG-1]+i2Tr[D-1(ϕ)G]+Γ2[ϕ,G],
where
(84)iDijab(t,t′)-1=δS[ϕ]δϕia(t)δϕjb(t′),
and Γ2 is the two-particle irreducible generating functional. Γ2 can be described as an infinite sum of diagrams from the expansion of interaction terms in S[ϕ+φ] of cubic and higher order in φia=Φia-ϕia [162].In practice, various approximation schemes keep some subset of these diagrams and can be classified based on which terms are kept. The simplest approximation is to assume thatΓ[ϕ,G]=S[ϕ], corresponding to the mean field approximation. The Bogoliubov (one-loop) approximation can be obtained if the effective action is taken to be (83) with Γ2 ignored. Approximations that take into account interactions more carefully include the time-dependent Hartree-Fock-Bogoliubov (HFB) approximation, in which Γ2 is truncated to keep only the first-order diagram in U. This approximation conserves energy and particle number but violates Goldstone’s theorem. The expansion can be viewed as one in Ut/J, where t is time, and hence will be most accurate for short times or small U/J. It is possible to keep higher-order terms in Γ2—second-order expansions were considered by Rey et al. [162] and Temme and Gasenzer [165], who both considered a 1/𝒩 expansion, which simplifies the second-order expansion of Γ2. From each of these various effective actions, saddle point equations of motion can be derived, which may be quite lengthy; the full expressions for the second order expansion are listed in Appendix B of [162]. These equations of motion may also be used to derive a quantum kinetic equation, as was done by Rey et al. in [163].An additional issue discussed in [162] is that of conservation laws: for a closed system, the evolution is unitary in time, so energy must be conserved. However, particular approximation schemes for Γ may not conserve energy. On the other hand, schemes which are ϕ derivable, that is, the self-energy may be written as the derivative of a functional of G [170] do conserve particle number, energy, and momentum, which is the satisfied for the cases considered in [162]. A comparison of the results obtained for the approximation schemes outlined in [162] are shown in Figure 10 and compared against the results of exact dynamical evolution in small systems. The condensate population as a function of time is shown for initial conditions in which the condensate starts all in a single well, for N=20, 40, and 80 particles with J=1/2 and NU/J=4. The Hartree-Fock Bogoliubov approximation is seen to be fairly poor at describing the time development of the condensate, with the 1/N approach giving the most accurate account of the time evolution for the condensate population per well as a function of time as N becomes larger.Figure 10
Condensate population per well in each well of a two-well system and total condensate population as a function of time (in units ofℏ/J) for parameters J=1/2, NU/J=4, and N=20,40,80 as calculated exactly, using second-order expansion, 1/N expansion, and the Hatree-Fock Bogoliubov approximation [162]. Adapted figure with permission from A. M. Rey, B. L. Hu, E. Calzetta, A. Roura, and C. W. Clark, Phys. Rev. A 69, 033610 (2004). Copyright 2004 by the American Physical Society.One can see that the range of couplings considered are very much in the weakly interacting limit—this suggests the need for other methods to tackle questions relating to the crossing of the superfluid Mott insulator phase transition.


#### 3.8.2. Strong Coupling

Sengupta and Dupuis [119] developed a strong coupling nonperturbative approach to the BHM with features that make it particularly attractive for generalization to real time. In particular, their approach is exact in the limit of both small (U=0) and large (J=0) interactions and obtains the expected features in the excitation spectrum in both the superfluid and the Mott insulating phases. It also naturally allows for the calculation of nonlocal correlations. I first briefly recap the equilibrium calculation of Sengupta and Dupuis and then discuss generalizations of this approach using the Schwinger-Keldysh method [112, 120].Sengupta and Dupuis started with a path integral form for the BHM and then performed two Hubbard Stratonovich transformations to write the generating functional as(85)Z=Z0∫‍𝒟[ψ*,ψ,ϕ*,ϕ]e(ψ∣Jψ)-[(ψ∣ϕ)+c.c.]+W[ϕ*,ϕ],
where Z0 is the generating functional when J=0, with the notation (ϕ∣ψ)=∫0β‍dτ∑r‍ϕ*(r)ψ(r), where τ is imaginary time and β is inverse temperature. The term
(86)W[ϕ*,ϕ]=∑R=1∞(-1)R(R!)2∑a1⋯aRG{ai,bi}Rcϕa1*⋯ϕaR*ϕbR⋯ϕb1,
that enters in the action is evaluated with all the fields at the same site, and
(87)G{ai,bi}Rc=(-1)R〈ψa1⋯ψaRψbR*⋯ψb1*〉,
is the connected local 2R point Green’s function. By integrating out the ϕ field in the generating functional, Sengupta and Dupuis obtained an effective action that to quartic order in the field ψ takes the form
(88)S[ψ*,ψ]=-∑a,bψa*(Gab-1+Jab)ψb+14∑a1,a2,b1,b2Γa1,a2,b1,b2IIψa1*ψa2*ψb2ψb1,
where G-1 is the Green’s function in the local limit and ΓII is the two particle vertex in the local limit. Sengupta and Dupuis approximated ΓII by its static value and introduced g=(1/2)ΓII|static. This leads to the effective action
(89)S=-∫0β‍dτdτ′∑r,r′‍ψr*(τ)[G-1(r,τ,r′,τ′)+Jr,r′δ(τ-τ′)]ψr′(τ′)+g2∫0β‍dτ∑rψr*ψr*ψrψr,
which has the attractive feature that it leads to the correct Green’s function -〈ψr(τ)ψr′*(τ′)〉 in both the strong interaction (J=0) and noninteracting (U=0) limits. This feature makes this approach particularly appealing for the study of out-of-equilibrium dynamics, since it gives the hope that one can accurately represent the behaviour of the system over a wide range of parameters that cover both the superfluid and Mott insulating phases.By taking a saddle point approximation to the action, and then expanding to quadratic order in fluctuations, Sengupta and Dupuis obtained the following results for the excitations in the Mott and superfluid phases, respectively. In the Mott phase, they found that the Green’s function𝒢(k,iω)=-〈ψ(k,iω)ψ*(k,iω)〉 has the form
(90)𝒢(k,iω)=1-zkiω-Ek-+zkiω-Ek+,
where
(91)Ek±=-δμ+ϵk2±12[ϵk2+4ϵkUx+U2]1/2,
with ϵk the Fourier transform of the hopping -Jr,r′,  x=n0+1/2 and
(92)zk=Ek++δμ+UxEk+-Ek-,
with δμ=μ-U(n0-1/2) where n0 is the number of bosons per site for chemical potential μ when J=0. In the Mott insulating phase, the modes Ek± are gapped and the system becomes unstable to superfluidity when the energy of one or other of the modes vanishes for k=0. In the superfluid phase, the expressions for the spectrum are more complicated, and four modes are found ±Ek±, two of which are gapless as k→0 and have the form Ek-=c|k|, and the other two are gapped, as illustrated in Figure 11.Figure 11
Excitation energies±Ek± for n=1 and U/J=20 [119]. Adapted figure with permission from K. Sengupta and N. Dupuis, Phys. Rev. A 71, 033629 (2005). Copyright 2005 by the American Physical Society.Both  Grass et al.  [120] and Kennett and Dalidovich [112] used the Schwinger-Keldysh technique to construct theories for the real-time dynamics of the BHM.  Grass et al.  developed a real-time Ginzburg-Landau theory based on the field theoretic approach developed in [171, 172] and used this theory to calculate collective modes of the BHM, finding amplitude and sound modes in the superfluid and particle and hole excitations in the Mott insulator, as expected.  Grass et al.  also illustrated the power of this approach by expanding the Green’s functions in the equations of motion in the superfluid regime, they were able to obtain the equation of motion
(93)i∂Ψi∂t=-∑jJijΨj-μΨi+U2Ψi|Ψi|2,
which is the lattice version of the Gross-Pitaevski equation (GPE) [102].A similar approach developed by Kennett and Dalidovich [112] generalizes the theory of Sengupta and Dupuis to real time. We determined the effective action to quartic order for the BHM within the strong-coupling Schwinger-Keldysh approach and then used this to derive equations of motion to study the out of equilibrium dynamics.The approach followed in Kennett and Dalidovich [112] was to write a path integral for the BHM
(94)𝒵=∫‍[𝒟a*][𝒟a]eiSBHM[a*,a],
where aia is the field at site i on contour a, with a=1 or 2, and with the notation τi representing the ith Pauli matrix that acts in Keldysh rather than spin space. We performed a Keldysh rotation from the 1, 2 basis to the quantum (q), classical (c) basis [168, 173–175]:
(95)(a1(t)a2(t))→(a~q(t)a~c(t))=L^(a1(t)a2(t)),
where  aq and ac are the quantum and classical components of the field and
(96)L^=12(1-111).
The effect of this on the action is that τ3 in the 1, 2 basis becomes τ1 in the q, c basis. Similarly to Sengupta and Dupuis, we performed two Hubbard-Stratonovich transformations, the first of which decouples the hopping term and the second leads to the effective action to quartic order in terms of z fields (one can show that the Green’s functions for z are the same as those for the original field a [112, 119]) as
(97)Seff[z*,z]=∫‍dt∑ijzia*(t)(2Jij)τab1zjb(t)+∫‍dt1dt2∑izia1*(t1)[Gia2a1(t2,t1)]-1zia2(t2)+14∫‍dt1dt2dt3dt4∑iua1a2a3a4(t1,t2,t3,t4)×zia1*(t1)zia2*(t2)zia3(t3)zia4(t4),
where
(98)ua1a2a3a4(t1,t2,t3,t4)=14∫‍dt5dt6dt5′dt6′Gia5a6a6′a5′2c(t5,t6,t6′,t5′)×{[Gia5a1(t5,t1)]-1[Gia6a2(t6,t2)]-1×[Gia3a5′(t3,t5′)]-1[Gia4a6′(t4,t6′)]-1+((a4,t4)⟷(a3,t3))+((a6,t6)⟷(a5,t5))+((a6,t6)⟷(a5,t5);(a4,t4)⟷(a3,t3))[Gia4a6′(t4,t6′)]-1}.
The matrix Green’s functions G contains the retarded, advanced, and Keldysh Green’s functions. These Green’s functions, along with the two-particle connected Green’s function G2c are evaluated with respect to the single-site Hamiltonian H^0: full expressions are presented in [112]. The mean field phase boundary can be determined from the effective action equation (97) from the vanishing of the coefficient of zq*zc. The diagram corresponding to the interaction term u is illustrated schematically in Figure 12.Figure 12
Schematic diagram for the coefficient of the quartic term in the effective actionua1a2a3a4(t1,t2,t3,t4).The following symmetry relations hold for the interaction kernelu from the definition above:
(99)uabcd(t1,t2,t3,t4)=ubacd(t2,t1,t3,t4)=uabdc(t1,t2,t4,t3),
(this corrects [112]) and it may also be shown [112] that
(100)Gia5a6a6′a5′2c(t5,t6,t6′,t5′)=Gia6a5a6′a5′2c(t6,t5,t6′,t5′)=Gia5a6a5′a6′2c(t5,t6,t5′,t6′).
Grass et al. [120, 168] also noted similar symmetry relations for this correlation function under interchange of Keldysh indices. These symmetries imply that G2c has only 8 independent components: Gqqqq2c,Gcqqq2c,Gqqqc2c,Gqqcc2c,Gccqq2c,Gcqcq2c,Gqccc2c, and Gcccq2c. The remaining four-point function Gcccc2c=0 by causality [161]. Explicit expressions for each of the nontrivial components were written down in [112].A simplified low frequency equation of motion for the superfluid order parameterz can be derived from the action that applies away from the degeneracy points of the Mott lobes (in which case only the component Gcqqq2c is required). Equations of motion that allow for high frequencies and spatial variations of the order parameter and its correlations can also be derived. Kennett and Dalidovich focused on k=0 and defined δ(t)=-2dj(t), where the hopping is taken to be J(t)=J0+j(t) with J0 chosen so that
(101)2dJ0+ν=0,
(i.e., J0 lies on the mean field phase boundary for a given μ). Hence, the approximate equation of motion is
(102)κ2∂2z∂t2+iλ∂z∂t+δ(t)z+u|z|2z=0,
where the coefficients κ2, λ, and u can all be calculated in terms of U and μ [112]. We studied (102) for fixed μ and time-varying J.For fixedμ there are two possibilities for the dynamics: (a) the particle-hole symmetric case (the transition at the tip of the Mott lobe), in which case λ=0, and (b) the generic case, in which λ≠0. The function δ(t) controls the traversal of the quantum critical point. We chose
(103)limt→-∞δ(t)=-δ0;limt→∞δ(t)=δ1.
In our numerics, we used
(104)δ(t)=(δ0+δ12)tanh(tτQ)+δ1-δ02,
where we introduced a timescale τQ which is the characteristic time for δ(t) to cross from -δ0 to δ1 [111]. In the vicinity of the transition t=0 is linear in t. Functional forms of δ(t) which are not linear in t in the vicinity of t=0 may lead to differing behaviour in the long-time limit [79].We solved (102) numerically with λ=0 for a variety of values of τQ, and found that for t≫τQ the form of the solution is that z(t) oscillates in a periodic manner with a magnitude that decreases with increasing τQ. When averaged over a period T at times t≫τQ,
(105)〈z〉T=1T∫tt+T‍dt~z(t~)=0,
as would be expected in the Mott insulating state. Defining zmax(τQ)=limt→∞|z(t)| we found that zmax(τQ) decreases with increasing τQ without any indication of saturation.In the generic case in whichλ≠0, we solved (102) numerically for a variety of values of τQ as is displayed with time scaled by τQ in Figure 13 for μ/U=0.25 (well away from both degeneracy and the particle hole symmetric case). Both the real and imaginary parts of z(t) are oscillatory functions, each of which has a vanishing average over a period for times greater than zero. The nonzero amplitude of z at long times indicates that there is memory of the nonzero value of ρ in the superfluid phase and that the final state is a metastable state rather than the Mott insulator in which ρ=0. The collapse of traces with different τQ after scaling the time with τQ indicates that |z| has a finite value even when τQ becomes large.Figure 13
Dynamics of|z(t)| normalized to unity in the generic case, for a variety of τQ with time rescaled by τQ for a quench from 2dJ/U=2.0 to 2dJ/U=0.0. Figure adapted from [112].To determine the fate of the system at long quench times, we definedzmax=limτQ→∞zmax(τQ) and calculated this as a function of chemical potential. Apart from the particle-hole symmetric point, it appears that there is a transition to a metastable state in which zmax≠0. Whilst it is true that the spatial dependence of z and higher frequency components were ignored in (102), it seems reasonable that the metastable state identified at the mean field level may persist even when these are taken into account.The calculations discussed above are for fixed chemical potential rather than fixed particle number, as is found in a cold atom system. However, within the LDA fixed chemical potential would correspond to studying dynamics at fixed radius. The results found for the long-time limit ofzmax are in accordance with this idea. There appears to be equilibration as zmax→0 in the particle-hole symmetric case, where there is no change in the number of bosons per site, in which only local equilibration takes place. For other values of the average number of bosons per site changes in crossing from the superfluid to the Mott insulator in which case global mass transport is required for equilibration (which is not captured within the equation of motion studied in [112]).


### 3.9. Exact Results

#### 3.9.1. Spectral Moments

Perturbative expansions inJ/U [41–43] have allowed for very accurate determination of the Mott insulator phase boundary and the momentum distribution in the Mott insulator phase [176]. Very recently Freericks et al. [177] calculated exact sum rules for the spectral functions of the retarded Green’s function GR and self-energy ΣR allowing for systems that are out of equilibrium and spatially inhomogeneous. The retarded Green’s function is defined as
(106)GijR(t1,t2)=GijR(T,t)=-iθ(t)〈[a^i(T+t2),a^j†(T-t2)]〉,
where operators are in the Heisenberg picture, T=(t1+t2)/2, t=t1-t2, and the average is taken with respect to the state operator
(107)ρ=e-βH(t→-∞)Z,
where Z=Tr[e-βH(t→-∞)]. The spectral moments of the retarded Green’s function are
(108)μnijR=-Im[in∂n∂tnGijR(T,t)]t=0+,
(note that terms involving derivatives of θ(t) do not contribute to these moments). Freericks et al. [177] calculated the general expressions for these moments and obtained specific analytic results in the strong coupling regime when U≫J. They also compared the moments as calculated within several different methods numerically: VCA, QMC, DMRG, and strong-coupling perturbation theory for the 1 dimensional BHM in the Mott insulating state with n=1. They found that VCA was the best performing method and that both QMC and DMRG had errors on the order of a few percent for the spectral moments. They argued that the reason for the differences in performance is a result of QMC and DMRG having a cutoff that means that they focus on low energy states, whereas VCA, a method that works best at strong coupling, retains information about higher energy states that can be important for higher order spectral moments. These sum rules are an important step towards benchmarking analytical and numerical methods, both in and out of equilibrium. They may prove particularly important for testing methods that try to span the quantum critical point.


#### 3.9.2. Quench fromJ=0 to U=0

One of the few examples of an exact result in this field was provided by Cramer et al. [104], who showed that for a quench from J=0 to U=0 with an initial Mott insulator ground state in 1 dimension, the final state, whilst not the superfluid state, is one that locally maximizes the entropy, in the sense that the entanglement between a particular site and the rest of the chain becomes maximal. In the infinite system size limit, whilst the system is locally relaxed, it never fully relaxes to an equilibrium state, and for large but finite lattices with N sites, there will be recurrences on times of order N. The physical picture behind this result is that for a lattice site (or sites) within the “lightcone” defined by excitations in the system, there is local mixing of the state, whereas contributions from outside the “lightcone” are exponentially suppressed similar to the Lieb-Robinson bound for spin systems [137]. In quenches from the Mott insulator to the superfluid phase, the idea of a “lightcone” reoccurs in both theoretical and experimental results. Cramer et al. argued that their results are quite general and also apply to initial states of the form ∏i=1N‍|mi〉, where the occupations on different sites mi may differ. They also argued that the lightcone physics underlying the arguments in one dimension should apply in higher dimensions, so that there will be local relaxation inside the causal cone.


### 3.10. Summary

Through the variety of techniques that have been applied to quantum quenches in the BHM and the particular calculations discussed in previous sections, several themes become clear. Both analytic and numerical techniques suggest that when there is a quench from superfluid to the Mott insulator, there can be equilibration if the final value ofU/J is not too large [97], but for larger values of U/J a frozen nonequilibrium state results [87, 97, 108, 112], and in the presence of a trap, the Mott insulating regions can impede relaxation through mass transport [94, 99]. These results leave a number of questions that can be pursued both theoretically and experimentally. Kollath et al. [97] asked about the eventual relaxation after a deep quench from superfluid to Mott insulator and whether there may be aging-like phenomena. Natural quantities to look at in the context of out-of-equilibrium physics are two-time correlation functions, such as the two-time density density correlation function
(109)C(r,r′,t,t′)=〈n(r,t)n(r′,t′)〉-〈n(r,t)〉〈n(r′,t′)〉.
If there is a separation of fast and slow timescales for relaxation then one might expect a quantity such as
(110)C(t,t′)=∑rC(r,r,t,t′),
to display behaviour similar to that seen for aging in glassy systems, in which case it could be separated into a fast time-translation invariant piece, and a slower piece that depends on the two times, t and t′:
(111)C(t,t′)=Cfast(t-t′)+Cslow(t,t′),
with Cfast decaying to zero with a timescale that is much shorter than the typical timescale for decay of Cslow, so that Cslow is essentially constant whilst Cfast decays. If aging is present in the Mott insulators for large U/J, then the simplest decay one might expect for Cslow is of the form (with t>t′)
(112)Cslow~f(h(t)h(t′)),
where h(t) is a monotonically increasing function of t, possibly h(t)∝t and f(x)~x-α for some α [178].In going from the Mott insulator to superfluid, there are a variety of results to reconcile. Several groups predict oscillations in the amplitude of the superfluid order parameter [91, 96], which are yet to be seen experimentally [138]. Scaling theory predictions of the residual energy and defect formation probability were also not found to match calculations by Trefzger and Sengupta [91] using a projection operator formalism. Navez and Schützhold [93] found a light-cone-like behaviour of the growth of the amplitude of the superfluid order parameter after a quench with diffusive growth of the phase coherence. Several groups found numerical and analytical evidence of a light-cone defined by a maximum velocity for the growth of density correlations both in quenches from superfluid to the Mott insulator and within the Mott insulating state [98, 99, 152] in 1 dimension. Regarding the quenches from the Mott insulator to superfluid, it would be desirable to see whether signatures of the KZM survive as expected in other methods than the projection operator approach, and whether the results of Navez and Schützhold can be found with other methods as well. It would also be very interesting to see whether the Lieb-Robinson-like bound of a maximal velocity found in the BHM in 1 dimension survives to higher dimensions than one—the results of the 1/Z expansion suggest this may be possible [93] although the velocity found in [93] had quite a different form as a function of J and U than found in one dimension.On the technical side, most of the methods considered are accurate when eitherU is small or J is small, and in some cases in both limits (e.g., the closed time path approaches based on strong coupling theory). Numerical methods such as QMC and DMRG have provided very accurate determinations of the phase diagram in equilibrium but are not easily generalizable either to out of equilibrium situations (QMC) or dimensions higher than 1 (DMRG and TEBD). Other methods such as the Gutzwiller mean field are reasonably straightforward to calculate with and can provide qualitatively correct behaviour for single-particle properties but in the simplest formulation do not contain information about nonlocal correlations. Formulations such as the generalization of Sengupta and Dupuis’ nonperturbative approach to real time in the context of closed time path methods seem to be an avenue that is worth exploring further for studying dimensions higher than 1. For all methods, however, the out-of-equilibrium sum rules introduced by Freericks et al. [177] are a particularly welcome development that can help to benchmark methods and provide consistency checks between methods.


## 3.1. Weak Interactions: The Bogoliubov Approach

The Bogoliubov approach to the BHM is a mean field approach based on the weak interaction limit [47], similar to the Bogoliubov approach to a weakly interacting superfluid in the absence of a lattice. This suggests working in momentum space, reflecting the extended nature of the state and the expectation of a condensate with k=0. Writing the boson annihilation operator in momentum space:
(17)a^i=1Ns∑ka^ke-ik·ri,
with Ns the number of sites, and similarly for the creation operator, one may rewrite the BHM Hamiltonian as [47]
(18)H=∑k(ϵk-μ)a^k†a^k+U2Ns∑k∑k′a^k†a^-k†a^k′a^-k′,
with ϵk=-2J∑j=1dcos(kja) on a d-dimensional cubic lattice. The approximation in the Bogoliubov approach comes from noting that the number of condensate atoms N0≫1, in which case N0=〈a^0†a^0〉≃〈a^0a^0†〉, and hence one can write N0=〈a^0†〉〈a^0〉. Thus, 〈a^0†〉=〈a^0〉=N0, where the phase is chosen so that the expectation values are real. This observation is used to replace the annihilation and creation operators with their average value plus a fluctuation:
(19)a^0†→N0+a^0†;a^0→N0+a^0.
Keeping only terms to quadratic order in the boson operators in the Hamiltonian after this substitution leads to a Hamiltonian that can be diagonalized with a Bogoliubov transformation to
(20)Heff=-12Un0N0+12∑k[ℏωk-(ϵk+Un0)]+∑kℏωkb^k†b^k,
where n0=N0/Ns, and quasiparticle excitations b^k† are related to the a^k† by
(21)(b^kb^-k†)=(ukvkvk*uk*)(a^ka^-k†),
with |uk|2-|vk|2=1 and their energies are given by
(22)ℏωk=ϵk2+2Un0ϵk.
Limitations of the Bogoliubov theory are that it is only valid when the number of atoms in the condensate N0≫1 and hence does not indicate the presence of a Mott insulator phase when the condensate vanishes [47]. This suggests that it will not be useful for giving a quantitative description of out-of-equilibrium dynamics which involve crossing the quantum phase transition from superfluid to Mott insulator. However, as pointed out by Fischer et al. [86], it is still possible to infer some features of these out-of-equilibrium dynamics using a Bogoliubov approach. Fischer et al. [86] took the hopping term in the BHM to have the form
(23)Hhop=J∑ijMija^i†a^j,
where Mij takes into account both the connectivity of the lattice and the possibility of spatially nonuniform hopping. They studied the Heisenberg equation of motion for the operator a^i (setting ℏ=1):
(24)i∂ta^i=J∑jMija^j+Un^ia^i,
focusing on large fillings n≫1, in which case 1/n may be used as a small expansion parameter. They wrote the operator a^ as
(25)a^i=(ψ0+χ^i+ζ^i)A^N^,
where A^=a^Σ(a^Σ†a^Σ)-1/2N^1/2, with a^Σ=∑ia^i, and N^=A^†A^=∑in^i. The leading term in the expansion is ψ0, which describes the condensate and is determined from the solution of the Gross-Pitaevskii equation. The quantum corrections to this mean field are split into the linear corrections, which correspond to quasiparticle excitations of the superfluid, χ^i and nonlinear corrections ζ^i. In the superfluid phase well away from the Mott insulator phase, the non-linear corrections ζ^i may be ignored, provided Utn≪1. Fischer et al. [86] decomposed the operator
(26)a^i=eiϕ^in^i,
into phase ϕ^ and number n^ operators and then related the fluctuations χ^ to phase δϕ^ and number δn^ fluctuations, respectively. The dynamical equations for these fluctuations when expanded in the basis of eigenvectors of the hopping matrix (with eigenvalues λκ where κ are generalized momenta) are
(27)∂tδn^κ=2n0Jλκδϕ^κ,∂tδϕ^κ=-(Jλκ2n0+U)δn^κ,
where U and J may both be time dependent. If U and J are both independent of time, these equations lead to the usual Bogoliubov spectrum
(28)ωκ2=J2λκ2+2n0UJλκ,
found above in (22). When U is taken to be constant and J is allowed to vary as a function of time, Fischer et al. defined a new time variable via dτκ=λκJdt, leading to an equation of motion for number fluctuations
(29)(∂2∂τκ2+1+UJ2n0λκ)δn^κ=0.
The solutions of this equation depend sensitively on the way in which J is taken to zero as a function of time. If J(t) is of the form J(t)=J0(γt)-α, then for α<2 the number fluctuations oscillate in time with decreasing amplitude and frequency; for α=2 the fluctuations decay to zero, and for faster decay of J(t), either as a power law or exponential form J(t)=J0e-γt, the fluctuations end up being frozen in at a finite value. This can be understood as α=2 providing the limit for horizon formation, in analogy to cosmological horizons. For the exponential sweep, the equation of motion for number fluctuations becomes scale invariant (i.e., independent of κ) and the dimensionless parameter governing the physics is ν=Un0/γ, the ratio of an internal energy scale Un0, to an external one, the sweep rate, γ. This allows a notion of slow (ν≫1) and fast (ν≪1) sweeps. The leading order behaviour of the number fluctuations and phase fluctuations was determined to be [86]
(30)〈δn^κ2〉≃n0(1-e-2πν2πν),〈δϕ^κ2〉≃ν(1-e-2πν2πn0)γ2t2.
A weakness of these results is that they come from extrapolation of the mean field theory for the weak coupling superfluid phase, where ψ0 is dominant over linear and non-linear fluctuations, into the strong coupling phase—they also apply only in the limit n0≫1. The approach taken by Fischer et al. to address this issue is to note that if the quench takes place sufficiently quickly such that the strong coupling regime is reached for times t≪1/nU and t≪1/γ, then the state reached via Bogoliubov dynamics can be viewed as an initial state for evolution after the quench with the Hamiltonian H^0. From these considerations, they found that the correlation function usually considered to determine the presence of off-diagonal long range order (ODLRO) takes the form
(31)〈a^i†a^j〉≃n0e-U2t2Δ2(n0),Δ2(n0)=n0(1-e-2πν2πν).
This prediction for the decay of the ODLRO with time could be tested in the momentum distribution function in time of flight experiments either in the peak at k=0, or perhaps in the visibility as a function of time.As mentioned above, Fischer et al. [86] noted the analogy between the dynamics of number fluctuations under a quantum quench and dynamics in the early universe. In the long wavelength limit, where the lattice structure is ignored, that is, for wavelengths λ≫a, the lattice spacing, the equation of motion for number fluctuations may be written as
(32)(∂∂t1J(t)∂∂t-Un0a2∇2)δn^(x,t)=0,
where the speed of sound cs2(t)=Un0a2J(t) decreases as a function of time. If this decrease is sufficiently quick, then there is a horizon beyond which excitations of a given wavelength cannot propagate, given by
(33)Δn(t)=∫t∞dt′cs(t′),
where for J(t)=J0e-γt, the integral converges and one has Δn(t)=(2J0Un0/γ)e-γt/2. This horizon shrinks as t→∞, so at a given time t, modes with wavelengths shorter than Δ will oscillate, but longer wavelength modes with λ≫Δ will be frozen out. Identifying J˙/J with the Hubble constant, the frozen fluctuations of the scalar field are analagous to the quantum fluctuations of the inflaton field [132]. The presence of such a horizon would imply that after a quench to the Mott insulator, there will be frozen in number fluctuations.


## 3.2. Infinite Dimensional Limit

The infinite dimensional limit of the BHM is somewhat artificial, but it is a situation where an exact solution can be obtained, which may have relevance to the finite dimensional case. Sciolla and Biroli [87, 88] took this approach and studied quantum quenches in the infinite dimensional BHM. They took advantage of the fact that the ground state of the BHM is site permutation symmetric for any value of U/J to note that after a quench the state is also site permutation symmetric. This allows one to parametrize the states with the parameters xi, where xi is the fraction of sites with i bosons, i=0,1,2,…. Considering the particular example where the maximum number of bosons per site is 2 (a constraint that is easily relaxed), one has the conditions,
(34)x0+x1+x2=1;N=V(x1+2x2),
where N is the total number of bosons and V is the volume of the system. From these relations, one may specify the state in terms of x=x1. Sciolla and Biroli considered quantum quenches in which the interaction strength is taken from Ui to Uf at time t=0 and studied the subsequent time evolution of x and the superfluid order parameter |Ψ0|2. They found the interesting result that there is a dynamical transition at a value of Ufd=(Ui+Uc)/2, where Uc=3+22 (in units of J), at which there is exponential relaxation to the Mott state. This is illustrated in Figure 4 where it can be seen that in addition to the dynamical transition noted above, the value of 〈|Ψ0|2〉 tends to a finite value at long times, as opposed to the value of 0 expected if the system is in equilibrium. Sciolla and Biroli noted that this behaviour is an artifact of the infinite dimensional limit but suggest that the mean field approximation may capture aspects of the short-time dynamics in finite dimensional systems, so for large quenches one may expect to see freezing in superfluid correlations, and a failure to relax to the Mott state. A study by Snoek [89] using Gutzwiller mean field dynamics also found a dynamical transition of the sort identified by Sciolla and Biroli.Figure 4
Evolution of the microcanonical average of the superfluid order parameter〈|Ψ0|2〉 as a function of Uf for a quench with Ui=0. The mean field value of |Ψ0|2 is shown for comparison [87]. Adapted figure with permission from B. Sciolla and G. Biroli, Phys. Rev. Lett. 105, 220401 (2010). Copyright 2010 by the American Physical Society.


## 3.3. Gutzwiller Mean Field

The Gutzwiller mean field approach [133, 134] is a widely used way to look at the equilibrium phase diagram and out-of-equilibrium dynamics in the BHM in part because it is relatively straightforward to implement numerically (dynamics in systems with as many as 653 sites have been investigated [90]). The approach involves approximating the many-body state of the system as a product of single-site states in the form
(35)|Ψ〉=∏i‍∑n=0nmaxfn(i)|n〉(i),
with the constraint that ∑n‍|fn(i)|2=1 [12, 134]. The state |Ψ〉 is to be viewed as a variational ansatz, and so the fn(i) are to be found by minimizing
(36)〈Ψ|H^|Ψ〉.
For an average of one boson per site, the transition from superfluid to the Mott insulator is found for zJ/U=1/5.8, as in regular mean field theory.The Gutzwiller technique generalizes to time-dependent problems fairly straightforwardly [90, 135] by requiring that
(37)〈Ψ|iddt-H^|Ψ〉=0,
which leads to the equations [12]
(38)iddtfn(i)=U2n(n-1)fn(i)-JΦi*n+1fn+1(i)-JΦinfn-1(i),
where
(39)Φi=∑〈i,j〉〈aj〉=∑〈i,j〉∑nnfn-1(j)*fn(j).As with any mean field theory, the estimate of the critical coupling for the phase transition is not particularly accurate within the Gutzwiller approach. A more important failing is that due to the assumption of a variational state which has a product form, two site correlations factorize into single-site quantities and so the method does not capture correlations involving different sites. This is most problematic for intermediate values ofU/J where the ground state is not close to product form, unlike the two limits U=0 and J=0. Improvements to the Gutzwiller mean field have been made using a variety of approaches [91–93, 136] and can allow for perturbative corrections to short-range correlations [136]. Nevertheless, because of its simplicity, it is a very useful method for gaining an understanding of physics for regimes where exact numerical results are not easily obtained. For example, for out-of-equilibrium dynamics in dimensions higher than 1, for example, [94] (Section 4.2.3), or in the presence of complicated space and time-dependent potentials such as the “Gaussian spoon” consisting of a usual harmonic trap V0 with a Gaussian with velocity v:
(40)V(r,t)=V0(r)+V1e-((x-vt)2+y2)/W2,
suggested by Lundh [95] as a means to excite Mott insulating states.


### 3.3.1. Projection Operator Approach

One suggestion to improve upon Gutzwiller mean field theory that allows the effects of quantum fluctuations to be included both in equilibrium and in dynamics is to use a projection operator approach [91, 92]. The particular projection operator introduced by Trefzger and Sengupta [91] was
(41)Pl=|n-〉〈n-|r⊗|n-〉〈n-|r′,
which lives on the link l between r and r′ and projects on to the manifold of states for which the occupation on sites r and r′ is n-. With this projection operator, the hopping term in the BHM may be written as (introducing the hopping operator Tl and replacing the sum over sites with a sum over links)
(42)-J∑〈rr′〉br†br′=∑lTl=∑l[PlTl+TlPl+Pl⊥TlPl⊥],
where Pl⊥=1-Pl. The terms PlTl+TlPl take the system out of the low energy manifold in Mott state where there are n- bosons on each site. The strategy employed in [91] to eliminate these terms from the Hamiltonian was to use a canonical transformation
(43)S=i∑l[Pl,Tl]U,
which leads to a transformed Hamiltonain
(44)H*=eiSHe-iS,
in which terms are kept up to order z2J2/U:
(45)H*=H0+∑lPl⊥TlPl⊥-1U∑l[PlTl2+Tl2Pl-PlTl2Pl-TlPlTl]-1U∑〈ll′〉[12PlTlTl′-TlPlTl′+12(TlPlPl′Tl′-PlTlTl′Pl′)+H.c.].
A Gutzwiller ansatz was then applied for the canonically transformed state |ψ′〉=eiS|ψ〉,
(46)|ψ′〉=∏r‍∑nfn(r)|n〉.
The phase diagram obtained using this state in [91] is significantly improved from mean field theory—in 2 dimensions the Mott insulator superfluid phase boundary is closer to QMC results than the mean field result is, and in 3 dimensions, the projection operator result is in very close agreement.As with Gutzwiller mean field theory, it is straightforward to generalize the method to allow time dependence. In particular, ifJ(t) is time dependent then the canonical transformation becomes time dependent and one has the following evolution equation:
(47)(iℏ∂t+∂S∂t)|ψ′〉=H*[J(t)]|ψ′〉,
provided J(t)/U≪1. This restricts the initial and final values of J to satisfy Ji/U,Jf/U≪1. On the other hand, it does not place restrictions on the ramp rate τ-1, which may be either fast or slow.Trefzger and Sengupta considered both quenches from the superfluid phase to the Mott insulator and from the Mott phase to the superfluid phase. In the quench from superfluid to the Mott insulator, they did not observe the expected universal scaling of either the residual energy(48)Q=〈ψf|H[Jf]|ψf〉-EG[ψf],
or the defect formation probability (which formed a plateau at large ramp times)
(49)P=1-|〈ψfΨG|ψf〉|2,
where |ψf〉 is the state at time tf after the ramp and |ΨG〉 is the ground state for the final value of J/U. EG[ψf] is the ground state energy determined by minimizing the energy functional for J=Jf. They observed that the timescale on which one might expect to see universality in P and Q can be estimated as the time for a boson to cross the system of size L, which may well exceed the system lifetime, so that the observed behaviour of P and Q is governed by local physics.On quenching from the Mott insulator to the superfluid phase, they found oscillations in the amplitude of the order parameter, with a period that scales likeT~(δJ)-α, where (δJ)=|Jf-Ji| and α=0.35±0.05 for quenches with a final value of J close to the critical J (higher values of J led to multiple frequency components being present). They also found P~(δJ)0.89 and Q~(δJ)1.90, both in disagreement with the exponents expected for a QCP with z=1 [13].The projection operator and canonical transformation used in this method have an advantage over Gutzwiller mean field theory in the Mott state in that they retain low energy excitations while removing high energy excitations via canonical transformation. In the superfluid phase close to the Mott insulator where number fluctuations are not too large, it is likely that this method may still be an improvement over regular Gutzwiller mean field theory, but deeper into the superfluid phase, where number fluctuations are much larger than in the Mott phase, it seems unlikely that this method will provide any improvement.


## 3.3.1. Projection Operator Approach

One suggestion to improve upon Gutzwiller mean field theory that allows the effects of quantum fluctuations to be included both in equilibrium and in dynamics is to use a projection operator approach [91, 92]. The particular projection operator introduced by Trefzger and Sengupta [91] was
(41)Pl=|n-〉〈n-|r⊗|n-〉〈n-|r′,
which lives on the link l between r and r′ and projects on to the manifold of states for which the occupation on sites r and r′ is n-. With this projection operator, the hopping term in the BHM may be written as (introducing the hopping operator Tl and replacing the sum over sites with a sum over links)
(42)-J∑〈rr′〉br†br′=∑lTl=∑l[PlTl+TlPl+Pl⊥TlPl⊥],
where Pl⊥=1-Pl. The terms PlTl+TlPl take the system out of the low energy manifold in Mott state where there are n- bosons on each site. The strategy employed in [91] to eliminate these terms from the Hamiltonian was to use a canonical transformation
(43)S=i∑l[Pl,Tl]U,
which leads to a transformed Hamiltonain
(44)H*=eiSHe-iS,
in which terms are kept up to order z2J2/U:
(45)H*=H0+∑lPl⊥TlPl⊥-1U∑l[PlTl2+Tl2Pl-PlTl2Pl-TlPlTl]-1U∑〈ll′〉[12PlTlTl′-TlPlTl′+12(TlPlPl′Tl′-PlTlTl′Pl′)+H.c.].
A Gutzwiller ansatz was then applied for the canonically transformed state |ψ′〉=eiS|ψ〉,
(46)|ψ′〉=∏r‍∑nfn(r)|n〉.
The phase diagram obtained using this state in [91] is significantly improved from mean field theory—in 2 dimensions the Mott insulator superfluid phase boundary is closer to QMC results than the mean field result is, and in 3 dimensions, the projection operator result is in very close agreement.As with Gutzwiller mean field theory, it is straightforward to generalize the method to allow time dependence. In particular, ifJ(t) is time dependent then the canonical transformation becomes time dependent and one has the following evolution equation:
(47)(iℏ∂t+∂S∂t)|ψ′〉=H*[J(t)]|ψ′〉,
provided J(t)/U≪1. This restricts the initial and final values of J to satisfy Ji/U,Jf/U≪1. On the other hand, it does not place restrictions on the ramp rate τ-1, which may be either fast or slow.Trefzger and Sengupta considered both quenches from the superfluid phase to the Mott insulator and from the Mott phase to the superfluid phase. In the quench from superfluid to the Mott insulator, they did not observe the expected universal scaling of either the residual energy(48)Q=〈ψf|H[Jf]|ψf〉-EG[ψf],
or the defect formation probability (which formed a plateau at large ramp times)
(49)P=1-|〈ψfΨG|ψf〉|2,
where |ψf〉 is the state at time tf after the ramp and |ΨG〉 is the ground state for the final value of J/U. EG[ψf] is the ground state energy determined by minimizing the energy functional for J=Jf. They observed that the timescale on which one might expect to see universality in P and Q can be estimated as the time for a boson to cross the system of size L, which may well exceed the system lifetime, so that the observed behaviour of P and Q is governed by local physics.On quenching from the Mott insulator to the superfluid phase, they found oscillations in the amplitude of the order parameter, with a period that scales likeT~(δJ)-α, where (δJ)=|Jf-Ji| and α=0.35±0.05 for quenches with a final value of J close to the critical J (higher values of J led to multiple frequency components being present). They also found P~(δJ)0.89 and Q~(δJ)1.90, both in disagreement with the exponents expected for a QCP with z=1 [13].The projection operator and canonical transformation used in this method have an advantage over Gutzwiller mean field theory in the Mott state in that they retain low energy excitations while removing high energy excitations via canonical transformation. In the superfluid phase close to the Mott insulator where number fluctuations are not too large, it is likely that this method may still be an improvement over regular Gutzwiller mean field theory, but deeper into the superfluid phase, where number fluctuations are much larger than in the Mott phase, it seems unlikely that this method will provide any improvement.


## 3.4.1/Z Corrections to Mean Field Theory

Mean field theory should be exact in the limit of infinite co-ordinationZ, but for finite dimensions, it is known to be inaccurate. Navez and Schützhold [93] developed a systematic expansion in 1/Z about mean field theory and used this expansion to explore a quench from a Mott insulator to superfluid. They started with reduced density matrices for 1 site (with site index μ):
(50)ρ^μ=Trμ{ρ^},
and 2 sites (site indices μ and ν)
(51)ρ^μν=Trμν{ρ^},
and similarly for higher numbers of sites. The correlated parts of these reduced density matrices, ρ^μνc, may be obtained from
(52)ρ^μν=ρ^μνc+ρ^μρ^ν.
Navez and Schützhold’s approach is based on the observation that the correlations in the Mott insulator have the scaling hierarchy
(53)ρ^𝒮c=O(Z1-|𝒮|),
where |𝒮| is the number of lattice sites in the set of lattice sets 𝒮. They showed that if this hierarchy is satisfied by the reduced density matrices ρ^𝒮c, then it is also satisfied by ∂tρ^𝒮c, so the hierarchy will remain valid for some finite time (of O(lnZ)) [93]. In the Mott phase with one particle per site ρ^μ=|1〉μ〈1|+O(1/Z). By making use of the hierarchy in (53) to neglect higher order reduced density matrices, they obtained an equation of motion for ρ^μνc, which leads to the following equations (with ℏ set to unity):
(54)(i∂t-U+3JTk)fk12=-2JTk(fk11+fk22+1),i∂tfk11=i∂tfk22=2JTk(fk12-fk21).
In (54), the fkij are the spatial Fourier transforms of the particle and hole correlation functions fμν11=Tr{ρ^h^μ†h^ν}, fμν12=Tr{ρ^h^μ†p^ν}, fμν21=Tr{ρ^p^μ†h^ν}, and fμν22=Tr{ρ^p^μ†p^ν}, where the local particle and hole operators are p^μ=|1〉μ〈2| and h^μ=|0〉μ〈1|, respectively, and Tk is (1/z)∑μ‍Tμνeik·(rμ-rν)≃1-k2/2m*, where Tμν is the connectivity matrix of the lattice.Navez and Schützhold considered a quench fromJ/U≪1 to a final much larger J, and hence assumed initial conditions fμν11=δμν, and all other fμνij=0. Writing the boson operators of the BHM as
(55)a^μ=h^μ+2p^μ+r^μ,
where  r^μ includes terms involving occupation numbers larger than 2, which decouple due to the scaling hierarchy assumption, they calculated 〈a^μ†(t)a^ν(t)〉 for an instantaneous quench from J≃0 to a final value of J=Jc(1+ϵ), where 0<ϵ≪1 and Jc is the critical value of J. In this limit, they estimated
(56)〈a^μ†(t)a^ν(t)〉≃𝒩(t)eγt2-(rμ-rν)2/c2,
where 𝒩(t) is weakly time dependent compared to the exponential, γ~J-Jc, and c2=3J(U-J)/m* with c a velocity scale. This form gives a constant propagation speed of the correlations, as found in the Lieb-Robinson bound [137]. The exponent in (56) also shows scaling behaviour in the vicinity of the critical point: as γ→γ′=λγ, the form is invariant under t→t′=t/λ and r→r′=r/λ. In addition to the growth of off-diagonal long range order, Navez and Schützhold also considered the growth of phase coherence by looking at the condensate fraction within a region 𝒮 with |𝒮|≪1. Defining the coarse-grained operator A^𝒮=∑μ∈𝒮a^μ/|𝒮|, which is the homogeneous mode which will give the largest eigenvalue of 〈a^μ†a^ν〉, this mode will have a macroscopic occupation N𝒮=〈A^𝒮†A^𝒮〉≫1 corresponding to the condensate fraction in the region 𝒮 after time t. For c2t2≪|𝒮|,  N𝒮/|𝒮|~1/|𝒮| and for |𝒮|≪c2t2, N𝒮/|𝒮|~eγt. This suggests that for distinct regions 𝒮 and 𝒮′ with 1≪|𝒮|,|𝒮′|≪c2t2, the relative phase between the regions correlates as
(57)〈ei(φ^𝒮-φ^𝒮′)〉≃e-γr2/2c2t,
where r is the distance between the regions and φ^𝒮 is defined via A^𝒮=N^𝒮eiφ^𝒮, and since 〈N^𝒮〉≫1, one can safely write N^𝒮≃N𝒮. This gives a growth of phase correlations which is diffusive, as opposed to ballistic growth of ODLRO found in (56).


## 3.5. Slave Particle Approach

Another approach in which the Hilbert space of the BHM is truncated to three states per site was put forward by Altman and Auerbach [96] who considered the BHM in the form
(58)H=U2∑i(n^i-n-)2-J∑〈ij〉(a^i†a^j+a^j†a^i)-μ∑i(n^i-n-),
and restricted the Hilbert space to states with n-,  n--1 or n-+1 bosons per site, allowing for a pseudospin 1 representation of the problem. This approach is most appropriate in and near the Mott insulating phase where number fluctuations are small. Altman and Auerbach introduced creation operators t^αi†, where |n-+α〉i=t^αi†|0〉i, with α=1,0,-1, subject to the constraint ∑αt^αi†t^αi=1. The original bosons of the BHM are represented as a^i†=n-t^0i†t^-1i+n-+1t^1i†t^0i. They also introduced modified coherent states of the form |Ω(θ,η,ϕ,χ)〉, with
(59)|Ω〉={cos(θ2)t^0†+eiηsin(θ2)×[eiϕsin(χ2)t^1†+e-iϕcos(χ2)t^-1†]cos(θ2)}|0〉,
and then chose a variational wave function of the form
(60)|ΦMF〉=∏i‍|Ω〉i.
The Mott phase is described by θ=0. When θ>0 there is a superfluid order parameter Ψ=n-sinθe-iϕ. Minimizing the variational energy leads to the usual mean field boundaries for the Mott lobes. To study excitations above the ground state, Altman and Auerbach performed the canonical transformation
(61)b^0i=cos(θ2)t^0+12sin(θ2)(t^1i+t^-1i),b^αi=sin(θ2)t^0-12cos(θ2)(t^1i+t^-1i),b^φi=12(t^1i-t^-1i),
subject to the constraint ∑mb^mi†b^mi=1. The mean field variational state may be written as |ΦMF〉=∏i‍b0i†|0〉, with bαi† and bφi† creating excitations about this state. They then used the constraint to eliminate b0 from the Hamiltonian and truncated the resulting Hamiltonian at quadratic order in the operators bα† and bφ†. After a Bogoliuibov transformation one obtains
(62)Hfluc=∑mkωmkβmk†βmk.
The modes associated with the two types of fluctuations in the superfluid phase are a massive amplitude (Higgs) mode and a phase (sound) mode [96]
(63)ωα(k)≃c2k2+Δ2,ωφ(k)≃c|k|,
where c and Δ may be expressed in terms of J,  U,  z, and n-. In the Mott phase, the two modes are degenerate and gapped and represent particle and hole excitations.Altman and Auerbach also constructed a path integral using the modified coherent states, and writing the action in terms ofΨ=n-sinθe-iφ, they obtained the action
(64)S=18zJn-2∫0t‍dt′∫‍ddr{|Ψ˙|2-(2Jn-)2z|∇Ψ|2-(2Jn-z)2(u-1)|Ψ|2-(zJ)2n-u|Ψ|4|Ψ˙|2}.
For large, n- the critical point is at u=U/4zJn-=1. From this action they obtained an equation of motion for the amplitude of the superfluid order parameter (after rescaling)
(65)Ψ¨=c2∇2Ψ+12Δ2Ψ(1-|Ψ|2).
The solution of this equation when Ψ is assumed to be spatially uniform shows oscillations as a function of time for a system prepared in the Mott state, as illustrated in Figure 5 (although these oscillations were not observed in a recent quench experiment [138])(a) Mean field energy in the superfluid state as a function ofΨ with the equiilbrium ground state indicated at (2); (b) solution of (65) as a function of time, compared to the equilibrium value of Ψ (dashed line) [96]. Adapted figure with permission from A. Altman and A. Auerbach, Phys. Rev. Lett. 89, 250404 (2002). Copyright 2002 by the American Physical Society.
(a)
(b)Altman and Auerbach noted some caveats to their results: the assumption of a uniform system may be spoiled by topological defects that get trapped by the Kibble-Zurek mechanism [63, 64]. For initial states near the transition, the correlation length should be large enough that such defects are widely spaced. However, for a quench from deep in the Mott phase, there will be a much higher initial density of defects. Altman and Auerbach argued that the k=0 mode will be the fastest growing, so that as J is increased to near the transition, defects will generally be widely separated by ξ=c/Δ. The question of the effects of the remaining vortices on the evolution of Ψ is an open question. The oscillations shown in Figure 5 are undamped—in general, one would expect that such oscillations would damp as a function of time. Altman and Auerbach argued that such damping could come from the coupling of the amplitude mode to low-energy phasons, which would lead to phason pair emission.The theory discussed in [96] is primarily for n-≫1, motivated by the experiments of Orzel et al. [139], for which n-~50. Further work by Huber et al. [118] extended the slave particle approach to allow for n- of order unity. Using a similar procedure in which they constructed the low energy excitations above the mean field variational ground state, they found qualitatively similar features to Altman and Auerbach, namely, a massive Higgs mode and a gapless sound mode in the superfluid and gapped particle hole excitations in the Mott insulator. The theory has similar caveats to the large n- theory that J/U not be too large and the density of excitations not be too large (which will certainly be satisfied in the Mott phase).


## 3.6. Numerical Methods

Numerical simulations give the opportunity to obtain quantitative results for both equilibrium and out-of-equilibrium dynamics both for the homogeneous BHM and for cold bosons in an optical lattice. Methods that have been used extensively for out-of-equilibrium dynamics include exact diagonalization and the two closely related methods: time-dependent density matrix renormalization group (t-DMRG), and time-evolving block decimation (TEBD). The t-DMRG [140, 141] and TEBD [142] methods are primarily limited to use in one dimension, due to the exponential increase in computational resources required in two dimensions [143], whereas exact diagonalization may be used in any dimension, but the Hilbert space grows exponentially with the number of bosons, limiting the number of sites that may be simulated. In both t-DMRG and TEBD, the Hilbert space is truncated, but low lying energy states are retained, allowing for the simulation of larger systems in 1 dimension than with exact diagonalization. Reviews by Schollwöck [143, 144] give a comprehensive discussion of details of these methods—my focus here is on the physics that has been obtained from numerical calculations for quenches from (i) the superfluid phase to the Mott insulating phase and (ii) the Mott insulator phase to the superfluid phase.


### 3.6.1. Superfluid to Mott Insulator Quench

Kollath et al. [97] used a combination of exact diagonalization and t-DMRG to study the behaviour of the homogeneous BHM with unit filling after a sudden quantum quench from the superfluid phase with an initial value of the interaction strength Ui to a final value Uf≫Ui. Their simulations encompassed up to 18 sites in two dimensions and up to 64 sites in 1 dimension. They calculated the time-and space-dependent correlation functions 〈a^x†(t)a^0(t)〉 (note that Kollath et al. denote the boson creation operators of the BHM by b† rather than a† as used here) to characterize the dynamics after such a quench. This correlation function is seen not to decay to zero, but after an initial relaxation on a timescale of order 1/J is seen to have small amplitude oscillations, with period 2π/Uf about a finite value after a quench, as illustrated in Figure 6.Time dependence of the nearest neighbour correlations forUi=2J and Uf=40J obtained using (a) exact diagonalization and t-DMRG methods in 1 dimension and (b) exact diagonalization in 2 dimensions [97]. Adapted figure with permission from C. Kollath, A. M. Lauchli, and E. Altman, Phys. Rev. Lett. 98, 180601 (2007). Copyright 2007 by the American Physical Society.
(a)
1D chain,Ui=2J,  Uf=40J
(b)
2D square,Ui=2J, Uf=40JKollath et al. investigated the decay of the time-dependent correlation functions with distance after the decay of the initial transient and found that there were two different behaviours, depending on the final value ofUf. For Uf/J not too large, the correlations fit reasonably well to thermal correlations determined with QMC and suggest that the final state can be viewed as thermalized. However, for large Uf, there was a strong memory of the initial state, and the decay of correlations with distance was much slower than found either for a thermal state or in the ground state at U=Uf, as illustrated in Figure 7, which they characterized as a nonthermal steady state.Figure 7
Average spatial correlations as a function of site separation after a quench fromUi=2J to Uf=40J in (a) 1 dimension and (b) 2 dimensions. Squares indicate the results from exact diagonalization and triangles indicate results from t-DMRG. The diamonds show the result of equilibrium QMC calculations, and the circles indicate the zero temperature correlations in the ground state for Ui and Uf [97]. Adapted figure with permission from C. Kollath, A. M. Läuchli, and E. Altman, Phys. Rev. Lett. 98, 180601 (2007). Copyright 2007 by the American Physical Society.Kollath et al. [97] rationalized the slow relaxation they found for large Uf by focusing on the particle and hole excitations of the Mott insulator. If quasiparticle interactions are ignored, then the effective Hamiltonian in the Mott regime may be written in the form H0=∑k,αωkβαk†βαk, with βαk† a creation operator for a quasiparticle in the Mott insulator. Thermalization then proceeds through the relaxation of the quasiparticle distribution given by the initial conditions. This requires processes in which quasiparticle number is not conserved, which are strongly suppressed if the gap Δ is larger than W/2 where W~4zJ is the quasiparticle bandwidth. Higher order processes could still lead to thermalization, but on much longer timescales than accessible numerically. This picture relies on having a dilute quasiparticle population, that is, an initial state close to the transition, but at least in their numerics, Ui=2J, appear to have been sufficiently close to the transition (Uc=3.37J in 1 dimension [44] and Uc=16.7J in 2 dimensions [41]) for this physics to be observable.The single particle correlations studied in [97] and also density-density correlations 〈n^jn^j+r〉 were calculated as a function of time and position with similar Ui and Uf=40J for average fillings 〈n〉=1 and 0.5 by two of the authors of [97] in [98]. They observed that the single particle and density-density correlations showed a “light-cone” effect in that there was a propagating front in the correlation with velocity vs. This behaviour is similar to that demonstrated by Lieb and Robinson [137] for one dimensional spin systems with finite range interactions, that at long distances there is a maximal velocity at which correlations can spread. The values of vs for 〈n〉=1, for which Uf corresponds to the Mott phase, were on the order of 5-6J, whereas the values of vs for 〈n〉=0.53 (for which the system remains superfluid) after the change in U were of the order of 4J. Additionally, Lauchli and Kollath [98] calculated the change in von Neumann entropy for blocks of length l as a function of time: S(l,t)-S(l,0), and found that this showed linear growth for short times and saturated at a time t*(l)≃l/2ve for periodic boundary conditions and t*(l)≃l/ve for open boundary conditions with ve≃vs in accord with predictions of Calabrese and Cardy [145].The dynamics of local observables both in homogeneous systems and in the presence of a parabolic trapping potential were studied by Bernier et al. [99] using a combination of exact diagonalization calculations and t-DMRG to investigate quenches of the form U(t)=Ui+(t/τ)(Uf-Ui) for the 1d BHM. By focusing on local observables, such as the local compressibility κj=〈nj2〉-〈nj〉2, and the probability of occupation of a site by n bosons, Pn, they observed that the dynamics split into two regimes: short ramp times and long ramp times as illustrated in Figure 8. For short ramp times (τ≲ℏ/J), the behaviour near the centre of a homogeneous system and a trap is essentially identical, and indicative of local dynamics, since the ramp timescale is less than the hopping timescale. At longer times, the inhomogeneity of the trap becomes apparent as density redistribution takes place, as seen in different values for the local quantities in the presence and absence of the trap. This behaviour was seen both for a quench from superfluid to the Mott insulator (Ui=2J and Uf=4J) and for a quench within the Mott insulator (Ui=4J and Uf=6J).Figure 8
Dependence of local observables on ramp timeτ for both homogeneous systems and at the centre of a parabolic trap, n: average density per site; κ: compressibility; and Pn: the probability that a site has an occupation of n particles. (a) For a protocol with Ui=2J and Uf=4J, (b) for the same protocol with Ui=4J and Uf=6J. Adapted figure with permission from J. S. Bernier, G. Roux, and C. Kollath, Phys. Rev. Lett. Lett. 106, 200601 (2011). Copyright 2011 by the American Physical Society.In subsequent work investigating local observables and local currents in a trap, Bernier et al. [100] observed that regions with density of order 1 boson per site (even if not in a Mott insulator state) can act as “Mott barriers” that impede particle transport, slowing down dynamics associated with density redistribution in the trap. These ideas are discussed further in the context of quench experiments by Hung et al. [146] in Sections 4.2.2 and 4.2.3.


### 3.6.2. Mott Insulator to Superfluid Quench

Clark and Jaksch [101] used TEBD for the one dimensional BHM to study the effect of varying the lattice depth V0(t) with time on the coherence in ramping from the Mott insulator to the superfluid phase. They considered slow ramps in which the lattice depth was evolved from the superfluid state with U/2J=2 to the Mott insulator with U/2J=40. Comparing the exact numerical solution with TEBD for a system of M=7 bosons, they established that their slow ramp had an infidelity of 1-F<10-4 so that their slow ramp could be viewed as adiabatic. They then considered larger systems with M=49 prepared with a slow ramp and studied how coherence was reestablished as the lattice depth was reduced. To do this, they studied a time-dependent correlation length ξc(t). They found that for a ramp in which V0 was reduced slowly in a similar manner to how it was increased, that the restoration of coherence was on a similar timescale either for a system that had been prepared with a slow ramp or for a system that had been prepared in the ground state for U/2J=20, with the timescale for single atom hopping governing the establishment of coherence. In contrast, for a system prepared with a slow ramp and then subjected to a fast ramp, Clark and Jaksch found that ξc increased quickly for short ramp times and reached the value expected in the superfluid on much shorter timescales than that of single-particle hopping, leading them to suggest that higher-order correlations must play a role in the establishment of coherence during a fast quench.


### 3.6.3. Entanglement Growth

Numerical simulation of the time dependence of quantum systems using approaches based on matrix product states, such as t-DMRG and TEBD, is limited in evolution times due to the growth of entanglement [145, 147]. Motivated by an experiment in which there was a quantum quench for which relaxation dynamics persisted to timescales not accessible with DMRG [148], Daley et al. [149] focused on the growth of entanglement during a quench from a Mott insulator to a superfluid and proposed a protocol to measure this entanglement growth by obtaining the arbitrary order Renyi entropies
(66)Sn(ρ)=11-nlogTr{ρn}.
They tested the scheme with t-DMRG to show its robustness against measurement imperfections. Implementation of this scheme experimentally could lead to quantitative insights into the performance of quantum simulators [149].


## 3.6.1. Superfluid to Mott Insulator Quench

Kollath et al. [97] used a combination of exact diagonalization and t-DMRG to study the behaviour of the homogeneous BHM with unit filling after a sudden quantum quench from the superfluid phase with an initial value of the interaction strength Ui to a final value Uf≫Ui. Their simulations encompassed up to 18 sites in two dimensions and up to 64 sites in 1 dimension. They calculated the time-and space-dependent correlation functions 〈a^x†(t)a^0(t)〉 (note that Kollath et al. denote the boson creation operators of the BHM by b† rather than a† as used here) to characterize the dynamics after such a quench. This correlation function is seen not to decay to zero, but after an initial relaxation on a timescale of order 1/J is seen to have small amplitude oscillations, with period 2π/Uf about a finite value after a quench, as illustrated in Figure 6.Time dependence of the nearest neighbour correlations forUi=2J and Uf=40J obtained using (a) exact diagonalization and t-DMRG methods in 1 dimension and (b) exact diagonalization in 2 dimensions [97]. Adapted figure with permission from C. Kollath, A. M. Lauchli, and E. Altman, Phys. Rev. Lett. 98, 180601 (2007). Copyright 2007 by the American Physical Society.
(a)
1D chain,Ui=2J,  Uf=40J
(b)
2D square,Ui=2J, Uf=40JKollath et al. investigated the decay of the time-dependent correlation functions with distance after the decay of the initial transient and found that there were two different behaviours, depending on the final value ofUf. For Uf/J not too large, the correlations fit reasonably well to thermal correlations determined with QMC and suggest that the final state can be viewed as thermalized. However, for large Uf, there was a strong memory of the initial state, and the decay of correlations with distance was much slower than found either for a thermal state or in the ground state at U=Uf, as illustrated in Figure 7, which they characterized as a nonthermal steady state.Figure 7
Average spatial correlations as a function of site separation after a quench fromUi=2J to Uf=40J in (a) 1 dimension and (b) 2 dimensions. Squares indicate the results from exact diagonalization and triangles indicate results from t-DMRG. The diamonds show the result of equilibrium QMC calculations, and the circles indicate the zero temperature correlations in the ground state for Ui and Uf [97]. Adapted figure with permission from C. Kollath, A. M. Läuchli, and E. Altman, Phys. Rev. Lett. 98, 180601 (2007). Copyright 2007 by the American Physical Society.Kollath et al. [97] rationalized the slow relaxation they found for large Uf by focusing on the particle and hole excitations of the Mott insulator. If quasiparticle interactions are ignored, then the effective Hamiltonian in the Mott regime may be written in the form H0=∑k,αωkβαk†βαk, with βαk† a creation operator for a quasiparticle in the Mott insulator. Thermalization then proceeds through the relaxation of the quasiparticle distribution given by the initial conditions. This requires processes in which quasiparticle number is not conserved, which are strongly suppressed if the gap Δ is larger than W/2 where W~4zJ is the quasiparticle bandwidth. Higher order processes could still lead to thermalization, but on much longer timescales than accessible numerically. This picture relies on having a dilute quasiparticle population, that is, an initial state close to the transition, but at least in their numerics, Ui=2J, appear to have been sufficiently close to the transition (Uc=3.37J in 1 dimension [44] and Uc=16.7J in 2 dimensions [41]) for this physics to be observable.The single particle correlations studied in [97] and also density-density correlations 〈n^jn^j+r〉 were calculated as a function of time and position with similar Ui and Uf=40J for average fillings 〈n〉=1 and 0.5 by two of the authors of [97] in [98]. They observed that the single particle and density-density correlations showed a “light-cone” effect in that there was a propagating front in the correlation with velocity vs. This behaviour is similar to that demonstrated by Lieb and Robinson [137] for one dimensional spin systems with finite range interactions, that at long distances there is a maximal velocity at which correlations can spread. The values of vs for 〈n〉=1, for which Uf corresponds to the Mott phase, were on the order of 5-6J, whereas the values of vs for 〈n〉=0.53 (for which the system remains superfluid) after the change in U were of the order of 4J. Additionally, Lauchli and Kollath [98] calculated the change in von Neumann entropy for blocks of length l as a function of time: S(l,t)-S(l,0), and found that this showed linear growth for short times and saturated at a time t*(l)≃l/2ve for periodic boundary conditions and t*(l)≃l/ve for open boundary conditions with ve≃vs in accord with predictions of Calabrese and Cardy [145].The dynamics of local observables both in homogeneous systems and in the presence of a parabolic trapping potential were studied by Bernier et al. [99] using a combination of exact diagonalization calculations and t-DMRG to investigate quenches of the form U(t)=Ui+(t/τ)(Uf-Ui) for the 1d BHM. By focusing on local observables, such as the local compressibility κj=〈nj2〉-〈nj〉2, and the probability of occupation of a site by n bosons, Pn, they observed that the dynamics split into two regimes: short ramp times and long ramp times as illustrated in Figure 8. For short ramp times (τ≲ℏ/J), the behaviour near the centre of a homogeneous system and a trap is essentially identical, and indicative of local dynamics, since the ramp timescale is less than the hopping timescale. At longer times, the inhomogeneity of the trap becomes apparent as density redistribution takes place, as seen in different values for the local quantities in the presence and absence of the trap. This behaviour was seen both for a quench from superfluid to the Mott insulator (Ui=2J and Uf=4J) and for a quench within the Mott insulator (Ui=4J and Uf=6J).Figure 8
Dependence of local observables on ramp timeτ for both homogeneous systems and at the centre of a parabolic trap, n: average density per site; κ: compressibility; and Pn: the probability that a site has an occupation of n particles. (a) For a protocol with Ui=2J and Uf=4J, (b) for the same protocol with Ui=4J and Uf=6J. Adapted figure with permission from J. S. Bernier, G. Roux, and C. Kollath, Phys. Rev. Lett. Lett. 106, 200601 (2011). Copyright 2011 by the American Physical Society.In subsequent work investigating local observables and local currents in a trap, Bernier et al. [100] observed that regions with density of order 1 boson per site (even if not in a Mott insulator state) can act as “Mott barriers” that impede particle transport, slowing down dynamics associated with density redistribution in the trap. These ideas are discussed further in the context of quench experiments by Hung et al. [146] in Sections 4.2.2 and 4.2.3.


## 3.6.2. Mott Insulator to Superfluid Quench

Clark and Jaksch [101] used TEBD for the one dimensional BHM to study the effect of varying the lattice depth V0(t) with time on the coherence in ramping from the Mott insulator to the superfluid phase. They considered slow ramps in which the lattice depth was evolved from the superfluid state with U/2J=2 to the Mott insulator with U/2J=40. Comparing the exact numerical solution with TEBD for a system of M=7 bosons, they established that their slow ramp had an infidelity of 1-F<10-4 so that their slow ramp could be viewed as adiabatic. They then considered larger systems with M=49 prepared with a slow ramp and studied how coherence was reestablished as the lattice depth was reduced. To do this, they studied a time-dependent correlation length ξc(t). They found that for a ramp in which V0 was reduced slowly in a similar manner to how it was increased, that the restoration of coherence was on a similar timescale either for a system that had been prepared with a slow ramp or for a system that had been prepared in the ground state for U/2J=20, with the timescale for single atom hopping governing the establishment of coherence. In contrast, for a system prepared with a slow ramp and then subjected to a fast ramp, Clark and Jaksch found that ξc increased quickly for short ramp times and reached the value expected in the superfluid on much shorter timescales than that of single-particle hopping, leading them to suggest that higher-order correlations must play a role in the establishment of coherence during a fast quench.


## 3.6.3. Entanglement Growth

Numerical simulation of the time dependence of quantum systems using approaches based on matrix product states, such as t-DMRG and TEBD, is limited in evolution times due to the growth of entanglement [145, 147]. Motivated by an experiment in which there was a quantum quench for which relaxation dynamics persisted to timescales not accessible with DMRG [148], Daley et al. [149] focused on the growth of entanglement during a quench from a Mott insulator to a superfluid and proposed a protocol to measure this entanglement growth by obtaining the arbitrary order Renyi entropies
(66)Sn(ρ)=11-nlogTr{ρn}.
They tested the scheme with t-DMRG to show its robustness against measurement imperfections. Implementation of this scheme experimentally could lead to quantitative insights into the performance of quantum simulators [149].


## 3.7. Fermionization

In the limit thatU/J→∞, there can be no double occupancy of bosons on a site, which leads to the notion of hard-core bosons, which obey bosonic commutation relations if they are not on the same site, that is, [a^i,a^j†]=[a^i,a^j]=[a^i†,a^j†]=0 if i≠j, but have fermionic commutation relations {a^i,a^i†}=1, a^i2=[a^i†]2=0 on the same site. In one dimension, the hard-core boson model can be mapped onto fermions by the use of the Jordan-Wigner transformation [150]
(67)a^i†=f^i†∏j=1i-1‍e-iπn^jf,a^i=∏j=1i-1‍eiπn^jff^i,
where the f^i† and f^i are fermionic creation and annihilation operators on site i, respectively, and n^if=f^i†f^i is the fermionic number operator on site i. This mapping removes the hard-core constraint and allows for the efficient simulation of out-of-equilibrium dynamics [150, 151] of hard-core bosons.An alternative fermionization scheme which allows for a treatment of out-of-equilibrium dynamics of the BHM whenU is finite was recently introduced by Barmettler et al. [152] for the 1-dimensional BHM in the strongly interacting regime, where there are small fluctuations around integer filling. They truncated the Hilbert space to three states on each site, similarly to the slave particle approach, but introduce different auxiliary bosons b^j±†, such that
(68)a^j†=n-+1b^j+†+n-b^j-,
and for the three states on site j: |n-+m〉j with m=±1,0, the operators satisfy
(69)b^j+†|n-〉j=|n-+1〉,b^j-†|n-〉j=|n--1〉j,bj±|n-〉j=0,
with additional hard core constraints to ensure no double occupancy: (b^j±†)2=b^j±2=0. In addition, there is not allowed to be double occupancy of different species of bosons, which implies the constraint
(70)b^j,+†b^j,+b^j-†b^j-=n^j+n^j-=0.
Barmettler et al. used a Jordan-Wigner transformation similar to (67), to introduce two species of fermions (reflecting two species of bosons):
(71)b^j,σ=Zj,σc^j,σ,
where
(72)Zj,+=eiπ∑σ,j′<jn^j′,σ,Zj,-=Zj,+eiπn^j,+,
are nonlocal string operators which satisfy Zj,σ†=Zj,σ and Zj,σ2=1. This form removes the hard-core constraints on b^j+ and b^j- individually, but not on the product n^j+n^j-. This condition may be satisfied with the projector 𝒫=∏j‍Pj, where Pj=1-n^j+n^j-, in which case the BHM Hamiltonian may be represented by
(73)H=∑j‍𝒫{U2-J(n-+1)c^j+†c^j+1,+-Jn-c^j+1,-†c^j--Jn-(n-+1)(c^j+†c^j+1,-†-c^j-c^j+1,+)+H.c.+U2(n^j++n^j-)}𝒫.
In order to circumvent the complications of calculating with the projector in the Hamiltonian, Barmettler et al. studied the problem with the approximation of unconstrained fermions (UF), that is, 𝒫=1. With this assumption, the Hamiltonian may be diagonalized with a Bogoliubov transformation to the form
(74)HUF=∑k,σϵσ(k)γ^k,σ†γ^k,σ,
where ϵσ(k)=-σJcos(k)+ℏω(k) and
(75)ℏω(k)=12[E+(k)+E-(k)]2+4|Δ(k)|2,
with
(76)E+(k)=-2J(n-+1)cos(k)+U2,E-(k)=-2Jn-cos(k)+U2,Δ(k)=2iJn-(n-+1)sin(k).
Barmettler et al. found that for U/J<4(n-+1), the UF approximation breaks down.Quench in the Mott Phase. In order to study the out-of-equilibrium dynamics, Barmettler et al. chose an initial Fock state with filling n- and considered an instantaneous quench from this J/U=0 state to a final J/U which remains within the Mott phase. They calculated the time-dependent density-density correlation function
(77)Cd(t)=〈n^j(t)n^j+d(t)〉-〈n^j(t)〉〈nj+d(t)〉,
for sites separated by d and found that for d>1,
(78)Cd(t)≃-(n-(n-+1)JdU)2(Jd(J~t)J~t)2,
where J~=2J(2n-+1)/ℏ and Jd is a Bessel function.Focusing in particular on the case ofn-=1, they studied density correlations Cd(t) with both DMRG and the UF approximation for U/J=18 and U/J=9. These two methods both showed close agreement and a propagating front of correlations, corresponding to doublon-holon pairs, with the same constant velocity in both methods. By approximating the expression (78) in the d≫1 limit as
(79)Cd(t)~-(2d2/321/3ℏ3Ut)2Ai2[-(2d)1/3(6Jtℏ-d)],
Barmettler et al. made use of the fact that the Airy function Ai(-21/3z) has a peak for z0≃1.02 to infer that there is a propagation front with correlations exponentially suppressed for t<tpeak where
(80)Jtpeakℏ≃16[d+z0(d2)1/3],
from which it can be inferred that v∞=6J/ℏ is the velocity at which correlations spread as d→∞ for large U/J. Barmettler et al. provided a similar analysis for noninteracting bosons to deduce v∞=4J/ℏ when U/J=0. In the regime of validity of the UF approximation, they considered the relative velocity of quasiparticle pairs analytically, which takes a maximum value
(81)vmax=maxk|v(k)|=maxk|2ddkω(k)|≃2J(2n-+1)ℏ(1-8n-(n-+1)J2(2n-+1)2U2)+O(J4U3),
which corresponds to vmax=6J/ℏ for n-=1 and decreases monotonically with increasing J/U from the J/U=0 limit.The fermionization procedure introduced by Barmettler et al. gives very nice results that go beyond previous analytical work on the one-dimensional BHM and, as discussed in Section4.4, also give good agreement with experiment. The derivation of a velocity which governs the spreading of correlations at large distances is an important result, as it gives a concrete value for a Lieb-Robinson-like bound on the speed at which correlations spread. The Lieb-Robinson bound was originally formulated for interacting spins on a lattice [137, 153–155], but these results indicate that it also appears to hold in the Mott insulating phase of the one-dimensional BHM.


## 3.8. Closed Time Path Approaches

The closed time path (CTP) or Schwinger-Keldysh technique [156–161] is an approach that allows for the calculation of both equilibrium and out-of-equilibrium quantum dynamics within the same formalism. In CTP methods the problem is formulated in real time (rather than imaginary time, as for the Matsubara formalism), allowing out-of-equilibrium problems to be tackled without analytic continuation. The price to pay is that the number of fields in the theory doubles, with a second copy of each field propagating backwards in time. Consequently in calculations of Green’s functions, the notion of time ordering needs to be replaced by that of contour ordering [159] on a contour such as illustrated in Figure 9. There are several reviews that provide a detailed introduction to CTP techniques [158–161], so few details will be given here.Figure 9
Contour for the Schwinger-Keldysh technique.Several authors have used these approaches to study the BHM both from the weakly interacting limit (largeJ/U) [162–166] and the strongly interacting limit (small J/U) [112, 120, 167, 168]. The Schwinger-Keldysh approach has also been used to study out-of-equilibrium dynamics for the BHM with periodic driving [20], including allowing for the possibility of ohmic dissipation [20, 169]. I first discuss the two-particle irreducible formalism which has a weak interaction starting point and then discuss the strongly interacting approach.


### 3.8.1. Two Particle Irreducible Formalism

The two particle irreducible (2PI) closed-time-path (CTP) formalism for the BHM has been developed in detail by Rey et al. [162] and Temme and Gasenzer [165] who studied the 1/𝒩 expansion.The 2PI formalism can be obtained by starting from the action for the complex fieldsΦ and Φ* in the Bose Hubbard model in the form
(82)S[Φ]=∫‍dt∑iΦia(t)τabyℏ∂tΦib(t)+∫‍dt{U4𝒩(JΦi+1a(t)τabxΦ1b(t))-U4𝒩[τabxΦia(t)Φib(t)]2},
where a compact notation has been introduced, with Φi1=Φi and Φi2=Φi*, and the matrices τx and τy have the same form as the corresponding Pauli matrices but are labelled so as not to be confused with spin. The quantities studied are ϕia(t)=〈Φia(t)〉 and Gijab(t,t′)=〈TCΦia(t)Φib(t′)〉-ϕia(t)ϕib(t′), where TC  indicates contour ordering. The theory is specified by the effective action Γ[ϕ,G], which takes the form
(83)Γ[ϕ,G]=S[ϕ]+i2Tr[lnG-1]+i2Tr[D-1(ϕ)G]+Γ2[ϕ,G],
where
(84)iDijab(t,t′)-1=δS[ϕ]δϕia(t)δϕjb(t′),
and Γ2 is the two-particle irreducible generating functional. Γ2 can be described as an infinite sum of diagrams from the expansion of interaction terms in S[ϕ+φ] of cubic and higher order in φia=Φia-ϕia [162].In practice, various approximation schemes keep some subset of these diagrams and can be classified based on which terms are kept. The simplest approximation is to assume thatΓ[ϕ,G]=S[ϕ], corresponding to the mean field approximation. The Bogoliubov (one-loop) approximation can be obtained if the effective action is taken to be (83) with Γ2 ignored. Approximations that take into account interactions more carefully include the time-dependent Hartree-Fock-Bogoliubov (HFB) approximation, in which Γ2 is truncated to keep only the first-order diagram in U. This approximation conserves energy and particle number but violates Goldstone’s theorem. The expansion can be viewed as one in Ut/J, where t is time, and hence will be most accurate for short times or small U/J. It is possible to keep higher-order terms in Γ2—second-order expansions were considered by Rey et al. [162] and Temme and Gasenzer [165], who both considered a 1/𝒩 expansion, which simplifies the second-order expansion of Γ2. From each of these various effective actions, saddle point equations of motion can be derived, which may be quite lengthy; the full expressions for the second order expansion are listed in Appendix B of [162]. These equations of motion may also be used to derive a quantum kinetic equation, as was done by Rey et al. in [163].An additional issue discussed in [162] is that of conservation laws: for a closed system, the evolution is unitary in time, so energy must be conserved. However, particular approximation schemes for Γ may not conserve energy. On the other hand, schemes which are ϕ derivable, that is, the self-energy may be written as the derivative of a functional of G [170] do conserve particle number, energy, and momentum, which is the satisfied for the cases considered in [162]. A comparison of the results obtained for the approximation schemes outlined in [162] are shown in Figure 10 and compared against the results of exact dynamical evolution in small systems. The condensate population as a function of time is shown for initial conditions in which the condensate starts all in a single well, for N=20, 40, and 80 particles with J=1/2 and NU/J=4. The Hartree-Fock Bogoliubov approximation is seen to be fairly poor at describing the time development of the condensate, with the 1/N approach giving the most accurate account of the time evolution for the condensate population per well as a function of time as N becomes larger.Figure 10
Condensate population per well in each well of a two-well system and total condensate population as a function of time (in units ofℏ/J) for parameters J=1/2, NU/J=4, and N=20,40,80 as calculated exactly, using second-order expansion, 1/N expansion, and the Hatree-Fock Bogoliubov approximation [162]. Adapted figure with permission from A. M. Rey, B. L. Hu, E. Calzetta, A. Roura, and C. W. Clark, Phys. Rev. A 69, 033610 (2004). Copyright 2004 by the American Physical Society.One can see that the range of couplings considered are very much in the weakly interacting limit—this suggests the need for other methods to tackle questions relating to the crossing of the superfluid Mott insulator phase transition.


### 3.8.2. Strong Coupling

Sengupta and Dupuis [119] developed a strong coupling nonperturbative approach to the BHM with features that make it particularly attractive for generalization to real time. In particular, their approach is exact in the limit of both small (U=0) and large (J=0) interactions and obtains the expected features in the excitation spectrum in both the superfluid and the Mott insulating phases. It also naturally allows for the calculation of nonlocal correlations. I first briefly recap the equilibrium calculation of Sengupta and Dupuis and then discuss generalizations of this approach using the Schwinger-Keldysh method [112, 120].Sengupta and Dupuis started with a path integral form for the BHM and then performed two Hubbard Stratonovich transformations to write the generating functional as(85)Z=Z0∫‍𝒟[ψ*,ψ,ϕ*,ϕ]e(ψ∣Jψ)-[(ψ∣ϕ)+c.c.]+W[ϕ*,ϕ],
where Z0 is the generating functional when J=0, with the notation (ϕ∣ψ)=∫0β‍dτ∑r‍ϕ*(r)ψ(r), where τ is imaginary time and β is inverse temperature. The term
(86)W[ϕ*,ϕ]=∑R=1∞(-1)R(R!)2∑a1⋯aRG{ai,bi}Rcϕa1*⋯ϕaR*ϕbR⋯ϕb1,
that enters in the action is evaluated with all the fields at the same site, and
(87)G{ai,bi}Rc=(-1)R〈ψa1⋯ψaRψbR*⋯ψb1*〉,
is the connected local 2R point Green’s function. By integrating out the ϕ field in the generating functional, Sengupta and Dupuis obtained an effective action that to quartic order in the field ψ takes the form
(88)S[ψ*,ψ]=-∑a,bψa*(Gab-1+Jab)ψb+14∑a1,a2,b1,b2Γa1,a2,b1,b2IIψa1*ψa2*ψb2ψb1,
where G-1 is the Green’s function in the local limit and ΓII is the two particle vertex in the local limit. Sengupta and Dupuis approximated ΓII by its static value and introduced g=(1/2)ΓII|static. This leads to the effective action
(89)S=-∫0β‍dτdτ′∑r,r′‍ψr*(τ)[G-1(r,τ,r′,τ′)+Jr,r′δ(τ-τ′)]ψr′(τ′)+g2∫0β‍dτ∑rψr*ψr*ψrψr,
which has the attractive feature that it leads to the correct Green’s function -〈ψr(τ)ψr′*(τ′)〉 in both the strong interaction (J=0) and noninteracting (U=0) limits. This feature makes this approach particularly appealing for the study of out-of-equilibrium dynamics, since it gives the hope that one can accurately represent the behaviour of the system over a wide range of parameters that cover both the superfluid and Mott insulating phases.By taking a saddle point approximation to the action, and then expanding to quadratic order in fluctuations, Sengupta and Dupuis obtained the following results for the excitations in the Mott and superfluid phases, respectively. In the Mott phase, they found that the Green’s function𝒢(k,iω)=-〈ψ(k,iω)ψ*(k,iω)〉 has the form
(90)𝒢(k,iω)=1-zkiω-Ek-+zkiω-Ek+,
where
(91)Ek±=-δμ+ϵk2±12[ϵk2+4ϵkUx+U2]1/2,
with ϵk the Fourier transform of the hopping -Jr,r′,  x=n0+1/2 and
(92)zk=Ek++δμ+UxEk+-Ek-,
with δμ=μ-U(n0-1/2) where n0 is the number of bosons per site for chemical potential μ when J=0. In the Mott insulating phase, the modes Ek± are gapped and the system becomes unstable to superfluidity when the energy of one or other of the modes vanishes for k=0. In the superfluid phase, the expressions for the spectrum are more complicated, and four modes are found ±Ek±, two of which are gapless as k→0 and have the form Ek-=c|k|, and the other two are gapped, as illustrated in Figure 11.Figure 11
Excitation energies±Ek± for n=1 and U/J=20 [119]. Adapted figure with permission from K. Sengupta and N. Dupuis, Phys. Rev. A 71, 033629 (2005). Copyright 2005 by the American Physical Society.Both  Grass et al.  [120] and Kennett and Dalidovich [112] used the Schwinger-Keldysh technique to construct theories for the real-time dynamics of the BHM.  Grass et al.  developed a real-time Ginzburg-Landau theory based on the field theoretic approach developed in [171, 172] and used this theory to calculate collective modes of the BHM, finding amplitude and sound modes in the superfluid and particle and hole excitations in the Mott insulator, as expected.  Grass et al.  also illustrated the power of this approach by expanding the Green’s functions in the equations of motion in the superfluid regime, they were able to obtain the equation of motion
(93)i∂Ψi∂t=-∑jJijΨj-μΨi+U2Ψi|Ψi|2,
which is the lattice version of the Gross-Pitaevski equation (GPE) [102].A similar approach developed by Kennett and Dalidovich [112] generalizes the theory of Sengupta and Dupuis to real time. We determined the effective action to quartic order for the BHM within the strong-coupling Schwinger-Keldysh approach and then used this to derive equations of motion to study the out of equilibrium dynamics.The approach followed in Kennett and Dalidovich [112] was to write a path integral for the BHM
(94)𝒵=∫‍[𝒟a*][𝒟a]eiSBHM[a*,a],
where aia is the field at site i on contour a, with a=1 or 2, and with the notation τi representing the ith Pauli matrix that acts in Keldysh rather than spin space. We performed a Keldysh rotation from the 1, 2 basis to the quantum (q), classical (c) basis [168, 173–175]:
(95)(a1(t)a2(t))→(a~q(t)a~c(t))=L^(a1(t)a2(t)),
where  aq and ac are the quantum and classical components of the field and
(96)L^=12(1-111).
The effect of this on the action is that τ3 in the 1, 2 basis becomes τ1 in the q, c basis. Similarly to Sengupta and Dupuis, we performed two Hubbard-Stratonovich transformations, the first of which decouples the hopping term and the second leads to the effective action to quartic order in terms of z fields (one can show that the Green’s functions for z are the same as those for the original field a [112, 119]) as
(97)Seff[z*,z]=∫‍dt∑ijzia*(t)(2Jij)τab1zjb(t)+∫‍dt1dt2∑izia1*(t1)[Gia2a1(t2,t1)]-1zia2(t2)+14∫‍dt1dt2dt3dt4∑iua1a2a3a4(t1,t2,t3,t4)×zia1*(t1)zia2*(t2)zia3(t3)zia4(t4),
where
(98)ua1a2a3a4(t1,t2,t3,t4)=14∫‍dt5dt6dt5′dt6′Gia5a6a6′a5′2c(t5,t6,t6′,t5′)×{[Gia5a1(t5,t1)]-1[Gia6a2(t6,t2)]-1×[Gia3a5′(t3,t5′)]-1[Gia4a6′(t4,t6′)]-1+((a4,t4)⟷(a3,t3))+((a6,t6)⟷(a5,t5))+((a6,t6)⟷(a5,t5);(a4,t4)⟷(a3,t3))[Gia4a6′(t4,t6′)]-1}.
The matrix Green’s functions G contains the retarded, advanced, and Keldysh Green’s functions. These Green’s functions, along with the two-particle connected Green’s function G2c are evaluated with respect to the single-site Hamiltonian H^0: full expressions are presented in [112]. The mean field phase boundary can be determined from the effective action equation (97) from the vanishing of the coefficient of zq*zc. The diagram corresponding to the interaction term u is illustrated schematically in Figure 12.Figure 12
Schematic diagram for the coefficient of the quartic term in the effective actionua1a2a3a4(t1,t2,t3,t4).The following symmetry relations hold for the interaction kernelu from the definition above:
(99)uabcd(t1,t2,t3,t4)=ubacd(t2,t1,t3,t4)=uabdc(t1,t2,t4,t3),
(this corrects [112]) and it may also be shown [112] that
(100)Gia5a6a6′a5′2c(t5,t6,t6′,t5′)=Gia6a5a6′a5′2c(t6,t5,t6′,t5′)=Gia5a6a5′a6′2c(t5,t6,t5′,t6′).
Grass et al. [120, 168] also noted similar symmetry relations for this correlation function under interchange of Keldysh indices. These symmetries imply that G2c has only 8 independent components: Gqqqq2c,Gcqqq2c,Gqqqc2c,Gqqcc2c,Gccqq2c,Gcqcq2c,Gqccc2c, and Gcccq2c. The remaining four-point function Gcccc2c=0 by causality [161]. Explicit expressions for each of the nontrivial components were written down in [112].A simplified low frequency equation of motion for the superfluid order parameterz can be derived from the action that applies away from the degeneracy points of the Mott lobes (in which case only the component Gcqqq2c is required). Equations of motion that allow for high frequencies and spatial variations of the order parameter and its correlations can also be derived. Kennett and Dalidovich focused on k=0 and defined δ(t)=-2dj(t), where the hopping is taken to be J(t)=J0+j(t) with J0 chosen so that
(101)2dJ0+ν=0,
(i.e., J0 lies on the mean field phase boundary for a given μ). Hence, the approximate equation of motion is
(102)κ2∂2z∂t2+iλ∂z∂t+δ(t)z+u|z|2z=0,
where the coefficients κ2, λ, and u can all be calculated in terms of U and μ [112]. We studied (102) for fixed μ and time-varying J.For fixedμ there are two possibilities for the dynamics: (a) the particle-hole symmetric case (the transition at the tip of the Mott lobe), in which case λ=0, and (b) the generic case, in which λ≠0. The function δ(t) controls the traversal of the quantum critical point. We chose
(103)limt→-∞δ(t)=-δ0;limt→∞δ(t)=δ1.
In our numerics, we used
(104)δ(t)=(δ0+δ12)tanh(tτQ)+δ1-δ02,
where we introduced a timescale τQ which is the characteristic time for δ(t) to cross from -δ0 to δ1 [111]. In the vicinity of the transition t=0 is linear in t. Functional forms of δ(t) which are not linear in t in the vicinity of t=0 may lead to differing behaviour in the long-time limit [79].We solved (102) numerically with λ=0 for a variety of values of τQ, and found that for t≫τQ the form of the solution is that z(t) oscillates in a periodic manner with a magnitude that decreases with increasing τQ. When averaged over a period T at times t≫τQ,
(105)〈z〉T=1T∫tt+T‍dt~z(t~)=0,
as would be expected in the Mott insulating state. Defining zmax(τQ)=limt→∞|z(t)| we found that zmax(τQ) decreases with increasing τQ without any indication of saturation.In the generic case in whichλ≠0, we solved (102) numerically for a variety of values of τQ as is displayed with time scaled by τQ in Figure 13 for μ/U=0.25 (well away from both degeneracy and the particle hole symmetric case). Both the real and imaginary parts of z(t) are oscillatory functions, each of which has a vanishing average over a period for times greater than zero. The nonzero amplitude of z at long times indicates that there is memory of the nonzero value of ρ in the superfluid phase and that the final state is a metastable state rather than the Mott insulator in which ρ=0. The collapse of traces with different τQ after scaling the time with τQ indicates that |z| has a finite value even when τQ becomes large.Figure 13
Dynamics of|z(t)| normalized to unity in the generic case, for a variety of τQ with time rescaled by τQ for a quench from 2dJ/U=2.0 to 2dJ/U=0.0. Figure adapted from [112].To determine the fate of the system at long quench times, we definedzmax=limτQ→∞zmax(τQ) and calculated this as a function of chemical potential. Apart from the particle-hole symmetric point, it appears that there is a transition to a metastable state in which zmax≠0. Whilst it is true that the spatial dependence of z and higher frequency components were ignored in (102), it seems reasonable that the metastable state identified at the mean field level may persist even when these are taken into account.The calculations discussed above are for fixed chemical potential rather than fixed particle number, as is found in a cold atom system. However, within the LDA fixed chemical potential would correspond to studying dynamics at fixed radius. The results found for the long-time limit ofzmax are in accordance with this idea. There appears to be equilibration as zmax→0 in the particle-hole symmetric case, where there is no change in the number of bosons per site, in which only local equilibration takes place. For other values of the average number of bosons per site changes in crossing from the superfluid to the Mott insulator in which case global mass transport is required for equilibration (which is not captured within the equation of motion studied in [112]).


## 3.8.1. Two Particle Irreducible Formalism

The two particle irreducible (2PI) closed-time-path (CTP) formalism for the BHM has been developed in detail by Rey et al. [162] and Temme and Gasenzer [165] who studied the 1/𝒩 expansion.The 2PI formalism can be obtained by starting from the action for the complex fieldsΦ and Φ* in the Bose Hubbard model in the form
(82)S[Φ]=∫‍dt∑iΦia(t)τabyℏ∂tΦib(t)+∫‍dt{U4𝒩(JΦi+1a(t)τabxΦ1b(t))-U4𝒩[τabxΦia(t)Φib(t)]2},
where a compact notation has been introduced, with Φi1=Φi and Φi2=Φi*, and the matrices τx and τy have the same form as the corresponding Pauli matrices but are labelled so as not to be confused with spin. The quantities studied are ϕia(t)=〈Φia(t)〉 and Gijab(t,t′)=〈TCΦia(t)Φib(t′)〉-ϕia(t)ϕib(t′), where TC  indicates contour ordering. The theory is specified by the effective action Γ[ϕ,G], which takes the form
(83)Γ[ϕ,G]=S[ϕ]+i2Tr[lnG-1]+i2Tr[D-1(ϕ)G]+Γ2[ϕ,G],
where
(84)iDijab(t,t′)-1=δS[ϕ]δϕia(t)δϕjb(t′),
and Γ2 is the two-particle irreducible generating functional. Γ2 can be described as an infinite sum of diagrams from the expansion of interaction terms in S[ϕ+φ] of cubic and higher order in φia=Φia-ϕia [162].In practice, various approximation schemes keep some subset of these diagrams and can be classified based on which terms are kept. The simplest approximation is to assume thatΓ[ϕ,G]=S[ϕ], corresponding to the mean field approximation. The Bogoliubov (one-loop) approximation can be obtained if the effective action is taken to be (83) with Γ2 ignored. Approximations that take into account interactions more carefully include the time-dependent Hartree-Fock-Bogoliubov (HFB) approximation, in which Γ2 is truncated to keep only the first-order diagram in U. This approximation conserves energy and particle number but violates Goldstone’s theorem. The expansion can be viewed as one in Ut/J, where t is time, and hence will be most accurate for short times or small U/J. It is possible to keep higher-order terms in Γ2—second-order expansions were considered by Rey et al. [162] and Temme and Gasenzer [165], who both considered a 1/𝒩 expansion, which simplifies the second-order expansion of Γ2. From each of these various effective actions, saddle point equations of motion can be derived, which may be quite lengthy; the full expressions for the second order expansion are listed in Appendix B of [162]. These equations of motion may also be used to derive a quantum kinetic equation, as was done by Rey et al. in [163].An additional issue discussed in [162] is that of conservation laws: for a closed system, the evolution is unitary in time, so energy must be conserved. However, particular approximation schemes for Γ may not conserve energy. On the other hand, schemes which are ϕ derivable, that is, the self-energy may be written as the derivative of a functional of G [170] do conserve particle number, energy, and momentum, which is the satisfied for the cases considered in [162]. A comparison of the results obtained for the approximation schemes outlined in [162] are shown in Figure 10 and compared against the results of exact dynamical evolution in small systems. The condensate population as a function of time is shown for initial conditions in which the condensate starts all in a single well, for N=20, 40, and 80 particles with J=1/2 and NU/J=4. The Hartree-Fock Bogoliubov approximation is seen to be fairly poor at describing the time development of the condensate, with the 1/N approach giving the most accurate account of the time evolution for the condensate population per well as a function of time as N becomes larger.Figure 10
Condensate population per well in each well of a two-well system and total condensate population as a function of time (in units ofℏ/J) for parameters J=1/2, NU/J=4, and N=20,40,80 as calculated exactly, using second-order expansion, 1/N expansion, and the Hatree-Fock Bogoliubov approximation [162]. Adapted figure with permission from A. M. Rey, B. L. Hu, E. Calzetta, A. Roura, and C. W. Clark, Phys. Rev. A 69, 033610 (2004). Copyright 2004 by the American Physical Society.One can see that the range of couplings considered are very much in the weakly interacting limit—this suggests the need for other methods to tackle questions relating to the crossing of the superfluid Mott insulator phase transition.


## 3.8.2. Strong Coupling

Sengupta and Dupuis [119] developed a strong coupling nonperturbative approach to the BHM with features that make it particularly attractive for generalization to real time. In particular, their approach is exact in the limit of both small (U=0) and large (J=0) interactions and obtains the expected features in the excitation spectrum in both the superfluid and the Mott insulating phases. It also naturally allows for the calculation of nonlocal correlations. I first briefly recap the equilibrium calculation of Sengupta and Dupuis and then discuss generalizations of this approach using the Schwinger-Keldysh method [112, 120].Sengupta and Dupuis started with a path integral form for the BHM and then performed two Hubbard Stratonovich transformations to write the generating functional as(85)Z=Z0∫‍𝒟[ψ*,ψ,ϕ*,ϕ]e(ψ∣Jψ)-[(ψ∣ϕ)+c.c.]+W[ϕ*,ϕ],
where Z0 is the generating functional when J=0, with the notation (ϕ∣ψ)=∫0β‍dτ∑r‍ϕ*(r)ψ(r), where τ is imaginary time and β is inverse temperature. The term
(86)W[ϕ*,ϕ]=∑R=1∞(-1)R(R!)2∑a1⋯aRG{ai,bi}Rcϕa1*⋯ϕaR*ϕbR⋯ϕb1,
that enters in the action is evaluated with all the fields at the same site, and
(87)G{ai,bi}Rc=(-1)R〈ψa1⋯ψaRψbR*⋯ψb1*〉,
is the connected local 2R point Green’s function. By integrating out the ϕ field in the generating functional, Sengupta and Dupuis obtained an effective action that to quartic order in the field ψ takes the form
(88)S[ψ*,ψ]=-∑a,bψa*(Gab-1+Jab)ψb+14∑a1,a2,b1,b2Γa1,a2,b1,b2IIψa1*ψa2*ψb2ψb1,
where G-1 is the Green’s function in the local limit and ΓII is the two particle vertex in the local limit. Sengupta and Dupuis approximated ΓII by its static value and introduced g=(1/2)ΓII|static. This leads to the effective action
(89)S=-∫0β‍dτdτ′∑r,r′‍ψr*(τ)[G-1(r,τ,r′,τ′)+Jr,r′δ(τ-τ′)]ψr′(τ′)+g2∫0β‍dτ∑rψr*ψr*ψrψr,
which has the attractive feature that it leads to the correct Green’s function -〈ψr(τ)ψr′*(τ′)〉 in both the strong interaction (J=0) and noninteracting (U=0) limits. This feature makes this approach particularly appealing for the study of out-of-equilibrium dynamics, since it gives the hope that one can accurately represent the behaviour of the system over a wide range of parameters that cover both the superfluid and Mott insulating phases.By taking a saddle point approximation to the action, and then expanding to quadratic order in fluctuations, Sengupta and Dupuis obtained the following results for the excitations in the Mott and superfluid phases, respectively. In the Mott phase, they found that the Green’s function𝒢(k,iω)=-〈ψ(k,iω)ψ*(k,iω)〉 has the form
(90)𝒢(k,iω)=1-zkiω-Ek-+zkiω-Ek+,
where
(91)Ek±=-δμ+ϵk2±12[ϵk2+4ϵkUx+U2]1/2,
with ϵk the Fourier transform of the hopping -Jr,r′,  x=n0+1/2 and
(92)zk=Ek++δμ+UxEk+-Ek-,
with δμ=μ-U(n0-1/2) where n0 is the number of bosons per site for chemical potential μ when J=0. In the Mott insulating phase, the modes Ek± are gapped and the system becomes unstable to superfluidity when the energy of one or other of the modes vanishes for k=0. In the superfluid phase, the expressions for the spectrum are more complicated, and four modes are found ±Ek±, two of which are gapless as k→0 and have the form Ek-=c|k|, and the other two are gapped, as illustrated in Figure 11.Figure 11
Excitation energies±Ek± for n=1 and U/J=20 [119]. Adapted figure with permission from K. Sengupta and N. Dupuis, Phys. Rev. A 71, 033629 (2005). Copyright 2005 by the American Physical Society.Both  Grass et al.  [120] and Kennett and Dalidovich [112] used the Schwinger-Keldysh technique to construct theories for the real-time dynamics of the BHM.  Grass et al.  developed a real-time Ginzburg-Landau theory based on the field theoretic approach developed in [171, 172] and used this theory to calculate collective modes of the BHM, finding amplitude and sound modes in the superfluid and particle and hole excitations in the Mott insulator, as expected.  Grass et al.  also illustrated the power of this approach by expanding the Green’s functions in the equations of motion in the superfluid regime, they were able to obtain the equation of motion
(93)i∂Ψi∂t=-∑jJijΨj-μΨi+U2Ψi|Ψi|2,
which is the lattice version of the Gross-Pitaevski equation (GPE) [102].A similar approach developed by Kennett and Dalidovich [112] generalizes the theory of Sengupta and Dupuis to real time. We determined the effective action to quartic order for the BHM within the strong-coupling Schwinger-Keldysh approach and then used this to derive equations of motion to study the out of equilibrium dynamics.The approach followed in Kennett and Dalidovich [112] was to write a path integral for the BHM
(94)𝒵=∫‍[𝒟a*][𝒟a]eiSBHM[a*,a],
where aia is the field at site i on contour a, with a=1 or 2, and with the notation τi representing the ith Pauli matrix that acts in Keldysh rather than spin space. We performed a Keldysh rotation from the 1, 2 basis to the quantum (q), classical (c) basis [168, 173–175]:
(95)(a1(t)a2(t))→(a~q(t)a~c(t))=L^(a1(t)a2(t)),
where  aq and ac are the quantum and classical components of the field and
(96)L^=12(1-111).
The effect of this on the action is that τ3 in the 1, 2 basis becomes τ1 in the q, c basis. Similarly to Sengupta and Dupuis, we performed two Hubbard-Stratonovich transformations, the first of which decouples the hopping term and the second leads to the effective action to quartic order in terms of z fields (one can show that the Green’s functions for z are the same as those for the original field a [112, 119]) as
(97)Seff[z*,z]=∫‍dt∑ijzia*(t)(2Jij)τab1zjb(t)+∫‍dt1dt2∑izia1*(t1)[Gia2a1(t2,t1)]-1zia2(t2)+14∫‍dt1dt2dt3dt4∑iua1a2a3a4(t1,t2,t3,t4)×zia1*(t1)zia2*(t2)zia3(t3)zia4(t4),
where
(98)ua1a2a3a4(t1,t2,t3,t4)=14∫‍dt5dt6dt5′dt6′Gia5a6a6′a5′2c(t5,t6,t6′,t5′)×{[Gia5a1(t5,t1)]-1[Gia6a2(t6,t2)]-1×[Gia3a5′(t3,t5′)]-1[Gia4a6′(t4,t6′)]-1+((a4,t4)⟷(a3,t3))+((a6,t6)⟷(a5,t5))+((a6,t6)⟷(a5,t5);(a4,t4)⟷(a3,t3))[Gia4a6′(t4,t6′)]-1}.
The matrix Green’s functions G contains the retarded, advanced, and Keldysh Green’s functions. These Green’s functions, along with the two-particle connected Green’s function G2c are evaluated with respect to the single-site Hamiltonian H^0: full expressions are presented in [112]. The mean field phase boundary can be determined from the effective action equation (97) from the vanishing of the coefficient of zq*zc. The diagram corresponding to the interaction term u is illustrated schematically in Figure 12.Figure 12
Schematic diagram for the coefficient of the quartic term in the effective actionua1a2a3a4(t1,t2,t3,t4).The following symmetry relations hold for the interaction kernelu from the definition above:
(99)uabcd(t1,t2,t3,t4)=ubacd(t2,t1,t3,t4)=uabdc(t1,t2,t4,t3),
(this corrects [112]) and it may also be shown [112] that
(100)Gia5a6a6′a5′2c(t5,t6,t6′,t5′)=Gia6a5a6′a5′2c(t6,t5,t6′,t5′)=Gia5a6a5′a6′2c(t5,t6,t5′,t6′).
Grass et al. [120, 168] also noted similar symmetry relations for this correlation function under interchange of Keldysh indices. These symmetries imply that G2c has only 8 independent components: Gqqqq2c,Gcqqq2c,Gqqqc2c,Gqqcc2c,Gccqq2c,Gcqcq2c,Gqccc2c, and Gcccq2c. The remaining four-point function Gcccc2c=0 by causality [161]. Explicit expressions for each of the nontrivial components were written down in [112].A simplified low frequency equation of motion for the superfluid order parameterz can be derived from the action that applies away from the degeneracy points of the Mott lobes (in which case only the component Gcqqq2c is required). Equations of motion that allow for high frequencies and spatial variations of the order parameter and its correlations can also be derived. Kennett and Dalidovich focused on k=0 and defined δ(t)=-2dj(t), where the hopping is taken to be J(t)=J0+j(t) with J0 chosen so that
(101)2dJ0+ν=0,
(i.e., J0 lies on the mean field phase boundary for a given μ). Hence, the approximate equation of motion is
(102)κ2∂2z∂t2+iλ∂z∂t+δ(t)z+u|z|2z=0,
where the coefficients κ2, λ, and u can all be calculated in terms of U and μ [112]. We studied (102) for fixed μ and time-varying J.For fixedμ there are two possibilities for the dynamics: (a) the particle-hole symmetric case (the transition at the tip of the Mott lobe), in which case λ=0, and (b) the generic case, in which λ≠0. The function δ(t) controls the traversal of the quantum critical point. We chose
(103)limt→-∞δ(t)=-δ0;limt→∞δ(t)=δ1.
In our numerics, we used
(104)δ(t)=(δ0+δ12)tanh(tτQ)+δ1-δ02,
where we introduced a timescale τQ which is the characteristic time for δ(t) to cross from -δ0 to δ1 [111]. In the vicinity of the transition t=0 is linear in t. Functional forms of δ(t) which are not linear in t in the vicinity of t=0 may lead to differing behaviour in the long-time limit [79].We solved (102) numerically with λ=0 for a variety of values of τQ, and found that for t≫τQ the form of the solution is that z(t) oscillates in a periodic manner with a magnitude that decreases with increasing τQ. When averaged over a period T at times t≫τQ,
(105)〈z〉T=1T∫tt+T‍dt~z(t~)=0,
as would be expected in the Mott insulating state. Defining zmax(τQ)=limt→∞|z(t)| we found that zmax(τQ) decreases with increasing τQ without any indication of saturation.In the generic case in whichλ≠0, we solved (102) numerically for a variety of values of τQ as is displayed with time scaled by τQ in Figure 13 for μ/U=0.25 (well away from both degeneracy and the particle hole symmetric case). Both the real and imaginary parts of z(t) are oscillatory functions, each of which has a vanishing average over a period for times greater than zero. The nonzero amplitude of z at long times indicates that there is memory of the nonzero value of ρ in the superfluid phase and that the final state is a metastable state rather than the Mott insulator in which ρ=0. The collapse of traces with different τQ after scaling the time with τQ indicates that |z| has a finite value even when τQ becomes large.Figure 13
Dynamics of|z(t)| normalized to unity in the generic case, for a variety of τQ with time rescaled by τQ for a quench from 2dJ/U=2.0 to 2dJ/U=0.0. Figure adapted from [112].To determine the fate of the system at long quench times, we definedzmax=limτQ→∞zmax(τQ) and calculated this as a function of chemical potential. Apart from the particle-hole symmetric point, it appears that there is a transition to a metastable state in which zmax≠0. Whilst it is true that the spatial dependence of z and higher frequency components were ignored in (102), it seems reasonable that the metastable state identified at the mean field level may persist even when these are taken into account.The calculations discussed above are for fixed chemical potential rather than fixed particle number, as is found in a cold atom system. However, within the LDA fixed chemical potential would correspond to studying dynamics at fixed radius. The results found for the long-time limit ofzmax are in accordance with this idea. There appears to be equilibration as zmax→0 in the particle-hole symmetric case, where there is no change in the number of bosons per site, in which only local equilibration takes place. For other values of the average number of bosons per site changes in crossing from the superfluid to the Mott insulator in which case global mass transport is required for equilibration (which is not captured within the equation of motion studied in [112]).


## 3.9. Exact Results

### 3.9.1. Spectral Moments

Perturbative expansions inJ/U [41–43] have allowed for very accurate determination of the Mott insulator phase boundary and the momentum distribution in the Mott insulator phase [176]. Very recently Freericks et al. [177] calculated exact sum rules for the spectral functions of the retarded Green’s function GR and self-energy ΣR allowing for systems that are out of equilibrium and spatially inhomogeneous. The retarded Green’s function is defined as
(106)GijR(t1,t2)=GijR(T,t)=-iθ(t)〈[a^i(T+t2),a^j†(T-t2)]〉,
where operators are in the Heisenberg picture, T=(t1+t2)/2, t=t1-t2, and the average is taken with respect to the state operator
(107)ρ=e-βH(t→-∞)Z,
where Z=Tr[e-βH(t→-∞)]. The spectral moments of the retarded Green’s function are
(108)μnijR=-Im[in∂n∂tnGijR(T,t)]t=0+,
(note that terms involving derivatives of θ(t) do not contribute to these moments). Freericks et al. [177] calculated the general expressions for these moments and obtained specific analytic results in the strong coupling regime when U≫J. They also compared the moments as calculated within several different methods numerically: VCA, QMC, DMRG, and strong-coupling perturbation theory for the 1 dimensional BHM in the Mott insulating state with n=1. They found that VCA was the best performing method and that both QMC and DMRG had errors on the order of a few percent for the spectral moments. They argued that the reason for the differences in performance is a result of QMC and DMRG having a cutoff that means that they focus on low energy states, whereas VCA, a method that works best at strong coupling, retains information about higher energy states that can be important for higher order spectral moments. These sum rules are an important step towards benchmarking analytical and numerical methods, both in and out of equilibrium. They may prove particularly important for testing methods that try to span the quantum critical point.


### 3.9.2. Quench fromJ=0 to U=0

One of the few examples of an exact result in this field was provided by Cramer et al. [104], who showed that for a quench from J=0 to U=0 with an initial Mott insulator ground state in 1 dimension, the final state, whilst not the superfluid state, is one that locally maximizes the entropy, in the sense that the entanglement between a particular site and the rest of the chain becomes maximal. In the infinite system size limit, whilst the system is locally relaxed, it never fully relaxes to an equilibrium state, and for large but finite lattices with N sites, there will be recurrences on times of order N. The physical picture behind this result is that for a lattice site (or sites) within the “lightcone” defined by excitations in the system, there is local mixing of the state, whereas contributions from outside the “lightcone” are exponentially suppressed similar to the Lieb-Robinson bound for spin systems [137]. In quenches from the Mott insulator to the superfluid phase, the idea of a “lightcone” reoccurs in both theoretical and experimental results. Cramer et al. argued that their results are quite general and also apply to initial states of the form ∏i=1N‍|mi〉, where the occupations on different sites mi may differ. They also argued that the lightcone physics underlying the arguments in one dimension should apply in higher dimensions, so that there will be local relaxation inside the causal cone.


## 3.9.1. Spectral Moments

Perturbative expansions inJ/U [41–43] have allowed for very accurate determination of the Mott insulator phase boundary and the momentum distribution in the Mott insulator phase [176]. Very recently Freericks et al. [177] calculated exact sum rules for the spectral functions of the retarded Green’s function GR and self-energy ΣR allowing for systems that are out of equilibrium and spatially inhomogeneous. The retarded Green’s function is defined as
(106)GijR(t1,t2)=GijR(T,t)=-iθ(t)〈[a^i(T+t2),a^j†(T-t2)]〉,
where operators are in the Heisenberg picture, T=(t1+t2)/2, t=t1-t2, and the average is taken with respect to the state operator
(107)ρ=e-βH(t→-∞)Z,
where Z=Tr[e-βH(t→-∞)]. The spectral moments of the retarded Green’s function are
(108)μnijR=-Im[in∂n∂tnGijR(T,t)]t=0+,
(note that terms involving derivatives of θ(t) do not contribute to these moments). Freericks et al. [177] calculated the general expressions for these moments and obtained specific analytic results in the strong coupling regime when U≫J. They also compared the moments as calculated within several different methods numerically: VCA, QMC, DMRG, and strong-coupling perturbation theory for the 1 dimensional BHM in the Mott insulating state with n=1. They found that VCA was the best performing method and that both QMC and DMRG had errors on the order of a few percent for the spectral moments. They argued that the reason for the differences in performance is a result of QMC and DMRG having a cutoff that means that they focus on low energy states, whereas VCA, a method that works best at strong coupling, retains information about higher energy states that can be important for higher order spectral moments. These sum rules are an important step towards benchmarking analytical and numerical methods, both in and out of equilibrium. They may prove particularly important for testing methods that try to span the quantum critical point.


## 3.9.2. Quench fromJ=0 to U=0

One of the few examples of an exact result in this field was provided by Cramer et al. [104], who showed that for a quench from J=0 to U=0 with an initial Mott insulator ground state in 1 dimension, the final state, whilst not the superfluid state, is one that locally maximizes the entropy, in the sense that the entanglement between a particular site and the rest of the chain becomes maximal. In the infinite system size limit, whilst the system is locally relaxed, it never fully relaxes to an equilibrium state, and for large but finite lattices with N sites, there will be recurrences on times of order N. The physical picture behind this result is that for a lattice site (or sites) within the “lightcone” defined by excitations in the system, there is local mixing of the state, whereas contributions from outside the “lightcone” are exponentially suppressed similar to the Lieb-Robinson bound for spin systems [137]. In quenches from the Mott insulator to the superfluid phase, the idea of a “lightcone” reoccurs in both theoretical and experimental results. Cramer et al. argued that their results are quite general and also apply to initial states of the form ∏i=1N‍|mi〉, where the occupations on different sites mi may differ. They also argued that the lightcone physics underlying the arguments in one dimension should apply in higher dimensions, so that there will be local relaxation inside the causal cone.


## 3.10. Summary

Through the variety of techniques that have been applied to quantum quenches in the BHM and the particular calculations discussed in previous sections, several themes become clear. Both analytic and numerical techniques suggest that when there is a quench from superfluid to the Mott insulator, there can be equilibration if the final value ofU/J is not too large [97], but for larger values of U/J a frozen nonequilibrium state results [87, 97, 108, 112], and in the presence of a trap, the Mott insulating regions can impede relaxation through mass transport [94, 99]. These results leave a number of questions that can be pursued both theoretically and experimentally. Kollath et al. [97] asked about the eventual relaxation after a deep quench from superfluid to Mott insulator and whether there may be aging-like phenomena. Natural quantities to look at in the context of out-of-equilibrium physics are two-time correlation functions, such as the two-time density density correlation function
(109)C(r,r′,t,t′)=〈n(r,t)n(r′,t′)〉-〈n(r,t)〉〈n(r′,t′)〉.
If there is a separation of fast and slow timescales for relaxation then one might expect a quantity such as
(110)C(t,t′)=∑rC(r,r,t,t′),
to display behaviour similar to that seen for aging in glassy systems, in which case it could be separated into a fast time-translation invariant piece, and a slower piece that depends on the two times, t and t′:
(111)C(t,t′)=Cfast(t-t′)+Cslow(t,t′),
with Cfast decaying to zero with a timescale that is much shorter than the typical timescale for decay of Cslow, so that Cslow is essentially constant whilst Cfast decays. If aging is present in the Mott insulators for large U/J, then the simplest decay one might expect for Cslow is of the form (with t>t′)
(112)Cslow~f(h(t)h(t′)),
where h(t) is a monotonically increasing function of t, possibly h(t)∝t and f(x)~x-α for some α [178].In going from the Mott insulator to superfluid, there are a variety of results to reconcile. Several groups predict oscillations in the amplitude of the superfluid order parameter [91, 96], which are yet to be seen experimentally [138]. Scaling theory predictions of the residual energy and defect formation probability were also not found to match calculations by Trefzger and Sengupta [91] using a projection operator formalism. Navez and Schützhold [93] found a light-cone-like behaviour of the growth of the amplitude of the superfluid order parameter after a quench with diffusive growth of the phase coherence. Several groups found numerical and analytical evidence of a light-cone defined by a maximum velocity for the growth of density correlations both in quenches from superfluid to the Mott insulator and within the Mott insulating state [98, 99, 152] in 1 dimension. Regarding the quenches from the Mott insulator to superfluid, it would be desirable to see whether signatures of the KZM survive as expected in other methods than the projection operator approach, and whether the results of Navez and Schützhold can be found with other methods as well. It would also be very interesting to see whether the Lieb-Robinson-like bound of a maximal velocity found in the BHM in 1 dimension survives to higher dimensions than one—the results of the 1/Z expansion suggest this may be possible [93] although the velocity found in [93] had quite a different form as a function of J and U than found in one dimension.On the technical side, most of the methods considered are accurate when eitherU is small or J is small, and in some cases in both limits (e.g., the closed time path approaches based on strong coupling theory). Numerical methods such as QMC and DMRG have provided very accurate determinations of the phase diagram in equilibrium but are not easily generalizable either to out of equilibrium situations (QMC) or dimensions higher than 1 (DMRG and TEBD). Other methods such as the Gutzwiller mean field are reasonably straightforward to calculate with and can provide qualitatively correct behaviour for single-particle properties but in the simplest formulation do not contain information about nonlocal correlations. Formulations such as the generalization of Sengupta and Dupuis’ nonperturbative approach to real time in the context of closed time path methods seem to be an avenue that is worth exploring further for studying dimensions higher than 1. For all methods, however, the out-of-equilibrium sum rules introduced by Freericks et al. [177] are a particularly welcome development that can help to benchmark methods and provide consistency checks between methods.


## 4. Experiment

The traversal of the quantum phase transition between the Mott insulating and superfluid phases achieved by Greiner et al. [6] provided a dramatic experimental demonstration of the realization of the BHM in a three-dimensional optical lattice. This followed earlier work by Orzel et al. [139] showing number squeezing in a two-well system and was quickly followed by a number of other realizations of the BHM and quantum phase transitions of bosons in optical lattices [57, 58, 179–185]. These results and more recent experiments on out-of-equilibrium dynamics in the BHM [138, 146, 186–188] are discussed in this section.Greiner et al. [6] were able to access the Mott insulator to superfluid transition by changing the strength of their optical lattice V0, which has the effect of changing the ratio J/U in the BHM as discussed in Section 2.3. As indicated in (14), the hopping, J, is much more sensitive to changes in V0 than the interaction strength U; hence, increasing the depth of the lattice increases the strength of interactions and makes the Mott phase more accessible. Greiner et al. obtained time of flight images of the atoms released from the trap, which showed a central peak and satellite peaks for weak lattice strengths, related to the quasimomentum distribution in the superfluid phase. As the depth of the optical lattice was increased, these peaks gained in strength up to about V0≃13ER and then decreased in strength with increasing V0, and no interference pattern was visible for V0≃20ER.To demonstrate that this change in the time of flight pattern was due to the formation of a Mott insulator rather than decoherence, Greiner et al. performed two experiments in which the depth of the lattice was lowered toV0=22ER, the system held for 20 ms, and then V0 was changed to 9ER (a value expected to correspond to the superfluid regime). The peaks in the time of flight images became sharp within 2 ms (on the order of the tunnelling time τtunnel≃ℏ/J between two neighbouring lattice sites) when this protocol was followed. On the other hand, when phase coherence was deliberately destroyed using a magnetic field gradient in the superfluid phase, no phase coherence was seen following the same protocol in V0 even 400 ms after the ramp down as illustrated in Figure 14. The results of extensive experiments that reveal the boundary of the Mott insulating phase as a function of lattice depth and characteristic density for a two-dimensional optical lattice [57] are shown in Figure 15.Figure 14
(a) Experimental sequence of lattice depth used in [6]. After a time t, the interference pattern was measured by releasing the atoms from the trapping potential. (b) The width of the central interference peak in time of flight images for different ramp down times t. The open circles correspond to a phase incoherent state and the closed circles correspond to a state reached following the protocol shown in (a). The solid line is a fit with a double exponential decay with time constants τ1=0.94(7) ms and τ2=10(5) ms. Also shown are absorption images of the interference patterns after ramp-down times of (c) t=0.1ms, (d) 4 ms, and (e) 14 ms. Adapted with permission from Macmillan Publishers Ltd: Nature, M. Greiner, O. Mandel, T. Esslinger, T. W. Hansch, and I. Bloch, Nature 415, 39 (2002), copyright 2002.Figure 15
Experimentally determined phase diagram of bosons in a two-dimensional optical lattice [57] as a function of lattice depth (equivalently interaction strength) and characteristic density. Each circle on the phase diagram corresponds to an individual measurement. The blue region was predicted to support a Mott insulating phase in QMC calculations [192]. The red ovals mark the experimentally determined transition boundary to be compared with the blue line of theory. The size of the ovals represent the respective uncertainties in the phase boundary in density and chemical potential. The yellow dashed line is a fit to the experimentally determined phase boundary, and the green dashed line shows a set of data at fixed atom number N2D≃3500. Adapted figure with permission of I. B. Spielman from K. Jiménez-Garcia, R. L. Compton, Y.-J. Lin, W. D. Phillips, J. V. Porto, and I. B. Spielman, Phys. Rev. Lett. 105, 110401 (2010). Copyright 2010 by the American Physical Society.The experimental protocols used for subsequent investigations of the transition from superfluid to Mott insulator or vice versa have mainly followed similar approaches, in which the depth of the optical lattice has been varied as a function of time to change the strength of interactions. Another possibility is the use of the Feshbach resonances to increase the strength of interactionsU, via changing the s-wave scattering length as, recently realized experimentally by Mark et al. [189]. In changing parameters in the Hamiltonian with time, if the goal is to study the BHM, then the changes must not be so fast as to cause significant excitation of bosons out of the lowest Bloch band. For a ramp time tramp, this implies that the characteristic frequency f≃1/tramp for the ramp should be such that f≪ν, where ν=4ERV0/2π, the energy for the lowest excitation to the first excited Bloch band [101, 190, 191].In this section, I focus on experiments investigating out-of-equilibrium behaviour in the Bose Hubbard model. In Section4.1, I discuss quantum phase revivals, and in Sections 4.2 and 4.3, I discuss quantum quenches from superfluid to Mott insulating states and from Mott insulating to superfluid states, respectively. In Section 4.4, I review recent experiments in which a quench within the Mott phase was studied. In Section 4.5, I discuss issues that may arise in trying to characterize out-of-equilibrium states, both in imaging and thermometry.


### 4.1. Quantum Phase Revivals

A related experiment to the crossing of the superfluid Mott insulator transition demonstrated in [6] is that of phase revivals [183, 193]. In such an experiment, an initially superfluid state is quickly quenched into the Mott phase by increasing the lattice depth nonadiabatically, so that tunnelling between sites is strongly suppressed. The state is held in the Mott phase for a time t, and then the trapping potential and the optical lattice potential are turned off and a time of flight image is taken. Greiner et al. [183] monitored the collapse and revival of coherent peaks in the interference pattern as illustrated in Figure 16 to determine the ratio of coherent atoms as a function of t.Figure 16
Ratio of coherent atoms as a function of timet after a quench from the superfluid phase to the Mott insulating phase [183]. The initial lattice depth was VA=8ER and the final lattice depth VB=35ER. The quench from VA to VB was made in 50 μs. Adapted with permission from Macmillan Publishers Ltd: Nature, M. Greiner, O. Mandel, T. Esslinger, T. W. Hänsch, and I. Bloch, Nature 419, 51 (2002), copyright 2002.The simplest picture to explain this was outlined in [183]. In the infinite system size limit, one can view an initial state which is a noninteracting Bose Einstein Condensate (BEC) as a product of local coherent states [194] for N bosons on M sites:
(113)|φ0〉=∏i=1M‍e-|αi|2/2eαia^i†|vac〉=∏i=1M‍e-|αi|2/2∑ni‍αinini!|ni〉,
where the coherent state on site i is parametrized by αi,  |ni〉 is a state with ni bosons on site i, and |vac〉 is the vacuum state. After a quench J/U→0, the time evolution arises solely from H^0 and so the time evolution operator factorizes and for a Fock state with occupation ni on site i, takes the form [97, 183] (with ℏ=1)
(114)U(t)=∏i=1M‍e(i/2)Uni(ni-1)t.
This leads to the time dependence
(115)|φ(t)〉=∏i=1M‍e-|αi|2/2∑niαinini!e(i/2)Uni(ni-1)t|ni〉,
and since ni(ni-1)/2 is an integer, the state will have a periodic time dependence, with a  revival time of trev=2π/U. Hence, the coherent fraction
(116)NcohN=1N∑i|ψi(t)|2,
where ψi(t)=〈φ(t)|a^i|φ(t)〉 will also have a periodic time dependence with a revival time trev.However, it was pointed out by Schachenmayer et al. [194] that this picture is too naive. This is because the state φ0 implies a coherent superposition of states with different numbers of particles. This is to be contrasted with the situation in a quantum revival experiment in which there are a fixed number of particles N. The authors of [194] used a number conserving initial state
(117)|BECN〉=1N!NN(∑i=1M‍βia^i†)N|vac〉,
which may be shown to have the form |BECN〉=PN|φ0〉, where PN projects states onto a subspace of the Hilbert space with fixed total particle number N. Schachenmayer et al. calculated
(118)nq(t)=1M∑i≠jM‍〈a^i†a^j〉teiq(xi-xj),
for both a homogeneous and a trapped system, making use of the projected state with fixed N. In a homogeneous system, they found that their results converged to the naive calculation for M≃5, but for bosons in a trap, they found that for N=M, the relative error in the naive calculation scales as ~1/N, corresponding to about a 1% error for a system with 50 sites. This is quantitatively small, but nevertheless, conceptually important.The simplest analysis of quantum phase revivals ignores tunnelling during the quench, or nonzero interactions before the quench that modify the initial state from a purely coherent state. Schachenmayer et al. considered the effect of this physics using TEBD and found that it could lead to damping of revivals with time, which can be seen in Figure16. Previous investigations of the effect of tunnelling on quantum phase revivals [108, 115] found similar results.The time dependence of quantum phase revivals has also been used to show evidence of multiparticle interactions that arise from virtual transitions of bosons from the lowest Bloch band to higher Bloch bands, up to six-body interactions [188]. Tiesinga and Johnson [106] showed (using a mean field approximation) that these multiband effects in quantum phase revivals can be used to infer the number squeezing in the initial superfluid state.


### 4.2. Quench from Superfluid to Mott Insulator

Two recent experiments [146, 186] explored relaxation after a quench from a superfluid to a Mott insulator, yet found rather different results. Bakr et al. [186] found fast relaxation, whereas Hung et al. [146] found slow relaxation, relative to natural timescales in the dynamics. I review these experiments and discuss the time-dependent Gutzwiller mean field calculations by Natu et al. [94] that suggest a picture of two timescales for equilibration in a trap.


#### 4.2.1. Fast Dynamics

In an important step in the imaging of cold atom systems, Bakr et al. [186] demonstrated single atom-single lattice site imaging that allows for site-resolved observations of dynamics in both space and time. They applied this imaging to follow the out-of-equilibrium dynamics after a quench from a superfluid to a Mott insulator. In an experiment using a few thousand Rb87 atoms in a two-dimensional optical lattice, they used the following protocol. After loading atoms into the trap, they adiabatically increased the depth of their trap to V0=11ER (corresponding to a superfluid state). They then ramped up the depth of their trap to V0=16ER, in a time tramp that ranged between 0.2 and 20 ms. They observed the probability that there were an odd or even number of atoms on each site (corresponding to n-=1 or n-=2, respectively, as a function of tramp as illustrated in Figure 17). The time constant for n=1 was found to be 3.5±0.5 ms and for n=2 it was found to be 3.9±1.3ms. These timescales were noted to be short compared to the timescale h/Jc~68ms (where Jc is the critical value of the hopping) and comparable to the timescale h/Uc~4.1ms.Figure 17
The probability of having an odd number of atoms on each site as a function of ramp time averaged over 19 data sets. The squares correspond to sites withn=1 and the circles to sites with n=2 and the curves are exponential fits to the data. The inset shows the ramp used in the experiment, with a fixed adiabatic exponential ramp to lattice depth V0=11ER followed by a linear ramp with variable duration  [186]. From W. S. Bakr, A. Peng, M. E. Tai, R. Ma, J. Simon, J. I. Gillen, S. Fölling, L. Pollet, and M. Greiner, Science 329, 547 (2010). Adapted with permission from AAAS.


#### 4.2.2.  Slow Dynamics

An example of very slow dynamics after a quench from superfluid to the Mott insulator was realised in the experiment of Hung et al. [146]. Studying the dynamics of trapped Cs133 atoms using in situ imaging, they prepared a condensate which was then loaded into a pancake trap, and then a two dimensional optical lattice was slowly turned on over 30 ms, to reach a lattice depth of V= 0.4ER. The lattice was then ramped to a final depth of Vf up to about 13ER (corresponding to U/J=41, implying a Mott insulator in a homogeneous system) using the protocol [101]
(119)V(t)=Vf(1+γ)1+γe4(t-τc)2/τc2,
so that V=Vf at time t=τc (chosen to be 20 ms for the data shown in Figure 18), and γ is chosen so that V(0)=0.4ER. At values of Vf higher than ~13ER, the timescale of the dynamics started to reach that of three-body recombination, limiting the resolution of mass transport. They characterized the relaxation to equilibrium by investigating the mean square deviation of the density after hold time thold from its long time (thold of order 500–800 ms) limit as
(120)Δ(thold)={∑i[n-i(thold)-n-longtimes]2}.
They found that even for timescales longer than the microscopic timescale t~h/J~ 15–75ms, there was still decay of Δ, which could be fitted with a single exponential decay. These results are illustrated in Figure 18.The evolution of the density profile following a 20 ms lattice ramp to (a)Vf=10ER and (b) Vf=12ER. The density profiles at three different hold times are illustrated: 0 (squares), 200 ms (circles), and 500 ms (triangles). The insets show the time evolution of Δ for each value of Vf, each fitted to a single exponential. The fits of the decay times against Vf are shown in (c), along with the appropriate phase of the system [146]. Adapted figure with permission from C.-L. Hung, X. Zhang, N. Gemelke, and C. Chin, Phys. Rev. Lett. 104, 160403 (2010). Copyright 2010 by the American Physical Society.
(a)
(b)
(c)In the deep trap limit, kinetic energy can be essentially ignored, and the average number of particles per site can be calculated fromn-=∑nnPn, where Pn=e-β(H-μn)/Ξ, with H≃(U/2)n(n-1) and the grand canonical partition function Ξ=∑ne-β(H-μn). Hung et al. induced three-body recombination losses after different hold times. By comparing these losses with the expectation for thermally occupied sites, for a hold time of thold~800ms they calculated an effective temperature for the insulating centre of the cloud of Teff~6 nK, considerably lower than the Teff~20 nK that they found in the superfluid wings of the trap. These results suggest that although there may have been local equilibration, even on the large timescales considered, global equilibrium was not established.


#### 4.2.3. Two Time Scales

At first glance the results obtained by Bakr et al. [186] and Hung et al. [146] do not seem to be particularly compatible in that in one case relaxation times much less than h/J were observed and in the other relaxation times larger than h/J were observed. However, calculations by Natu et al. [94] illustrate how they can be reconciled within a picture of two timescales for relaxation, specifically, very quick local equilibration and then slower mass transport. They noted that the initial state before the quench can strongly influence equilibration depending on whether it leads to impeded mass transport or not.Natu et al. [94] used the time-dependent Gutzwiller approach to study the Bose Hubbard model with an external trapping potential. By making use of the explicit lattice depth dependence of J and U (14), they tracked the effects of changing the depth of the lattice with time in the manner
(121)V(t)=Vi+(Vf-Vi)tτr,
where Vi,f are the initial and final values of the lattice potential and τr is the ramp time. By comparing a homogeneous system with 1 atom per site with different initial configurations of bosons in a trap as in experiments in [146, 186, 195], Natu et al. argued that local equilibration will be fast, but the speed of global relaxation is strongly affected by whether there needs to be mass transport from one region of the trap to another, particularly if there is a Mott insulating phase, which can act as a barrier to particle diffusion (this is the “Mott barrier” idea discussed by Bernier et al. [99]).Evidence of slow relaxation and local, but not global equilibration is illustrated in Figure19 for an initially superfluid configuration. The figure shows the initial and final density distribution, μ(n(r))+Vtrap, and the coherences Ci(t). In equilibrium, μ(n(r))+Vtrap=μ, but this is not the case in Figure 19, where it can be seen to have local plateaux, but is not globally flat. The plot of coherences Ci=-〈ai〉∑j〈aj*〉 illustrates that for parameters appropriate to Rb87, relaxation timescales can be on the order of 200 ms, much longer than h/Ui~3.3ms, as seen by Hung et al.Slow dynamics after a quench in a trap. (a) Density of particles in the trap as a function of radius—the initial superfluid state is indicated with a solid line, and the final density profile with a dashed line, which differs strongly from the equilibrium configuration (dots). (b)μ(n(r))+Vtrap as a function of radius after the quench—for large values of r there are no particles, and the points represent V(r). (c) Density plot of the time-dependent coherences shows the growing Mott insulating regions that separate the brighter superfluid regions [94]. Adapted figure with permission from S. S. Natu, K. R. A. Hazzard, and E. J. Mueller, Phys. Rev. Lett. 106, 125301 (2011). Copyright 2011 by the American Physical Society.
(a)
(b)
(c)Further confirming the two-timescales picture, Natu et al. investigated other initial conditions in which particle transport across a Mott insulating region was not needed for global equilibration and found that this can lead to fast global equilibration as seen by Bakr et al.


### 4.3. Quench from Mott Insulator to Superfluid

A quench from a Mott insulator to a superfluid phase may lead to a proliferation of defects through the Kibble-Zurek mechanism. An experiment that investigated such issues was carried out by Chen et al. [138]. The authors investigated Rb87 in a three-dimensional optical lattice in which the Mott insulating layers emerge for lattice depths greater than V0=13ER. They considered quenches at rates 1/τQ=(d/dt)(J/U) which were sufficiently slow so as not to excite atoms into higher vibrational states: 1/τQ~ (0.001–0.2)Uc/h. The quenches considered were from variable initial lattice depth to a final lattice depth of ~4ER corresponding to J~U.In order to obtain a measure of the number of excitations produced during the quench, Chen et al. obtained time of flight (TOF) images after a relatively long expansion of about 50 ms, to enhance the visibility of vortices and other excitations through large density fluctuations. They fitted the TOF images with a smooth profilef, which was the combination of a Thomas-Fermi profile with a Gaussian and then measured the amount of excitation χ~2 as the deviation from the smooth profile via
(122)χ~2=∑ijχ~ij2=α∑ij((𝒪ij-fij)2/fij)∑ij𝒪ij,
where i and j are pixel indices, 𝒪 is the measured optical depth, and α  is  a  constant that they determined via numerical simulation. They found in each case that a condensate was present, with a condensate fraction that varied between 0.35 and 0.6. Based on numerical simulations of the three dimensional time-dependent Gross-Pitaevskii equation, Chen et al. argued that in the noninteracting limit, χ~2 accurately represents the fraction of the Bogoliubov excitations for a trapped condensate. They also pointed out that in the non-interacting limit, χ~2 should be proportional to the number of excited atoms.In the Mott insulating phase, there are large fluctuations in the relative phases between different sites, and in the Kibble-Zurek mechanism, these are frozen in after a quench is started due to the diverging relaxation time near the critical point. After the system reaches the superfluid phase, these excitations can start to relax. The relaxation of these excitations as measured byχ~2-χ~02 is illustrated in Figure 20 (note that χ~02≃0.06 was determined by an average over all images without a Mott insulating shell). The authors also considered kinetic energy as a function of quench rate to obtain an estimate of the quench-induced heating.Dependence of excitations on quench rate1/τQ. (a), χ~2 with an offset subtracted is shown against quench rate, and (b), the kinetic energy (KE) per particle (with the quench KE0 subtracted) is plotted against quench rate. The experimental protocol for the lattice depth is shown in the inset [138]. Adapted figure with permission from D. Chen, M. White, C. Borries, and B. DeMarco, Phys. Rev. Lett. 106, 235304 (2011). Copyright 2011 by the American Physical Society.
(a)
(b)In the Kibble-Zurek mechanism, the quench rate1/τQ controls defect formation and the expectation is that the number of excitations in 3-dimension scales as [66]
(123)nex~1τQ3ν/(νz+1).
Figure 20 illustrates that (χ~2-χ~02)~1/τQr, with r=0.31±0.03. The experimental value differs substantially for the expectations for either a generic phase transition nex~1/τQ3/4 or for the transition in the XY universality class, for which nex~1/τQ. There are a number of possibilities as to why the predictions of scaling theory do not match with experiments. First, the finite length scales for the phases in the trap may cut off the quantum critical scaling, and thermal effects may also play a role in changing the density of excitations. Additionally, Chen et al. pointed out that the nature of the excitation may affect the scaling, for instance, for vortex excitations nex~1/τQ1/2. Further theoretical investigation of the Mott insulator to superfluid quench in a trap is needed in order to understand these experimental results for the number of excitations.


### 4.4. Quench within a Mott Insulator

Light-cone-like spreading of correlations was recently observed in experiments by Cheneau et al. [196] on the one-dimensional BHM. A two-dimensional gas of Rb87 was loaded into 10 decoupled 1-dimensional chains, with 10–18 particles in each chain, with V0=20Er(U/J=40), corresponding to an n-=1 Mott insulator in each chain. From this initial state, the lattice was reduced in depth in a quench to give values of U/J of roughly 5.0, 7.0 and 9.0. Cheneau et al. measured the two-point parity correlation function [124]
(124)Cd(t)=〈s^j(t)s^j+d(t)〉-〈s^j(t)〉〈s^j+d(t)〉,
where
(125)sj(t)=eiπ[n^j(t)-n-],
and d is the distance between sites. Note that sj gives +1 for no quasiparticles on a site and −1 for the presence of a quasiparticle. The results obtained from measuring Cd(t) and averaging over more than 1000 chains are shown in Figure 21. The correlations showed a front-like spreading, with a constant velocity for d = 2–6, in good agreement with t-DMRG calculations and an analytical model developed by Barmettler et al. [152] discussed in Section 3.7.Figure 21
Time dependence of the parity correlation functionCd(t) measured for bosons in one dimensional chains prepared with U/J=40 and then quenched to U/J=9. The circles are experimental data averaged over more than 1000 chains, the green line is the result of t-DMRG calculations for a homogeneous system, and the dashed-line is the result of calculations with the analytical model developed in [152]. The inset shows the experimental data as a colour map [196]. Adapted with permission from Macmillan Publishers Ltd: Nature, M. Cheneau, P. Barmettler, D. Poletti, M. Endres, P. Schauß, T. Fukuhara, C. Gross, I. Bloch, C. Kollath, and S. Kuhr, Nature 481, 484 (2012), copyright 2012.


#### 4.4.1. Superlattice Initial State

Another example of a quench within the Mott state was studied by Trotzky et al. [148] for small one-dimensional systems. They prepared a collection of isolated 1d Bose Hubbard chains, with an average of N=31 particles, each with an initial state |ψ0〉=|…,1,0,1,0,1,…〉 in which all atoms were on “even” sites, prepared using a superlattice potential. They then removed the superlattice in 200 μs and followed the subsequent ensemble average of the filling of bosons on odd sites as a function of time, nodd, which relaxed from 0 to ~0.5 on a timescale of a few ms. Ensemble averaged DMRG calculations were fairly successful in fitting the time dependence of the nodd(t) but gave slightly more accurate results for U/J≤6 when a next-nearest neighbour hopping term was allowed in the BHM, as calculated from the single-particle band structure. The amplitude of the oscillations in nodd was observed to decay as t-α with α=0.86±0.07 over a range of intermediate U/J≲10, which is considerably faster than the value of α=0.5 expected in the limits of no interactions or hardcore bosons (infinite interactions) [148]. This behaviour is captured well by time-dependent DMRG, but is currently a challenge for analytical models. For U/J≥10, the agreement between experiment and theory was less strong—Trotzky et al. suggested this might be due to non-adiabatic heating or interchain tunnelling. This work presented a nice comparison of theoretical and experimental work but highlighted a limitation of t-DMRG, in that the experiment ran for longer times than the calculations could be performed for—the limiting factor being the growth of entanglement entropy [145].


### 4.5. Interpreting Experiments

In order to interpret and understand experiments on cold atoms in Bose Hubbard systems both imaging and thermometry are key. I discuss both of these in the context of studying out-of-equilibrium dynamics.


#### 4.5.1. Imaging

The imaging techniques for cold atoms in optical lattices can roughly be classified as those that measure properties of the system in momentum space, and are most suited to the superfluid phase, and those that measure properties in real space and are more suited to the Mott insulating phase. I will focus first on momentum space.The experiments by Greiner et al. [6] demonstrated the traversal of the quantum critical region separating superfluid and Mott insulator by measuring the momentum distribution of atoms in the trap, taking advantage of the mapping of momentum onto position after atoms are released from the trapping potential. The resulting interference pattern in space takes the form [185]
(126)n(r)=(mℏt)3|w~(k=mrℏt)|2S(k=mrℏt)
and depends on S(k)=∑i,jeik·(ri-rj)〈a^i†a^j〉, multiplied by a Wannier envelope |w~|2. If there is long-range phase coherence, then the correlation function 〈a^i†a^j〉 has a correlation length of the order of the lattice, leading to a sharp peak at k=0 along with satellite peaks due to the presence of a lattice. The disappearance of these peaks was argued to indicate the disappearance of superfluidity and the entry into a Mott phase. However, peaks are also present in both the normal and Mott phases, albeit less sharp [38]. It has been suggested that a better quantity to measure is the visibility [38, 185], which, if measured so that the Wannier envelope cancels, takes the form [185]
(127)𝒱=nmax-nminnmax+nmin=Smax-SminSmax+Smin,
where nmin and nmax are the minimum and the maximum of the density distribution. At small U, the visibility tends to 1, indicating the presence of a superfluid and then drops rapidly for U/zJ larger than the critical value [179, 181, 185]. At lowest order in perturbation theory in J/U, it can be estimated that [185]
(128)𝒱≃43(n0+1)zJU
which is in accordance with calculations using QMC [197] and an analytic strong coupling approach [198] both of which suggest 𝒱~(U/zJ)-1. Gerbier et al. [185] found good experimental agreement with this form up to a lattice depth of V0=30ER corresponding to U/zJ≃200, where they suggested that adiabaticity in the preparation of the state may have broken down.The issue of adiabaticity also arises in imaging as was recently discussed in detail by Natu et al. [199]. In the band mapping protocol lattice strength is reduced over some timescale τ, during which the quasimomentum distribution is mapped to momentum, followed by a ballistic expansion in which the momentum distribution is mapped to position and then imaged as n(r). This can be viewed as an accurate representation of the quasimomentum distribution provided that τ is slow enough (τ≫h/Ebg, where Ebg is the band gap) that the quasimomentum states can adiabatically evolve into momentum states and fast enough (τ≪1/ν, where ν is the trap frequency) that occupation numbers of quasimomentum states do not change. These conditions can be met for a sufficiently weakly interacting Bose gas [199]. However, interactions introduce an additional timescale h/U. For band mapping to be successful one should have τ≪h/U,1/ν as well as τ≫1/Ebg which may become difficult to satisfy for large interactions. Natu et al. suggest that turning off interactions prior to band mapping might be a way to avoid the problem of quasimomentum redistribution that can result from interactions during band mapping.It was pointed out in [200] that noise correlations can be used to characterize the Mott insulating state. Spielman et al. [181] used noise correlations determined from the autocorrelation function of n(k) determined from time of flight images
(129)S(δkx,δky)=〈∫‍dkxdkyn(kx,ky)n(kx+δkx,ky+δky)〉,
where 〈⋯〉 indicates an average over images to locate the transition between superfluid and the Mott insulator in a 2D BHM. This was done by tracking the growth of the area under the noise correlation curve with increasing J/U. Whilst not necessarily an easy measurement, the equilibrium baseline established in [181] might allow for the determination of whether superfluid order persists to U/J values larger than the equilibrium transition in a quench from superfluid to Mott insulator, by measuring autocorrelations of n(kx,ky,t). If there was persistence of superfluid order, this would imply a divergence in the area under the noise correlation curve at a smaller value of J/U than seen in the equilibrium case for times much less than the global equilibration time.Real space images of bosons in optical lattices have become available in the past few years due to the work of several groups on in situ observation of coarse-grained density [146, 201, 202] of atoms in trapped two-dimensional optical lattices and single-site imaging [186, 195, 203, 204]. These techniques allow for detailed characterization of the properties of the Mott insulating phase. By measuring the local compressibility κ(r)=∂n/∂μ(r), Gemelke et al. [201] were able to observe the vanishing of the compressibility in the centre of a trap for which the density was n-=1 as expected for the Mott insulator. In the single site fluorescence imaging used in [186, 195, 203], the technique allows the determination of whether there is an even or an odd number of bosons on a site. This has allowed for the imaging of the “wedding cake” structure of the BHM in a trap with a Mott insulator with four bosons per site in the centre [186]. It has also allowed the measurement of correlation functions such as the time-dependent parity correlation function illustrated in Figure 21 and the prospect of further out-of-equilibrium measurements of number statistics and correlations to come.In addition to the techniques discussed above, methods such as lattice modulation [179, 180, 205, 206] or Bragg spectroscopy [207–209] allow for the measurement of collective excitations although it is less clear that such finite frequency measurements can be useful for out-of-equilibrium dynamics where time translation invariance is broken.


#### 4.5.2. Thermometry

There are two crucial challenges relating to thermometry in cold atom experiments—the first is to achieve cooling to very low temperatures, so that the quantum phenomena of interest are accessible, the second to develop thermometry to establish that the desired temperatures have in fact been reached [210].In order for the physics of quantum quenches to be observable in experiment, it is essential that atoms be cooled to very low temperatures and that changes of parameters in time (e.g., changing the depth of the optical lattice) do not lead to significant heating [210–212]. Whilst out-of-equilibrium experiments are not carried out adiabatically, it is useful to compare to an adiabatic baseline. Pollet et al. [197] considered the temperature changes for isentropic trajectories in which U/J was increased for both the 1- and 2-dimensional BHM in a trap. In both cases, they found that βJ decreased, and βU increased, with β the inverse temperature. These results would suggest that for low initial temperature, it is possible to reach the Mott insulator phase in a quench from a superfluid even given concerns about the possibility of heating [211, 212]. A complicating factor for experiments involving the BHM is that the zero temperature phase diagram is considerably simpler than that at finite temperature. The Mott insulating phase is only truly incompressible at zero temperature, but signatures of strongly suppressed compressibility are calculated to persist to temperatures of order kBT*~0.2U [39, 213].   Additionally, at finite temperature, the portion of the phase diagram with a superfluid state shrinks with increasing temperature leading to normal as well as superfluid regions—all of these regimes may be present simultaneously in a parabolic trap. The possibility of overlapping quantum critical regions [60] at higher temperatures futher complicates the picture.As stated succinctly in the review by McKay and DeMarco [210], this leads to the challenge “to develop thermometry methods that do not rely on unverified theoretical results and that can be experimentally validated”. Techniques that have been used or suggested for thermometry include using the quasimomentum distribution [58, 181, 214] or the in-trap size of the gas [214] or probes that make use of real space imaging. In the Mott phase, either fitting the density distribution [146, 201] or site occupancy statistics have been used for thermometry [186]. A conceptually attractive approach is to use the fluctuation dissipation theorem, which can be expressed within the LDA as [215, 216]
(130)kBT∂〈n^(r)〉∂μ(r)=∫‍dr′[〈n^(r)n^(r′)〉-〈n^(r)〉〈n^(r′)〉],
where the 〈n^(r)〉 are to be understood as coarse-grained over some small region in space. Gemelke et al. [201] implemented this approach and found results that were in qualitative agreement with the fluctuation dissipation theorem. The advantage of the fluctuation-dissipation approach to thermometry is that it is model independent, but it does rely on the assumption that the LDA holds [216].The validity of the LDA for cold atoms in a trap was recently questioned based on the phase diagram obtained in experiments in [57] (Figure 15). Zero temperature QMC calculations [192] along with the LDA give a fair account of the experimental data, but there are systematic deviations in the data.  Jiménez-García et al.  [57] argued that these indicated the breakdown of the LDA, a point of view that was also expressed by Mahmud et al. [39], who performed finite temperature QMC calculations for cold atoms in a trap with T/J=1 (as compared to the reported Tinitial/J=0.9 in the experiment). These calculations were also not in complete agreement with experiment but captured the possibility that the could be points reported as the Mott insulator (i.e., nonsuperfluid) below the line indicating the edge of the Mott insulator phase in Figure 15. Obtaining quantitative agreement between experiment and QMC simulations for the phase boundaries of trapped bosons in an optical lattice is still work in progress, especially for stronger interactions [56].From the theory side, there has been considerable work aside from the QMC simulations mentioned above on the properties of the BHM at finite temperature [48, 56, 112, 125, 126, 197, 198, 211, 213, 217–220], but most has focused on the equilibrium properties of the model. There is still plenty of scope to consider the effect of thermal initial states on out-of-equilibrium dynamics.


## 4.1. Quantum Phase Revivals

A related experiment to the crossing of the superfluid Mott insulator transition demonstrated in [6] is that of phase revivals [183, 193]. In such an experiment, an initially superfluid state is quickly quenched into the Mott phase by increasing the lattice depth nonadiabatically, so that tunnelling between sites is strongly suppressed. The state is held in the Mott phase for a time t, and then the trapping potential and the optical lattice potential are turned off and a time of flight image is taken. Greiner et al. [183] monitored the collapse and revival of coherent peaks in the interference pattern as illustrated in Figure 16 to determine the ratio of coherent atoms as a function of t.Figure 16
Ratio of coherent atoms as a function of timet after a quench from the superfluid phase to the Mott insulating phase [183]. The initial lattice depth was VA=8ER and the final lattice depth VB=35ER. The quench from VA to VB was made in 50 μs. Adapted with permission from Macmillan Publishers Ltd: Nature, M. Greiner, O. Mandel, T. Esslinger, T. W. Hänsch, and I. Bloch, Nature 419, 51 (2002), copyright 2002.The simplest picture to explain this was outlined in [183]. In the infinite system size limit, one can view an initial state which is a noninteracting Bose Einstein Condensate (BEC) as a product of local coherent states [194] for N bosons on M sites:
(113)|φ0〉=∏i=1M‍e-|αi|2/2eαia^i†|vac〉=∏i=1M‍e-|αi|2/2∑ni‍αinini!|ni〉,
where the coherent state on site i is parametrized by αi,  |ni〉 is a state with ni bosons on site i, and |vac〉 is the vacuum state. After a quench J/U→0, the time evolution arises solely from H^0 and so the time evolution operator factorizes and for a Fock state with occupation ni on site i, takes the form [97, 183] (with ℏ=1)
(114)U(t)=∏i=1M‍e(i/2)Uni(ni-1)t.
This leads to the time dependence
(115)|φ(t)〉=∏i=1M‍e-|αi|2/2∑niαinini!e(i/2)Uni(ni-1)t|ni〉,
and since ni(ni-1)/2 is an integer, the state will have a periodic time dependence, with a  revival time of trev=2π/U. Hence, the coherent fraction
(116)NcohN=1N∑i|ψi(t)|2,
where ψi(t)=〈φ(t)|a^i|φ(t)〉 will also have a periodic time dependence with a revival time trev.However, it was pointed out by Schachenmayer et al. [194] that this picture is too naive. This is because the state φ0 implies a coherent superposition of states with different numbers of particles. This is to be contrasted with the situation in a quantum revival experiment in which there are a fixed number of particles N. The authors of [194] used a number conserving initial state
(117)|BECN〉=1N!NN(∑i=1M‍βia^i†)N|vac〉,
which may be shown to have the form |BECN〉=PN|φ0〉, where PN projects states onto a subspace of the Hilbert space with fixed total particle number N. Schachenmayer et al. calculated
(118)nq(t)=1M∑i≠jM‍〈a^i†a^j〉teiq(xi-xj),
for both a homogeneous and a trapped system, making use of the projected state with fixed N. In a homogeneous system, they found that their results converged to the naive calculation for M≃5, but for bosons in a trap, they found that for N=M, the relative error in the naive calculation scales as ~1/N, corresponding to about a 1% error for a system with 50 sites. This is quantitatively small, but nevertheless, conceptually important.The simplest analysis of quantum phase revivals ignores tunnelling during the quench, or nonzero interactions before the quench that modify the initial state from a purely coherent state. Schachenmayer et al. considered the effect of this physics using TEBD and found that it could lead to damping of revivals with time, which can be seen in Figure16. Previous investigations of the effect of tunnelling on quantum phase revivals [108, 115] found similar results.The time dependence of quantum phase revivals has also been used to show evidence of multiparticle interactions that arise from virtual transitions of bosons from the lowest Bloch band to higher Bloch bands, up to six-body interactions [188]. Tiesinga and Johnson [106] showed (using a mean field approximation) that these multiband effects in quantum phase revivals can be used to infer the number squeezing in the initial superfluid state.


## 4.2. Quench from Superfluid to Mott Insulator

Two recent experiments [146, 186] explored relaxation after a quench from a superfluid to a Mott insulator, yet found rather different results. Bakr et al. [186] found fast relaxation, whereas Hung et al. [146] found slow relaxation, relative to natural timescales in the dynamics. I review these experiments and discuss the time-dependent Gutzwiller mean field calculations by Natu et al. [94] that suggest a picture of two timescales for equilibration in a trap.


### 4.2.1. Fast Dynamics

In an important step in the imaging of cold atom systems, Bakr et al. [186] demonstrated single atom-single lattice site imaging that allows for site-resolved observations of dynamics in both space and time. They applied this imaging to follow the out-of-equilibrium dynamics after a quench from a superfluid to a Mott insulator. In an experiment using a few thousand Rb87 atoms in a two-dimensional optical lattice, they used the following protocol. After loading atoms into the trap, they adiabatically increased the depth of their trap to V0=11ER (corresponding to a superfluid state). They then ramped up the depth of their trap to V0=16ER, in a time tramp that ranged between 0.2 and 20 ms. They observed the probability that there were an odd or even number of atoms on each site (corresponding to n-=1 or n-=2, respectively, as a function of tramp as illustrated in Figure 17). The time constant for n=1 was found to be 3.5±0.5 ms and for n=2 it was found to be 3.9±1.3ms. These timescales were noted to be short compared to the timescale h/Jc~68ms (where Jc is the critical value of the hopping) and comparable to the timescale h/Uc~4.1ms.Figure 17
The probability of having an odd number of atoms on each site as a function of ramp time averaged over 19 data sets. The squares correspond to sites withn=1 and the circles to sites with n=2 and the curves are exponential fits to the data. The inset shows the ramp used in the experiment, with a fixed adiabatic exponential ramp to lattice depth V0=11ER followed by a linear ramp with variable duration  [186]. From W. S. Bakr, A. Peng, M. E. Tai, R. Ma, J. Simon, J. I. Gillen, S. Fölling, L. Pollet, and M. Greiner, Science 329, 547 (2010). Adapted with permission from AAAS.


### 4.2.2.  Slow Dynamics

An example of very slow dynamics after a quench from superfluid to the Mott insulator was realised in the experiment of Hung et al. [146]. Studying the dynamics of trapped Cs133 atoms using in situ imaging, they prepared a condensate which was then loaded into a pancake trap, and then a two dimensional optical lattice was slowly turned on over 30 ms, to reach a lattice depth of V= 0.4ER. The lattice was then ramped to a final depth of Vf up to about 13ER (corresponding to U/J=41, implying a Mott insulator in a homogeneous system) using the protocol [101]
(119)V(t)=Vf(1+γ)1+γe4(t-τc)2/τc2,
so that V=Vf at time t=τc (chosen to be 20 ms for the data shown in Figure 18), and γ is chosen so that V(0)=0.4ER. At values of Vf higher than ~13ER, the timescale of the dynamics started to reach that of three-body recombination, limiting the resolution of mass transport. They characterized the relaxation to equilibrium by investigating the mean square deviation of the density after hold time thold from its long time (thold of order 500–800 ms) limit as
(120)Δ(thold)={∑i[n-i(thold)-n-longtimes]2}.
They found that even for timescales longer than the microscopic timescale t~h/J~ 15–75ms, there was still decay of Δ, which could be fitted with a single exponential decay. These results are illustrated in Figure 18.The evolution of the density profile following a 20 ms lattice ramp to (a)Vf=10ER and (b) Vf=12ER. The density profiles at three different hold times are illustrated: 0 (squares), 200 ms (circles), and 500 ms (triangles). The insets show the time evolution of Δ for each value of Vf, each fitted to a single exponential. The fits of the decay times against Vf are shown in (c), along with the appropriate phase of the system [146]. Adapted figure with permission from C.-L. Hung, X. Zhang, N. Gemelke, and C. Chin, Phys. Rev. Lett. 104, 160403 (2010). Copyright 2010 by the American Physical Society.
(a)
(b)
(c)In the deep trap limit, kinetic energy can be essentially ignored, and the average number of particles per site can be calculated fromn-=∑nnPn, where Pn=e-β(H-μn)/Ξ, with H≃(U/2)n(n-1) and the grand canonical partition function Ξ=∑ne-β(H-μn). Hung et al. induced three-body recombination losses after different hold times. By comparing these losses with the expectation for thermally occupied sites, for a hold time of thold~800ms they calculated an effective temperature for the insulating centre of the cloud of Teff~6 nK, considerably lower than the Teff~20 nK that they found in the superfluid wings of the trap. These results suggest that although there may have been local equilibration, even on the large timescales considered, global equilibrium was not established.


### 4.2.3. Two Time Scales

At first glance the results obtained by Bakr et al. [186] and Hung et al. [146] do not seem to be particularly compatible in that in one case relaxation times much less than h/J were observed and in the other relaxation times larger than h/J were observed. However, calculations by Natu et al. [94] illustrate how they can be reconciled within a picture of two timescales for relaxation, specifically, very quick local equilibration and then slower mass transport. They noted that the initial state before the quench can strongly influence equilibration depending on whether it leads to impeded mass transport or not.Natu et al. [94] used the time-dependent Gutzwiller approach to study the Bose Hubbard model with an external trapping potential. By making use of the explicit lattice depth dependence of J and U (14), they tracked the effects of changing the depth of the lattice with time in the manner
(121)V(t)=Vi+(Vf-Vi)tτr,
where Vi,f are the initial and final values of the lattice potential and τr is the ramp time. By comparing a homogeneous system with 1 atom per site with different initial configurations of bosons in a trap as in experiments in [146, 186, 195], Natu et al. argued that local equilibration will be fast, but the speed of global relaxation is strongly affected by whether there needs to be mass transport from one region of the trap to another, particularly if there is a Mott insulating phase, which can act as a barrier to particle diffusion (this is the “Mott barrier” idea discussed by Bernier et al. [99]).Evidence of slow relaxation and local, but not global equilibration is illustrated in Figure19 for an initially superfluid configuration. The figure shows the initial and final density distribution, μ(n(r))+Vtrap, and the coherences Ci(t). In equilibrium, μ(n(r))+Vtrap=μ, but this is not the case in Figure 19, where it can be seen to have local plateaux, but is not globally flat. The plot of coherences Ci=-〈ai〉∑j〈aj*〉 illustrates that for parameters appropriate to Rb87, relaxation timescales can be on the order of 200 ms, much longer than h/Ui~3.3ms, as seen by Hung et al.Slow dynamics after a quench in a trap. (a) Density of particles in the trap as a function of radius—the initial superfluid state is indicated with a solid line, and the final density profile with a dashed line, which differs strongly from the equilibrium configuration (dots). (b)μ(n(r))+Vtrap as a function of radius after the quench—for large values of r there are no particles, and the points represent V(r). (c) Density plot of the time-dependent coherences shows the growing Mott insulating regions that separate the brighter superfluid regions [94]. Adapted figure with permission from S. S. Natu, K. R. A. Hazzard, and E. J. Mueller, Phys. Rev. Lett. 106, 125301 (2011). Copyright 2011 by the American Physical Society.
(a)
(b)
(c)Further confirming the two-timescales picture, Natu et al. investigated other initial conditions in which particle transport across a Mott insulating region was not needed for global equilibration and found that this can lead to fast global equilibration as seen by Bakr et al.


## 4.2.1. Fast Dynamics

In an important step in the imaging of cold atom systems, Bakr et al. [186] demonstrated single atom-single lattice site imaging that allows for site-resolved observations of dynamics in both space and time. They applied this imaging to follow the out-of-equilibrium dynamics after a quench from a superfluid to a Mott insulator. In an experiment using a few thousand Rb87 atoms in a two-dimensional optical lattice, they used the following protocol. After loading atoms into the trap, they adiabatically increased the depth of their trap to V0=11ER (corresponding to a superfluid state). They then ramped up the depth of their trap to V0=16ER, in a time tramp that ranged between 0.2 and 20 ms. They observed the probability that there were an odd or even number of atoms on each site (corresponding to n-=1 or n-=2, respectively, as a function of tramp as illustrated in Figure 17). The time constant for n=1 was found to be 3.5±0.5 ms and for n=2 it was found to be 3.9±1.3ms. These timescales were noted to be short compared to the timescale h/Jc~68ms (where Jc is the critical value of the hopping) and comparable to the timescale h/Uc~4.1ms.Figure 17
The probability of having an odd number of atoms on each site as a function of ramp time averaged over 19 data sets. The squares correspond to sites withn=1 and the circles to sites with n=2 and the curves are exponential fits to the data. The inset shows the ramp used in the experiment, with a fixed adiabatic exponential ramp to lattice depth V0=11ER followed by a linear ramp with variable duration  [186]. From W. S. Bakr, A. Peng, M. E. Tai, R. Ma, J. Simon, J. I. Gillen, S. Fölling, L. Pollet, and M. Greiner, Science 329, 547 (2010). Adapted with permission from AAAS.


## 4.2.2.  Slow Dynamics

An example of very slow dynamics after a quench from superfluid to the Mott insulator was realised in the experiment of Hung et al. [146]. Studying the dynamics of trapped Cs133 atoms using in situ imaging, they prepared a condensate which was then loaded into a pancake trap, and then a two dimensional optical lattice was slowly turned on over 30 ms, to reach a lattice depth of V= 0.4ER. The lattice was then ramped to a final depth of Vf up to about 13ER (corresponding to U/J=41, implying a Mott insulator in a homogeneous system) using the protocol [101]
(119)V(t)=Vf(1+γ)1+γe4(t-τc)2/τc2,
so that V=Vf at time t=τc (chosen to be 20 ms for the data shown in Figure 18), and γ is chosen so that V(0)=0.4ER. At values of Vf higher than ~13ER, the timescale of the dynamics started to reach that of three-body recombination, limiting the resolution of mass transport. They characterized the relaxation to equilibrium by investigating the mean square deviation of the density after hold time thold from its long time (thold of order 500–800 ms) limit as
(120)Δ(thold)={∑i[n-i(thold)-n-longtimes]2}.
They found that even for timescales longer than the microscopic timescale t~h/J~ 15–75ms, there was still decay of Δ, which could be fitted with a single exponential decay. These results are illustrated in Figure 18.The evolution of the density profile following a 20 ms lattice ramp to (a)Vf=10ER and (b) Vf=12ER. The density profiles at three different hold times are illustrated: 0 (squares), 200 ms (circles), and 500 ms (triangles). The insets show the time evolution of Δ for each value of Vf, each fitted to a single exponential. The fits of the decay times against Vf are shown in (c), along with the appropriate phase of the system [146]. Adapted figure with permission from C.-L. Hung, X. Zhang, N. Gemelke, and C. Chin, Phys. Rev. Lett. 104, 160403 (2010). Copyright 2010 by the American Physical Society.
(a)
(b)
(c)In the deep trap limit, kinetic energy can be essentially ignored, and the average number of particles per site can be calculated fromn-=∑nnPn, where Pn=e-β(H-μn)/Ξ, with H≃(U/2)n(n-1) and the grand canonical partition function Ξ=∑ne-β(H-μn). Hung et al. induced three-body recombination losses after different hold times. By comparing these losses with the expectation for thermally occupied sites, for a hold time of thold~800ms they calculated an effective temperature for the insulating centre of the cloud of Teff~6 nK, considerably lower than the Teff~20 nK that they found in the superfluid wings of the trap. These results suggest that although there may have been local equilibration, even on the large timescales considered, global equilibrium was not established.


## 4.2.3. Two Time Scales

At first glance the results obtained by Bakr et al. [186] and Hung et al. [146] do not seem to be particularly compatible in that in one case relaxation times much less than h/J were observed and in the other relaxation times larger than h/J were observed. However, calculations by Natu et al. [94] illustrate how they can be reconciled within a picture of two timescales for relaxation, specifically, very quick local equilibration and then slower mass transport. They noted that the initial state before the quench can strongly influence equilibration depending on whether it leads to impeded mass transport or not.Natu et al. [94] used the time-dependent Gutzwiller approach to study the Bose Hubbard model with an external trapping potential. By making use of the explicit lattice depth dependence of J and U (14), they tracked the effects of changing the depth of the lattice with time in the manner
(121)V(t)=Vi+(Vf-Vi)tτr,
where Vi,f are the initial and final values of the lattice potential and τr is the ramp time. By comparing a homogeneous system with 1 atom per site with different initial configurations of bosons in a trap as in experiments in [146, 186, 195], Natu et al. argued that local equilibration will be fast, but the speed of global relaxation is strongly affected by whether there needs to be mass transport from one region of the trap to another, particularly if there is a Mott insulating phase, which can act as a barrier to particle diffusion (this is the “Mott barrier” idea discussed by Bernier et al. [99]).Evidence of slow relaxation and local, but not global equilibration is illustrated in Figure19 for an initially superfluid configuration. The figure shows the initial and final density distribution, μ(n(r))+Vtrap, and the coherences Ci(t). In equilibrium, μ(n(r))+Vtrap=μ, but this is not the case in Figure 19, where it can be seen to have local plateaux, but is not globally flat. The plot of coherences Ci=-〈ai〉∑j〈aj*〉 illustrates that for parameters appropriate to Rb87, relaxation timescales can be on the order of 200 ms, much longer than h/Ui~3.3ms, as seen by Hung et al.Slow dynamics after a quench in a trap. (a) Density of particles in the trap as a function of radius—the initial superfluid state is indicated with a solid line, and the final density profile with a dashed line, which differs strongly from the equilibrium configuration (dots). (b)μ(n(r))+Vtrap as a function of radius after the quench—for large values of r there are no particles, and the points represent V(r). (c) Density plot of the time-dependent coherences shows the growing Mott insulating regions that separate the brighter superfluid regions [94]. Adapted figure with permission from S. S. Natu, K. R. A. Hazzard, and E. J. Mueller, Phys. Rev. Lett. 106, 125301 (2011). Copyright 2011 by the American Physical Society.
(a)
(b)
(c)Further confirming the two-timescales picture, Natu et al. investigated other initial conditions in which particle transport across a Mott insulating region was not needed for global equilibration and found that this can lead to fast global equilibration as seen by Bakr et al.


## 4.3. Quench from Mott Insulator to Superfluid

A quench from a Mott insulator to a superfluid phase may lead to a proliferation of defects through the Kibble-Zurek mechanism. An experiment that investigated such issues was carried out by Chen et al. [138]. The authors investigated Rb87 in a three-dimensional optical lattice in which the Mott insulating layers emerge for lattice depths greater than V0=13ER. They considered quenches at rates 1/τQ=(d/dt)(J/U) which were sufficiently slow so as not to excite atoms into higher vibrational states: 1/τQ~ (0.001–0.2)Uc/h. The quenches considered were from variable initial lattice depth to a final lattice depth of ~4ER corresponding to J~U.In order to obtain a measure of the number of excitations produced during the quench, Chen et al. obtained time of flight (TOF) images after a relatively long expansion of about 50 ms, to enhance the visibility of vortices and other excitations through large density fluctuations. They fitted the TOF images with a smooth profilef, which was the combination of a Thomas-Fermi profile with a Gaussian and then measured the amount of excitation χ~2 as the deviation from the smooth profile via
(122)χ~2=∑ijχ~ij2=α∑ij((𝒪ij-fij)2/fij)∑ij𝒪ij,
where i and j are pixel indices, 𝒪 is the measured optical depth, and α  is  a  constant that they determined via numerical simulation. They found in each case that a condensate was present, with a condensate fraction that varied between 0.35 and 0.6. Based on numerical simulations of the three dimensional time-dependent Gross-Pitaevskii equation, Chen et al. argued that in the noninteracting limit, χ~2 accurately represents the fraction of the Bogoliubov excitations for a trapped condensate. They also pointed out that in the non-interacting limit, χ~2 should be proportional to the number of excited atoms.In the Mott insulating phase, there are large fluctuations in the relative phases between different sites, and in the Kibble-Zurek mechanism, these are frozen in after a quench is started due to the diverging relaxation time near the critical point. After the system reaches the superfluid phase, these excitations can start to relax. The relaxation of these excitations as measured byχ~2-χ~02 is illustrated in Figure 20 (note that χ~02≃0.06 was determined by an average over all images without a Mott insulating shell). The authors also considered kinetic energy as a function of quench rate to obtain an estimate of the quench-induced heating.Dependence of excitations on quench rate1/τQ. (a), χ~2 with an offset subtracted is shown against quench rate, and (b), the kinetic energy (KE) per particle (with the quench KE0 subtracted) is plotted against quench rate. The experimental protocol for the lattice depth is shown in the inset [138]. Adapted figure with permission from D. Chen, M. White, C. Borries, and B. DeMarco, Phys. Rev. Lett. 106, 235304 (2011). Copyright 2011 by the American Physical Society.
(a)
(b)In the Kibble-Zurek mechanism, the quench rate1/τQ controls defect formation and the expectation is that the number of excitations in 3-dimension scales as [66]
(123)nex~1τQ3ν/(νz+1).
Figure 20 illustrates that (χ~2-χ~02)~1/τQr, with r=0.31±0.03. The experimental value differs substantially for the expectations for either a generic phase transition nex~1/τQ3/4 or for the transition in the XY universality class, for which nex~1/τQ. There are a number of possibilities as to why the predictions of scaling theory do not match with experiments. First, the finite length scales for the phases in the trap may cut off the quantum critical scaling, and thermal effects may also play a role in changing the density of excitations. Additionally, Chen et al. pointed out that the nature of the excitation may affect the scaling, for instance, for vortex excitations nex~1/τQ1/2. Further theoretical investigation of the Mott insulator to superfluid quench in a trap is needed in order to understand these experimental results for the number of excitations.


## 4.4. Quench within a Mott Insulator

Light-cone-like spreading of correlations was recently observed in experiments by Cheneau et al. [196] on the one-dimensional BHM. A two-dimensional gas of Rb87 was loaded into 10 decoupled 1-dimensional chains, with 10–18 particles in each chain, with V0=20Er(U/J=40), corresponding to an n-=1 Mott insulator in each chain. From this initial state, the lattice was reduced in depth in a quench to give values of U/J of roughly 5.0, 7.0 and 9.0. Cheneau et al. measured the two-point parity correlation function [124]
(124)Cd(t)=〈s^j(t)s^j+d(t)〉-〈s^j(t)〉〈s^j+d(t)〉,
where
(125)sj(t)=eiπ[n^j(t)-n-],
and d is the distance between sites. Note that sj gives +1 for no quasiparticles on a site and −1 for the presence of a quasiparticle. The results obtained from measuring Cd(t) and averaging over more than 1000 chains are shown in Figure 21. The correlations showed a front-like spreading, with a constant velocity for d = 2–6, in good agreement with t-DMRG calculations and an analytical model developed by Barmettler et al. [152] discussed in Section 3.7.Figure 21
Time dependence of the parity correlation functionCd(t) measured for bosons in one dimensional chains prepared with U/J=40 and then quenched to U/J=9. The circles are experimental data averaged over more than 1000 chains, the green line is the result of t-DMRG calculations for a homogeneous system, and the dashed-line is the result of calculations with the analytical model developed in [152]. The inset shows the experimental data as a colour map [196]. Adapted with permission from Macmillan Publishers Ltd: Nature, M. Cheneau, P. Barmettler, D. Poletti, M. Endres, P. Schauß, T. Fukuhara, C. Gross, I. Bloch, C. Kollath, and S. Kuhr, Nature 481, 484 (2012), copyright 2012.


### 4.4.1. Superlattice Initial State

Another example of a quench within the Mott state was studied by Trotzky et al. [148] for small one-dimensional systems. They prepared a collection of isolated 1d Bose Hubbard chains, with an average of N=31 particles, each with an initial state |ψ0〉=|…,1,0,1,0,1,…〉 in which all atoms were on “even” sites, prepared using a superlattice potential. They then removed the superlattice in 200 μs and followed the subsequent ensemble average of the filling of bosons on odd sites as a function of time, nodd, which relaxed from 0 to ~0.5 on a timescale of a few ms. Ensemble averaged DMRG calculations were fairly successful in fitting the time dependence of the nodd(t) but gave slightly more accurate results for U/J≤6 when a next-nearest neighbour hopping term was allowed in the BHM, as calculated from the single-particle band structure. The amplitude of the oscillations in nodd was observed to decay as t-α with α=0.86±0.07 over a range of intermediate U/J≲10, which is considerably faster than the value of α=0.5 expected in the limits of no interactions or hardcore bosons (infinite interactions) [148]. This behaviour is captured well by time-dependent DMRG, but is currently a challenge for analytical models. For U/J≥10, the agreement between experiment and theory was less strong—Trotzky et al. suggested this might be due to non-adiabatic heating or interchain tunnelling. This work presented a nice comparison of theoretical and experimental work but highlighted a limitation of t-DMRG, in that the experiment ran for longer times than the calculations could be performed for—the limiting factor being the growth of entanglement entropy [145].


## 4.4.1. Superlattice Initial State

Another example of a quench within the Mott state was studied by Trotzky et al. [148] for small one-dimensional systems. They prepared a collection of isolated 1d Bose Hubbard chains, with an average of N=31 particles, each with an initial state |ψ0〉=|…,1,0,1,0,1,…〉 in which all atoms were on “even” sites, prepared using a superlattice potential. They then removed the superlattice in 200 μs and followed the subsequent ensemble average of the filling of bosons on odd sites as a function of time, nodd, which relaxed from 0 to ~0.5 on a timescale of a few ms. Ensemble averaged DMRG calculations were fairly successful in fitting the time dependence of the nodd(t) but gave slightly more accurate results for U/J≤6 when a next-nearest neighbour hopping term was allowed in the BHM, as calculated from the single-particle band structure. The amplitude of the oscillations in nodd was observed to decay as t-α with α=0.86±0.07 over a range of intermediate U/J≲10, which is considerably faster than the value of α=0.5 expected in the limits of no interactions or hardcore bosons (infinite interactions) [148]. This behaviour is captured well by time-dependent DMRG, but is currently a challenge for analytical models. For U/J≥10, the agreement between experiment and theory was less strong—Trotzky et al. suggested this might be due to non-adiabatic heating or interchain tunnelling. This work presented a nice comparison of theoretical and experimental work but highlighted a limitation of t-DMRG, in that the experiment ran for longer times than the calculations could be performed for—the limiting factor being the growth of entanglement entropy [145].


## 4.5. Interpreting Experiments

In order to interpret and understand experiments on cold atoms in Bose Hubbard systems both imaging and thermometry are key. I discuss both of these in the context of studying out-of-equilibrium dynamics.


### 4.5.1. Imaging

The imaging techniques for cold atoms in optical lattices can roughly be classified as those that measure properties of the system in momentum space, and are most suited to the superfluid phase, and those that measure properties in real space and are more suited to the Mott insulating phase. I will focus first on momentum space.The experiments by Greiner et al. [6] demonstrated the traversal of the quantum critical region separating superfluid and Mott insulator by measuring the momentum distribution of atoms in the trap, taking advantage of the mapping of momentum onto position after atoms are released from the trapping potential. The resulting interference pattern in space takes the form [185]
(126)n(r)=(mℏt)3|w~(k=mrℏt)|2S(k=mrℏt)
and depends on S(k)=∑i,jeik·(ri-rj)〈a^i†a^j〉, multiplied by a Wannier envelope |w~|2. If there is long-range phase coherence, then the correlation function 〈a^i†a^j〉 has a correlation length of the order of the lattice, leading to a sharp peak at k=0 along with satellite peaks due to the presence of a lattice. The disappearance of these peaks was argued to indicate the disappearance of superfluidity and the entry into a Mott phase. However, peaks are also present in both the normal and Mott phases, albeit less sharp [38]. It has been suggested that a better quantity to measure is the visibility [38, 185], which, if measured so that the Wannier envelope cancels, takes the form [185]
(127)𝒱=nmax-nminnmax+nmin=Smax-SminSmax+Smin,
where nmin and nmax are the minimum and the maximum of the density distribution. At small U, the visibility tends to 1, indicating the presence of a superfluid and then drops rapidly for U/zJ larger than the critical value [179, 181, 185]. At lowest order in perturbation theory in J/U, it can be estimated that [185]
(128)𝒱≃43(n0+1)zJU
which is in accordance with calculations using QMC [197] and an analytic strong coupling approach [198] both of which suggest 𝒱~(U/zJ)-1. Gerbier et al. [185] found good experimental agreement with this form up to a lattice depth of V0=30ER corresponding to U/zJ≃200, where they suggested that adiabaticity in the preparation of the state may have broken down.The issue of adiabaticity also arises in imaging as was recently discussed in detail by Natu et al. [199]. In the band mapping protocol lattice strength is reduced over some timescale τ, during which the quasimomentum distribution is mapped to momentum, followed by a ballistic expansion in which the momentum distribution is mapped to position and then imaged as n(r). This can be viewed as an accurate representation of the quasimomentum distribution provided that τ is slow enough (τ≫h/Ebg, where Ebg is the band gap) that the quasimomentum states can adiabatically evolve into momentum states and fast enough (τ≪1/ν, where ν is the trap frequency) that occupation numbers of quasimomentum states do not change. These conditions can be met for a sufficiently weakly interacting Bose gas [199]. However, interactions introduce an additional timescale h/U. For band mapping to be successful one should have τ≪h/U,1/ν as well as τ≫1/Ebg which may become difficult to satisfy for large interactions. Natu et al. suggest that turning off interactions prior to band mapping might be a way to avoid the problem of quasimomentum redistribution that can result from interactions during band mapping.It was pointed out in [200] that noise correlations can be used to characterize the Mott insulating state. Spielman et al. [181] used noise correlations determined from the autocorrelation function of n(k) determined from time of flight images
(129)S(δkx,δky)=〈∫‍dkxdkyn(kx,ky)n(kx+δkx,ky+δky)〉,
where 〈⋯〉 indicates an average over images to locate the transition between superfluid and the Mott insulator in a 2D BHM. This was done by tracking the growth of the area under the noise correlation curve with increasing J/U. Whilst not necessarily an easy measurement, the equilibrium baseline established in [181] might allow for the determination of whether superfluid order persists to U/J values larger than the equilibrium transition in a quench from superfluid to Mott insulator, by measuring autocorrelations of n(kx,ky,t). If there was persistence of superfluid order, this would imply a divergence in the area under the noise correlation curve at a smaller value of J/U than seen in the equilibrium case for times much less than the global equilibration time.Real space images of bosons in optical lattices have become available in the past few years due to the work of several groups on in situ observation of coarse-grained density [146, 201, 202] of atoms in trapped two-dimensional optical lattices and single-site imaging [186, 195, 203, 204]. These techniques allow for detailed characterization of the properties of the Mott insulating phase. By measuring the local compressibility κ(r)=∂n/∂μ(r), Gemelke et al. [201] were able to observe the vanishing of the compressibility in the centre of a trap for which the density was n-=1 as expected for the Mott insulator. In the single site fluorescence imaging used in [186, 195, 203], the technique allows the determination of whether there is an even or an odd number of bosons on a site. This has allowed for the imaging of the “wedding cake” structure of the BHM in a trap with a Mott insulator with four bosons per site in the centre [186]. It has also allowed the measurement of correlation functions such as the time-dependent parity correlation function illustrated in Figure 21 and the prospect of further out-of-equilibrium measurements of number statistics and correlations to come.In addition to the techniques discussed above, methods such as lattice modulation [179, 180, 205, 206] or Bragg spectroscopy [207–209] allow for the measurement of collective excitations although it is less clear that such finite frequency measurements can be useful for out-of-equilibrium dynamics where time translation invariance is broken.


### 4.5.2. Thermometry

There are two crucial challenges relating to thermometry in cold atom experiments—the first is to achieve cooling to very low temperatures, so that the quantum phenomena of interest are accessible, the second to develop thermometry to establish that the desired temperatures have in fact been reached [210].In order for the physics of quantum quenches to be observable in experiment, it is essential that atoms be cooled to very low temperatures and that changes of parameters in time (e.g., changing the depth of the optical lattice) do not lead to significant heating [210–212]. Whilst out-of-equilibrium experiments are not carried out adiabatically, it is useful to compare to an adiabatic baseline. Pollet et al. [197] considered the temperature changes for isentropic trajectories in which U/J was increased for both the 1- and 2-dimensional BHM in a trap. In both cases, they found that βJ decreased, and βU increased, with β the inverse temperature. These results would suggest that for low initial temperature, it is possible to reach the Mott insulator phase in a quench from a superfluid even given concerns about the possibility of heating [211, 212]. A complicating factor for experiments involving the BHM is that the zero temperature phase diagram is considerably simpler than that at finite temperature. The Mott insulating phase is only truly incompressible at zero temperature, but signatures of strongly suppressed compressibility are calculated to persist to temperatures of order kBT*~0.2U [39, 213].   Additionally, at finite temperature, the portion of the phase diagram with a superfluid state shrinks with increasing temperature leading to normal as well as superfluid regions—all of these regimes may be present simultaneously in a parabolic trap. The possibility of overlapping quantum critical regions [60] at higher temperatures futher complicates the picture.As stated succinctly in the review by McKay and DeMarco [210], this leads to the challenge “to develop thermometry methods that do not rely on unverified theoretical results and that can be experimentally validated”. Techniques that have been used or suggested for thermometry include using the quasimomentum distribution [58, 181, 214] or the in-trap size of the gas [214] or probes that make use of real space imaging. In the Mott phase, either fitting the density distribution [146, 201] or site occupancy statistics have been used for thermometry [186]. A conceptually attractive approach is to use the fluctuation dissipation theorem, which can be expressed within the LDA as [215, 216]
(130)kBT∂〈n^(r)〉∂μ(r)=∫‍dr′[〈n^(r)n^(r′)〉-〈n^(r)〉〈n^(r′)〉],
where the 〈n^(r)〉 are to be understood as coarse-grained over some small region in space. Gemelke et al. [201] implemented this approach and found results that were in qualitative agreement with the fluctuation dissipation theorem. The advantage of the fluctuation-dissipation approach to thermometry is that it is model independent, but it does rely on the assumption that the LDA holds [216].The validity of the LDA for cold atoms in a trap was recently questioned based on the phase diagram obtained in experiments in [57] (Figure 15). Zero temperature QMC calculations [192] along with the LDA give a fair account of the experimental data, but there are systematic deviations in the data.  Jiménez-García et al.  [57] argued that these indicated the breakdown of the LDA, a point of view that was also expressed by Mahmud et al. [39], who performed finite temperature QMC calculations for cold atoms in a trap with T/J=1 (as compared to the reported Tinitial/J=0.9 in the experiment). These calculations were also not in complete agreement with experiment but captured the possibility that the could be points reported as the Mott insulator (i.e., nonsuperfluid) below the line indicating the edge of the Mott insulator phase in Figure 15. Obtaining quantitative agreement between experiment and QMC simulations for the phase boundaries of trapped bosons in an optical lattice is still work in progress, especially for stronger interactions [56].From the theory side, there has been considerable work aside from the QMC simulations mentioned above on the properties of the BHM at finite temperature [48, 56, 112, 125, 126, 197, 198, 211, 213, 217–220], but most has focused on the equilibrium properties of the model. There is still plenty of scope to consider the effect of thermal initial states on out-of-equilibrium dynamics.


## 4.5.1. Imaging

The imaging techniques for cold atoms in optical lattices can roughly be classified as those that measure properties of the system in momentum space, and are most suited to the superfluid phase, and those that measure properties in real space and are more suited to the Mott insulating phase. I will focus first on momentum space.The experiments by Greiner et al. [6] demonstrated the traversal of the quantum critical region separating superfluid and Mott insulator by measuring the momentum distribution of atoms in the trap, taking advantage of the mapping of momentum onto position after atoms are released from the trapping potential. The resulting interference pattern in space takes the form [185]
(126)n(r)=(mℏt)3|w~(k=mrℏt)|2S(k=mrℏt)
and depends on S(k)=∑i,jeik·(ri-rj)〈a^i†a^j〉, multiplied by a Wannier envelope |w~|2. If there is long-range phase coherence, then the correlation function 〈a^i†a^j〉 has a correlation length of the order of the lattice, leading to a sharp peak at k=0 along with satellite peaks due to the presence of a lattice. The disappearance of these peaks was argued to indicate the disappearance of superfluidity and the entry into a Mott phase. However, peaks are also present in both the normal and Mott phases, albeit less sharp [38]. It has been suggested that a better quantity to measure is the visibility [38, 185], which, if measured so that the Wannier envelope cancels, takes the form [185]
(127)𝒱=nmax-nminnmax+nmin=Smax-SminSmax+Smin,
where nmin and nmax are the minimum and the maximum of the density distribution. At small U, the visibility tends to 1, indicating the presence of a superfluid and then drops rapidly for U/zJ larger than the critical value [179, 181, 185]. At lowest order in perturbation theory in J/U, it can be estimated that [185]
(128)𝒱≃43(n0+1)zJU
which is in accordance with calculations using QMC [197] and an analytic strong coupling approach [198] both of which suggest 𝒱~(U/zJ)-1. Gerbier et al. [185] found good experimental agreement with this form up to a lattice depth of V0=30ER corresponding to U/zJ≃200, where they suggested that adiabaticity in the preparation of the state may have broken down.The issue of adiabaticity also arises in imaging as was recently discussed in detail by Natu et al. [199]. In the band mapping protocol lattice strength is reduced over some timescale τ, during which the quasimomentum distribution is mapped to momentum, followed by a ballistic expansion in which the momentum distribution is mapped to position and then imaged as n(r). This can be viewed as an accurate representation of the quasimomentum distribution provided that τ is slow enough (τ≫h/Ebg, where Ebg is the band gap) that the quasimomentum states can adiabatically evolve into momentum states and fast enough (τ≪1/ν, where ν is the trap frequency) that occupation numbers of quasimomentum states do not change. These conditions can be met for a sufficiently weakly interacting Bose gas [199]. However, interactions introduce an additional timescale h/U. For band mapping to be successful one should have τ≪h/U,1/ν as well as τ≫1/Ebg which may become difficult to satisfy for large interactions. Natu et al. suggest that turning off interactions prior to band mapping might be a way to avoid the problem of quasimomentum redistribution that can result from interactions during band mapping.It was pointed out in [200] that noise correlations can be used to characterize the Mott insulating state. Spielman et al. [181] used noise correlations determined from the autocorrelation function of n(k) determined from time of flight images
(129)S(δkx,δky)=〈∫‍dkxdkyn(kx,ky)n(kx+δkx,ky+δky)〉,
where 〈⋯〉 indicates an average over images to locate the transition between superfluid and the Mott insulator in a 2D BHM. This was done by tracking the growth of the area under the noise correlation curve with increasing J/U. Whilst not necessarily an easy measurement, the equilibrium baseline established in [181] might allow for the determination of whether superfluid order persists to U/J values larger than the equilibrium transition in a quench from superfluid to Mott insulator, by measuring autocorrelations of n(kx,ky,t). If there was persistence of superfluid order, this would imply a divergence in the area under the noise correlation curve at a smaller value of J/U than seen in the equilibrium case for times much less than the global equilibration time.Real space images of bosons in optical lattices have become available in the past few years due to the work of several groups on in situ observation of coarse-grained density [146, 201, 202] of atoms in trapped two-dimensional optical lattices and single-site imaging [186, 195, 203, 204]. These techniques allow for detailed characterization of the properties of the Mott insulating phase. By measuring the local compressibility κ(r)=∂n/∂μ(r), Gemelke et al. [201] were able to observe the vanishing of the compressibility in the centre of a trap for which the density was n-=1 as expected for the Mott insulator. In the single site fluorescence imaging used in [186, 195, 203], the technique allows the determination of whether there is an even or an odd number of bosons on a site. This has allowed for the imaging of the “wedding cake” structure of the BHM in a trap with a Mott insulator with four bosons per site in the centre [186]. It has also allowed the measurement of correlation functions such as the time-dependent parity correlation function illustrated in Figure 21 and the prospect of further out-of-equilibrium measurements of number statistics and correlations to come.In addition to the techniques discussed above, methods such as lattice modulation [179, 180, 205, 206] or Bragg spectroscopy [207–209] allow for the measurement of collective excitations although it is less clear that such finite frequency measurements can be useful for out-of-equilibrium dynamics where time translation invariance is broken.


## 4.5.2. Thermometry

There are two crucial challenges relating to thermometry in cold atom experiments—the first is to achieve cooling to very low temperatures, so that the quantum phenomena of interest are accessible, the second to develop thermometry to establish that the desired temperatures have in fact been reached [210].In order for the physics of quantum quenches to be observable in experiment, it is essential that atoms be cooled to very low temperatures and that changes of parameters in time (e.g., changing the depth of the optical lattice) do not lead to significant heating [210–212]. Whilst out-of-equilibrium experiments are not carried out adiabatically, it is useful to compare to an adiabatic baseline. Pollet et al. [197] considered the temperature changes for isentropic trajectories in which U/J was increased for both the 1- and 2-dimensional BHM in a trap. In both cases, they found that βJ decreased, and βU increased, with β the inverse temperature. These results would suggest that for low initial temperature, it is possible to reach the Mott insulator phase in a quench from a superfluid even given concerns about the possibility of heating [211, 212]. A complicating factor for experiments involving the BHM is that the zero temperature phase diagram is considerably simpler than that at finite temperature. The Mott insulating phase is only truly incompressible at zero temperature, but signatures of strongly suppressed compressibility are calculated to persist to temperatures of order kBT*~0.2U [39, 213].   Additionally, at finite temperature, the portion of the phase diagram with a superfluid state shrinks with increasing temperature leading to normal as well as superfluid regions—all of these regimes may be present simultaneously in a parabolic trap. The possibility of overlapping quantum critical regions [60] at higher temperatures futher complicates the picture.As stated succinctly in the review by McKay and DeMarco [210], this leads to the challenge “to develop thermometry methods that do not rely on unverified theoretical results and that can be experimentally validated”. Techniques that have been used or suggested for thermometry include using the quasimomentum distribution [58, 181, 214] or the in-trap size of the gas [214] or probes that make use of real space imaging. In the Mott phase, either fitting the density distribution [146, 201] or site occupancy statistics have been used for thermometry [186]. A conceptually attractive approach is to use the fluctuation dissipation theorem, which can be expressed within the LDA as [215, 216]
(130)kBT∂〈n^(r)〉∂μ(r)=∫‍dr′[〈n^(r)n^(r′)〉-〈n^(r)〉〈n^(r′)〉],
where the 〈n^(r)〉 are to be understood as coarse-grained over some small region in space. Gemelke et al. [201] implemented this approach and found results that were in qualitative agreement with the fluctuation dissipation theorem. The advantage of the fluctuation-dissipation approach to thermometry is that it is model independent, but it does rely on the assumption that the LDA holds [216].The validity of the LDA for cold atoms in a trap was recently questioned based on the phase diagram obtained in experiments in [57] (Figure 15). Zero temperature QMC calculations [192] along with the LDA give a fair account of the experimental data, but there are systematic deviations in the data.  Jiménez-García et al.  [57] argued that these indicated the breakdown of the LDA, a point of view that was also expressed by Mahmud et al. [39], who performed finite temperature QMC calculations for cold atoms in a trap with T/J=1 (as compared to the reported Tinitial/J=0.9 in the experiment). These calculations were also not in complete agreement with experiment but captured the possibility that the could be points reported as the Mott insulator (i.e., nonsuperfluid) below the line indicating the edge of the Mott insulator phase in Figure 15. Obtaining quantitative agreement between experiment and QMC simulations for the phase boundaries of trapped bosons in an optical lattice is still work in progress, especially for stronger interactions [56].From the theory side, there has been considerable work aside from the QMC simulations mentioned above on the properties of the BHM at finite temperature [48, 56, 112, 125, 126, 197, 198, 211, 213, 217–220], but most has focused on the equilibrium properties of the model. There is still plenty of scope to consider the effect of thermal initial states on out-of-equilibrium dynamics.


## 5. Conclusions and Future Directions

I briefly highlight some possible future directions that might give opportunities for studying out-of-equilibrium dynamics both using quantum quenches and with protocols in the BHM beyond those mentioned already, and then provide some concluding comments relating the material reviewed in this paper.


### 5.1. Future Directions

The following directions are not meant to be an exhaustive list of possible interesting extensions of current research, but each seems to have the potential for new results both in out-of-equilibrium dynamics and broader contexts, such as quantum simulation.


#### 5.1.1. Two and Higher Component Bose Hubbard Models

Two (and higher) component Bose-Hubbard models have started to attract theoretical [5, 221–227] and experimental attention [207, 228] as they allow the prospect of studying bosons with spin (strictly pseudospin). The interplay of spin degrees of freedom with all of the other physics discussed in this paper suggests a smorgasbord of out-of-equilibrium physics to be explored in many varied parameter regimes [229–231] and in exotic quantum phases such as spin liquids [232]. Such models also open the door to studying spin-orbit coupling of bosons [233, 234].


#### 5.1.2. Designer Initial States

Weitenberg et al. [203] recently demonstrated the ability to selectively remove atoms with site resolved accuracy with a fidelity of around 95% from a Mott insulating state. As they point out, this leads to the very exciting possibility of observing quantum dynamics at the single atom level with designer initial states. It has also recently been shown theoretically [235] that for special initial states with clusters of doubly occupied sites in the 1-dimensional BHM, there can be freezing of the dynamical evolution. These results may help to shed light on the nature of dynamical freezing seen after superfluid to Mott insulator quenches.


#### 5.1.3. Novel Geometries

Beyond the cubic lattice geometry considered here, several other geometries have been used to consider the BHM, and in some cases to study out-of-equilibrium dynamics. The creation of toroidal confining potentials, through the intersection of two different red detuned laser beams [236, 237] and the engineering of weak links to control superflow of a condensate in such a trap [237], leads naturally to the thought that it might be possible to create a Bose Hubbard model with a circular geometry. Such a situation was considered in [238], where QMC simulations of such a 1 dimensional BHM were used to determine the existence of a compressible, non-superfluid state in addition to the superfluid and Mott states dubbed the “local Mott” phase. Experimental work on out-of-equilibrium dynamics of the BHM in coupled one dimensional systems [239] has stimulated theoretical investigations of the BHM on ladders [240, 241]. Finally, more exotically, Hamma et al. [242] have studied the BHM with an evolving graph as a toy model for emerging spacetime.


### 5.2. Conclusions

Much has been learned about out-of-equilibrium dynamics in the BHM, both experimentally and theoretically, but there remain many opportunities to improve our understanding of this system. On the theory side, the most pressing need is the ability to have theories or numerical approaches that can go beyond the mean field level and address questions of out-of-equilibrium dynamics in dimensions larger than one accurately in both the superfluid and the Mott insulating regimes. This should allow clearer investigations of the slow dynamics seen in quenches from the superfluid to the Mott insulator and whether (or how) the Lieb-Robinson bounds seen in density and parity correlation functions in one dimension persist to higher dimensions. On the experimental side, investigations in two or higher dimensions of the light-cone physics observed in [196] would be welcomed to investigate the question of Lieb-Robinson bounds in dimensions larger than one in the BHM. In both theory and experiment, there also seem to be opportunities to explore situations where there is a lack of thermalization and possible frozen states in the Mott insulators in more detail through the study of quantities that have proven very useful in glassy systems, such as two time correlations and out-of-equilibrium effective temperatures extracted from the fluctuation dissipation relation (FDR). The separation of fast and slow timescales after a quench into the Mott phase seen in some experiments [146] and theoretical works [94, 97, 112] suggests the possibility of scaling in two time correlation functions, as seen in glassy systems [178] and outlined in Section 3.10. In addition, the study of the two space and two time correlation functions C(r,r′,t,t′) defined in (109) should allow for the study of a growing dynamic correlation length when there is slow dynamics. The FDR expression (130) should hold in equilibrium, but it can also be inverted to define a temporally and spatially varying effective temperature that could be used to monitor the equilibration in space and time after a quantum quench via
(131)kBTeff(r,t)=∫‍dr′[〈n^(r,t)n^(r′,t)〉-〈n^(r,t)〉〈n^(r′,t)〉](∂〈n^(r,t)〉/∂μ(r)).
Work in this direction was already undertaken in [146]: it was noted that the temperature in the centre of the trap was 6 nK after a quench, whereas the temperature in the superfluid wings was 20 nK even after long hold times of 800 ms, showing evidence of a lack of global equilibration. Finally, further study of quenches from the Mott insulating to superfluid phase seems warranted, given the intriguing scalings found by Chen et al. [138] that do not seem to be in full accord with the predictions of the Kibble-Zurek mechanism. Investigation of a wider range of quench rates, especially very slow quench rates might also allow access to some of the predictions of [68, 69] in which it is noted that for slow quench rates, there can continue to be evolution as the system passes through the critical region, contrary to the assumptions in the usual Kibble-Zurek mechanism.Certainly, there appears to be much room for exploration and hopefully surprises in this field in the years ahead.


## 5.1. Future Directions

The following directions are not meant to be an exhaustive list of possible interesting extensions of current research, but each seems to have the potential for new results both in out-of-equilibrium dynamics and broader contexts, such as quantum simulation.


### 5.1.1. Two and Higher Component Bose Hubbard Models

Two (and higher) component Bose-Hubbard models have started to attract theoretical [5, 221–227] and experimental attention [207, 228] as they allow the prospect of studying bosons with spin (strictly pseudospin). The interplay of spin degrees of freedom with all of the other physics discussed in this paper suggests a smorgasbord of out-of-equilibrium physics to be explored in many varied parameter regimes [229–231] and in exotic quantum phases such as spin liquids [232]. Such models also open the door to studying spin-orbit coupling of bosons [233, 234].


### 5.1.2. Designer Initial States

Weitenberg et al. [203] recently demonstrated the ability to selectively remove atoms with site resolved accuracy with a fidelity of around 95% from a Mott insulating state. As they point out, this leads to the very exciting possibility of observing quantum dynamics at the single atom level with designer initial states. It has also recently been shown theoretically [235] that for special initial states with clusters of doubly occupied sites in the 1-dimensional BHM, there can be freezing of the dynamical evolution. These results may help to shed light on the nature of dynamical freezing seen after superfluid to Mott insulator quenches.


### 5.1.3. Novel Geometries

Beyond the cubic lattice geometry considered here, several other geometries have been used to consider the BHM, and in some cases to study out-of-equilibrium dynamics. The creation of toroidal confining potentials, through the intersection of two different red detuned laser beams [236, 237] and the engineering of weak links to control superflow of a condensate in such a trap [237], leads naturally to the thought that it might be possible to create a Bose Hubbard model with a circular geometry. Such a situation was considered in [238], where QMC simulations of such a 1 dimensional BHM were used to determine the existence of a compressible, non-superfluid state in addition to the superfluid and Mott states dubbed the “local Mott” phase. Experimental work on out-of-equilibrium dynamics of the BHM in coupled one dimensional systems [239] has stimulated theoretical investigations of the BHM on ladders [240, 241]. Finally, more exotically, Hamma et al. [242] have studied the BHM with an evolving graph as a toy model for emerging spacetime.


## 5.1.1. Two and Higher Component Bose Hubbard Models

Two (and higher) component Bose-Hubbard models have started to attract theoretical [5, 221–227] and experimental attention [207, 228] as they allow the prospect of studying bosons with spin (strictly pseudospin). The interplay of spin degrees of freedom with all of the other physics discussed in this paper suggests a smorgasbord of out-of-equilibrium physics to be explored in many varied parameter regimes [229–231] and in exotic quantum phases such as spin liquids [232]. Such models also open the door to studying spin-orbit coupling of bosons [233, 234].


## 5.1.2. Designer Initial States

Weitenberg et al. [203] recently demonstrated the ability to selectively remove atoms with site resolved accuracy with a fidelity of around 95% from a Mott insulating state. As they point out, this leads to the very exciting possibility of observing quantum dynamics at the single atom level with designer initial states. It has also recently been shown theoretically [235] that for special initial states with clusters of doubly occupied sites in the 1-dimensional BHM, there can be freezing of the dynamical evolution. These results may help to shed light on the nature of dynamical freezing seen after superfluid to Mott insulator quenches.


## 5.1.3. Novel Geometries

Beyond the cubic lattice geometry considered here, several other geometries have been used to consider the BHM, and in some cases to study out-of-equilibrium dynamics. The creation of toroidal confining potentials, through the intersection of two different red detuned laser beams [236, 237] and the engineering of weak links to control superflow of a condensate in such a trap [237], leads naturally to the thought that it might be possible to create a Bose Hubbard model with a circular geometry. Such a situation was considered in [238], where QMC simulations of such a 1 dimensional BHM were used to determine the existence of a compressible, non-superfluid state in addition to the superfluid and Mott states dubbed the “local Mott” phase. Experimental work on out-of-equilibrium dynamics of the BHM in coupled one dimensional systems [239] has stimulated theoretical investigations of the BHM on ladders [240, 241]. Finally, more exotically, Hamma et al. [242] have studied the BHM with an evolving graph as a toy model for emerging spacetime.


## 5.2. Conclusions

Much has been learned about out-of-equilibrium dynamics in the BHM, both experimentally and theoretically, but there remain many opportunities to improve our understanding of this system. On the theory side, the most pressing need is the ability to have theories or numerical approaches that can go beyond the mean field level and address questions of out-of-equilibrium dynamics in dimensions larger than one accurately in both the superfluid and the Mott insulating regimes. This should allow clearer investigations of the slow dynamics seen in quenches from the superfluid to the Mott insulator and whether (or how) the Lieb-Robinson bounds seen in density and parity correlation functions in one dimension persist to higher dimensions. On the experimental side, investigations in two or higher dimensions of the light-cone physics observed in [196] would be welcomed to investigate the question of Lieb-Robinson bounds in dimensions larger than one in the BHM. In both theory and experiment, there also seem to be opportunities to explore situations where there is a lack of thermalization and possible frozen states in the Mott insulators in more detail through the study of quantities that have proven very useful in glassy systems, such as two time correlations and out-of-equilibrium effective temperatures extracted from the fluctuation dissipation relation (FDR). The separation of fast and slow timescales after a quench into the Mott phase seen in some experiments [146] and theoretical works [94, 97, 112] suggests the possibility of scaling in two time correlation functions, as seen in glassy systems [178] and outlined in Section 3.10. In addition, the study of the two space and two time correlation functions C(r,r′,t,t′) defined in (109) should allow for the study of a growing dynamic correlation length when there is slow dynamics. The FDR expression (130) should hold in equilibrium, but it can also be inverted to define a temporally and spatially varying effective temperature that could be used to monitor the equilibration in space and time after a quantum quench via
(131)kBTeff(r,t)=∫‍dr′[〈n^(r,t)n^(r′,t)〉-〈n^(r,t)〉〈n^(r′,t)〉](∂〈n^(r,t)〉/∂μ(r)).
Work in this direction was already undertaken in [146]: it was noted that the temperature in the centre of the trap was 6 nK after a quench, whereas the temperature in the superfluid wings was 20 nK even after long hold times of 800 ms, showing evidence of a lack of global equilibration. Finally, further study of quenches from the Mott insulating to superfluid phase seems warranted, given the intriguing scalings found by Chen et al. [138] that do not seem to be in full accord with the predictions of the Kibble-Zurek mechanism. Investigation of a wider range of quench rates, especially very slow quench rates might also allow access to some of the predictions of [68, 69] in which it is noted that for slow quench rates, there can continue to be evolution as the system passes through the critical region, contrary to the assumptions in the usual Kibble-Zurek mechanism.Certainly, there appears to be much room for exploration and hopefully surprises in this field in the years ahead.



---

*Source: 393616-2013-06-12.xml*