Contour-time approach to the disordered Bose-Hubbard model in the strong coupling regime

There has been considerable interest in the disordered Bose Hubbard model (BHM) in recent years, particularly in the context of thermalization and many-body localization. We develop a two-particle irreducible (2PI) strong-coupling approach to the disordered BHM that allows us to treat both equilibrium and out-of-equilibrium situations. We obtain equations of motion for spatio-temporal correlations and explore their equilibrium solutions. We study the equilibrium phase diagram as a function of disorder strength and discuss applications of the formalism to out-of-equilibrium situations. We also note that the disorder strengths where the emergence of non-ergodic dynamics was observed in a recent experiment [Choi $et \,al.$, Science $\bf{352}$, 1547 (2016)] appear to correspond to the Mott insulator -- Bose glass phase boundary.


I. INTRODUCTION
Isolated strongly interacting quantum systems with quenched disorder may fail to thermalize and enter a phase in which the entire spectrum is composed of localized states [1,2].Such many-body localized (MBL) states have been studied intensively in recent years .In addition to their interest from a fundamental point of view, MBL states have also been suggested as having potential for use as quantum memories.
Experimentally, there are indications of localization in disordered, interacting many body cold atom systems in optical lattices in two and higher dimensions [59][60][61][62][63].Of particular interest for our work is the experiment by Choi et al. [61], in which the relaxation dynamics of disordered bosons in a two dimensional optical lattice were studied.Starting from an initial condition in which all of the atoms were localized on one side of a trap, Choi et al. observed the imbalance as a function of time and found that beyond a critical disorder strength, their system failed to thermalize in the time window of their experiment.Yan et al. [64] applied Gutzwiller mean-field theory (GMFT) to the two-dimensional disordered Bose Hubbard model and were able to reproduce the main experimental results, even though GMFT is unable to capture MBL, raising the possibility that the experiments probe glassy dynamics rather than MBL.This highlights the need to develop theoretical methods to investigate the out-of-equilibrium dynamics of the disordered Bose Hubbard model in dimensions greater than one.
An alternative approach has been developed by two of us that is based on a two particle irreducible (2PI) out-of-equilibrium strong coupling approach to the BHM (2PISC) [91][92][93][94][95].This approach allows the treatment of the dynamics of the order parameter and correlation functions on an equal footing and we have previously used it to demonstrate excellent agreement [96] with experiments investigating the spreading of correlations for bosons in optical lattices in one and two dimensions [78,97].It also has the attractive feature that it allows for the inclusion of disorder averaging, which we make use of to study the disordered Bose-Hubbard model.
The presence of disorder in the BHM can lead to an additional phase in between the superfluid and Mott insulator, the Bose glass [98], and can reduce the size of the Mott lobes [98,99].The experiments by Choi et al. [61] have focused attention on the out-of-equilibrium dynamics of the disordered Bose Hubbard model and both many-body localization and glassiness for bosons in one [100][101][102][103][104] and two or more dimensions [53,57,64,102,[105][106][107][108][109][110][111].This activity motivates our extension of the 2PISC formalism for the BHM to include disorder arXiv:2304.11260v1[cond-mat.quant-gas]21 Apr 2023 to provide an additional route to investigate the out-ofequilibrium dynamics of the disordered BHM.
The main result of this paper is that we develop a 2PI framework that allows us to treat both the equilibrium and out-of-equilibrium behaviour of the disordered Bose Hubbard model.This allows us to obtain equations of motion for the superfluid order parameter and spatio-temporal correlations.We obtain solutions of these equations in the equilibrium case and investigate the Mott insulator phase boundary as a function of disorder strength and calculate the collective excitation spectrum both in and outside the Mott phase.We find that our results compare favourably with quantum Monte Carlo (QMC) simulations in two [112] and three [113] dimensions.We also note that the disorder strengths at which Ref. [61] found the emergence of non-ergodic dynamics appear to correspond to the Mott insulator -Bose glass phase boundary.
This paper is structured as follows: in Sec.II we introduce the disordered Bose Hubbard model and formalism, deriving an effective theory.We use our effective theory to obtain 2PI equations of motion which we then solve for equilibrium properties of the model in Sec.III.We conclude and discuss our results in Sec.IV

II. MODEL AND FORMALISM
In this section, we introduce the disordered Bose Hubbard model and discuss the generalization of the strongcoupling approach developed in Refs.[91,92,114] for the standard BHM to the disordered case allowing for both equilibrium and out-of-equilibrium behaviour.The Hamiltonian for the disordered BHM is where 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, V r a harmonic trapping potential, µ the chemical potential, and r an on-site disorder potential.The disorder potential is drawn from a Gaussian distribution 1. Contour for a system initially prepared at time ti.t f is the maximum real-time considered in the problem, which may be set to t f → ∞ without loss of generality.
with ∆ being the full-width at half maximum for the distribution.The notation r 1 , r 2 indicates a sum over nearest neighbours only.

A. Contour-time formalism
The general formalism that we discuss and adopt in this paper was developed in a previous paper by two of us; we refer the reader to Ref. [92] for further details on the formalism.We use the contour-time formalism [115][116][117][118][119][120], which replaces the notion of real time along the real line with contour time, a complex valued time on a contour in the complex plane.Furthermore, an appropriate choice of contour is particularly attractive for studying disordered systems as it eliminates the need to use replicas in carrying out the average over the quenched disorder [121,122].For systems initially prepared in outof-equilibrium states, one can work with a contour C of the form illustrated in Fig. 1.A popular alternative to this contour is the Schwinger-Keldysh (SK) closed-time path [115,116] which is also suitable for certain out-ofequilibrium problems.However, unlike contour C, the SK contour ignores transient phenomena and, more importantly, information about the initial state.Given that we are interested in comparing long-time density profiles with that of the initial state, contour C is a more appropriate choice.

B. Contour-ordered Green's functions
In deriving our effective theory of the disordered BHM, we calculate various contour-ordered Green's functions (COGFs).We define the n-point COGF as [120] 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 representing the initial state of the system and the a i upper indices are defined such that and âa r (τ ; ) are the bosonic fields in the Heisenberg picture with respect to ĤBHM (τ ; ) [Eq. (1)] 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.Given the somewhat cumbersome notation in expressions such as that in Eq. ( 6), we make extensive use of a compact notation where we write an arbitrary function X as X a1...an r1... rn,τ1...τn; ≡ X a1...an r1... rn (τ 1 . . .τ n ; ) , (10) and introduce the following implicit summation convention where σ i is the i th Pauli matrix, 1 = 2 and 2 = 1.Note that we only include the parameter in Eq. (10) if the function X depends on the disorder configuration.

C. Generating functional Z [f ; ]
The COGFs above can be derived from a generating functional Z [f ; ], which can be cast in the following path integral form [91,92,94,118]: where S [a; ] is the action for the disordered BHM with with and where S f [a] is the source term where contour ordering is now implicit in the path integral representation [123].Occasionally, we drop the action subscript . . .S → . . .for brevity.

D. Disorder averaging
We are ultimately interested in calculating disorder averaged COGFs where for a quantity θ we denote the disorder average θ with a carat.Using Eq. ( 23) we can determine an expression for calculating the disorder-averaged COGFs: where E. Effective theory of the disordered BHM We develop an effective theory that is suitable for studying the dynamics of the disordered BHM in the strong coupling regime.The approach can be outlined as follows: first we calculate the disorder average of Z [f ; ], which gives us an effective theory S dis eff in terms of the original a-fields, then we apply various Hubbard-Stratonovich transformations such that we can obtain a strong coupling expansion of the theory.The resulting effective strong coupling theory introduces two auxiliary fields z and Q.We derive identities relating the correlators of these two auxiliary fields to those of the original a-fields.One can then apply a two-particle irreducible effective action approach [124] to the effective theory to obtain equations of motion for the correlation functions.
We begin by performing the disorder average of where Next, using Eqs.( 5) and ( 16) we calculate e iS∆ [a] : where and To decouple the quartic a-field term in Eq. ( 30), we perform a Hubbard-Stratonovich transformation (similarly to e.g.Ref. [125]) where and Q is an auxiliary field introduced by the transformation.At this point, the generating functional Z [f, K] can be written as Next, following Refs.[91,92,114,126], we decouple the hopping term by performing another Hubbard-Stratonovich transformation where with and ψ is another auxiliary field.By making a field substitution, ψ a r,τ → −ψ a r,τ + f a r,τ , and rearranging terms in Eq. ( 36) we get where In this context ψ and Q take the same form as f and K in the 2PI generating functionals introduced in Ref. [92], i.e. by taking functional derivatives of with respect to ψ and Q one can generate all n-point COGFs (or CCOGFs).In this case the generating functionals Z 0 [ψ, Q] and W 0 [ψ, Q] are governed by a different theory than that introduced in Ref. [92].
W 0 [ψ, Q] generates all the n-point CCOGFs in the limit of zero disorder and hopping for a system prepared in the initial state ρi G a1...an,c r1... rn,τ1...τn = (−1) as well as a set of generalized CCOGFs defined by: These functions are connected in a particular sense: indices that are paired inside a pair of brackets . . .should be thought of as indices belonging to a single field.If we assume an initial state of the form the CCOGFs defined in Eqs. ( 41) and ( 42) vanish unless all site indices are equal.Moreover, when Eq. ( 43) holds, correlators of the form a a1 r1,τ1 . . .a an rn,τn S0 vanish unless the number of a * -fields equals the number of a-fields.By inverting Eq. ( 42) we may rewrite W 0 as where Θ (x) is the Heaviside function.Once again, following Refs.[91,92,114,126], we perform another Hubbard-Stratonovich transformation to decouple the inverse hopping term such that where with The effective theory is obtained by adding all of the action terms excluding sources We next perform a cumulant expansion of W 0 , similar to that found in Refs.[91,92,114,126] although the calculation is more complicated in the disordered case.Even in the compact notation we introduced in Sec.II B, the resulting expression for the effective theory is quite cumbersome to write out.We therefore condense the notation further such that Using the above shorthand notation, the effective theory can be expressed as follows where The couplings for quadratic terms in the theory are: with and the vertices g χ1...χn are combinations of the CCOGFs generated from W 0 .The presence of the contour ordering operator T C in the CCOGFs leads to symmetry under permutations {p 1 , . . ., p n } of the sequence {1, . . ., n}: The action in Eq. ( 51) contains an infinite sum, therefore for practical calculations we truncate the action keeping only terms to order O [ψ n Q m ] where 4 − 2m ≥ n ≥ 0 and 2 ≥ m ≥ 0. Out of these terms, the only non-vanishing terms (modulo index-permutations) are where and In this work, we do not consider a trapping potential for simplicity.Therefore, local quantities such as n-point CCOGFs in the atomic limit have no spatial dependency, and henceforth we will drop the r index for these quantities.

III. EQUATIONS OF MOTION
We now apply a 2PI approach to our effective theory to obtain equations of motion for the mean-field and the full two-point CCOGF for the Φ-fields (the "full propagator" from now on).In a related work [92], we applied this approach to the homogeneous BHM in the strong coupling limit.We follow the same general procedure in this work with a few modifications to account for the additional vertices and auxiliary fields that appear in our effective theory of the disordered BHM as compared to that of the homogeneous BHM.Here we only briefly outline the 2PI calculation.For a more detailed exposition of the 2PI approach, see Ref. [92].
We define the mean-field V χ1 and full propagator V χ1χ2,c as follows Next, we consider the 2PI Dyson's equation where Σ (1) χ1χ2 is the "1-loop" self-energy Σ (2PI) χ1χ2 is the 2PI self-energy and Γ 2 V (1) , V (2) is the sum of all 2PI connected vacuum diagrams in the theory with vertices determined by the action To first order in the vertices we have If we define the full self-energy Σ χ1χ2 to be then we may rearrange the Dyson's equation as follows Additionally, the equation of motion for the mean field As will be clear shortly, we are mostly interested in finding V z1 and V z1z2,c , that can be calculated from Eqs. ( 78) and ( 77) respectively, for the equilibrium and out-ofequilibrium scenarios .First, we start with V z1z2,c .From Dyson's equation ( 77) we can write where To make further progress, we cast Eq. ( 79) in terms of the correlation functions of the original a-fields.By inspection, one can see from Eq. ( 45) that Z [f, K] is the generator of COGFs of the z-fields in addition to the afields.This implies that where φ denotes the superfluid order parameter.Since we consider hopping strengths below the critical value of the Mott insulator to superfluid transition and assume that our initial state is of the form given in Eq. ( 43), we may safely assume that V z1 = 0. Note that this assumption is not valid when the system is in the superfluid phase as will be discussed in Sec.III B 2. For V Q12 , we follow a similar calculation to that in Ref. [127]: using Eqs.( 25) and ( 39) we can write then integrate by parts to get where Continuing with Eq. ( 85), we have hence Substituting Eqs. ( 53), ( 61)-( 64), ( 82), (83), and (88) into Eq.( 79) gives Ǧa4a2,c r4 r2,τ4τ2 , (89) where and the contributions to the self energy are and The next step is to calculate the mean field V z1 .Using Eqs. ( 51), ( 71) and ( 72) we can rewrite Eqs. ( 78) as where we used V Q12z3,c = V z1Q23,c .Next we need to calculate V Q12z3,c .Using Eq. ( 77) again, we get In order to make progress we treat the disorder strength perturbatively.Our results here are based on keeping terms to O ∆ 2 .Since [g 0 ] Q12Q45 is of order ∆ 2 , it forces the last term in Eq. ( 95) to be of order ∆ 4 .Therefore we approximate Eq. ( 95) by where we also used the fact that [g 0 ] Q12z3 = 0. Now, using Eqs.( 56), ( 61), ( 80), ( 81), ( 82), ( 83), ( 88) and (96) we can rewrite Eq. ( 94) as follows We will also be interested in calculating particle number to obtain the Mott insulator phase boundary, which can be calculated from G. To calculate G, we solve Eq. ( 89), however, the form shown here is still not particularly amenable to solution.We now discuss simplifications that allow us to obtain more tractable equations of motion.

A. Low-frequency approximation
Equation ( 89), whilst having a compact form in our notation, contains as many as four time-integrals, making it computationally expensive to solve the equations numerically.This suggests that some level of approximation beyond simply truncating the self-energy is required in order to obtain physical insight from the equations above.Following Refs.[91,92], we focus on the lowfrequency components of the equations of motion, specifically the self-energy terms Σ zz 2 and Σ zz 3 .Σ zz 1 + Σ zz 2 is almost identical in form to the self-energy obtained in Ref. [93], the only difference being that the u-vertices have a trivial spatial dependency.Therefore the low-frequency calculation of Σ zz 2 is almost identical to that in Ref. [93] and can be written as where ň r (τ ) is the disorder-averaged particle number at site r and contour time τ , and u 1 comes from taking the low-frequency approximation of u a1a2a3a4 τ1τ2τ3τ4 , the expression for which is given in Appendix B. Provided µ/U is not close to an integer, then u 1 /U 1 [91,92]; and we also focus on the small ∆ 2 limit, we keep terms of either order u 1 or ∆ 2 but not terms of order u 1 ∆ 2 or higher.Now, to calculate the low-frequency approximation to By casting u a1a2 a3a4 τ1τ2 τ3τ4 in this form, one can then perform a similar calculation to that for Σ zz 2 to obtain After applying the low-frequency approximation, we obtain a self-energy that is identical in Keldysh structure to the low-frequency self-energy obtained in Ref. [93] (for a discussion about Keldysh structure, see Ref. [92]).Therefore the equations of motion of the full propagator for the disordered system are identical in structure to those obtained in Ref. [93].With this in mind, it is straightforward to show that the equations of motion can be written as follows: where Ǧ(R) r r (t, t ) and Ǧ(A) r r (t, t ) are the disorder averaged full retarded and advanced Green's functions respectively obtained from performing a disorder average (as in Eq. ( 24)) of the quantities and Ǧ(K) r (t, t ) is the disorder averaged full kinetic Green's function obtained from disorder averaging 101) and ( 102) are the W 0 -generated retarded, advanced and kinetic Green's functions, respectively, which are all time-translational invariant (expressions for each are presented in Appendix A).It is important to note that ň r (t) can be obtained from Ǧ(K) r r (t, t ) as follows: which introduces an element of nonlinearity into the equations of motion above.Finally, in the low-frequency limit, we can approximate Eq. ( 97) as where is the low-frequency approximation of the inverse retarded Green's function obtained from W 0 , the expression for which is given in Appendix B.

B. Equilibrium solution
In studying the equilibrium solution to the equations of motion derived in Sec.III A we consider the system to be at zero temperature.We also work in k-space rather than real space.Whilst disordered systems are not homogeneous due to the random potential in the Hamiltonian, the disorder-averaged COGFs respect translation invariance.Therefore we follow the same procedure as Ref. [92] in order to obtain the Mott insulator phase boundary in the presence of disorder.The only difference being that now the COGFs are replaced by disorder-averaged COGFs.In k-space Eq. ( 89) becomes Ǧa1a2,(R) where where ň is the average local particle number for J = 0 Now, from Eq. ( 108) we can write Ǧ12,(R) Ǧ22,(R) In the following, we will discuss equilibrium solutions for the Mott insulator and superfluid phases.For simplicity we only consider the static limit, i.e., ω = 0.As we demonstrate later, this assumption is acceptable at least in equilibrium.

Mott Insulator
In the Mott insulator phase φ and Therefore Eq. ( 113) reduces to Ǧ12,(R) which we can rewrite as [92] Ǧ12,(R) where z(±) and ∆E (±) are the excitation energies in the atomic limit (i.e J = 0) As explained in Ref. [92], the local particle number can be expressed as for different disorder strengths ∆ in the Mott insulator phase where the relation between ∆ and ∆ is given in Eq. (30) and from now on we drop the index in ∆ to avoid confusion.In addition, we compare the quasi-momentum distribution ň k for two different disorder strengths: ∆/U = 0.1 and ∆/U = 0.3.
Here we have a 1000 × 1000 square lattice, the chemical potential is µ/U = 0.42, the chosen hopping term is J/U = 0.02, and βU = ∞.As can be seen from Fig. 2, in the Mott phase, the excitation energies decrease with increasing disorder strength.In addition, with increasing disorder strength, the quasi-momentum distribution ň k becomes more localized, which implies that for the same hopping term, the system is closer to the transition point.
(125) Now, in the superfluid phase, using Eqs.( 109), (110), (111) and (125) we have that the self-energy is and Next, we calculate Ǧ(R) : starting from Eq. ( 114), one can show that Ǧ12,(R) where and It is important to note that the expressions for B and C differ from those for the Mott insulator in that B has units of energy for the MI but units of [energy] 2 for the SF.Following the same approach as Ref. [92], the quasi-momentum where In Fig. 3 we show the collective mode spectra and quasi-particle spectral weight in the superfluid phase for different disorder strengths as calculated from Eqs. ( 130) and (134).To perform the numerical calculations we used a 1000×1000 square lattice, and set the chemical potential µ/U = 0.36, the hopping to J/U = 0.03, and βU = ∞.

Mott insulator phase boundary
To obtain the Mott Insulator -Bose Glass (MI-BG) phase boundary, we calculate the critical hopping J c at which φ = 0.This can be done numerically using Eq.(125).In Fig. 4 we show the phase boundary for different disorder strengths for dimensions one, two, and three.
In Refs.[112,113], the MI-BG transition for the chemical potential at the tip of the Mott lobe was calculated for two and three dimensional cubic lattices with random disorder uniformly distributed on the interval [−∆, ∆].In Fig. 5, using the disordered BHM in the strong coupling regime, we calculated the MI-BG transition and compared our results keeping terms to O ∆ 2 with those obtained using QMC simulations [112,113].We find good agreement between our results and QMC simulations.It should be noted that we consider a Gaussian distribution of disorder, while Refs.[112,113] used a box distribution.order parameter and two, three and four-point correlations.These equations apply both in and out of equilibrium.A strength of the 2PISC is that it is applicable in one, two, and three spatial dimensions.This is particularly advantageous for out-of-equilibrium dynamics, where numerical methods that are essentially exact, such as exact diagonalization or DMRG, are limited to one  [112] and for d = 3 from Ref. [113].
dimension or very small system sizes.A weakness of the 2PISC method is that in order to make progress, one needs to truncate the effective action which was done at quartic order, and this is an uncontrolled approximation.However, previous results in the clean case for both the phase boundary in equilibrium [92] and agreement with exact diagonalization for out-of-equilibrium dynamics [93] give confidence in its usefulness.It should be noted that the accuracy appears to be greatest for larger U/J [96].The main result of this paper is the derivation of the effective theory and the 2PI equations of motion, but as a check on the theory, we solved the disorder-averaged equations of motion in the equilibrium limit.We obtained the collective excitation spectra for the disordered BHM in both the Mott and superfluid phases and also obtained the Mott insulator phase boundary at a variety of disorder strengths in one, two, and three dimensions.We compared our results with QMC simulations performed in Refs.[112,113] and found very good agreement with the exact phase bound- ary.
Previous comparison of Mott insulator phase boundaries with QMC calculations in the clean case [92] found the 2PISC method gave a big improvement over the mean-field theory but was not in complete quantitative agreement with QMC calculation.We find a similar situation in the disordered case.However, there are caveats, in that we consider a Gaussian distribution of disorder rather than the box distribution used in Refs.[112] and [113].We also treat the disorder perturbatively and keep the highest-order term only (calculations to order ∆2 ).Given that there are uncontrolled approximations in the 2PISC method, these results give confidence in the results here and future applications to the out-ofequilibrium dynamics of the disordered Bose-Hubbard model.One limitation of our method is that we have not been able to determine the Bose glass -Superfuid phase boundary, which corresponds to the vanishing of the superfluid stiffness.
We noted that a motivation for our work was the experiments by Choi et al. [61] on thermalization in the disordered two-dimensional BHM.In that work, there was an identification of a thermal to MBL transition at a critical disorder value of ∆/U = 5.3(2) when U/J = 24.4.Using these same parameter values, in Fig. 6, we show that this point appears to lie essentially at the Mott insulator -Bose glass phase transition identified in QMC simulations.While our calculations and the QMC calculations in Ref. [112] do not include a trap, this result adds further to the questions raised in Ref. [64] as to whether the experiments in Ref. [61] probe an MBL transition or a glass transition.We intend to explore this question further in future work on out-of-equilibrium dynamics of the disordered BHM.

FIG. 5 .
FIG. 5. Comparison of the results from effective theory to O ∆ 2 and QMC for the Mott insulator -Bose glass transition for (a) d = 2 and (b) d = 3. QMC data for both the Mott insulator -Bose glass and Bose glass -superfluid transition taken for d = 2 from Ref.[112] and for d = 3 from Ref.[113].

FIG. 6 .
FIG.6.Comparison of the experimentally identified thermal-MBL transition point (Red dot) at U/J = 24.4 for unit filling with the Mott insulator -Bose glass transition curves for ∆/J = 5.1 and ∆/J = 5.5 (corresponding to ∆/J = 5.3(2)) obtained using the effective theory to O ∆2 .Blue and green stars are the corresponding QMC Mott insulator -Bose glass transition points for ∆/J = 5.1 and ∆/J = 5.5, respectively taken from Ref.[112].