Nonequilibrium dynamical mean-field theory for bosonic lattice models

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.


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 [1,2]. They enable an almost ideal realization of the lowenergy effective Hamiltonians (the fermionic and bosonic Hubbard models [3,4]) 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 [5]. For bosonic atoms, the Mott insulating and superfluid regime can easily be accessed [1] and the experimental control is so precise that the use of cold atoms as "quantum simulators" becomes a realistic option [6] (for a recent review see Ref. 7).
A particularly interesting aspect of cold atom experiments is the possibility to study the time-evolution of interacting many-body systems [8][9][10][11][12][13][14][15][16]. This was beautifully demonstrated in the seminal work by Greiner et al. [8], 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 [17] can be computed accurately using Monte Carlo simulations [18], the real-time evolution of interacting bosonic lattice systems is a big computational challenge.
In one dimension, density matrix renormalization group (DMRG) methods [19] 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 [15]. Still, DMRG calculations have provided important insights into the short time dynamics, as measured in 1D optical lattices [14][15][16]. Kollath et al. [20] used non-local correlators to study relaxation and ther-malization. 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 [21,22] and debated [23,24] 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 [25]. It has also been applied to 2D systems and is not inherently limited to any dimensionality [26]. 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 [27].
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 [16], out-of-equilibrium phenomena in the three dimensional Bose-Hubbard model remain largely unexplored [14][15][16] 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 collapseand-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 [8,10].
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 ther-modynamical limit, and thus enables the study of relaxation and thermalization phenomena in infinite systems [38]. The equilibrium version of BDMFT [39][40][41][42] produces phase diagrams, condensate fractions, and correlation functions with remarkable accuracy [42]. While the extension of this formalism to nonequilibrium situations is analogous to the fermionic case [38], 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. 41 cannot easily be applied to nonequilibrium problems, because of a dynamical sign problem [43], while exact diagonalization based solvers are even more limited than in the fermionic case [44], 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) [45]. 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 A], the NCA impurity solver [Sec. II B], the energy calculations [Sec.
II C], and our numerical implementation [Sec. II D]. In Sec. III we first present benchmark calculations showing density and energy conservation and discuss the lowest order spectral moments [Sec. III A]. The dynamical regimes in the normal phase are mapped out in Sec. III B. In Sec. III C we consider superfluid initial states, and after an overview of the relaxation regimes in Sec. III C 1, we study the dynamics for short times in Sec. III C 2, and long times in Sec. III C 3. The findings are summarized in Sec. III C 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 [4,5] where b † i (b i ) andn i are the bosonic creation (annihilation) and number operators acting on site i, µ is the chemical potential, and U the local pair interaction which competes with the nearest neighbor hopping J that we take as our unit of energy.

A. Nonequilibrium bosonic dynamical mean-field theory
By extending the equilibrium bosonic dynamical mean field theory (BDMFT) [41,42] to the three-branch Kadanoff-Baym contour C (0 → t max → 0 → −iβ) [38,46], we obtain the bosonic impurity action , ∆(t, t ) the hybridization function, and Φ † eff the effective symmetry breaking field, which is defined in terms of ∆, the local condensate fraction Φ † = (φ * , φ), and the lattice coordination number z as 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 [47], which would be obtained by taking the infinite dimensional limit z → ∞ at fixed zJ (or analogously ∆ → 0) [42]. The solution of the impurity model yields the connected impurity Green's function where T C is the time-ordering operator on the contour C, and the local condensate fraction is The BDMFT self-consistency loop is closed by computing the lattice Green's function G k from G, and then expressing the hybridization function ∆ in terms of the local lattice Green's function G L = 1 N k k G k (at self consistency G L = G) [38]. In the present study we employ the simplified self-consistency relation and set z = 6, which corresponds to a non-interacting semi-circular density of states (DOS) with the same bandwidth W = 12J and lattice coordination number z as the 3D cubic lattice with nearest neighbor hopping J.
NCA diagram representations of a) the pseudoparticle self-energyΣ, and b) the single-particle Green's function Gγν .

B. Non-crossing approximation impurity solver
The previous BDMFT equilibrium studies employed a hybridization expansion CT-QMC impurity solver [41,42]. 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 [48] 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. [49]. 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 p Γ and p † Γ for each local occupation number manybody state |Γ . This maps the local Hamiltonian to a quadratic term ΓΓ Ĥ (t) ΓΓ p † Γ p Γ , while the hybridization ∆ turns into a pseudo-particle interaction. Expanding to first order in ∆ gives the NCA of Ref. 49 generalized to Nambu formalism.
Within NCA,Ĝ andΣ are calculated self-consistently, and local observables are determined from the reduced local density matrixρ(t) = iĜ < (t, t), yielding the local condensate as while response functions must be determined diagrammatically. In particular, the connected single-particle impurity Green's function G is obtained from the bubble diagram without hybridization insertions (see Fig. 1 and

C. Total energy components
The total energy E t of the system, is the sum of the (connected) kinetic energy E k , the condensate energy (or disconnected kinetic energy) E c , and the local interaction and E i can be written in terms of n 2 (t) = Tr Γ [n 2ρ (t)] and n (t) = Tr Γ [nρ(t)] as E i (t) = U (t)( n 2 (t) − n (t))/2.

D. Numerical implementation
We solve the pseudo-particle Dyson equation using a fifth order multi-step method [49,50] 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 [49] (the gauge property µ → µ + δµ(t) ⇒ b → be −i t 0 dt δµ(t) ensures ∂ t n = 0). In principle the local Fock space is unbounded, but for U > 0 it can safely be truncated, keeping only N max states. The cut-off error is controlled by monitoring the drift in Tr Γ [ρ] away from unity. Close to the n = 1 superfluid transition at U J, the results are converged for N max = 5 to 11.
The computational limitations of our real-time BDMFT+NCA implementation are very similar to the real-time fermionic DMFT+NCA case [49]. Memory is the limiting factor when working with two time response (color online) BDMFT superfluid phase boundary for CT-QMC [42] (blue) and NCA (red) on the 3D cubic lattice, NCA with a semi circular DOS (green), and meanfield theory (black). Panel a) shows the (J/U , µ/U ) plane at T = 1.5, and panel b) the (U , T ) plane for n = 1.
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 E Γ ∼ Γ|Un 2 |Γ = U n 2 Γ . This induces a pseudo-particle time dependenceĜ ΓΓ (t, t ) ∼ e −iU n 2 Γ (t−t ) , which means that including higher occupation number states, by increasing N max , also requires a finer time discretization.

A. Benchmark calculations
Even though BDMFT neglects spatial fluctuations, the equilibrium results for the 3D Bose-Hubbard model are in good quantitative agreement [41,42] with high precision lattice QMC calculations [17] and high-order perturbation theory [51], for both the phase diagram and local correlation functions. 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 (J/U ) c ≈ 0.0340 (at T = 1.5), as expected, considering the success of the linked cluster expansion [52]. The simplified self-consistency based on the semicircular DOS leads to a shift in the phase boundaries ( Fig. 2) with (J/U ) c ≈ 0.0378, 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 n = 1. 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. 17 and define the Mott regime as the whole region U > U c (T = 0), where the low temperature superfluid phase is absent.
For instantaneous interaction quenches the final total energy E (f ) t is given by the initial equilibrium total energy E (i) t and an additional interaction energy contribution E and U f the effective temperature T eff of the system after thermalization can be determined using separate equilibrium calculations. The resulting non-equilibrium (U f , T eff ) pair of a quench can be used to determine the final state after eventual thermalization by direct comparison with the equilibrium (U, T ) phase boundaries. This will be used throughout this study in order to produce combined equilibrium (U, T ) and non-equilibrium (U f , T eff ) "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 E t and the particle number n 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 [53]. To assess these errors it is useful to check the accuracy to which spectral sum rules (valid also in a nonequilibrium setting [54]) are fulfilled. The moments µ R n (T ) of the spectral function A R (T, ω), , are given by the higher order derivatives of the retarded Green's function where T and t are the absolute and relative time respectively, see Ref. 54. The moments can also be determined using the equation of motion, in terms of operator expectation values, , where n denotes the nth moment of the non-interacting density of states, see Ref. 42. 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 µ R 1 in the Mott insulating phase (U = 96, T = 6), and a 10% error in the vicinity of the superfluid phase boundary (U = 40, T = 6). For the second moment µ R 2 , the relative errors are 11% and 36%, respectively.

B. Quenches from the Mott insulator
As a first application, we study quenches within the symmetric Mott and normal phases (|φ| = 0), i.e. suppressing symmetry breaking (superfluid states). Since the Gutzwiller mean-field description only contains the condensate energy E c and the interaction energy E i , the symmetric state in this approximation is simply the atomic limit with E c ∝ |φ| 2 = 0. 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 E i into kinetic energy E k , and vice versa, which leads to a non-trivial quench dynamics.
To search for thermalization we study the relative change κ in the double occupancy n 2 , defined so that κ(t = 0) = 0 and κ = 1 for a thermalized state. Using this quantity, we identify an intermediate region of rapid thermalization in the (U, T ) plane [ Fig. 4a] in the following way: For a given initial state (U i , T i ) we locate the maximum of |1 − κ(t = t max )| −1 as a function of U f at the longest accessible time t max = 2.66, as shown explicitly for U i = 30 and T i = 6 in Fig. 4b. The values of U f and the corresponding effective temperatures T eff are shown in Fig. 4a, and the result turns out to be insensitive to the initial interaction (U i = 30, 45). In both the weak and strong-coupling U f regimes the system is trapped in a long-lived "prethermalized state", reminiscent of the relaxation dynamics in the paramagnetic Fermi-Hubbard model [55]. 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 U = 0 and U = ∞ [22]. Interestingly the relaxation behavior in the two regimes differ. In the strong-coupling regime the exponential decay of |1 − κ| slows down as U f increases (green, red, cyan lines in Fig. 4d). In the low U f regime κ very rapidly reaches a plateau value (blue line in Fig. 4d), which increases roughly exponentially as U f is decreased. The crossover between these two disparate behaviors is hard to pinpoint, and the indicated regions in Fig. 4a are only qualitative as |1 − κ(t = t max )| −1 is t max dependent (the region becomes narrower and shifts to slightly higher U f with increasing t max ).
C. Quenches from the superfluid

Relaxation regimes
The non-equilibrium dynamics after a quench from the superfluid (|φ| > 0) with weak interaction U i = 6 to larger interactions U f > U i generates a variety of dynamical behaviors depending on U f . Apart from U f the system has two other characteristic energies (or inverse timescales), namely the bandwidth W = 12J and the condensate coupling zJ = 6J (where J = 1 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 U f W, zJ, the condensate oscillates with the frequency ω of the final interaction strength, ω ≈ U f , while relaxing exponentially (green line). The relaxation rate strongly depends on the initial temperature T i . For high T i (as in Fig. 5b) the system displays relaxation to the Mott phase, while for low T i the system is trapped in a non-equilibrium superfluid state for long times, see Sec. III C 3.
ii) In the intermediate coupling regime U f W > zJ the interaction driven oscillations compete with the kinetic time scale and only a few oscillations can be observed, and after the condensate time scale 2π/(zJ) the system displays exponential relaxation (red line).
iii) For W U f > zJ 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 U f 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 ∆U = U f − U i , (W > U f ≈ zJ) the initial transient is weak as the quench-energy scales with ∆U . 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 equi-  For U f = 6.6 the (equilibrium) thermal reference state at (U f , T eff ) is also superfluid and the non-equilibrium dynamics displays a rapid transient growth of the condensate (cyan). Close to the phase boundary in the normal phase the system is trapped in a superfluid for long times and the quench induced excitation is transferred to an amplitude mode at longer times (magenta). For intermediate U f a constant trapped superfluid persists (blue) and after passing the dynamical transition U dyn c the system undergoes exponential relaxation to the normal phase (red). For large final interactions U f W, zJ the condensate displays "collapse-and-revival" oscillations with frequency ω ≈ U f (green). The evolution from U f = 6.6 to U f = 11.4 is also shown (gray lines). For the growing condensate at U f = 6.6 panel c) shows the accuracy in total energy Et and density n and panel d) shows the energy conversion during the time of rapid condensate growth. librium 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 N max = 11 and ∆t = 0.005 the drift in total energy ∆E t = |E t (t) − E t (0)| and density ∆ n = | n(t) − n(0) | is of the order 10 −6 , 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 E t is conserved, see Fig. 5d. The interaction energy E i and the normal component of the kinetic energy E 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 [56]. 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 [57], also observed experimentally [58]. In the Bose-Hubbard model the exponential condensate growth of symmetric initial states has previously been studied in mean-field [36].
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 U f W , zJ, i.e. in regime (i), the superfluid quenches from U i = 6 display oscillations with a frequency ω scaling with U f , ω ≈ U f , 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 (t > 2π/W ≈ 0.52) to be dominated by hopping. In order to study the short time dynamics we perform a series of quenches to U f = 48 for several initial temperatures T i = 3.00, . . . , 5.25, see Fig. 6. While ω scales with U f there are important contributions from other frequency components. Pure 2π/U f -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 oc-cupation number states in the initial state. The first revival maximum occurs at the period of the final interaction 2π/U f , the second revival has a pronounced two peak structure with the first peak occurring at 2 · 2π/U f , and in the third revival the 3 · 2π/U f 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 T i , 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. 59. To investigate this we fit the simple exponential model |φ|(t) ≈ |φ|(0)e −t/τ to the real time data, where the relaxation rate τ is approximated using |φ| at the first revival maximum t r = 2π/U f as τ = t r / log(|φ|(t r )/|φ|(0)). 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 2π/(zJ), 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 [29], where in the particle-hole symmetric limit n = 1 the condensate phase θ and amplitude |φ| (where φ = |φ|e iθ ) are constrained by (∂ t θ)|φ| 2 = C for some constant C. 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 [32][33][34][35][36][37]. The most prominent nonequilibrium effect is a dynamical transi- [34]. It is important to note, however, that for low U i only the mean-field calculation using a constrained basis including the lowest three bosonic occupation number states (N max = 3) produces a sharp transition. If the physically important states with higher occupations are also considered, the transition turns into a crossover [60]. In a broader perspective the occurrence of a dynamical transition is not specific to the Bose-Hubbard model, but has also been observed in the Fermi-Hubbard model and other systems on the mean-field level [35,61,62].
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 [55], the Gutzwiller approxi-mation including Gaussian fluctuations [63,64], and 1/N expansions [65]. 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 U i and temperatures T i close to the superfluid phase boundary, we fix U i = 6, far away from the zero temperature transition, U c (T = 0) ≈ 26.4, see Fig. 2. To see the enhanced relaxation in the vicinity of the phase boundary, located at T c (U = 6) ≈ 5.49, we consider the two initial temperatures T i = 4.5 and 5.1 with relatively weak superfluidity |φ| 2 0.5, see Fig. 7a.
To analyze the dynamics of the condensate amplitude |φ| we first look at windowed time averages| φ|(t), using a Gaussian window with width t w = 1/3 to filter out oscillations. In Fig. 7b we plot the window aver-age| φ|(t = t max ) at the longest time as a function of U f , 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| φ|(t = t max ) exhibits a crossover for T i = 4.5, from high values at low U f close to U i = 6, through a minimum at intermediate U f , and increasing again for U f 30, in qualitative agreement with the mean-field dynamical transition [34]. 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 U f plateau are drastically reduced, going from T i = 4.5 to T i = 5.1. As we will show, this reduction is due to a rapid condensate relaxation rate τ −1 |φ|c emerging close to the phase boundary (Fig. 7g).
The BDMFT real-time evolution in the three regimes is shown in Fig. 7c. For small U f 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 collapseand-revival oscillations, in qualitative agreement with the experimental results of Greiner et al. [8]. Interestingly, for larger U f 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| φ|(t = t max ) indicates a crossover, the dynamical transition U dyn c can be accurately located by studying the condensate phase θ(t), where φ(t) = |φ|(t) e iθ(t) . Its linear component ∂ t θ(t) ≈ ω θ exhibits a kink at U f = U dyn c , see Fig. 7e, in direct analogy with mean-field predictions where θ (mapped to a conjugate momentum p) also has a slope discontinuity at U dyn c [34]. Similar to the fast thermalization region in the symmetric phase, the double occupancy thermalizes rapidly for U f ≈ (U dyn c ) + as can be seen in Fig. 7f from the drastic increase in the relaxation τ −1 To get a qualitative understanding of the condensate amplitude |φ| relaxation dynamics, we fit the late timeevolution (t > 1.3) to a damped two component model with a non-oscillatory component (|φ| c ) and a coherent amplitude-mode (|φ| AM ), and relaxation τ −1 |φ|c and damping τ −1 |φ| AM respectively, see Fig. 7f and 7g. The amplitude-mode frequency ω has the same general behavior as ω θ (not shown), and ω, ω θ → U f in the large U f limit. Analogous to the rapid relaxation of the double occupancy, the amplitude mode is strongly damped for U f ≈ (U dyn c ) + , but it retains a finite lifetime τ −1 |φ| AM > 0 for large U f , see Fig. 7f. The relaxation of the nonoscillatory component shows two distinct behaviors: For U f < U dyn c the system is trapped in a superfluid state and the condensate relaxation is almost zero, τ −1 |φ|c ≈ 0, while for U f U dyn c it becomes finite, reaching a maximum at intermediate U f , see Fig. 7g. For large U f and T i = 5.1, τ −1 |φ|c stays finite and the system eventually thermalizes to the Mott state, while for T i = 4.5, τ −1 |φ|c becomes small as U f → ∞, which means that the system is trapped for  [22]. 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 [66]. 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 U i = 6 and the series T i = 3.50, 4.00, 4.50, 5.10, and 5.25 of initial temperatures, we locate the boundaries of the three dynamical regimes in the equilibrium normal phase region, regime (i), the high-U 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 nonequilibrium "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 [20], and for longer times using time-dependent variational Monte Carlo [25]. 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 [20]. At unity filling this behavior can be understood in terms of a reduced effective scattering of holon and doublon excitations at strong interactions [67]. 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 [20].

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 [34,35], and the strong coupling based real-time approach [29] is limited to zero temperature and slow dynamics. The hopping perturbation expansion [28] 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 [68,69] and inhomogeneous systems (e.g. with a trapping potential) [70] 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 [71]. While calculations based on unitary-time evolution [72] suffice to understand experiments [10] in the J → 0 limit, BDMFT can extend the theoretical treatment to finite J.
An inhomogeneous extension of BDMFT will require a more advanced parallelization scheme than that applied by Dirks et al. [73] 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 [31], but without incorporating thermal and retarded correlation effects.
In a broader perspective it should be productive to apply nonequilibrium BDMFT or variants of this formalism to nonequilibrium Bose-condensation in, e.g., polaritonic systems and field theories [56,74,75].
on a lattice with nearest neighbor hopping J and a local pair interaction U , is a pure two particle interaction counting the number of pairs on site i, and i, j denotes the sum over all nearest neighbor pairs i and j. Using the Nambu-spinor notation b † = (b † , b) and collecting the local terms on site i into H i = Un i (n i −1)/2−µn i , the Hamiltonian H can be expressed as where we have used that

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 C (0 → t max → 0 → −iβ) [38,46]. The partition function Z of the initial state can be expressed as Z = Tr[T C e −iS ] where S is the action defined on the contour C, S = C dz H(z), T C is the time-ordering operator on C, and the trace Tr[·] runs over the Hilbert space of H. Time-dependent operator expectation values can be expressed in the trace formalism as and the single-particle Green's function on the contour, G(t, t ), is given by Nambu generalization of the single-particle Green's function is a 2 × 2 matrix, which can be expressed in spinor notation as For a general introduction to the Kadanoff-Baym contour formalism, see Ref. [46] and for a DMFT specific introduction see Ref. [38].

Real-time generating functional
To construct the generating functional on the contour C we introduce source fields η i on each site i and the source action Using S η and the action S of the system the generating functional Z[η] can be defined as It can be used to compute any connected response function by taking derivatives with respect to the source fields

Cavity construction
To derive a local effective action we use the standard cavity construction [77] and separate the Hamiltonian into three parts, where H 0 acts on the site i = 0, ∆H connects the zeroth site to its neighbors, and H (0) is the lattice with a cavity at the zeroth site, i.e.
which in turn separates the action S into Analogously the source term can be decomposed into according to the same protocol, with which yields the corresponding terms of the source action Using this separation of the zeroth site's degrees of freedom the generating functional can be written as where Tr 0 [·] 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.

Cumulant expansion
We are now ready to perform a cumulant expansion [78] of the expectation value e −i∆S+S (0) η S (0) in Eq. (A17). Formally this corresponds to expanding ln e −i∆S+S (0) η S (0) in an infinite sum of response functions with respect to S (0) . 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. 78 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 z → ∞ (using a J → J/ √ z scaling of the hopping) [77]. This yields the usual hybridization function term For Bosons, however, anomalous contributions due to symmetry breaking scale linearly with the coordination number z, requiring a 1/z scaling of the hopping to obtain a finite z → ∞ limit [41,42]. This procedure results in the mean field effective action [47] 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 1/z corrections in the effective action [76].
We write the second order approximation of the cumulant expansion as collecting the source-free terms in the effective action and the terms containing source fields η in the effective source action Hence, by truncating the expansion at second order we obtain an effective action S eff and generating functional Z eff [η] according to

Explicit 2nd order form
To obtain the explicit form for the local effective action we need to look into the details of the expansion giving the actions S (0) eff and S (0) eff,η . At a later stage, we will also make use of the effective generating functional in order to arrive at the contour generalization of the (selfconsistent) B-DMFT effective action, previously derived for equilibrium in [42,79].
The operators appearing in the expansion of S where in the last steps we have used the hermitian property of Nambu creation-annihilation operator products [Eq. (A2)]. Hence the first order expectation values take the form The second order terms can be obtained using the two different ways of expressing the operators in Eqs. (A23) and (A24) in order to arrive at Nambu response function expressions as in Eq. (A4). The second order term of S A27) where we have introduced the connected single particle Green's function G S (0) of the lattice with cavity and the total hybridization function ∆ of the zeroth lattice site Hence, the source-free action S (0) eff can be written as a. Local effective source action Next we consider the second order terms of the effective source action S and the first mixed term becomes By an interchange of integration variables it is possible to show that the other mixed term gives an equal contribution.
Collecting all the terms we arrive at the final expression for the local effective source action

Local anomalous term
We see in Eq. (A29) that the symmetry breaking of the infinite lattice system induces a local symmetry breaking term on the zeroth lattice site. The strength of the symmetry breaking field is however determined by the anomalous expectation values b † i (z) For finite coordination numbers z the removal of the cavity site affects the neighboring sites, hence this expectation value is not equal to that of the original homogeneous system [76] b † i (t) (c) S eff we calculate the latter using the effective generating functional Z eff and the Nambu generalization of Thus the local anomalous term of S (0) eff in Eq. (A29) can be rewritten, using only expectation values with respect to S eff , as where in the last step, we have introduced the local anomalous amplitude Φ †

Local effective action
Substituting Eqs. (A29) and (A32) into Eq. (A22), rewriting the local symmetry breaking using Eq. (A35) and setting the sources to zero η = 0, we obtain the BDMFT local effective action for the Bose-Hubbard model where we have dropped site indices and the complex field a. Equilibrium form On the imaginary time branch the field is constant Φ † (τ ) = Φ † , the hybridization function is time translational invariant ∆(τ, τ ) = ∆(τ − τ ), and the action simplifies to in agreement with Refs. 42 and 79, up to a minus sign on the hybridization function due to a different notation.

One-loop correction in 1/z
To see that the BDMFT effective action in Eq. (A36) is a one-loop correction in the inverse coordination number 1/z one must study the scaling of its terms. The only non-trivial contribution comes from the hybridization function ∆ = J 2 0,i , 0,j G  ii (t, t ) . (A38) On more general lattices, the power counting in z gives the same leading order result, but is more elaborate [77]. Substituting this result into Eq. (A36) and making a J → J/z rescaling of the hopping gives the rescaled actioñ S eff = C dt −µn(t) + U 2 n(t)(n(t) − 1) which makes it evident that the terms containing G ii , i.e. the hybridization terms in the BDMFT effective action [Eq. (A36)] corresponds to a 1/z correction of the mean-field action (obtained by setting ∆ = 0).

Second order fluctuation expansion
While BDMFT can bee seen as a one-loop expansion in the inverse coordination number it is also a second order expansion in the condensate fluctuations, as discussed in Ref. 42. This can be made explicit by rewriting the terms containing the hybridization in the effective action using the fluctuation operators δb, defined as δb ≡ b − Φ with δb = 0. Inserting these in Eq. (A36) yields S eff = C dt −µn(t) + U 2 n(t)(n(t) − 1) − zJΦ † (t)b(t) where the hybridization term is the exact 2nd order contribution of the fluctuations. Hence BDMFT correctly describes the deep superfluid where fluctuations are suppressed, i.e. the weakly interacting Bose gas limit (WIBG) [80].
a vertex creates a "Nambu-particle" by insertion of b † γ giving the matrix element Γ|b † γ |Γ , and an interaction line leaving the vertex annihilates a "Nambu-particle" through b γ giving Γ|b γ |Γ . In the following we will use the operator symbols b † γ and b γ to represent these matrix elements as they act in the same Fock-space as the pseudo-particle propagatorĜ and the pseudo-particle self-energyΣ.

Pseudo particle self-energy
Following the diagram rules of Ref. 49 the pseudoparticle self-energyΣ at first order in ∆, corresponding to the non-crossing approximation (NCA), takes the form of shell-diagramŝ where ξ = ±1 for bosons and fermions respectively. Using the Nambu generalization of propagators and vertices gives the contour expression for the first diagram with a forward propagating hybridization line according tô with implicit matrix multiplications in the many-body state indices Λ and Λ . The second diagram is constructed analogouslŷ Suppressing many-body state indices in these diagrams yields Fig. 1a in Sec. II B. Collecting all terms, we obtain which corresponds to Eq. (4) in Sec. II B.
To perform actual calculations we work with a subset of Keldysh components [38], namely, the imaginary time Matsubara componentΣ M (τ ), the real-time greater com-ponentΣ > (t, t ), the real-time lesser componentΣ < (t, t ), and the right-mixing componentΣ (t, τ ). These components can be derived from the general contour expression forΣ using the Langreth product rules [46,81]. This is because the pseudo-particle self energyΣ is given by contour time products of the hybridization function ∆ and the pseudo-particle propagatorĜ,Σ ∝ ∆Ĝ. Note that the two contributionsΣ (1) andΣ (2) [Eqs. (B1) and (B1)] must be treated differently as the order of the time arguments in ∆ differs.

Single particle Green's function
The NCA approximation for the single particle Green's function is given by the pseudo-particle bubble equipped with two vertices [49].
The Nambu generalization amounts to adding Nambu indices to the vertices and gives the non-connected Green's function as where the hybridization line stubs denote the insertion of a Nambu creation or annihilation operator, and the trace corresponds to the summation over the Γ manybody state index. To obtain the connected Green's function G fromG, the symmetry broken contribution must be removed, i.e. G(t, t ) =G(t, t ) + iΦ(t)Φ † (t ), which corresponds to Eq. (5) and Fig. 1b in Sec. II B.