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.

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 [1012] 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 [1420], setting up states with finite superfluid current [2123], adding dissipation [24, 25], and setting up states with a linear potential (artificial electric field) [2628]. 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 Section 2, 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 where and are annihilation and creation operators for bosons on site , respectively, with the usual commutation relations , and where is the number operator on site . The are the hopping parameters (the choice that will generally be considered here is that for nearest neighbour sites and zero otherwise on a cubic lattice), 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 that contains only single site terms, with and corresponding to the kinetic energy. The term in 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 and .

A phase diagram which is qualitatively correct in dimensions higher than one can be determined from simple considerations and mean field theory. If , 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 , the ground state can be determined by considering the single site Hamiltonian . 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 bosons on each site when the chemical potential satisfies . This phase with an integer number of bosons per site is the incompressible Mott insulator. For finite , 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 and [1, 3032] 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 [3340], series expansions [4143], 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.

In general it is not simple to write down the ground state when both and are finite, but in the limits and , the ground state takes the form of a particularly simple product state. When , one may write the ground state as a product of identical single particle states [46]: where satisfies the Gross-Pitaevski equation, and when , the ground state is a product of Fock states: where 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 , then the hopping terms in the Hamiltonian simplify to and the Hamiltonian may be approximated as a sum of single site terms: Take the hopping term to be a perturbation to the Mott insulator phase, which for has the ground state given in (4) when there are bosons per site. One may calculate the appropriate Landau expansion of the ground state energy as a series in powers of to obtain [47] where , and . An estimate for the phase boundary of the Mott insulator may be found from requiring that which yields the expressions [47] 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 Figure 1 is qualitatively correct for dimensions 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 on both sides of the transition, mean field theory predicts [1, 30, 42, 43, 48], where is the coordination of the lattice ( for a cubic lattice). Quantum Monte calculations and series expansions have found the following positions for the transitions in dimensions up to 3: in [44], in [41], and in [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 a -dimensional BHM is in the universality class of the dimensional XY model [1, 31, 32, 49] for which the lower and upper critical dimensions are and , respectively [1]. Hence, for , 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] 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].

In and , one finds a scaling form for the gap in the Mott insulator of where for , , [50] and for , and [1, 31]. For the generic transition, where the density is not fixed, for example, at fixed chemical potential, then , in , 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 and , in good agreement with the scaling predictions of Fisher et al. [1] of and .

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 The strength of the optical lattice is usually expressed in units of the recoil energy , where is the wavevector of the laser. The lattice constant of the optical lattice is , 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 and a trap potential , then the Hamiltonian for a 3-dimensional system is where is a boson field operator for atoms in a given internal atomic state and is the -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 (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 () [46] Note that the hopping parameter is much more sensitive than the interaction parameter 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 and 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 to the BHM for a homogeneous system, where is the strength of the trap at site . 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 in the trap (assuming rotational symmetry) Within the local density approximation (LDA), the phase of the system at radius may be determined to be that of the homogeneous system found at . This implies that for a sufficiently deep trap and not too large a value of , there can be concentric shells of the Mott insulator and superfluid phases, as follows from the schematic line 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, 5658] 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 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 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 , and the compressibility as , where and are the zero temperature values of density and compressibility at given , , and for . Then a plot of against 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 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 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 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 [6567]. If one focuses on a time-dependent coupling which is close to the critical coupling , then close to the critical point the gap scales as . When gets sufficiently close to , then becomes too small to maintain adiabaticity at some time . At this timescale, the energy scale for excitations will be . The correlation length at time will be , which implies a density of excitations of [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 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 [7275].

3. Theory

In Section 2.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 of and being time-dependent parameters, usually such that the ground states for the initial and final values of are in different phases. In order for the system to be out of equilibrium, the change of 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, 7685] 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 in the BHM in the form of quantum quenches [66, 86117].

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 [118122] 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 (Section 3.1); an infinite dimensional mean field approach (Section 3.2); Gutzwiller mean field theory (Section 3.3); a method for 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) [125129], 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 . Writing the boson annihilation operator in momentum space: with the number of sites, and similarly for the creation operator, one may rewrite the BHM Hamiltonian as [47] with on a -dimensional cubic lattice. The approximation in the Bogoliubov approach comes from noting that the number of condensate atoms , in which case , and hence one can write . Thus, , 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: 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 where , and quasiparticle excitations are related to the by with and their energies are given by Limitations of the Bogoliubov theory are that it is only valid when the number of atoms in the condensate 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 where 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 (setting ): focusing on large fillings , in which case may be used as a small expansion parameter. They wrote the operator as where , with , and . The leading term in the expansion is , 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, and nonlinear corrections . In the superfluid phase well away from the Mott insulator phase, the non-linear corrections may be ignored, provided . Fischer et al. [86] decomposed the operator into phase and number operators and then related the fluctuations to phase and number 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 where and may both be time dependent. If and are both independent of time, these equations lead to the usual Bogoliubov spectrum found above in (22). When is taken to be constant and is allowed to vary as a function of time, Fischer et al. defined a new time variable via , leading to an equation of motion for number fluctuations The solutions of this equation depend sensitively on the way in which is taken to zero as a function of time. If is of the form , then for the number fluctuations oscillate in time with decreasing amplitude and frequency; for the fluctuations decay to zero, and for faster decay of , either as a power law or exponential form , the fluctuations end up being frozen in at a finite value. This can be understood as 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 , the ratio of an internal energy scale , to an external one, the sweep rate, . This allows a notion of slow and fast sweeps. The leading order behaviour of the number fluctuations and phase fluctuations was determined to be [86] A weakness of these results is that they come from extrapolation of the mean field theory for the weak coupling superfluid phase, where is dominant over linear and non-linear fluctuations, into the strong coupling phase—they also apply only in the limit . 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 and , then the state reached via Bogoliubov dynamics can be viewed as an initial state for evolution after the quench with the Hamiltonian . 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 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 , 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 , the lattice spacing, the equation of motion for number fluctuations may be written as where the speed of sound 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 where for , the integral converges and one has . This horizon shrinks as , so at a given time , modes with wavelengths shorter than will oscillate, but longer wavelength modes with will be frozen out. Identifying 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 to note that after a quench the state is also site permutation symmetric. This allows one to parametrize the states with the parameters , where is the fraction of sites with bosons, . Considering the particular example where the maximum number of bosons per site is 2 (a constraint that is easily relaxed), one has the conditions, where is the total number of bosons and is the volume of the system. From these relations, one may specify the state in terms of . Sciolla and Biroli considered quantum quenches in which the interaction strength is taken from to at time and studied the subsequent time evolution of and the superfluid order parameter . They found the interesting result that there is a dynamical transition at a value of , where (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 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.

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 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 with the constraint that [12, 134]. The state is to be viewed as a variational ansatz, and so the are to be found by minimizing For an average of one boson per site, the transition from superfluid to the Mott insulator is found for , as in regular mean field theory.

The Gutzwiller technique generalizes to time-dependent problems fairly straightforwardly [90, 135] by requiring that which leads to the equations [12] where

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 of where the ground state is not close to product form, unlike the two limits and . Improvements to the Gutzwiller mean field have been made using a variety of approaches [9193, 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 with a Gaussian with velocity : 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 which lives on the link between and and projects on to the manifold of states for which the occupation on sites and is . With this projection operator, the hopping term in the BHM may be written as (introducing the hopping operator and replacing the sum over sites with a sum over links) where . The terms take the system out of the low energy manifold in Mott state where there are bosons on each site. The strategy employed in [91] to eliminate these terms from the Hamiltonian was to use a canonical transformation which leads to a transformed Hamiltonain in which terms are kept up to order : A Gutzwiller ansatz was then applied for the canonically transformed state , 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, if is time dependent then the canonical transformation becomes time dependent and one has the following evolution equation: provided . This restricts the initial and final values of to satisfy . On the other hand, it does not place restrictions on the ramp rate , 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 or the defect formation probability (which formed a plateau at large ramp times) where is the state at time after the ramp and is the ground state for the final value of . is the ground state energy determined by minimizing the energy functional for . They observed that the timescale on which one might expect to see universality in and can be estimated as the time for a boson to cross the system of size , which may well exceed the system lifetime, so that the observed behaviour of and 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 like , where and for quenches with a final value of close to the critical (higher values of led to multiple frequency components being present). They also found and , both in disagreement with the exponents expected for a QCP with [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. Corrections to Mean Field Theory

Mean field theory should be exact in the limit of infinite co-ordination , but for finite dimensions, it is known to be inaccurate. Navez and Schützhold [93] developed a systematic expansion in 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 ): and 2 sites (site indices and ) and similarly for higher numbers of sites. The correlated parts of these reduced density matrices, , may be obtained from Navez and Schützhold’s approach is based on the observation that the correlations in the Mott insulator have the scaling hierarchy 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 , then it is also satisfied by , so the hierarchy will remain valid for some finite time (of ) [93]. In the Mott phase with one particle per site . By making use of the hierarchy in (53) to neglect higher order reduced density matrices, they obtained an equation of motion for , which leads to the following equations (with set to unity): In (54), the are the spatial Fourier transforms of the particle and hole correlation functions , , , and , where the local particle and hole operators are and , respectively, and is , where is the connectivity matrix of the lattice.

Navez and Schützhold considered a quench from to a final much larger , and hence assumed initial conditions , and all other . Writing the boson operators of the BHM as where   includes terms involving occupation numbers larger than 2, which decouple due to the scaling hierarchy assumption, they calculated for an instantaneous quench from to a final value of , where and is the critical value of . In this limit, they estimated where is weakly time dependent compared to the exponential, , and with 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 and . 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 . Defining the coarse-grained operator , which is the homogeneous mode which will give the largest eigenvalue of , this mode will have a macroscopic occupation corresponding to the condensate fraction in the region after time . For ,   and for , . This suggests that for distinct regions and with , the relative phase between the regions correlates as where is the distance between the regions and is defined via , and since , one can safely write . 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 and restricted the Hilbert space to states with ,   or 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 , where , with , subject to the constraint . The original bosons of the BHM are represented as . They also introduced modified coherent states of the form , with and then chose a variational wave function of the form The Mott phase is described by . When there is a superfluid order parameter . 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 subject to the constraint . The mean field variational state may be written as , with and creating excitations about this state. They then used the constraint to eliminate from the Hamiltonian and truncated the resulting Hamiltonian at quadratic order in the operators and . After a Bogoliuibov transformation one obtains 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] where and may be expressed in terms of ,  ,  , and . 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 , they obtained the action For large, the critical point is at . From this action they obtained an equation of motion for the amplitude of the superfluid order parameter (after rescaling) 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])

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 mode will be the fastest growing, so that as is increased to near the transition, defects will generally be widely separated by . 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 , motivated by the experiments of Orzel et al. [139], for which . Further work by Huber et al. [118] extended the slave particle approach to allow for 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 theory that 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 to a final value . 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 (note that Kollath et al. denote the boson creation operators of the BHM by rather than 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 is seen to have small amplitude oscillations, with period about a finite value after a quench, as illustrated in Figure 6.

Kollath 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 of . For 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 , 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 , as illustrated in Figure 7, which they characterized as a nonthermal steady state.

Kollath et al. [97] rationalized the slow relaxation they found for large 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 , with 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 where 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, , appear to have been sufficiently close to the transition ( in 1 dimension [44] and in 2 dimensions [41]) for this physics to be observable.

The single particle correlations studied in [97] and also density-density correlations were calculated as a function of time and position with similar and for average fillings 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 . 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 for , for which corresponds to the Mott phase, were on the order of 5-6, whereas the values of for (for which the system remains superfluid) after the change in were of the order of . Additionally, Lauchli and Kollath [98] calculated the change in von Neumann entropy for blocks of length as a function of time: , and found that this showed linear growth for short times and saturated at a time for periodic boundary conditions and for open boundary conditions with 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 for the 1 BHM. By focusing on local observables, such as the local compressibility , and the probability of occupation of a site by bosons, , 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 , 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 ( and ) and for a quench within the Mott insulator ( and ).

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 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 to the Mott insulator with . Comparing the exact numerical solution with TEBD for a system of bosons, they established that their slow ramp had an infidelity of so that their slow ramp could be viewed as adiabatic. They then considered larger systems with 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 . They found that for a ramp in which 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 , 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 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 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 that , 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, if , but have fermionic commutation relations , 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] where the and are fermionic creation and annihilation operators on site , respectively, and is the fermionic number operator on site . 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 when 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 , such that and for the three states on site :  with , the operators satisfy with additional hard core constraints to ensure no double occupancy: . In addition, there is not allowed to be double occupancy of different species of bosons, which implies the constraint Barmettler et al. used a Jordan-Wigner transformation similar to (67), to introduce two species of fermions (reflecting two species of bosons): where are nonlocal string operators which satisfy and . This form removes the hard-core constraints on and individually, but not on the product . This condition may be satisfied with the projector , where , in which case the BHM Hamiltonian may be represented by 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, . With this assumption, the Hamiltonian may be diagonalized with a Bogoliubov transformation to the form where and with Barmettler et al. found that for , 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 and considered an instantaneous quench from this state to a final which remains within the Mott phase. They calculated the time-dependent density-density correlation function for sites separated by and found that for , where and is a Bessel function.

Focusing in particular on the case of , they studied density correlations with both DMRG and the UF approximation for and . 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 limit as Barmettler et al. made use of the fact that the Airy function has a peak for to infer that there is a propagation front with correlations exponentially suppressed for where from which it can be inferred that is the velocity at which correlations spread as for large . Barmettler et al. provided a similar analysis for noninteracting bosons to deduce when . In the regime of validity of the UF approximation, they considered the relative velocity of quasiparticle pairs analytically, which takes a maximum value which corresponds to for and decreases monotonically with increasing from the 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 Section 4.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, 153155], 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 [156161] 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 [158161], so few details will be given here.

Several authors have used these approaches to study the BHM both from the weakly interacting limit (large ) [162166] and the strongly interacting limit (small ) [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 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 where a compact notation has been introduced, with and , and the matrices and have the same form as the corresponding Pauli matrices but are labelled so as not to be confused with spin. The quantities studied are and , where   indicates contour ordering. The theory is specified by the effective action , which takes the form where and is the two-particle irreducible generating functional. can be described as an infinite sum of diagrams from the expansion of interaction terms in of cubic and higher order in [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 , corresponding to the mean field approximation. The Bogoliubov (one-loop) approximation can be obtained if the effective action is taken to be (83) with ignored. Approximations that take into account interactions more carefully include the time-dependent Hartree-Fock-Bogoliubov (HFB) approximation, in which is truncated to keep only the first-order diagram in . This approximation conserves energy and particle number but violates Goldstone’s theorem. The expansion can be viewed as one in , where is time, and hence will be most accurate for short times or small . It is possible to keep higher-order terms in —second-order expansions were considered by Rey et al. [162] and Temme and Gasenzer [165], who both considered a expansion, which simplifies the second-order expansion of . 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 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 [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 , , and particles with and . The Hartree-Fock Bogoliubov approximation is seen to be fairly poor at describing the time development of the condensate, with the approach giving the most accurate account of the time evolution for the condensate population per well as a function of time as becomes larger.

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 and large 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 where is the generating functional when , with the notation , where is imaginary time and is inverse temperature. The term that enters in the action is evaluated with all the fields at the same site, and is the connected local 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 where is the Green’s function in the local limit and is the two particle vertex in the local limit. Sengupta and Dupuis approximated by its static value and introduced . This leads to the effective action which has the attractive feature that it leads to the correct Green’s function in both the strong interaction and noninteracting 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 has the form where with the Fourier transform of the hopping ,   and with where is the number of bosons per site for chemical potential when . In the Mott insulating phase, the modes are gapped and the system becomes unstable to superfluidity when the energy of one or other of the modes vanishes for . In the superfluid phase, the expressions for the spectrum are more complicated, and four modes are found , two of which are gapless as and have the form , and the other two are gapped, as illustrated in Figure 11.

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 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 where is the field at site on contour , with or 2, and with the notation representing the th Pauli matrix that acts in Keldysh rather than spin space. We performed a Keldysh rotation from the , basis to the quantum , classical basis [168, 173175]: where   and are the quantum and classical components of the field and The effect of this on the action is that in the , basis becomes in the , 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 fields (one can show that the Green’s functions for are the same as those for the original field [112, 119]) as where The matrix Green’s functions contains the retarded, advanced, and Keldysh Green’s functions. These Green’s functions, along with the two-particle connected Green’s function are evaluated with respect to the single-site Hamiltonian : 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 . The diagram corresponding to the interaction term is illustrated schematically in Figure 12.

The following symmetry relations hold for the interaction kernel from the definition above: (this corrects [112]) and it may also be shown [112] that Grass et al. [120, 168] also noted similar symmetry relations for this correlation function under interchange of Keldysh indices. These symmetries imply that has only 8 independent components: , and . The remaining four-point function 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 parameter can be derived from the action that applies away from the degeneracy points of the Mott lobes (in which case only the component 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 and defined , where the hopping is taken to be with chosen so that (i.e., lies on the mean field phase boundary for a given ). Hence, the approximate equation of motion is where the coefficients , , and can all be calculated in terms of and [112]. We studied (102) for fixed and time-varying .

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 , and (b) the generic case, in which . The function controls the traversal of the quantum critical point. We chose In our numerics, we used where we introduced a timescale which is the characteristic time for to cross from to [111]. In the vicinity of the transition is linear in . Functional forms of which are not linear in in the vicinity of may lead to differing behaviour in the long-time limit [79].

We solved (102) numerically with for a variety of values of , and found that for the form of the solution is that oscillates in a periodic manner with a magnitude that decreases with increasing . When averaged over a period at times , as would be expected in the Mott insulating state. Defining we found that decreases with increasing without any indication of saturation.

In the generic case in which , we solved (102) numerically for a variety of values of as is displayed with time scaled by in Figure 13 for (well away from both degeneracy and the particle hole symmetric case). Both the real and imaginary parts of are oscillatory functions, each of which has a vanishing average over a period for times greater than zero. The nonzero amplitude of 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 . The collapse of traces with different after scaling the time with indicates that has a finite value even when becomes large.

To determine the fate of the system at long quench times, we defined 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 . Whilst it is true that the spatial dependence of 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 of are in accordance with this idea. There appears to be equilibration as 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 in [4143] 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 and self-energy allowing for systems that are out of equilibrium and spatially inhomogeneous. The retarded Green’s function is defined as where operators are in the Heisenberg picture, , , and the average is taken with respect to the state operator where . The spectral moments of the retarded Green’s function are (note that terms involving derivatives of 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 . 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 . 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 from to

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 to 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 sites, there will be recurrences on times of order . 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 , where the occupations on different sites 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 of is not too large [97], but for larger values of 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 If there is a separation of fast and slow timescales for relaxation then one might expect a quantity such as 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, and : with decaying to zero with a timescale that is much shorter than the typical timescale for decay of , so that is essentially constant whilst decays. If aging is present in the Mott insulators for large , then the simplest decay one might expect for is of the form (with ) where is a monotonically increasing function of , possibly and 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 expansion suggest this may be possible [93] although the velocity found in [93] had quite a different form as a function of and than found in one dimension.

On the technical side, most of the methods considered are accurate when either is small or 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, 179185]. These results and more recent experiments on out-of-equilibrium dynamics in the BHM [138, 146, 186188] 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 , which has the effect of changing the ratio in the BHM as discussed in Section 2.3. As indicated in (14), the hopping, , is much more sensitive to changes in than the interaction strength ; 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 and then decreased in strength with increasing , and no interference pattern was visible for .

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 to , the system held for 20 ms, and then was changed to (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 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 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.

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 interactions , via changing the -wave scattering length , 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 , this implies that the characteristic frequency for the ramp should be such that , where , 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 Section 4.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 , 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 .

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 bosons on sites: where the coherent state on site is parametrized by ,   is a state with bosons on site , and is the vacuum state. After a quench , the time evolution arises solely from and so the time evolution operator factorizes and for a Fock state with occupation on site , takes the form [97, 183] (with ) This leads to the time dependence and since is an integer, the state will have a periodic time dependence, with a revival time of . Hence, the coherent fraction where will also have a periodic time dependence with a revival time .

However, it was pointed out by Schachenmayer et al. [194] that this picture is too naive. This is because the state 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 . The authors of [194] used a number conserving initial state which may be shown to have the form , where projects states onto a subspace of the Hilbert space with fixed total particle number . Schachenmayer et al. calculated for both a homogeneous and a trapped system, making use of the projected state with fixed . In a homogeneous system, they found that their results converged to the naive calculation for , but for bosons in a trap, they found that for , the relative error in the naive calculation scales as ~, 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 Figure 16. 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 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 (corresponding to a superfluid state). They then ramped up the depth of their trap to , in a time 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 or , respectively, as a function of as illustrated in Figure 17). The time constant for was found to be  ms and for it was found to be ms. These timescales were noted to be short compared to the timescale ms (where is the critical value of the hopping) and comparable to the timescale ms.

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 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  0.4. The lattice was then ramped to a final depth of up to about 13 (corresponding to , implying a Mott insulator in a homogeneous system) using the protocol [101] so that at time (chosen to be 20 ms for the data shown in Figure 18), and is chosen so that . At values of higher than ~13, 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 from its long time ( of order 500–800 ms) limit as They found that even for timescales longer than the microscopic timescale  15–75ms, there was still decay of , which could be fitted with a single exponential decay. These results are illustrated in Figure 18.

In the deep trap limit, kinetic energy can be essentially ignored, and the average number of particles per site can be calculated from , where , with and the grand canonical partition function . 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 ms they calculated an effective temperature for the insulating centre of the cloud of  nK, considerably lower than the  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 were observed and in the other relaxation times larger than 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 and (14), they tracked the effects of changing the depth of the lattice with time in the manner where are the initial and final values of the lattice potential and 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 Figure 19 for an initially superfluid configuration. The figure shows the initial and final density distribution, , and the coherences . In equilibrium, , 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 illustrates that for parameters appropriate to , relaxation timescales can be on the order of 200 ms, much longer than ms, as seen by Hung et al.

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 in a three-dimensional optical lattice in which the Mott insulating layers emerge for lattice depths greater than . They considered quenches at rates which were sufficiently slow so as not to excite atoms into higher vibrational states:  (0.001–0.2). The quenches considered were from variable initial lattice depth to a final lattice depth of ~4 corresponding to .

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 profile , which was the combination of a Thomas-Fermi profile with a Gaussian and then measured the amount of excitation as the deviation from the smooth profile via where and 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, accurately represents the fraction of the Bogoliubov excitations for a trapped condensate. They also pointed out that in the non-interacting limit, 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 is illustrated in Figure 20 (note that 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.

In the Kibble-Zurek mechanism, the quench rate controls defect formation and the expectation is that the number of excitations in 3-dimension scales as [66] Figure 20 illustrates that , with . The experimental value differs substantially for the expectations for either a generic phase transition or for the transition in the XY universality class, for which . 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 . 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 was loaded into 10 decoupled 1-dimensional chains, with 10–18 particles in each chain, with , corresponding to an Mott insulator in each chain. From this initial state, the lattice was reduced in depth in a quench to give values of of roughly 5.0, 7.0 and 9.0. Cheneau et al. measured the two-point parity correlation function [124] where and is the distance between sites. Note that gives +1 for no quasiparticles on a site and −1 for the presence of a quasiparticle. The results obtained from measuring and averaging over more than 1000 chains are shown in Figure 21. The correlations showed a front-like spreading, with a constant velocity for = 2–6, in good agreement with t-DMRG calculations and an analytical model developed by Barmettler et al. [152] discussed in Section 3.7.

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 Bose Hubbard chains, with an average of particles, each with an initial state 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, , 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 but gave slightly more accurate results for 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 was observed to decay as with over a range of intermediate , which is considerably faster than the value of 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 , 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] and depends on , multiplied by a Wannier envelope . If there is long-range phase coherence, then the correlation function has a correlation length of the order of the lattice, leading to a sharp peak at 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] where and are the minimum and the maximum of the density distribution. At small , the visibility tends to 1, indicating the presence of a superfluid and then drops rapidly for larger than the critical value [179, 181, 185]. At lowest order in perturbation theory in , it can be estimated that [185] which is in accordance with calculations using QMC [197] and an analytic strong coupling approach [198] both of which suggest . Gerbier et al. [185] found good experimental agreement with this form up to a lattice depth of corresponding to , 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 . This can be viewed as an accurate representation of the quasimomentum distribution provided that is slow enough (, where is the band gap) that the quasimomentum states can adiabatically evolve into momentum states and fast enough (, 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 . For band mapping to be successful one should have as well as 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 determined from time of flight images 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 . Whilst not necessarily an easy measurement, the equilibrium baseline established in [181] might allow for the determination of whether superfluid order persists to values larger than the equilibrium transition in a quench from superfluid to Mott insulator, by measuring autocorrelations of . 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 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 , Gemelke et al. [201] were able to observe the vanishing of the compressibility in the centre of a trap for which the density was 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 [207209] 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 [210212]. 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 was increased for both the 1- and 2-dimensional BHM in a trap. In both cases, they found that decreased, and 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 [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] where the 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 (as compared to the reported 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, 217220], 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, 221227] 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 [229231] 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 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 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.

Acknowledgments

The author thanks Andrew Daley, Tobias Grass, Jesko Sirker, and Joseph Thywissen for conversations relating to cold atoms and the BHM and in particular thanks Denis Dalidovich for extensive discussions and collaboration on the BHM and Jeff McGuirk for helpful comments on the paper. This work was supported by NSERC.