Contour-time approach to the Bose-Hubbard model in the strong coupling regime: Studying two-point spatio-temporal correlations at the Hartree-Fock-Bogoliubov level

We develop a formalism that allows the study of correlations in space and time in both the superfluid and Mott insulating phases of the Bose-Hubbard Model. Specifically, we obtain a two particle irreducible effective action within the contour-time formalism that allows for both equilibrium and out of equilibrium phenomena. We derive equations of motion for both the superfluid order parameter and two-point correlation functions. To assess the accuracy of this formalism, we study the equilibrium solution of the equations of motion and compare our results to existing strong coupling methods as well as exact methods where possible. We discuss applications of this formalism to out of equilibrium situations.


Introduction
The out of equilibrium dynamics of cold atoms trapped in optical lattices has received considerable attention in recent years [1][2][3][4][5][6]. The ability to tune experimental parameters over a wide range of values in real time makes these systems very versatile and gives the opportunity to study quantum systems out of equilibrium in a controlled fashion. Quantum quenches, in which parameters in the Hamiltonian are varied in time faster than the system can respond adiabatically, e.g. when a system is driven through a quantum critical point, are a protocol that is natural to study in this context and have been studied intensely both theoretically and experimentally.
The Bose-Hubbard model (BHM) [7] has been shown to describe interacting ultracold bosons in an optical lattice [8], allowing the opportunity for experiments to probe the out of equilibrium dynamics of the model [8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25]. The BHM is a particularly convenient context for studying quantum quenches as it displays a quantum phase transition between the superfluid and Mott-insulator phases (or vice versa) as the ratio of intersite hopping J to the on-site repulsion U is varied, as observed by Greiner et al. [9]. Theoretical studies of the BHM suggest that whether equilibration occurs or not after a quantum quench depends sensitively on the initial and final values of J/U and the chemical potential [26][27][28][29][30][31][32][33]. In the case of quenches from superfluid (large J/U ) to Mott insulator (small J/U ) there have been suggestions that there may be aging behaviour and glassiness that might be experimentally observable in two time correlations or in violations of the fluctuation dissipation theorem [6, 26-28, 31, 33]. In the alternative quench from Mott insulator to superfluid, it has been suggested that Kibble-Zurek [34][35][36] scaling of defects should be observed [37,38], which has recently been tested experimentally [10].
In experiments, the combination of a harmonic trap and small J/U leads to a wedding cake structure of the equilibrium density, with alternating Mott insulating and superfluid regions [39,40]. The presence of Mott insulating regions has been predicted to retard relaxation to equilibrium after a quench to small J/U by impeding mass transport of bosons through these regions [41,42] which has also been observed experimentally [43]. This gives a picture in which relaxation after a quench takes place in two steps -fast relaxation to local equilibrium followed by slower relaxation via mass transport [41,44].
In addition to slow dynamics, several analytical and numerical studies have also shown a Lieb-Robinsonlike [45] bound of a maximal velocity which leads to a light-cone like spreading of density correlations in one dimensional systems for quenches from the superfluid to Mott-insulating regime as well as quenches within the superfluid [46] or Mott-insulating phases [29,42,47,48]. The latter case was recently observed experimentally by Cheneau et al. [49]. Similar predictions have been made for higher dimensional systems [46,50,51]. The results summarized above motivate the study of the temporal and spatial correlations of the BHM after a quantum quench in order to elucidate the dynamics observed after quenches.
A generic problem in the theoretical description of quantum quenches is that it is necessary to have a formalism that is able to describe the physics in the phases on both sides of a quantum critical point. In the case of the Bose Hubbard model, numerical approaches such as exact diagonalization and the timedependent density matrix renormalization group (t-DMRG) [24,26,42,47,49,52,53] can be essentially exact in all parts of parameter space but are limited by system size and usually are practical only in one dimension. For dimensions higher than one, methods such as time-dependent Gutzwiller mean field theory [4,41,54,55] and dynamical mean field theory [32] have been used which can capture the presence of a quantum phase transition, but in their simplest form do not capture spatial correlations, although there has been work on including perturbative corrections [50,[56][57][58][59][60][61]. An analytical approach based on using two Hubbard Stratonovich transformations to capture both weak-coupling and strong-coupling physics in the same formalism was developed by Sengupta and Dupuis [62]. Within their effective theory, they performed a mean-field calculation of the superfluid order parameter and a Bogoliubov (1-loop) approximation to the two-point Green's function to study the excitation spectrum. Their work was generalized by one of us from an equilibrium theory to out of equilibrium by using the Schwinger-Keldysh formalism to obtain a one-particle irreducible (1PI) effective action which was then used to study the superfluid order parameter after a quench [31].
Here, we extend the approach developed in Ref. [31] to obtain a two-particle irreducible (2PI) effective action using the contour-time formalism, which is a generalisation of the Schwinger-Keldysh formalism. In the 2PI approach, the evolution of the order parameter and the two-point Green's functions are treated on the same footing [63] which allows us to describe correlations both in the broken symmetry (superfluid) phase and the Mott phase. Moreover, the method provides a systematic way to go beyond the mean-field or the 1-loop approximation. We obtain two main results. First, we develop the 2PI strong coupling formalism for the BHM. Second, we derive equations of motion within a Hartree-Fock-Bogoliubov-Popov approximation suitable for both equilibrium and out of equilibrium calculations. We obtain equilibrium solutions of these equations that allow us to obtain phase boundaries and excitation spectra that we compare to previous equilibrium results obtained in a 1-loop calculation [62] and numerically exact results where possible. This paper is structured as follows. In Section 2, we describe the model that we study and derive the 2PI effective action for the BHM. In Section 3, we obtain the equations of motion for both the order parameter and the two-particle Green's function by taking appropriate variations of the 2PI effective action. In Section 4, we study the equilibrium solution of the equations of motion at the HFBP level. Finally in Section 5 we discuss our results and present our conclusions.

Model and formalism
In this section we introduce the Bose Hubbard model and discuss the generalization of the 1PI approach developed in Ref. [31] to a 2PI effective action within the Schwinger-Keldysh formalism. The Hamiltonian for the BHM, allowing for a time dependent hopping term, iŝ where Figure 1: Contour for a system initially prepared at time t i in a thermal state with inverse temperature β. t f is the maximum real-time considered in the problem, which may be set to t f → ∞ without loss of generality.
withâ † r andâ r annihilation and creation operators for bosons on lattice site r respectively,n r ≡â † râ r the number operator, U the interaction strength, and µ the chemical potential. The notation r 1 , r 2 indicates a sum over nearest neighbours only. We allow J r1 r2 (t), the hopping amplitude between sites r 1 and r 2 , to be time dependent.

Contour-time formalism
We use the contour-time formalism [64][65][66][67][68][69], which treats time as a complex variable lying along a contour. For systems initially prepared in thermal states, which we consider here, one can work with a contour C of the form illustrated in Fig. 1. One obtains the imaginary-time Matsubara formalism, which is restricted to equilibrium problems, by setting t f = t i . If one does not work in the Matsubara formalism, t f can be set to ∞ without loss of generality [70]. Furthermore, if one were to set instead t i → −∞, then one can obtain the real-time Schwinger-Keldysh closed-time path, which is suitable for both equilibrium and out of equilibrium problems, as the imaginary part of the contour would not contribute anything to the dynamics of the system. By setting t i → −∞, one is effectively discarding transient effects. Since we are interested in studying transient phenomena, we do not set t i → −∞ and instead work with the general contour illustrated in Fig. 1. A number of authors have applied contour-time approaches to the BHM [31,63,[71][72][73][74][75][76][77][78][79] -our work differs from previous approaches in that we apply a 2PI approach within the contour formalism that is appropriate for strong coupling as well as weak coupling [63,76].

Green's functions and the 1PI generating functionals
To characterize spatio-temporal correlations in the BHM we calculate contour-ordered Green's functions (COGFs). We generalize the work in Ref. [31] to include Green's functions with unequal numbers of annihilation and creation operators to allow for the study of broken symmetry phases. We frequently use the compact notationâ a r for the bosonic fields, defined bŷ We define the n-point COGF as [69] G a1...an r1... rn (τ 1 , . . . , τ n ) ≡ (−i) n−1 Tr ρ i T C â a1 r1 (τ 1 ) . . .â an rn (τ n ) whereρ i is the state operator for a thermal state representing the initial state of our system andâ a r (τ ) are the bosonic fields in the Heisenberg picture with respect toĤ Here we have introduced explicitly the complex contour time argument τ , the sub-contour C (τ, τ ′ ) which goes from τ to τ ′ along the contour C, and the contour time ordering operator T C , which orders strings of operators according to their position on the contour, with operators at earlier contour times placed to the right. Note that the presence of T C in Eq. (5) leads to symmetry under permutations {p 1 , . . . , p n } of the sequence {1, . . . , n}: At times it will be useful to express the contour time τ in terms of a contour label α (commonly called a Keldysh index) indicating a contour time located on C α and a positive real parameter s such that e.g. we can rewrite the bosonic fieldsâ a r (τ ) aŝ a a r,α (s) ≡â a r (τ ) , In order for the Heisenberg fieldsâ a r (τ ) to be well-defined, we need to analytically continue the BHM Hamiltonian [Eq. (1)]. For the contour considered in this paper,Ĥ BHM (τ ) is analytically continued as followsĤ

4
The COGFs above can be derived from a generating functional Z [f ] defined as in the (+, −, T ) basis, the f s are source currents, the overscored index in f a r,α (s) is defined by

Path integral form of Z [f ]
We cast the generating functional Z [f ] in the path integral form [67], which for the case of the BHM is where S BHM is the action for the BHM, and´[Da a ] is the coherent-state measure. We absorb overall constants into the measure as they will cancel out in the calculation of the COGFs due to the factor of 1/Z [f = 0] in Eq. (18). Note that in the path-integral formalism a 1 r,α = a r,α and a 2 r,α = a * r,α . In this formalism, we can rewrite averages of the form T C [. . .] ρi as follows where contour ordering is now implicit in the path integral representation [80]. In addition to the generating functional, we make extensive use of the generator of connected COGFs (CCOGFs) defined by The where . . . c indicates that only connected diagrams are kept. Note that the CCOGFs satisfy the same symmetry property as the COGFs

Keldysh rotation
For the n-point CCOGF defined in Eq. (22) there are 3 n Keldysh components. However, as a consequence of causality, we can eliminate n−1 m=0 ( n m ) of these components by performing the following transformation on the bosonic fields [65]  whereã q andã c are the quantum and classical components of the field respectively [74,[81][82][83], andã T = a T . After the above basis transformation (+, −, T ) → (q, c, T ), the matrix τ 3 becomeŝ the limits of integration become Following the argument given in Ref. [74], multiplying out the products in the second T C [. . .] yields 2 n−m path-ordered terms. The key point to note is that within any one of these path-ordered products the position of the field with the largest s does not depend on its Keldysh index. This implies that for each path-ordered product there is another path-ordered product is which is identical except with opposite sign. Therefore every term cancels out. It immediately follows that the associated CCOGFs vanish as well: Moreover, any permutation of the Keldysh indices in Eq. (29) will also yield a vanishing CCOGF. Since there are ( n m ) distinct permutations for fixed n and m, there are n−1 m=0 ( n m ) components that will vanish in total. This completes the proof. Note that if we were working with a closed-time path, where there is no imaginary appendix to the contour, we recover the special case where only ( n 0 ) = 1 Keldysh component vanishes, namelyG a1...an,c r1... rn,q...q (s 1 , . . . , s n ) [65,74]. After performing the Keldysh transformation, the BHM action takes the form [31] (dropping tildes) where r,α1 (s) a a2 r,α2 (s) a a3 r,α3 (s) a a4 r,α4 (s) , In the (q, c, T ) basis, the source term becomes and the CCOGFs are

Effective theory for the Bose-Hubbard model
In order to study quench dynamics in the BHM, we make use of an effective theory that can describe both the weak and strong coupling limits of the model in the same formalism. Such an approach was developed in imaginary time by Sengupta and Dupuis [62] by using two Hubbard-Stratonovich transformations and generalized to real-time in Ref. [31]. A similar real-time theory was also obtained based on a Ginzburg-Landau approach using the Schwinger-Keldysh technique [72][73][74]. A brief discussion of the derivation of the effective theory along with minor corrections to several expressions presented in Ref. [31] is given in Appendix A. The effective theory obtained in Ref. [31] for the z fields (which are obtained after two Hubbard Stratonovich transformations and have the same correlations as the original a fields [62]) is 8 where (G c ) −1 is the inverse of the two-point CCOGF in the atomic limit (i.e. J = 0), u (4) is and the inverse of an arbitrary two-point function X satisfieŝ Both (G c ) −1 and u (4) are independent of site index r, hence we write them without site labels. However, throughout this paper we occasionally include the site labels when it serves to provide more clarity to the reader. One would have to include the site labels if for instance one considers the BHM with a harmonic potential as is realised experimentally.
Equation (40) is the key result from Ref. [31] that we use to develop the 2PI formalism in Section 2.6. However, before applying the 2PI formalism to this action, we need to include an additional correction term: whereũ (2) contains an infinite set of diagrams, although here we truncate it keeping only the lowest order term:ũ a1a2 α1α2 (s 1 , This correction term ensures that our equations of motion are accurate to first order in G (4),c (see Appendix A for further discussion). Moreover, it ensures that the equations of motion for the two-point CCOGF we derive in Section 3 are exact in the atomic (J = 0) limit, which is essential when considering quenches beginning in the atomic limit. This action also gives the exact two-point CCOGF in the noninteracting (U = 0) limit [62]. These features make this theory particularly appealing for the study of quench dynamics, since it gives the hope that one can accurately describe the behaviour of the system in both the superfluid and Mott-insulating regimes [6].

2PI Formalism and the effective action
In order to obtain the full two-point CCOGF (the "full propagator" from now on), which encodes nonlocal spatial and temporal correlations, we adopt a 2PI approach. Unlike 1PI approaches [31,[72][73][74], the 2PI formalism describes the evolution of the mean field (i.e. superfluid order parameter for the BHM) and the full propagator on equal footing [63]. Several authors [63,75,76] have applied the 2PI formalism to the BHM to derive equations of motion for the mean field and the full propagator for weak interactions.
Here, we develop a real-time 2PI approach based on the strong-coupling theory of Sengupta and Dupuis [31,62] to capture behaviour of correlations across a quantum quench. We adopt a compact notation where we write an arbitrary function X as We extend the Einstein summation convention to the τ subindices such that for two arbitrary functions X and Y we have We can rewrite Eq. (40) (with the correction term [Eq. (43)] included) in the condensed notation as where we have introduced the generalized inverse bare propagator g −1 with In the 2PI formalism [70,84], physical quantities are expressed in terms of the mean field φ and the full Note that G c is symmetric: G a1a2,c r1 r2,τ1τ2 = G a2a1,c r2 r1,τ2τ1 . The equations of motion for φ and G c are obtained by requiring the 2PI effective action Γ [φ, G c ] be stationary with respect to variations of φ and G c . This is similar to the 1PI case where the equations of motion for φ are obtained by requiring the 1PI effective action Γ [φ] to be stationary with respect to variations of φ. The full propagator from the 2PI effective action allows one to take into account broken symmetry states [70,84], which is necessary to describe quenches in the superfluid regime.
To obtain the effective action we define the 2PI generating functional for Green's functions Z where in addition to the single-particle source current f , we have included a (symmetric) two-particle source current K. Note that φ and G c are obtained by calculating the following functional derivatives of W [f, K]: These equations implicitly give f and K as functions of φ and G c : where f and K should be understood as being expressed in terms of φ and G c . The following identities can be derived [70,84] from Eq. (58) δΓ Defining the effective action can be shown to take the form [70,84] where is the sum of all 2PI connected vacuum diagrams in the theory with vertices determined by the action and the propagator lines determined by G c , i.e.
, and φ a 1 One can use Eq. (64) along with Wick's theorem to generate all the diagrams in Fig. 2 up to second-order in the four-point vertex u (4) . The solid dots represent the interaction vertices u (4) , the solid lines represent G c , and the dashed lines represent φ (as illustrated in Fig. 3). In this paper, we only consider the first diagram in Fig. 2, i.e. the double-bubble (D.B.) diagram, which was also considered (along with the remaining two diagrams) in Refs. [63,76] where the BHM was studied at weak coupling. However, there is an important distinction between the calculations here and those in Refs. [63,76], which is that the interaction vertices in Refs. [63,76] are local in both space and time, whereas the interaction vertices we consider are local in space but nonlocal in time -this leads to additional features in the equations of motion. The contribution from the D.B. diagram is

Equations of motion
To calculate the equations of motion, first we use Eqs. (59) and (60) and set the sources to zero, giving δS δφ a1 r1,τ1 and where the second equation is Dyson's equation with the 2PI self energy. Given the form of the bare propagator in our strong-coupling theory, the equations of motion Eq. (66) and (67) in their above formulations are not suitable for dynamical calculations. We begin by reformulating Eq. (66). First, we explicitly calculate the first term in Eq. (66) The second term in Eq. (66) can be written as We act on both sides of Eq. (66) with G c from the left and rearrange terms to get where we have introduced the quantity Eq. (71) is a much more suitable form for dynamical calculations.
Next we reformulate Eq. (67) into a more appropriate form. First, we separate D −1 a1a2 r1 r2,τ1τ2 as follows where is the 1-loop contribution to the total self energy. If we define the full self energy as then Eq. (67) becomes After rearranging a few terms, one obtains which is a more suitable form for dynamical calculations. That being said, the form shown here is still not particularly amenable to solution. We now discuss simplifications that allow us to obtain equations of motion that are easier to solve.

Low-frequency approximation
Equations (71) and (77), whilst having a compact form in our notation, contain as many as four timeintegrals, making it computationally expensive to solve the equations numerically. This suggests that some level of approximation is required in order to obtain more physical insight from the equations above. Following Refs. [31], we focus on the low frequency components of the equations of motion. In a quench protocol this would correspond to considering changes that are slow enough that the equations of motion are dominated by low frequency terms. The approximation also applies to equilibrium calculations where there is no quench at all. The low-frequency approximation we consider involves taking the static-limit of the four-point vertex u (4) . If we only consider values of the chemical potential away from the degeneracy points between adjacent Mott lobes, i.e. µ ≈ U r, with r an integer, then the static limit of u (4) can be expressed as [31,62,74] u a1a2a3a4 where u 1 and u 2 2 are defined in Appendix D, ζ a1a2a3a4 α1α2α3α4 is defined in Eq. (34) and Numerical evaluation of u 1 and u 2 2 for a homogeneous system, shown in Fig. 4 demonstrates that unless µ/U is close to an integer, the u 1 terms will dominate the u 2 2 terms. Moreover, for low temperatures, u 2 2 becomes negligible and goes to zero as β → ∞. Hence, to simplify the equations of motion, we further assume that the temperature is sufficiently low such that u 2 2 can be safely ignored. The end result is that the equations of motion contain single time-integrals only.
3.2. Keldysh structure of φ, G c , Ω, Σ Before presenting numerical results, it is worth discussing the explicit Keldysh structure of the mean field φ, full propagator G c , and their respective interaction terms Σ and Ω. Starting with the mean field φ, we have where φ a1 r1 (s 1 ) is the superfluid order parameter 15 Note that φ 2 r1 (s 1 ) = [φ r1 (s 1 )] * . Then, following Ref. [85], we can express G c as follows where G (R) and G (A) are the retarded and advanced Green's functions respectively, G (K) is the Keldysh or Kinetic Green's function, G (⌈) and G (⌉) are the left and right Green's functions respectively, and G (M) is the Matsubara Green's function.
Next we have Ω, which takes on the following Keldysh structure where to first order in u 1 we have The self energy Σ is similar in structure to G where we have where Σ (R) and Σ (A) have the same properties of causality as G (R) and G (A) respectively. To first order in and Lastly, we rewrite the equations of motion Eqs. (71) and (77) (90) together form one of the main results of this paper. These can be readily used to study out of equilibrium dynamics for strongly interacting systems. By considering only terms up to first order in u 1 , our approximation can be thought of in some sense as a Hartree-Fock-Bogoliubov (HFB) approximation in the strong-coupling regime. In future works we will study these equations of motion for various nonequilibrium scenarios. In the remainder of this paper however, we study the equilibrium solutions to the equations of motion above, going beyond the work in Ref. [62] in which only the equilibrium solutions at the one-loop level in the imaginary-time formalism were studied.

Equilibrium solution
In studying the equilibrium solution to the equations of motion derived in the previous section we consider a homogeneous system at zero temperature. As a result, it is easier to work in k-space rather than real space. In equilibrium, the mean field equation of motion Eq. (96) reduces to [85] where we used the fact that the superfluid order parameter is constant in time, φ 1 (s 1 ) = φ. Expressions for G 12,(R) (ω) and G 12,(R) (ω ′ = 0) are given by Eqs. (C.8) and (D.2) respectively. We also have that in equilibrium all the various real-time Green's functions may be expressed in terms of the spectral function G (ρ) One can calculate G (K) from G (ρ) via the fluctuation dissipation theorem (FDT) [70,85], which at zero temperature is hence one need only focus on the G (R) equation of motion directly. In equilibrium, it is easier to work in frequency space, hence we may rewrite the G (R) equation of motion as [85] G a1a2,(R) k where and n and n 0 are the average particle densities for J = 0 and J = 0 respectively. Note that With a bit of algebra, one can show that From here, the next step is to simplify G where s is the branch number, ∆E are the corresponding spectral weights. Once written in this form, we can simply read off the expressions for the desired quantities. We do this in the following by considering the Mott insulator and superfluid cases separately.

Mott insulator phase
In the Mott insulator phase, φ = Σ One can rewrite Eq. (114) as where ∆E (±) and ∆E (±) are the excitation energies in the atomic limit (i.e J = 0) Using Eq. (103) along with the Sokhotski-Plemelj theorem we obtain for the spectral function which in turn depends on G 12,(K) k (s ′ = 0) through in the Mott insulator phase. Using Eq. (104) we obtain for G and therefore Hence the self-consistent solution can be formulated as follows: 1. Make an initial guess for n. , and the quasi-momentum distribution n k for a square lattice system with µ/U = 0.42, J/U = 0.04, and βU = ∞. The 1-loop solution, which was studied in Ref. [62], amounts to approximating the self-energy by Σ 12,(R) k = ǫ k in the Mott-insulating phase. From Fig.  5 we see that there is little qualitative change in the excitation energies between the two approximations. The same can be said for the spectral weights for values of k well away from zero, however there are appreciable differences in the long-wavelength limit. These differences can be more clearly visualised in the quasi-momentum distribution where we see that the k = 0 peak is sharper in the 1-loop approximation than the HFB approximation.
One way to account for the differences in the spectral weights is to consider how well each solution scheme approximates the phase boundary between Mott insulating and superfluid phases. In Fig. 6 we compare the mean-field (MF) and HFB approximations of the phase boundary along with the exact calculation. Figure 6 clearly shows that there is significant quantitative improvement in the phase boundary calculation when going from the MF level to the HFB level. Moreover, in 1 dimension, where the MF approximation is expected to be poor, we have a clear qualitative improvement in the phase boundary calculation. The closer we are to the phase boundary (in the Mott-insulator phase), the sharper the k = 0 peak is in n k . Since the MF approximation always underestimates the location of the phase boundary more than the HFB approximation, the 1-loop approximation -which uses the MF approximation of φ -will wrongly predict a sharper peak as compared to that in the HFB case. Equivalently, the 1-loop approximation will always overestimate the values of the spectral weights in the neighbourhood of k = 0.
Another way to assess the accuracy of the two approximation schemes in the Mott-insulating phase is to look at the average particle density n [Eq. (125)]. In the Mott-insulating phase, n = ⌈µ/U ⌉. For the same parameter values mentioned above, we have n ≈ 1.08, (HFB), (130) n = 1.00, where we see that the HFB approximation yields a significant improvement as compared to the 1-loop approximation.

Superfluid phase
In the superfluid phase, φ and Σ

22,(R,A) k
are non-zero, hence we must use the full form of Eq. (111). We begin by calculating φ from Eqs. (102) and (90). Without loss of generality, we can assume that φ is real which further implies that the quantities iG (s ′ = 0) are real. Based on these assumptions we obtain As is clear from Eq. (132) the mean field φ needs to be solved self-consistently along with the full propagator G. We now calculate G (R) . Starting from Eq. (111), one can show that where ∆E (s) In a moment we will show that the ∆E (s) SF, k are the excitation energies in the SF phase. Before doing so, it is worth commenting on our approximation for the self energy in the superfluid phase. In Appendix E we show that in the full HFB approximation the excitation spectrum is not gapless, violating Goldstone's Theorem, whereas if we ignore contributions from the anomalous Keldysh Green's function iG there is a gapless spectrum. The latter scheme is called the HFB-Popov (HFBP) approximation [89]. Thus in the HFBP approximation we have The HFBP approximation is most accurate for values of the chemical potential away from integer values which is evident from the fact that G , which in turn is proportional to u 1 , which is small for values of the chemical potential away from integer values. Therefore iG

k = 0
When k = 0, we can derive the spectral function from the retarded Green's function as we did above in Sec. 4.1 using the Sokhotski-Plemelj formula as we did in the MI case [Eq. (123)] where z (s,±) are the excitation energies and spectral weights respectively. Moreover, for each branch the particle excitation energy is equal to the hole excitation energy. Using Eq. (104) we have for the Keldysh Green's function (141)

k = 0
In the zero-quasi-momentum case, G (142) One cannot use the same Sokhotski-Plemelj formula as we did above in deriving the spectral function, instead one must used a generalized version of the formula Doing so yields the following spectral function

25
In both k cases, G
In the case where k = 0, one needs to be careful when calculating G

12,(K) k=0
(ω) as the FDT [Eq. (104)] is ill-defined for ω = 0. Fortunately, G 12,(K) k (ω) = 0 (see Appendix F for a proof). Therefore we have for 4.2.3. Calculating n k and n One can calculate n k from where And lastly, the average particle density n is calculated using Eq. (125). Therefore, at the HFBP level, the system can be solved self-consistently as follows: 1. Make an initial guess for n. 7. Use n k to recalculate n via Eq. (125). 8. Repeat steps 2 to 7 until self-consistency is reached.
In Fig. 7, we compare the 1-loop and HFBP equilibrium solutions in the superfluid phase by calculating the excitation energies ∆E = u 1 (φ) 2 in the superfluid phase. We see that there is little qualitative change in the excitation energies between the two approximations. Moreover, the spectral weights in the second branch s = 2 change very little as well. We do observe appreciable differences in the spectral weights for the first branch s = 1 in the long-wavelength limit, similar to the Mott-insulator case. As was argued for in the Mott-insulator case, since the HFBP calculation yields a more accurate phase boundary, we believe this method will also yield a more accurate result for z for the second branch. Note that Γ = (0, 0), M = (π, π), and X = (π, 0).

Phase boundary
To calculate the phase boundary, we make a slight modification to our solution scheme for the MI phase. The modification comes from the extra step of calculating the critical hopping J c . Consider again the φ-equation Eq. (132). At the boundary, φ = G 22,(K) r=0 (s ′ = 0) = 0. Solving for J we get With this established, we can outline the phase boundary solution as follows 1. Make an initial guess for the average particle density n 2. Use n to calculate the hopping J c , see Eq. (150) 3. Use n and J c to calculate the self-energy Σ This calculation ends up reproducing the phase boundary found from the Mott insulating side since the anomalous Green's functions vanish at the phase boundary.

Discussion and Conclusions
The ability to address single sites in cold atom experiments [11] has allowed for experimental exploration of spatio-temporal correlations in the BHM [49]. This has led to theoretical investigations of these correlations in both one [48] and higher dimensions [46,51,59,61] in the presence of a quench. In dimensions higher than one, where numerical approaches are limited, a theoretical challenge has been to develop a framework which can treat correlations in both the superfluid and Mott insulating phases over the course of a quench. An important result in this paper is that we have developed a formalism that allows for the description of the space and time dependence of correlations in both phases during a quench. The specific approach we took was to derive a 2PI effective action for the BHM using the contour-time technique building on the 1PI real-time strong-coupling theory developed in Ref. [31] which generalized the imaginary-time theory developed in Ref. [62]. From this 2PI effective action we were able to derive equations of motion that treat the superfluid order parameter and the full two-point Green's functions on equal footing. We emphasise that our formalism is applicable even in the limit of low occupation number per site.
Even at the level of the 1PI real-time theory, the quartic coupling becomes non-local in time, which in the 2PI theory leads to complicated expressions in the equations of motion, involving up to four time integrals, even at the first order in the interaction vertices. We showed that by taking a low frequency approximation, this complexity can be reduced to at most a single time integral. The equations of motion obtained at this point are somewhat similar to previous 2PI studies of the out of equilibrium dynamics of interacting bosons [63,75,[90][91][92][93]. However, in contrast to these previous studies, the equations of motion we obtain are a series of integral equations rather than integro-differential equations.
We showed that taking a HFB(P) approximation of the 2PI effective action yields significant improvements to the calculation of the particle density and phase boundary when compared to the 1-loop approximation considered in Ref. [62]. Our results also suggest that the HFB(P) approximation gives a better account of the spectral weights in the long-wavelength limit. These improvements in the equilibrium case suggest that our formalism should be suitable for accurately describing spatio-temporal correlations in nonequilibrium scenarios.
The space and time dependence of correlations after a quantum quench give insight into the propagation of excitations generated by that quench, and hence we hope that the formalism we have developed here will allow further theoretical investigation of the excitations after quenches in the BHM, to complement experimental efforts in the same direction. In future work we plan to investigate a broad range of quench protocols, including quenches in the Mott phase where one can study the light-cone-like spreading of singleparticle correlations. Other quench protocols of interests are those beginning in the superfluid phase and then ending in the Mott phase. In such scenarios, one may be interested in studying for example the possibility of aging-like phenomena. Lastly, we plan to investigate generalizations such as the inclusion of a harmonic trap, coupling to a bath [71,94] Note that Eq. (A.8) corrects Eq. (6) in Ref. [31]. Moreover, note that for the uniform BHM as considered here, the atomic CCOGFs are independent of site index, and so we drop these indices when they do not affect the clarity of the exposition in this paper. Inverting Eq. (22), with G c → G c , we may rewrite W 0 as As pointed out in Ref. [62], the quadratic terms in the equilibrium action of the form in Eq. (A.12) allow one to calculate the mean-field phase boundary, however it yields an unphysical excitation spectrum in the superfluid regime. This issue is circumvented by performing a second Hubbard-Stratonovich transformation [31,62]. Starting from Eq. (A.5) (keeping the source currents f this time), we introduce a second field z such that  20) and theũ vertices contain an infinite set of "anomalous" diagrams, i.e. diagrams that contain internal inverse bare propagator lines. Such diagrams have no physical meaning and should not contribute to the physical quantities [95]. It should be noted that in addition to the physical diagrams, the u vertices also generate "anomalous" terms. In Appendix B, we show that these anomalous terms cancel one another out when calculating the superfluid order parameter φ and the full two-point CCOGF. That being said, the action in Eq. (A.19) contains an infinite sum, therefore one will eventually have to truncate said action which will ultimately lead to only certain subclasses of "anomalous" terms cancelling out.
In this paper, we truncate the action to quartic order in the z fields and neglect any contributions fromũ (4) . In Refs. [31,62], allũ terms were neglected. By including theũ term given in Eq. (A.22), one obtains equations of motion which are accurate to first order in G (4),c , which is not the case in Refs. [31,62]. Lastly, we stress that this approach leads to a strong-coupling theory that is not simply an expansion order by order in J/U .

Appendix B. Cancellation of anomalous diagrams
In this appendix, we show that the anomalous terms introduced in Appendix A do not contribute when calculating the mean field φ and the two-point CCOGF G c of the original field a. For the sake of economy in writing, we adopt the notation introduced in Section 2.6 and condense it even further such that X x1...xn ≡ X a1...an r1... rn,τ1...τn , (B.1) where we performed the field substitution ψ x → −ψ x . We first establish a relationship between the expectation values of the a-field, φ x , and of the ψ-field, V x . To do this, we start by calculating φ x1 = a x1 as follows and then integrate by parts to get which establishes a relation between φ x and V x . Note that where Φ is some arbitrary field. By similar calculation, one can show that where V c x1x2 is the two-point CCGOF for the field ψ. Taking the inverses of the above relations yields We now use the ψ theory to calculate the 2PI equations of motion for V x1 and V c x1x2 . The action S aux [ψ] for the auxiliary field ψ can be expressed as and hence using this action in Eqs. (66) and (67) and rearranging terms, we obtain the following relations where Ξ and Σ are obtained from the corresponding Γ 2 . Next, we apply Eqs. (B.8) and (B.9) to obtain recursive expressions for φ and G c (B.14) We now derive recursive relations for φ and G c by an alternative approach: we apply the 2PI approach to the theory of the z-fields, allowing for anomalous terms, which is given by Eq. As noted in Appendix A, the Green's functions for the z-fields are the same as those for the a-fields.
Similarly to the calculations leading to the recursive relations V x1 and V c x1x2 , we calculate the following recursive 2PI relations for φ and G c and n 0 and E n are given by Eqs. (110) and (122) respectively. Given that the Fourier transforms G 12,(R,K) (ω) are used throughout this paper, it is worth explicitly writing out the expressions for these particular Keldysh components Appendix D. Low frequency approximation to four-point vertex u (4) To calculate the low frequency approximation to the four-point vertex u a1a2a3a4 α1α2α3α4 (s 1 , s 2 , s 3 , s 4 ), we begin with Eq. (41). We make use of the time-translational invariance of the atomic two-point Green's function and take the low-frequency approximation, which gives (noting that there is no contribution from the Keldysh Green's function except at points where the Mott lobes are degenerate) [31]