# Nonequilibrium dynamical mean-field theory for bosonic lattice models

## Abstract

We develop the nonequilibrium extension of bosonic dynamical mean field theory (BDMFT) and a Nambu real-time strong-coupling perturbative impurity solver. In contrast to Gutzwiller mean-field theory and strong coupling perturbative approaches, nonequilibrium BDMFT captures not only dynamical transitions, but also damping and thermalization effects at finite temperature. We apply the formalism to quenches in the Bose-Hubbard model, starting both from the normal and Bose-condensed phases. Depending on the parameter regime, one observes qualitatively different dynamical properties, such as rapid thermalization, trapping in metastable superfluid or normal states, as well as long-lived or strongly damped amplitude oscillations. We summarize our results in non-equilibrium “phase diagrams” which map out the different dynamical regimes.

###### pacs:

71.10.Fd, 03.75.Kk, 05.70.Ln, 37.10.Jk## I Introduction

Cold atomic gases trapped in an optical lattice provide a unique play-ground to explore equilibrium and nonequilibrium properties of interacting many-particle systems Morsch and Oberthaler (2006); Bloch *et al.* (2008). They enable an almost ideal realization of the low-energy effective Hamiltonians (the fermionic and bosonic Hubbard models Hubbard (1963); Fisher *et al.* (1989)) which have been studied in the condensed matter context for a long time, and whose properties are still not yet fully understood. A big advantage of cold atoms, as compared to condensed matter systems, is that interaction parameters can be tuned almost arbitrarily, and that the lattice spacings and characteristic time-scales are much larger Jaksch *et al.* (1998). For bosonic atoms, the Mott insulating and superfluid regime can easily be accessed Morsch and Oberthaler (2006) and the experimental control is so precise that the use of cold atoms as “quantum simulators” becomes a realistic option Trotzky *et al.* (2010) (for a recent review see Ref. Kennett, 2013).

A particularly interesting aspect of cold atom experiments is the possibility to study the time-evolution of interacting many-body systems Greiner *et al.* (2002); Sebby-Strabley *et al.* (2007); Will *et al.* (2010); Bakr *et al.* (2010); Bissbort *et al.* (2011); Endres *et al.* (2012); Cheneau *et al.* (2012); Trotzky *et al.* (2012); Braun *et al.* (2014). This was beautifully demonstrated in the seminal work by Greiner et al. Greiner *et al.* (2002), who measured the condensate collapse-and-revival oscillations after a quench in a Bose-Hubbard system from the superfluid to the Mott regime. In contrast to equilibrium, where the phase diagram and correlation functions of the Bose-Hubbard model Capogrosso-Sansone *et al.* (2007) can be computed accurately using Monte Carlo simulations Pollet (2012), the real-time evolution of interacting bosonic lattice systems is a big computational challenge.

In one dimension, density matrix renormalization group (DMRG) methods Schollwöck (2005) can be used to simulate the time-evolution after a quench on relatively large lattices, but a rapid entanglement growth limits the accessible time-scale Trotzky *et al.* (2012). Still, DMRG calculations have provided important insights into the short time dynamics, as measured in 1D optical lattices Cheneau *et al.* (2012); Trotzky *et al.* (2012); Braun *et al.* (2014). Kollath et al. Kollath *et al.* (2007) used non-local correlators to study relaxation and thermalization. They showed that an initially superfluid system is trapped in a nonthermal steady state after quenching the interaction deep into the Mott regime, while thermalization occurs after quenches to intermediate interactions.
Also the eigenstate thermalization hypothesis has been explored Roux (2009); Sorg *et al.* (2014) and debated Rigol (2010); Roux (2010) in this context.
A more recent development is the time-dependent variational Monte Carlo (tVMC) approach that shows good agreement with DMRG in 1D without being limited in time Carleo *et al.* (2012). It has also been applied to 2D systems and is not inherently limited to any dimensionality Carleo *et al.* (2014). While tVMC is well suited for studying the spread of correlations, it is a method that treats finite systems,
which complicates the study of thermalization Biroli *et al.* (2010).

In three dimensions, perturbation theory Trefzger and Sengupta (2011); Kennett and Dalidovich (2011); Dutta *et al.* (2012, 2014), and Gutzwiller mean-field (MF) Huber *et al.* (2007); Wolf *et al.* (2010); Sciolla and Biroli (2010, 2011); Snoek (2011); Krutitsky and Navez (2011) calculations have been performed. Both work in specific regions of the phase diagram, but generally fail to describe finite temperature relaxation and thermalization phenomena.
Hence, while being accessible experimentally Braun *et al.* (2014), out-of-equilibrium phenomena in the three dimensional Bose-Hubbard model remain largely unexplored Trotzky *et al.* (2012); Cheneau *et al.* (2012); Braun *et al.* (2014) from the theoretical point of view.
Describing the generic relaxation phenomena and nonthermal transient states, as well as mapping out the different dynamical regimes of this model is fundamental to our understanding of nonequilibrium lattice bosons.
A clear picture of the nonequilibrium properties of the homogeneous bulk-system is also important for the interpretation of more complicated experimental set-ups.
For example, one open question is whether damped superfluid collapse-and-revival oscillations are a dynamical feature of the homogeneous system, or an effect of the trapping potential or other processes not considered in the Bose-Hubbard description Greiner *et al.* (2002); Will *et al.* (2010).

A computationally tractable and promising scheme, which allows to address such issues, is the nonequilibrium generalization of bosonic dynamical mean field theory (BDMFT).
This method is formulated in the thermodynamical limit, and thus enables the study of relaxation and thermalization phenomena in infinite systems Aoki *et al.* (2014).
The equilibrium version of BDMFT Byczuk and Vollhardt (2008); Hubener *et al.* (2009); Anders *et al.* (2010, 2011) produces phase diagrams, condensate fractions, and correlation functions with remarkable accuracy Anders *et al.* (2011).
While the extension of this formalism to nonequilibrium situations is analogous to the fermionic case Aoki *et al.* (2014), and essentially involves the replacement of the imaginary-time interval by a Kadanoff-Baym contour, there are a number of practical challenges.
The most important one is the development of a suitable bosonic impurity solver. The exact continuous-time quantum Monte Carlo (CT-QMC) impurity solver of Ref. Anders *et al.*, 2010 cannot easily be applied to nonequilibrium problems, because of a dynamical sign problem Werner *et al.* (2009), while exact diagonalization based solvers are even more limited than in the fermionic case Gramsch *et al.* (2013), due to the larger local Hilbert space. Weak-coupling perturbation theory is not an option if one is interested in Mott physics. Instead, we will develop and benchmark an impurity solver based on the lowest order strong-coupling perturbation theory, i.e. the non-crossing approximation (NCA) Keiter and Kimball (1971).
As a first application of this new scheme, we will map out the different dynamical regimes of both the symmetric and symmetry broken states, searching for thermalization and trapping phenomena after a quench of the interaction parameter.

This paper is organized as follows: In section II we give an overview of the Bose-Hubbard model, the nonequilibrium generalization of BDMFT [Sec. II.1], the NCA impurity solver [Sec. II.2], the energy calculations [Sec. II.3], and our numerical implementation [Sec. II.4]. In Sec. III we first present benchmark calculations showing density and energy conservation and discuss the lowest order spectral moments [Sec. III.1]. The dynamical regimes in the normal phase are mapped out in Sec. III.2. In Sec. III.3 we consider superfluid initial states, and after an overview of the relaxation regimes in Sec. III.3.1, we study the dynamics for short times in Sec. III.3.2, and long times in Sec. III.3.3. The findings are summarized in Sec. III.3.4 in the form of a nonequilibrium “phase diagram”. Sec. IV is devoted to conclusions. We also provide a derivation of nonequilibrium BDMFT in Appendix A, and discuss the details of the Nambu generalization of NCA in Appendix B.

## Ii Theory

We consider the simplest model for bosonic atoms in an optical lattice, namely the Bose-Hubbard model Fisher *et al.* (1989); Jaksch *et al.* (1998)

(1) |

where () and are the bosonic creation (annihilation) and number operators acting on site , is the chemical potential, and the local pair interaction which competes with the nearest neighbor hopping that we take as our unit of energy.

### ii.1 Nonequilibrium bosonic dynamical mean-field theory

By extending the equilibrium bosonic dynamical mean field theory (BDMFT) Anders *et al.* (2010, 2011) to the three-branch Kadanoff-Baym contour () Aoki *et al.* (2014); Stefanucci and van
Leeuwen (2013), we obtain the bosonic impurity action

(2) | ||||

where is the Nambu spinor , the hybridization function, and the effective symmetry breaking field, which is defined in terms of , the local condensate fraction , and the lattice coordination number as

(3) |

For a detailed derivation of this action, see App. A. Note that the single-particle fluctuations (the term in Eq. (2) and Eq. (3)) enter as a correction to the mean-field action Sachdev (1999), which would be obtained by taking the infinite dimensional limit at fixed (or analogously ) Anders *et al.* (2011).

The solution of the impurity model yields the connected impurity Green’s function

where is the time-ordering operator on the contour , and the local condensate fraction is

The BDMFT self-consistency loop is closed by computing the lattice Green’s function from , and then expressing the hybridization function in terms of the local lattice Green’s function (at self consistency ) Aoki *et al.* (2014).
In the present study we employ the simplified self-consistency relation

and set , which corresponds to a non-interacting semi-circular density of states (DOS) with the same bandwidth and lattice coordination number as the 3D cubic lattice with nearest neighbor hopping .

### ii.2 Non-crossing approximation impurity solver

The previous BDMFT equilibrium studies employed a hybridization expansion CT-QMC impurity solver Anders *et al.* (2010, 2011). However, the extension of this technique to the contour-action in Eq. (2) does not look promising, because the dynamical sign problem from the expansion along the real-time branches Werner *et al.* (2010) will add to the inherent sign problem of the hybridization expansion (in the superfluid regime).
We therefore solve the BDMFT effective impurity action using the first order self-consistent strong coupling expansion.
The generalization of strong coupling expansions to real-time impurity problems has been presented in Ref. Eckstein and Werner (2010). To treat the BDMFT effective action in Eq. (2) we have generalized the formalism to systems with symmetry breaking, as discussed in App. B.

In short we follow the standard procedure and introduce pseudo-particle second quantization operators and for each local occupation number many-body state . This maps the local Hamiltonian to a quadratic term , while the hybridization turns into a pseudo-particle interaction. Expanding to first order in gives the NCA of Ref. Eckstein and Werner, 2010 generalized to Nambu formalism.

The corresponding NCA pseudo-particle self-energy consists of the two shell diagrams with a directed hybridization line (see Fig. 1 and App. B.1)

(4) |

where is the pseudo-particle Green’s function, and are Nambu indices, and is the tensor (operator products are implicit matrix products). The pseudo-particle Dyson equation takes the form

where is the static part in Eq. (2), , and denotes cyclic convolution on , Eckstein and Werner (2010).

Within NCA, and are calculated self-consistently, and local observables are determined from the reduced local density matrix , yielding the local condensate as

while response functions must be determined diagrammatically. In particular, the connected single-particle impurity Green’s function is obtained from the bubble diagram without hybridization insertions (see Fig. 1 and App. B.2)

(5) |

### ii.3 Total energy components

The total energy of the system, is the sum of the (connected) kinetic energy , the condensate energy (or disconnected kinetic energy) , and the local interaction energy , .
Using and , is given by Aoki *et al.* (2014)

depends on as

and can be written in terms of and as

### ii.4 Numerical implementation

We solve the pseudo-particle Dyson equation using a fifth order multi-step method Brunner and van der Houwen (1986); Eckstein and Werner (2010) on an uniformly discretized time grid. To ensure negligible real-time discretization errors we monitor the total energy and density, which both are constants of motion of the conserving NCA Eckstein and Werner (2010) (the gauge property ensures ). In principle the local Fock space is unbounded, but for it can safely be truncated, keeping only states. The cut-off error is controlled by monitoring the drift in away from unity. Close to the superfluid transition at , the results are converged for to .

The computational limitations of our real-time BDMFT+NCA implementation are very similar to the real-time fermionic DMFT+NCA case Eckstein and Werner (2010). Memory is *the* limiting factor when working with two time response functions, whose storage size scales quadratically with the number of time steps. The local Fock space in the bosonic case adds one or two orders of magnitude in memory usage, compared to the single-band fermionic case. A further limitation is the quadratic energy dependence of the local occupation number states , scaling with . This induces a pseudo-particle time dependence , which means that including higher occupation number states, by increasing , also requires a finer time discretization.

## Iii Results

### iii.1 Benchmark calculations

Even though BDMFT neglects spatial fluctuations, the equilibrium results for the 3D Bose-Hubbard model are in good quantitative agreement Anders *et al.* (2010, 2011) with high precision lattice QMC calculations Capogrosso-Sansone *et al.* (2007) and high-order perturbation theory Teichmann *et al.* (2009), for both the phase diagram and local correlation functions.
For example, the critical couplings at the superfluid-Mott transition are (BDMFT at ), and (lattice QMC), see Ref. Anders *et al.*, 2011 for an explicit comparison of phase diagrams.

To assess the validity of the NCA approximation we compare its superfluid phase boundary for the 3D cubic lattice with the (within BDMFT) exact CT-QMC result, see Fig. 2. It is evident that already this lowest order strong coupling expansion provides a very good approximation with (at ), as expected, considering the success of the linked cluster expansion Kauch *et al.* (2012). The simplified self-consistency based on the semi-circular DOS leads to a shift in the phase boundaries (Fig. 2) with , but we expect that the qualitative features of the solution, both in and out of equilibrium, remain unchanged.

Note that the Mott phase is only present at integer fillings. Hence, in order to study quenches between the superfluid and Mott insulator, we limit our calculations to . Strictly speaking the Mott insulator exist only at zero temperature, with a smooth crossover to the normal phase, see Fig. 2. However we follow Ref. Capogrosso-Sansone *et al.*, 2007 and define the *Mott regime* as the whole region , where the low temperature superfluid phase is absent.

For instantaneous interaction quenches the final total energy is given by the initial equilibrium total energy and an additional interaction energy contribution (due to the sudden change of from to at ). Given and the effective temperature of the system after thermalization can be determined using separate equilibrium calculations. The resulting non-equilibrium pair of a quench can be used to determine the final state after eventual thermalization by direct comparison with the equilibrium phase boundaries. This will be used throughout this study in order to produce combined equilibrium and non-equilibrium “phase diagrams”.

BDMFT captures the conversion between interaction, kinetic, and condensate energy, as well as the relaxation to the predicted thermal values (Fig. 3). Despite a nontrivial time-evolution of the individual components the total energy and the particle number is conserved to high accuracy by our 5th order solver (right panel).

We should note, however, that the NCA solution yields an approximate spectral function, as for the Fermi-Hubbard model Pruschke *et al.* (1993). To assess these errors it is useful to check the accuracy to which spectral sum rules (valid also in a nonequilibrium setting Freericks *et al.* (2013)) are fulfilled.
The moments of the spectral function , , are given by the higher order derivatives of the retarded Green’s function at , , where and are the absolute and relative time respectively, see Ref. Freericks *et al.*, 2013. The moments can also be determined using the equation of motion, in terms of operator expectation values, , , and , where denotes the th moment of the non-interacting density of states, see Ref. Anders *et al.*, 2011. For an approximate solution of the BDMFT equations, these approaches do not yield the same result.
In equilibrium, BDMFT+NCA gives a 1.6% relative error of the first spectral moment in the Mott insulating phase (, ), and a 10% error in the vicinity of the superfluid phase boundary (, ). For the second moment , the relative errors are 11% and 36%, respectively.

### iii.2 Quenches from the Mott insulator

As a first application, we study quenches within the symmetric Mott and normal phases (), i.e. suppressing symmetry breaking (superfluid states). Since the Gutzwiller mean-field description only contains the condensate energy and the interaction energy , the symmetric state in this approximation is simply the atomic limit with . So, for these quenches, no energy conversion occurs in the Gutzwiller treatment, resulting in an unphysical constant time-evolution. BDMFT, however, retains temporal fluctuations and enables the conversion of interaction energy into kinetic energy , and vice versa, which leads to a non-trivial quench dynamics.

To search for thermalization we study the relative change in the double occupancy ,

(6) |

defined so that and for a thermalized state.
Using this quantity, we identify an intermediate region of rapid thermalization in the plane [Fig. 4a] in the following way: For a given initial state we locate the maximum of as a function of at the longest accessible time , as shown explicitly for and in Fig. 4b. The values of and the corresponding effective temperatures are shown in Fig. 4a, and the result turns out to be insensitive to the initial interaction (, ). In both the weak and strong-coupling regimes the system is trapped in a long-lived “prethermalized state”, reminiscent of the relaxation dynamics in the paramagnetic Fermi-Hubbard model Eckstein *et al.* (2009).
The observed absence of thermalization in these regimes can be understood in terms of proximity to an integrable point, since the Bose-Hubbard model is integrable both for and Sorg *et al.* (2014).
Interestingly the relaxation behavior in the two regimes differ. In the strong-coupling regime the exponential decay of slows down as increases (green, red, cyan lines in Fig. 4d). In the low regime very rapidly reaches a plateau value (blue line in Fig. 4d), which increases roughly exponentially as is decreased.
The crossover between these two disparate behaviors is hard to pinpoint, and the indicated regions in Fig. 4a are only qualitative as is dependent (the region becomes narrower and shifts to slightly higher with increasing ).

### iii.3 Quenches from the superfluid

#### Relaxation regimes

The non-equilibrium dynamics after a quench from the superfluid () with weak interaction to larger interactions generates a variety of dynamical behaviors depending on . Apart from the system has two other characteristic energies (or inverse timescales), namely the bandwidth and the condensate coupling (where is our unit of energy). In general the time evolution can be separated into five regimes, see Fig. 5a and 5b.

i) For quenches deep into the Mott regime, i.e. for large , the condensate oscillates with the frequency of the final interaction strength, , while relaxing exponentially (green line). The relaxation rate strongly depends on the initial temperature . For high (as in Fig. 5b) the system displays relaxation to the Mott phase, while for low the system is trapped in a non-equilibrium superfluid state for long times, see Sec. III.3.3.

ii) In the intermediate coupling regime the interaction driven oscillations compete with the kinetic time scale and only a few oscillations can be observed, and after the condensate time scale the system displays exponential relaxation (red line).

iii) For the thermal reference state is closer to the phase boundary in the normal phase. In this regime, after an initial transient undershoot in , the system becomes trapped in a non-equilibrium superfluid state with a constant non-zero condensate (blue line).

iv) In the same range of an amplitude mode is excited at longer times (magenta line) with a roughly constant frequency but growing amplitude as the phase boundary is approached from the normal phase side.

v) For small , () the initial transient is weak as the quench-energy scales with . First the condensate undergoes a weak oscillatory transient followed by a sudden rapid growth (cyan line). This growth occurs when the final state is in the equilibrium superfluid region. The rapidly growing condensate is a numerical challenge because of the occupation of high-energy (i.e. high occupation-number) states. On the one hand, the cutoff in the bosonic Fock space must be chosen large enough to accommodate this, and on the other hand, the fast oscillations of the high-energy modes require a small time-step. In Fig. 5b We plot the results up to the point to which they can be fully converged both in the size of the time-step and the size of the Hilbert space. For and the drift in total energy and density is of the order , see Fig. 5c. Hence, we conclude that the growth is a robust feature of our BDMFT+NCA calculations.

During the growth there is a rapid conversion between the different energy components of the system, while the total energy is conserved, see Fig. 5d. The interaction energy and the normal component of the kinetic energy rapidly increase, while the condensate energy and the anomalous component of the kinetic energy decrease by the same total amount.

The self-amplified transient growth of the condensate fraction resembles the quantum turbulence driven dual cascade with non-equilibrium Bose-Einstein condensation observed in scalar field theories Berges and Sexty (2012).
Also in other contexts, dynamical instabilities with diverging solutions have been observed in lattice boson systems. In the weak coupling limit, a Gross-Pitaevskii treatment yields dynamically unstable solutions Machholm *et al.* (2003), also observed experimentally De Sarlo *et al.* (2005). In the Bose-Hubbard model the exponential condensate growth of symmetric initial states has previously been studied in mean-field Snoek (2011).

However, for the weak interaction quenches within the superfluid we cannot rule out that the sudden condensate growth is an artifact of the NCA treatment. Higher-order implementations of the self-consistent strong coupling expansion may in the future help to clarify this issue.

#### Short time dynamics after quenches deep into the Mott regime

In the limit of large final interaction , , i.e. in regime (i), the superfluid quenches from display oscillations with a frequency scaling with , , and the short time behavior is dominated by the interaction, as it defines the shortest time scale of the system. The short time relaxation is expected to be driven by local decoherence and the long time relaxation () to be dominated by hopping. In order to study the short time dynamics we perform a series of quenches to for several initial temperatures , see Fig. 6. While scales with there are important contributions from other frequency components. Pure -oscillations can only be observed in the first few revivals, while they are at later times washed out by the off-diagonal mixing of local occupation number states in the initial state. The first revival maximum occurs at the period of the final interaction , the second revival has a pronounced two peak structure with the first peak occurring at , and in the third revival the peak only appears as a shoulder of the main peak. From Fig. 6 it is also evident that the short time decoherence strongly depends on the initial temperature , with higher temperature resulting in faster damping.

An interesting question is whether the long time relaxation rate can already be inferred from the short time decoherence, in the spirit of the strong coupling analysis of Ref. Fischer and Schützhold, 2008. To investigate this we fit the simple exponential model to the real time data, where the relaxation rate is approximated using at the first revival maximum as . Figure 6 shows that the relaxation rate is overestimated for low temperature initial states and underestimated for high temperature initial states. Hence in this regime the condensate relaxation can not be inferred from the first revival maximum. Infact the long time exponential relaxation rate is only established after the characteristic condensate time scale , see e.g. the green and red lines in Fig. 5b.

We also note that the BDMFT calculation does not involve any approximation concerning the timescales of the dynamics. This sets it apart from, for example, the low frequency approximation applied in the Schwinger-Keldysh generalization of the strong coupling approach Kennett and Dalidovich (2011), where in the particle-hole symmetric limit the condensate phase and amplitude (where ) are constrained by for some constant . As shown in the lower panel of Fig. 6, the BDMFT dynamics has a non-trivial time dependence in this quantity.

#### Long time dynamics

The long time dynamics for quenches from the superfluid has been investigated in a number of zero temperature Gutzwiller mean-field studies Huber *et al.* (2007); Wolf *et al.* (2010); Sciolla and Biroli (2010, 2011); Snoek (2011); Krutitsky and Navez (2011). The most prominent nonequilibrium effect is a dynamical transition at Sciolla and Biroli (2010).
It is important to note, however, that for low only the mean-field calculation using a constrained basis including the lowest three bosonic occupation number states () produces a sharp transition. If the physically important states with higher occupations are also considered, the transition turns into a crossover ^{1}

The quantum fluctuations missing in mean-field treatments are expected to heavily modify the dynamical transition, as previously shown for other systems, using dynamical mean field theory Eckstein *et al.* (2009), the Gutzwiller approximation including Gaussian fluctuations Schiró and Fabrizio (2011); Sandri *et al.* (2012), and expansions Sciolla and Biroli (2013).
Here we show how the dynamical transition in the Bose-Hubbard model is affected when we go beyond the simple mean-field treatment, starting from a thermal initial state, and include quantum fluctuations using BDMFT.

As the hopping induced relaxation is most prominent for small and temperatures close to the superfluid phase boundary, we fix , far away from the zero temperature transition, , see Fig. 2. To see the enhanced relaxation in the vicinity of the phase boundary, located at , we consider the two initial temperatures and with relatively weak superfluidity , see Fig. 7a.

To analyze the dynamics of the condensate amplitude we first look at windowed time averages , using a Gaussian window with width to filter out oscillations. In Fig. 7b we plot the window average at the longest time as a function of , thereby restricting ourselves to the regimes (i)-(iii) above, i.e., when the final equilibrium state is not in the superfluid phase and the order parameter does not show self-amplified growth. From Fig. 7b it is evident that exhibits a crossover for , from high values at low close to , through a minimum at intermediate , and increasing again for , in qualitative agreement with the mean-field dynamical transition Sciolla and Biroli (2010). Also the general temperature dependence, namely that a higher temperature leads to lower condensate averages, agrees with mean-field. However the thermal effects in BDMFT are much stronger: both the minimum and the large plateau are drastically reduced, going from to . As we will show, this reduction is due to a rapid condensate relaxation rate emerging close to the phase boundary (Fig. 7g).

The BDMFT real-time evolution in the three regimes is shown in Fig. 7c. For small in regime (iii), stabilizes at a finite value after an initial transient, even though the thermal reference state is in the normal phase. The intermediate regime (ii) shows fast thermalization with a rapidly decaying condensate and damped collapse-and-revival oscillations, in qualitative agreement with the experimental results of Greiner et al. Greiner *et al.* (2002). Interestingly, for larger in regime (iii) the system is again trapped in a nonthermal superfluid phase, now exhibiting coherent amplitude oscillations with finite life-time.
Note that none of these relaxation and thermalization effects are captured by the Gutzwiller mean-field approximation, which predicts an oscillatory behavior (Fig. 7d).

While the minimum of indicates a crossover, the dynamical transition can be accurately located by studying the condensate phase , where . Its linear component exhibits a kink at , see Fig. 7e, in direct analogy with mean-field predictions where (mapped to a conjugate momentum ) also has a slope discontinuity at Sciolla and Biroli (2010). Similar to the fast thermalization region in the symmetric phase, the double occupancy thermalizes rapidly for as can be seen in Fig. 7f from the drastic increase in the relaxation of

as .

To get a qualitative understanding of the condensate amplitude relaxation dynamics, we fit the late time-evolution to a damped two component model

with a non-oscillatory component () and a coherent amplitude-mode (), and relaxation and damping respectively, see Fig. 7f and 7g.
The amplitude-mode frequency has the same general behavior as (not shown), and in the large limit.
Analogous to the rapid relaxation of the double occupancy, the amplitude mode is strongly damped for , but it retains a finite lifetime for large , see Fig. 7f. The relaxation of the non-oscillatory component shows two distinct behaviors: For the system is trapped in a superfluid state and the condensate relaxation is almost zero, , while for it becomes finite, reaching a maximum at intermediate , see Fig. 7g.
For large and , stays finite and the system eventually thermalizes to the Mott state, while for , becomes small as , which means that the system is trapped for a very long time in a superfluid state.
The stability of the superfluid can be understood in terms of a simple two fluid model of doublons and hard-core bosons Sorg *et al.* (2014). In this picture the quench generates long-lived doublons and depletes the hard-core boson gas away from unity filling, where it can remain a superfluid for any local interaction Schneider (2014).
This case is particularly interesting as it opens up the possibility to study the Higgs amplitude mode in a metastable superfluid.

#### Nonequilibrium “phase diagram”

We summarize the results for the long-time dynamics in the non-equilibrium “phase diagram” shown in Fig. 8. By repeating the analysis for and the series , , , , and of initial temperatures, we locate the boundaries of the three dynamical regimes in the equilibrium normal phase region, regime (i), the high- region characterized by a trapped superfluid and amplitude mode (green), regime (ii), the crossover region with rapid thermalization (red), and regime (iii), the trapped superfluid in the vicinity of the equilibrium superfluid phase (blue). While this non-equilibrium “phase diagram” depends on the initial states and the quench protocol, it gives an overview of the different relaxation and trapping phenomena and their location in parameter space. An experimental verification of these different dynamical regimes of the Bose-Hubbard model would be very interesting and presumably possible.

In the one dimensional Bose-Hubbard model a similar behavior has been theoretically observed using DMRG Kollath *et al.* (2007), and for longer times using time-dependent variational Monte Carlo Carleo *et al.* (2012). Quenches from the zero temperature superfluid display a region of thermalization at intermediate final interactions and a trapping in nonthermal states for long times at strong final coupling Kollath *et al.* (2007). At unity filling this behavior can be understood in terms of a reduced effective scattering of holon and doublon excitations at strong interactions Altman and Auerbach (2002). Our results from BDMFT indicate that this phenomenon is also relevant in three dimensions (green region in Fig. 8). However, we also identified a transient trapping at low interactions (blue region in Fig. 8) which has not been reported for 1D. It is an open question whether this feature is specific to high-dimensional models.
We also note that while BDMFT allows us to compare nonequilibrium and equilibrium states within the same formalism, the DMRG studies involve comparisons between time-dependent correlators and finite-temperature QMC results Kollath *et al.* (2007).

## Iv Conclusions

We developed the nonequilibrium BDMFT formalism and its implementation in combination with a NCA type bosonic impurity solver. We have demonstrated its ability to capture nontrivial dynamical effects in quenched Bose-Hubbard systems, including dynamical transitions, fast-thermalization crossovers, and trapped superfluid phases with long-lived but damped amplitude oscillations. These results were collected into two nonequilibrium “phase diagrams” (Figs. 4 and 8), which illustrate the transitions and crossovers that occur as one varies the quench parameters. Particularly noteworthy results are the prediction of a very long-lived transient superfluid state with an amplitude mode after quenches from the superfluid phase into the Mott regime, our finding of a trapped superfluid state after quenches into the vicinity of the superfluid phase boundary, and the nonequilibrium Bose condensation (growing condensate) after small quenches within the superfluid phase.

The ability of BDMFT to describe hopping induced relaxation phenomena at finite temperature goes beyond all current competing theoretical approaches. The Gutzwiller mean-field formalism lacks all hopping induced phenomena Sciolla and Biroli (2010, 2011), and the strong coupling based real-time approach Kennett and Dalidovich (2011) is limited to zero temperature and slow dynamics. The hopping perturbation expansion Trefzger and Sengupta (2011) looks promising but has so far only been applied at zero temperature. A comparative study with the finite temperature extension of this approach would be very interesting.

Extensions of the nonequilibrium BDMFT formalism to multi-flavor Hamiltonians Söyler *et al.* (2009); Capogrosso-Sansone *et al.* (2010a) and inhomogeneous systems (e.g. with a trapping potential) Eckstein and Werner (2013) should enable direct comparisons with cold atom experiments.
Multi-orbital effects such as virtual excitations to higher orbitals can trivially be included in BDMFT in terms of effective three body interactions Johnson *et al.* (2009).
While calculations based on unitary-time evolution Tiesinga and Johnson (2011) suffice to understand experiments Will *et al.* (2010) in the limit, BDMFT can extend the theoretical treatment to finite .

An inhomogeneous extension of BDMFT will require a more advanced parallelization scheme than that applied by Dirks et al. Dirks *et al.* (2014) for the fermionic case, but it would enable studies of very important phenomena, such as mass transport and the effects of the trapping potential in general cold-atom systems out-of-equilibrium.
The big challenge in these systems is the inherent disparity of the hopping and mass transport time-scales.
Simpler approximations such as the hopping expansion has successfully been applied in this context Dutta *et al.* (2014),
but without incorporating thermal and retarded correlation effects.

###### Acknowledgements.

The authors would like to acknowledge fruitful discussions with H. Aoki, T. Ayral, J. Berges, N. Buchheim, D. Golez, F. Heidrich-Meisner, A. Herrmann, S. Hild, D. Hügel, M. P. Kennett, Y. Murakami, L. Pollet, U. Schneider, N. Tsuji, and L. Vidmar. The calculations have been performed on the UniFr cluster. HS and PW are supported by FP7/ERC starting grant No. 278023.## Appendix A Nonequilibrium bosonic dynamical mean field theory

The bosonic dynamical mean field theory (BDMFT) for the Bose-Hubbard model in *equilibrium* has been derived in Ref. Anders *et al.* (2011) in three alternative ways, using the kinetic energy functional, an effective medium approach, and the quantum cavity method, which is very similar to the cavity construction by Snoek and Hofstetter Snoek and Hofstetter (2010).
In this appendix we follow the latter approach, which combines the cavity construction with a generating functional formalism and a cumulant expansion to second order. By performing the derivation on the three branch Kadanoff-Baym contour we obtain the *nonequilibrium* generalization of BDMFT.
We will also show that BDMFT corresponds to the first order correction in the inverse coordination number , and the second order correction in the fluctuations, of the mean-field approximation for the Bose-Hubbard model.

### a.1 The Bose-Hubbard model

We consider the Bose-Hubbard model [Eq. (1)]

on a lattice with nearest neighbor hopping and a local pair interaction , where is a pure two particle interaction counting the number of pairs on site , and denotes the sum over all nearest neighbor pairs and . Using the Nambu-spinor notation and collecting the local terms on site into , the Hamiltonian can be expressed as

(7) |

where we have used that if . Note that in the Nambu notation is hermitian, i.e.

(8) |

### a.2 Kadanoff-Baym and Nambu formalism

To treat an arbitrary time evolution starting from a finite temperature equilibrium state we formulate the theory on the three-branch Kadanoff-Baym contour () Aoki *et al.* (2014); Stefanucci and van
Leeuwen (2013). The partition function of the initial state can be expressed as where is the action defined on the contour , , is the time-ordering operator on , and the trace runs over the Hilbert space of .
Time-dependent operator expectation values can be expressed in the trace formalism as

(9) |

and the single-particle Green’s function on the contour, , is given by . The Nambu generalization of the single-particle Green’s function is a matrix, which can be expressed in spinor notation as

(10) |

For a general introduction to the Kadanoff-Baym contour formalism, see Ref. Stefanucci and van
Leeuwen (2013) and for a DMFT specific introduction see Ref. Aoki *et al.* (2014).

### a.3 Real-time generating functional

To construct the generating functional on the contour we introduce source fields on each site and the source action

(11) |

Using and the action of the system

(12) |

the generating functional can be defined as

(13) |

It can be used to compute any *connected* response function by taking derivatives with respect to the source fields

(14) |

### a.4 Cavity construction

To derive a local effective action we use the standard cavity construction Georges *et al.* (1996) and separate the Hamiltonian into three parts,

(15) |

where acts on the site , connects the zeroth site to its neighbors, and is the lattice with a cavity at the zeroth site, i.e.

(16) | ||||

(17) | ||||

(18) |

which in turn separates the action into

(19) |

Analogously the source term can be decomposed into

(20) |

according to the same protocol, with

(21) |

which yields the corresponding terms of the source action

(22) |

Using this separation of the zeroth site’s degrees of freedom the generating functional can be written as

(23) |

where denotes the trace over the Fock-space of the zeroth site. In this form the generating functional can be approximated and/or taken to e.g. the infinite connectivity limit, which results in different types of dynamical mean field theory (DMFT) approximations.

### a.5 Cumulant expansion

We are now ready to perform a *cumulant expansion* Kubo (1962) of the expectation value in Eq. (23). Formally this corresponds to expanding in an infinite sum of response functions with respect to . The initial logarithm ensures that the series enters in the exponent, and for this reason the procedure is often referred to as “re-exponentiation”.

Following Ref. Kubo, 1962 the cumulant expansion becomes

In the derivation of the fermionic dynamical mean field effective action, the cumulant expansion terminates at second order in the limit of infinite dimensions (using a scaling of the hopping) Georges *et al.* (1996). This yields the usual hybridization function term

For Bosons, however, anomalous contributions due to symmetry breaking scale linearly with the coordination number , requiring a scaling of the hopping to obtain a finite limit Anders *et al.* (2010, 2011). This procedure results in the mean field effective action Sachdev (1999) which does not include quantum fluctuations of non-condensed Bosons.
In order to retain fluctuations we therefore avoid taking the infinite connectivity limit and instead truncate the cumulant expansion at second order, which (as we will see) yields corrections in the effective action Snoek and Hofstetter (2010).

We write the second order approximation of the cumulant expansion as

(24) |

collecting the source-free terms in the effective action

(25) |

and the terms containing source fields in the effective *source* action