{'fraction_non_alphanumeric': 0.047698342214328004, 'fraction_numerical': 0.012729425695677915, 'mean_word_length': 4.50295255548768, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A grand-canonical approach to the disordered Bose gas

We study the problem of disordered interacting bosons within grand-canonical thermodynamics and Bogoliubov theory. We compute the fractions of condensed and non-condensed particles and corrections to the compressibility and the speed of sound due to interaction and disorder. There are two small parameters, the disorder strength compared to the chemical potential and the dilute-gas parameter.


Introduction: grand canonical formalism
We approach the weakly interacting Bose gas with the grand-canonical Hamiltonian [1][2][3] The annihilation (creation) operatorsΨ ( †) =Ψ ( †) (r) obey bosonic canonical commutator relations, µ is the chemical potential, g 0 > 0 the repulsive s-wave interaction strength, and U (r) an external one-body potential that renders the gas inhomogeneous. As an application, we have in mind either a weak lattice or a random potential. In the latter case, meaningful quantities will involve the ensemble average (·). In order to describe the thermodynamic properties of the gas, one would like to know the ensemble-averaged grand potential (GP) Ω(β, µ), where β = 1/k B T is the inverse temperature: −βΩ = ln Ξ = ln{tr[exp(−βĤ gc )]}. (2) Ξ(β, µ) is known as the grand partition function. Other than on β and µ, the partition function and the Gibbs stateρ = Ξ −1 exp{−βĤ gc } depend also on all the parameters appearing in the Hamiltonian (1), such as the detailed configuration of the external potential U (r). The grand potential, on the other hand, only contains those properties that are relevant after the ensemble average. The advantage of this approach is that one obtains relevant physical quantities directly by differentiating the GP.
In particular, the average particle number is N (β, µ) = tr{ρN } = −∂Ω/∂µ. Often, one prefers to treat intensive quantities in the thermodynamic limit, such as the density N/V = n. The functional dependence n(β, µ) is known as the equation of state. By further differentiation, one has access to thermodynamic response functions, such as the inverse compressibility κ −1 = n 2 ∂µ/∂n.
Due to the interplay of interactions and external potential in (1), though, it is in general impossible to compute the partition function, let alone the GP, in closed form without further approximations. Here, we are interested in the thermodynamics of the Bose-condensed phase. Therefore, we resort to Bogoliubov's prescription Ψ (r) = Φ(r) + δΨ (r), where a macroscopically occupied condensate mode Φ(r) is separated from the quantum fluctuations δΨ (r). The condensate plays a role analogous to the classical trajectory in Feynman's pathintegral formulation of quantum mechanics: on the meanfield level, Φ(r) is an extremum, actually a minimum, of the functional (1) inside the trace (2). The minimization condition is known as the Gross-Pitaevskii or non-linear Schrödinger equation. Including contributions from the quantum fluctuations, one is later led to minimize more generally the grand-canonical energy, and thus finds beyond-mean-field corrections to the equation of state.
In this article, we explore the consequences brought about by quantum fluctuations in the presence of an external potential U (r). These corrections to the 'classical' mean-field solution can be computed by a quadratic expansion of the Hamiltonian (1) and subsequent Gaussian integration for the grand potential (2). We take advantage of the effective impurity-scattering Hamiltonian derived in [4] to take into account the external potential's effect on the condensate ("condensate deformation") as well as on the fluctuations. Specifically, we de-rive beyond-mean-field corrections to the particle density ("condensate depletion") in a disordered Bose fluid at finite temperature, thus complementing the zero-temperature results of Ref. [5]. Furthermore, we compare some of our findings on the mean-field level with recent measurements [6].
Before tackling the general, inhomogeneous case in Sec. 3, we find it instructive to first introduce our readers to the subtleties of grand-canonical Bogoliubov theory in the homogeneous case, treated in the following Sec. 2. Notably, it is shown how to recover the celebrated beyond-mean-field corrections to the equation of state first derived by Lee, Huang, and Yang [7].

Homogeneous case
In the homogeneous case U (r) = 0, condensation occurs in the k = 0 mode [8]. Therefore, we only need to determine the population N c of that mode, but not its shape. The Bogoliubov approximation consists in replacing the condensate field operator with a c-number,â 0 = N 1/2 c , and treating all other k-space modes as quantum fluctuations. In the following, we will first establish the effective Hamiltonian, then determine the GP, and analyze in detail the ground-state density. We close this section with a discussion of the condensate fraction at zero and finite temperature.

Hamiltonian
Expanding the HamiltonianĤ gc = H 0 +Ĥ 2 + . . . to second order in the fluctuations (Ĥ 1 vanishes by momentum conservation), one finds for the mean-field energy, where n c = |Φ c | 2 = N c /V is the condensate density. If one minimizes the meanfield energy alone, ∂H 0 /∂n c = 0, one finds g 0 n c = µ, and recovers canonical Bogoliubov theory [2,3]. Here, we postpone the minimization until the complete GP is known, in order to obtain beyond-mean-field corrections. The fluctuations are described by the Hamiltonian The primed sum indicates that k = 0 is omitted, and ε 0 k is the single-particle dispersion. 1 In order to avoid a UV divergence of the ground-state energy later on, 1 We note ε 0 k with a vector index to cover cases where the dispersion is anisotropic, e.g., in a tight-binding lattice [9].
For concrete examples within this paper, we consider only the free-space case, where ε 0 k = 2 k 2 /2m is isotropic. In all cases, we assume parity invariance, one renormalizes the interaction constant in (3) as g 0 = g + k g 2 /2ε 0 k V [3], which adds a c-number term under the sum in (4). The quadratic HamiltonianĤ 2 becomes diagonal after a transformation to the Bogoliubov quasi- defined in terms of All these quantities still depend separately on the chemical potential µ and the condensate population n c . Only with the choice gn c = µ, the energy (7) becomes purely real and gapless, as it should according to a theorem by Hugenholtz and Pines [10], and turns into the celebrated Bogoliubov dispersion relation, Yet, in order to be able to differentiate with respect to µ at fixed n c (or vice versa), we keep both quantities and remember to choose gn c = µ in all final expressions relating to the excitations. Thus, the grand-canonical Hamiltonian takes the form The first term is the grand-canonical candidate for the ground-state energy, to which fluctuations contribute with their commutators: (10) At this point, we still have the freedom to choose the condensate density n c by minimizing the energy (10) at fixed µ. Requiring ∂E 0 /∂n c | µ = 0 results in where we have introduced the characteristic length ξ via µ = 2 /(2mξ 2 ). Inserting this result in (10)-at mean-field precision inside the fluctuation sum-yields the ground-state energy density

Grand potential
The ground-state energy E 0 (µ) thus determined is a constant in the Hilbert space of Bogoliubov excitations and pulls out of the trace (2) for the GP, which evaluates in the thermodynamic limit to The density derives as where ν k = [e βε k − 1] −1 is the Bose-Einstein distribution function for the occupation of excitation modes. Remember that the second contribution, namely the thermal contribution of fluctuations, should be differentiated with respect to µ at fixed gn c , and then evaluated with gn c = µ at the end. Alternatively, it is also possible to derive the condensate density n c not from the ground-state energy (10) (i.e., at zero temperature), but by minimizing the full GP, eq. (13), as function of µ and arbitrary T . Thus, one is able to account for the thermal depletion of the condensate at fixed µ. The final results (as presented in the following) come out the same, as described in the Appendix A. But for technical reasons that will become apparent in Sec. 3 below, we prefer to use the 'semicanonical' prescription, in which n c is kept independent from µ when differentiating, and only substituted later at the required precision.

Zero-temperature equation of state
The density n = −V −1 ∂E 0 /∂µ derived from (12) thus determines the zero-temperature equation of state The difference between this total density and the condensate density (11) is the so-called quantum depletion, The depletion must be small compared to n (and thus n c ) in order for the Bogoliubov ansatz to hold. In this case we can express ξ = (8πna) −1/2 through the s-wave scattering length a and total density n c ≈ n itself, and recover the equivalent canonical expression [3, eq. (4.34)] Here, the dilute-gas parameter √ na 3 has come into play, which must be small for this correction to be meaningful. One can further derive the compressibility The corresponding speed of sound, determined by κ −1 = nmc 2 [11,12], includes the Lee-Huang-Yang correction. 2

Condensate fraction at zero and finite temperature
Often, one is interested in the condensate fraction n c /n as function of temperature and fixed total density. Equation (14) with the help of eq. (10) gives the well-known formula [3, eq. (4.42)] Figure 1 visualizes this result by showing (a) the integrand, or single-particle momentum distribution â † kâ k = v 2 k +(u 2 k +v 2 k )ν k , as function of reduced momentum kξ for different temperatures and (b) the resulting condensate fraction as function of temperature. Bogoliubov theory can be expected to give reasonably accurate results when the condensate fraction is large, i.e., for weak interaction and low temperatures.

Inhomogeneous (disordered) case
The presence of an external potential substantially complicates the situation, especially if U (r) is a random function. For a given realization, the bosons condense into a macroscopically populated eigenmode of the onebody density matrix [8], whose precise form is shaped by the interplay of kinetic, interaction, and potential energy.
In the spirit of Bogoliubov theory, one first needs to find the deformed condensate amplitude, given as a functional Φ(r) = Φ[U (r)] and depending of course also on µ and g. The total occupation number of this mode, is now larger than the occupation N 0 = |Φ 0 | 2 of the coherent mode k = 0 alone [16]. The inhomogeneous [5] or "glassy fraction" [17]. In a second step, one may then describe the quadratic fluctuations around this deformed condensate.
We are assured to find a well-defined set of elementary excitations whenever the external potential is weak enough not to fragment the condensate.
In this section, we first determine the condensate amplitudes for a given external potential on the mean-field level, and thus derive disorder corrections to the meanfield equation of state of the previous section. Our prediction for the resulting dependence of the compressibility on disorder strength compares rather well with recent measurements with ultracold molecules confined to 2D in presence of laser-speckle disorder [6].
In a second step, we put the quadratic Hamiltonian of the fluctuations to use and calculate their contribution to the GP. From there, we derive disorder corrections to the condensate depletion, recovering the zero-temperature results of [5] and extending them to finite temperatures.

Mean-field equation of state and compressibility
As before, we expand the grand-canonical Hamiltonian H gc = H 0 +Ĥ 2 + . . . up to second order in the fluctuations. On the mean-field level, the condensate amplitude minimizes the Gross-Pitaevskii functional, which reads in momentum representation and thus generalizes (3). For a weak external potential, whose smoothed Fourier components [18] are a set of small numbers, one can compute a pertur- Using this solution in (22), we find the ground-state GP By virtue of n c = −V −1 ∂Ω 0 /∂µ, or by inserting the perturbative solution (24)-(25) directly into (21), the ensemble-averaged equation of state becomes The first-order effect of the external potential is to shift the chemical potential by its mean value U (r) = U 0 . To second order, at fixed µ − U 0 , the external potential draws more particles into the condensate. We thus find the mean-field compressibility In all of the preceding expressions appears the paircorrelation function (with V = L d in d dimensions) which contains information about the strength of disorder, via the variance U 2 of on-site fluctuations. It also specifies spatial correlations, via the correlator C(σk). This function typically decays over a spatial correlation length σ, which is of the order of a micron in experiments involving laser speckle. Then, (28) can be written Thus, the compressibility is expected to decrease quadratically with increasing disorder strength U/µ at fixed correlation ratio σ/ξ. The compressibility can be measured with some precision in cold-atom experiments such as [6]. There, a quadratic decrease of the compressibility is measured for weak disorder, in quantitative agreement with (30), when evaluated in 2D with the correlation length σ ≈ ξ comparable to the healing length, as shown in Fig. 2. Since the 2D molecular BEC in the experiment is rather strongly interacting, with a depletion of order unity already without disorder, beyondmean-field corrections to the homogeneous compressibility κ κ (0) (1 − αU 2 ) therefore starts from the experimentally measured value forκ (0) . The agreement is satisfactory for disorder strengths U/µ not exceeding unity, as expected for a random potential whose correlation length is of the order of the healing length itself.
Comparing the two equations of state considered so far, (15) and (27), we see that there are two small parameters: the dilute-gas parameter √ na 3 , and the dimensionless disorder strength U 2 /µ 2 . Within the scope of this article, we are interested in their first-order effects. Therefore, we do not consider cross terms that would come from a higher-order solution of the ground state energy (the equivalent of (10) including contributions from the fluctuations), which is a somewhat ill-defined quantity in the presence of disorder anyway. 3 Also, we do not attempt to derive the condensate amplitude by minimizing the full GP including the fluctuations (to be described shortly), and thus forego a direct access to the thermal depletion of the condensate. But as the homogeneous case showed, we are allowed to use the 'semicanonical' method by keeping µ and condensate Φ(r) formally independent, and differentiating with respect to µ alone, inserting the mean-field solution Φ(r) to the required precision at the end. In the following, we take into account the first-order effect of disorder via a compensating shift in µ, and can assume without loss of generality that the potential has zero mean, U 0 = 0, such that only second-order corrections need to be discussed.

Quadratic Fluctuation Hamiltonian
Turning now to the quantum fluctuations, also affectionately called "bogolons", we first must ensure that they live in the space orthogonal to the condensate [21]. This constraint can be respected in the density-phase parametrization (see also [22] for a number-conserving approach, and [23,24] for the connection between both approaches), where the excitations are given by [4,5] The transformation matrices contain the Fourier coefficients Φ k of the condensate amplitude Φ(r) and its inverseΦ(r) = n c /Φ(r), which encode the dependence on the external potential U (r) or rather its Fourier components U k , as described in (24) and (25). In the absence of an external potential, the condensate Φ (0) k = φ 0 δ k0 renders this transformation diagonal in k, and by choosing a k = (ε k /ε k ) 1/2 , one recovers the homogeneous transformation (5). Now, we seek the effective quadratic Hamiltonian that describes these excitations. As the condensate mode Φ minimizes H 0 , the linear term in the expansion vanishes and the relevant term is the second-order fluctuation HamiltonianĤ whereĤ (0) 2 = k ε kγ † kγ k formally looks like the freespace contribution-but please be reminded that the excitations defined via (31) are not the plane-wave modes of the homogeneous case. Furthermore, we recall that µ = gn c has to be taken in expressions relating to the fluctuations at the end, thus ensuring a real, gapless excitation spectrum. And finally, we have discarded the zero-point contribution of commutators, which would result, as explained above, in a beyond-mean-field modification of the ground-state energy, which is not investigated here. More importantly, the spatial inhomogeneity leads to the appearance of the scattering potential The impurity-scattering matrices can be traced back to the terms h kk = (ε 0 k − µ)δ kk + U k−k + 2gn ck−k and gn ck−k in the inhomogeneous generalization of eq. (4): In contrast to the otherwise equivalent formulas in Refs. [4,5,9], we keep the chemical potential µ and the condensate mode [Φ(r) = n c (r) = n c /Φ(r)] separately here, in order to be able to take partial derivatives with respect to µ only. By virtue of the perturbative expansion (24) and (25), the effective scattering potentialÛ = ∞ n=1Û (n) can be expanded in powers of the external potential strength, U/µ, up to the desired order.

Fluctuation Grand Potential
The GP of fluctuations, described by the quadratic Hamiltonian (34), can be split into the sum of two terms: Indeed, the homogeneous contribution Ω can now be computed straightforwardly by perturbation theory in powers of the effective scattering potential (35). Taking the logarithm leaves us with the connected correlations, which up to second order read With the help of Wick's theorem, all correlations can be expressed by Matsubara Green functions connecting the matrix elements (36) and (37) ofÛ . Expanding these in turn to second order in the external potential, we find This expression for the disorder-induced correction to the GP of quantum fluctuations is the central result of this article. Let us emphasize that our approach, starting from the deformed condensate background, takes into account all contributions that are of second order in U (and thus goes beyond the method of Huang and Meng [26,27]). Furthermore, this result only involves elementary perturbation theory and does not rely on the replica method. Its obvious drawback is its perturbative nature. We therefore do not make any claims concerning the strong-disorder regime, nor do we cover high temperatures or strongly interacting regimes, where interactions between excitations become important (see [20,28]). Here, we rather wish to provide an account as complete as possible of disorder-effects up to second order in U/µ at low temperatures, where Bogoliubov theory applies.

Condensate depletion
We are now in the position to compute the particle number shift δN 2 = −∂δΩ (2) 2 /∂µ due to the disorder. When this number is compared to the number of particles in the condensate, N c , we find the additional condensate depletion caused by the inhomogeneous potential, or "potential depletion" for short. In the partial derivative of (44) with respect to −µ at fixed n c , we set µ ≈ gn c ≈ gn in the end because the expression is already of first order in the dilute-gas parameter and second order in the disorder strength. We find the following collection of identities helpful: −∂ε k /∂µ = 1, as well as withň q = [n 2 c /n c (r)] q . Applying these to Eqs. (36) and (37), one finds Via (36)-(39) and n k = V −1 k Φ k−k Φ k , we express the right hand sides in terms of the perturbative solution of the Gross-Pitaevskii equation (24)- (25). The relevant expressions up to second order read and withw (1) kp andỹ (1) kp as given in Ref. [5] and The term λ kq = (ε 0 k−q +ε 0 k+q −2ε 0 k −2ε 0 q )/2ε 0 k vanishes in the present case of a quadratic dispersion relation ε 0 k ∝ k 2 , but is nonzero in the case of a lattice potential [9].
We can write down the potential depletion as kk+q defined in terms of the envelopẽ where ν k = ∂ν k /∂ε k | µ=gnc = β(ν k + ν 2 k ). Equation (53) leaves a certain freedom to exchange p and k in the individual components ofM (2) kp . In the spirit of Ref. [5], we have used this freedom to write down eq. (54) in a way that allows the identification of δn with the momentum distribution of the condensate depletion induced by the disorder. Now we evaluate eq. (53) in the case of a disorder potential with strength U and correlation length σ: Results are shown in Fig. 3(b) as a function of the disorder correlation length for different temperature values. In all cases there is an increase of the depletion due to disorder. At first sight surprisingly, this increase diminishes with temperature in the regime of uncorrelated disorder σ ξ. However, the disorder correction as shown in Fig. 3(b) is expressed in units of the homogeneous quantum depletion at zero temperature (16). Thus, it goes on top of the thermal depletion discussed in Fig. 1, such that both temperature and disorder deplete the condensate. We have found [Figs. 1(b) and 3(b)] that for temperatures up to k B T µ and for not too strong disorder U µ, the depletion δn (0) + δn (2) remains of the same order of magnitude as the zero-temperature homogeneous depletion δn 0 , which a posteriori validates the Bogoliubov method.
Finally, we note that the bulky expression for the envelope function (54) can be simplified significantly in the Thomas-Fermi regime σ ξ, where the condensate profile can faithfully follow the variations of the disorder potential on the length scale σ. In this case, the disorder correlation (55) tends to a Dirac δ-function and the potential depletion δn (2) is dominated by the diagonal elementsM (2) kk , which are given as Summing over k then gives G(0) [left edge of Fig. 3(a)], which is proportional to the depletion in the Thomas-Fermi limit σ/ξ → ∞ [right edge of Fig. 3(b)].

Connection to the canonical frame
At a given value of the chemical potential, the disorder potential draws more particles into the condensate [eq. (27)]. In Refs. [4,5] the canonical frame is used, where this effect is compensated by a shift ∆µ = − q ε 0 q |Ũ q | 2 of the chemical potential [4], which results in different second-order expressions (49) and (52). In particular, the expression for the potential depletion given in Ref. [5,Eqs. (48) and (49)] differs slightly from the one given here in Eqs. (53) and (54), even at zero temperature. Equation (54) goes over to its canonical form at fixed n c by replacingw (2) kp andỹ (2) kp with their canonical expressions (given in Ref. [5]) and by dropping the term ε 0 k−p /µ in the second line of eq. (54). In fact, the difference between the two frames can be written as δn (2) canonical − δn (2) gc = ∆µ where is the homogeneous depletion at finite temperature. With this shift, the kernel functions shown in Fig. 3(a) would begin with the same values at q = 0 and then decrease monotonically as function of q without crossing each other. Likewise, the disorder-induced depletion of Fig. 3(b) would become monotonic without crossings; the zero temperature curve takes the same form as in Figure 4 of [5].

Conclusions
We have applied the grand-canonical formulation to the problem of disordered Bose-Einstein condensates, which brings conceptual advantages over the conventional canonical frame. Once the grand potential is determined, one obtains relevant physical quantities by differentiation. The condensate mode Φ(r) plays a special role in the grand-canonical Bogoliubov approach. In principle, it has to minimize the grand potential, which includes a back action of the excitations on the condensate. For the main work of this article, we have chosen the equivalent 'semicanonical' approach, where one keeps the condensate as a parameter that is inserted only at the end of the calculation (i.e., after taking derivatives). To the desired precision, it is then sufficient to determine the condensate mode by minimizing the ground-state energy.
Concerning physical results, we have mainly focused on the speed of sound, the compressibility, as well as on the particle fractions condensate fraction and condensate depletion. In particular, we have reproduced previous results [5] on the disorder-induced condensate depletion from the perspective of the grand-canonical picture and have extended them to the case of finite temperatures.
A Grand-canonical condensate density with beyond-mean-field corrections In Sec. 2, we have determined the condensate density n c by minimizing the ground-state energy E 0 at fixed µ. This amounts to determining the condensate density, once and for all, at zero temperature. When the temperature is raised, then of course thermal excitations will appear, which deplete the condensate. This effect can be explicitly accounted for by determining n c directly from the GP Ω at finite temperature. Then, using this (now µ and T dependent) solution, one has a GP Ω(µ, T ) that depends only on µ, and not separately on gn c anymore. The total density then derives by differentiation with respect to this µ alone. This proper grand-canonical procedure yields the same results than the 'semicanonical' method used in Sec. 2 above, as demonstrated in the following.
Requiring that the homogeneous GP, eq. (13), be stationary, ∂Ω/∂n c | µ = 0, yields the condensate density with now a T -dependent contribution. Inserting this solution into (10) yields actually the same ground-state energy (12) as before. The reason is that the beyondmean-field correction n c = (µ/g) + ∆n c does not contribute there to lowest order, since this correction is only used in the mean-field term H 0 , eq. (3), for which (gn c − 2µ)gn c = −(µ − g∆n c )(µ + g∆n c ) = −µ 2 (59) to the order considered. Thus, the ground-state energy E 0 (µ) is unchanged, just as the GP, eq. (13). The difference now is that the excitation energy inside the fluctuations is to be taken at the Bogoliubov dispersion ε B k from eq. (8) as function of µ alone. Thus, the total number of particles (for the same µ) is now different, namely, with ∂ε B k /∂µ = ε 0 k /ε B k . However, also the condensate density (58) is now temperature dependent. The difference between these two densities is the depleted density with the zero-temperature depletion δn 0 given by eq. (16). As for the thermal depletion, the dispersion relation is such that the difference of derivatives appearing there is precisely the result we had before: Thus, (20) still holds as before, and all zero-T results are identical anyway. It is largely a matter of taste whether one wants to have a T -dependent contribution to n c or not, and whether one wants to have the GP depend really on µ alone. Both approaches, strict grand canonical and 'semicanonical', are equivalent.