ISRN Condensed Matter Physics

Volume 2013 (2013), Article ID 393616, 39 pages

http://dx.doi.org/10.1155/2013/393616

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

Department of Physics, Simon Fraser University, 8888 University Drive, Burnaby, BC, Canada V5A 1S6

Received 14 November 2012; Accepted 6 December 2012

Academic Editors: K. Kasamatsu and A. A. Kordyuk

Copyright © 2013 Malcolm P. Kennett. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### 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 ^{4}He in porous media [4]. The proposal by Jaksch et al. [5] that the BHM could be realized by cold atoms in an optical lattice and the subsequent experimental demonstration of a superfluid to Mott insulator transition in this system by Greiner et al. [6] has lead the focus of work on this model to shift to cold atoms. The tunability of parameters in cold atom systems [7, 8], particularly as a function of time, naturally leads to interest in the out-of-equilibrium dynamics of the BHM especially in the vicinity of a quantum critical point [6].

The out-of-equilibrium dynamics of interacting quantum systems is an area of very active research, due to the challenge of understanding such a nontrivial problem. Unlike equilibrium statistical mechanics, where there are clear prescriptions for determining the state of a system, in out-of-equilibrium systems, not only do the current parameters of the Hamiltonian play a role, but also its history. Finding general principles and widely applicable calculational approaches to obtain greater insight into these problems are important goals of experimental and theoretical work in this field. In this context, the BHM offers up an example of a system in which there is strong interplay between theory and experiment and hope for advancing these goals, both for the BHM itself, and potentially in a broader context.

In addition to the interest in understanding the out-of-equilibrium dynamics of interacting bosons, there has been much excitement about the prospects for using cold atoms in optical lattices for quantum simulations [9]. Any quantum simulator should be characterized and reliable on well-understood problems before being used to simulate poorly understood problems. The temporal evolution of quantum systems might well be a future task for a quantum simulator, and hence it is worthwhile to try to characterize the out-of-equilibrium dynamics for one of the simplest models of interacting bosons.

There have been a number of extensive reviews on cold atoms in optical lattices, which include discussions of the BHM [10–12] and a recent review on the out-of-equilibrium dynamics of closed interacting quantum systems [13]; hence, one might ask whether a review of out-of-equilibrium dynamics in the BHM is required. The large number of publications on this topic in the past ten years and the current activity in the field suggest (at least to this author) that it may be valuable to provide a narrower focus to take stock of what has been learned and what directions may be fruitful to explore in the future.

On the theoretical side, there are a number of challenges to studying the out-of-equilibrium dynamics of the BHM. In particular, even though there are a wide array of theoretical approaches that have been used to study the BHM, all have their drawbacks. Of analytic approaches, most are accurate only in some portion of parameter space, either for weak interactions or strong interactions, or have limitations that there are a large average number of bosons per site. Numerical approaches can be essentially exact but are limited in the size of system that can be simulated due to finite computational resources and for the same reason may be restricted in the dimensions for which they are useful. As a consequence, especially in situations such as a quantum quench, which may require a description of both weakly and strongly interacting limits of the BHM within the same calculation, it may be challenging to find methods that are fully trustworthy for all parameters in the calculation. Nevertheless, this should not discourage such efforts but does imply that the results that should generally be most trusted are those that can be derived from a variety of different methods, with differing strengths and weaknesses.

On the experimental side, as mentioned above, the main focus of work on the BHM is in the context of cold atoms in optical lattices. The flexibility in engineering the dimensionality of the lattice and parameters in the BHM, as well as their temporal dependence through appropriate choices of laser configurations and parameters, makes these systems well suited to studying out-of-equilibrium dynamics. Additionally, to a very good approximation, cold atoms in optical lattices can be treated as isolated quantum systems, so the interaction of the system with an environment is not a primary consideration and the observed dynamics is intrinsic to the system. However, there is the complicating feature that cold atoms require a harmonic confining potential as well as the optical lattice potential, as opposed to the homogeneous systems often preferred by theorists. The trap introduces an inhomogeneous density profile, so that there can be multiple phases present in the same trap, which allows a broad range of parameter space to be probed in a single experiment but places constraints on the ability to observe quantum critical behaviour. Out-of-equilibrium effects in imaging methods and thermometry also complicate the interpretation of some experiments.

In addition to the quantum quench protocol mentioned above, several other protocols have been considered for studying out-of-equilibrium physics in the BHM, both theoretically and experimentally. These include periodic driving in time [14–20], setting up states with finite superfluid current [21–23], adding dissipation [24, 25], and setting up states with a linear potential (artificial electric field) [26–28]. I have chosen to focus primarily on out-of-equilibrium dynamics associated with quantum quenches, which have attracted the bulk of the attention of work on out-of-equilibrium dynamics in the BHM, but wish to highlight these other efforts in the interest of completeness.

This paper is structured as follows: in 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, 30–32] as described in Section 2.1 which leads to the well-known Mott insulator lobes illustrated in Figure 1. Considerable effort has been expended to obtain more accurate determinations of the phase boundaries, using quantum Monte Carlo [33–40], series expansions [41–43], and the density matrix renormalization group [44, 45]. The mean field phase diagram is qualitatively similar to more exact calculations in dimensions 2 and higher but differs greatly from accurate determinations of the phase boundary in 1 dimension, as illustrated in Figure 2.

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, 56–58] have indicated good agreement for equilibrium systems. In a particularly detailed recent study, Trotzky et al. [56] compared measurements of the critical temperature for the finite temperature phase transition between a normal fluid and the superfluid phase as a function of 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 [65–67]. 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 [72–75].

#### 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, 76–85] as the resulting systems give examples of out-of-equilibrium dynamics in interacting quantum systems, a class of problem that is still under active investigation. In the context of the Bose Hubbard model, quenches have generally taken one of two forms. Either a quench from the Mott insulating phase to the superfluid phase or a quench from the superfluid phase to the Mott insulating phase. Extensive theoretical effort has been expended on the effects of time-dependent in the BHM in the form of quantum quenches [66, 86–117].

There are several features that are desirable for any theory of equilibrium or out-of-equilibrium dynamics of the Bose-Hubbard model. First, it should allow for both superfluid and Mott insulating ground states (the Bogoliubov theory fails in this instance). Second, a variety of methods concur on the nature of the collective excitation spectrum in the two phases. In the superfluid phase, there is a gapless sound mode that arises from phase fluctuations of the superfluid order parameter. Near the fixed density quantum critical point (QCP) there is also a gapped (Higgs) mode corresponding to amplitude fluctuations of the superfluid order parameter that has been predicted in a variety of calculations [118–122] and observed experimentally for bosons in two-and-three dimensional optical lattices [123, 124]. In the Mott phase, the collective modes are particle-hole excitations. Third, the approach should be capable of describing nonlocal correlations in space and time.

The methods discussed here mostly have the features outlined above and are: the weak interactions Bogoliubov approach (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) [125–129], the variational cluster approximation (VCA) [130, 131] and Quantum Monte Carlo (QMC) have been used to study the equilibrium properties of the BHM, but to the best of my knowledge these have not been used to address out-of-equilibrium dynamics in the BHM.

##### 3.1. Weak Interactions: The Bogoliubov Approach

The Bogoliubov approach to the BHM is a mean field approach based on the weak interaction limit [47], similar to the Bogoliubov approach to a weakly interacting superfluid in the absence of a lattice. This suggests working in momentum space, reflecting the extended nature of the state and the expectation of a condensate with . 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 [91–93, 136] and can allow for perturbative corrections to short-range correlations [136]. Nevertheless, because of its simplicity, it is a very useful method for gaining an understanding of physics for regimes where exact numerical results are not easily obtained. For example, for out-of-equilibrium dynamics in dimensions higher than 1, for example, [94] (Section 4.2.3), or in the presence of complicated space and time-dependent potentials such as the “Gaussian spoon” consisting of a usual harmonic trap 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, 153–155], but these results indicate that it also appears to hold in the Mott insulating phase of the one-dimensional BHM.

##### 3.8. Closed Time Path Approaches

The closed time path (CTP) or Schwinger-Keldysh technique [156–161] is an approach that allows for the calculation of both equilibrium and out-of-equilibrium quantum dynamics within the same formalism. In CTP methods the problem is formulated in real time (rather than imaginary time, as for the Matsubara formalism), allowing out-of-equilibrium problems to be tackled without analytic continuation. The price to pay is that the number of fields in the theory doubles, with a second copy of each field propagating backwards in time. Consequently in calculations of Green’s functions, the notion of time ordering needs to be replaced by that of contour ordering [159] on a contour such as illustrated in Figure 9. There are several reviews that provide a detailed introduction to CTP techniques [158–161], so few details will be given here.

Several authors have used these approaches to study the BHM both from the weakly interacting limit (large ) [162–166] 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, 173–175]: 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