Condensation vs Ordering: From the Spherical Models to BEC in the Canonical and Grand Canonical Ensemble

In this paper we take a fresh look at the long standing issue of the nature of macroscopic density fluctuations in the grand canonical treatment of the Bose-Einstein condensation (BEC). Exploiting the close analogy between the spherical and mean-spherical models of magnetism with the canonical and grand canonical treatment of the ideal Bose gas, we show that BEC stands for different phenomena in the two ensembles: an ordering transition of the type familiar from ferromagnetism in the canonical ensemble and condensation of fluctuations, i.e. growth of macroscopic fluctuations in a single degree of freedom, without ordering, in the grand canonical case. We further clarify that this is a manifestation of nonequivalence of the ensembles, due to the existence of long range correlations in the grand canonical one. Our results shed new light on the recent experimental realization of BEC in a photon gas, suggesting that the observed BEC when prepared under grand canonical conditions is an instance of condensation of fluctuations.


I. INTRODUCTION
In general statistical ensembles are constructed to be equivalent in the thermodynamic limit, but there are exceptions to this rule. This paper deals with phenomena arising when this equivalence breaks down. Pairs of conjugate ensembles are obtained by controlling the system either by fixing the value of some extensive quantity, through appropriate isolating walls, or by putting it in contact with a reservoir of that same quantity. A familiar example, which will be of central interest in the following, is that of the canonical and grand canonical pair resulting from fixing either the density of particles or the chemical potential, while keeping the system thermalized.
Basically, equivalence holds in situations in which correlations are short ranged. Then, the central limit theorem guarantees that fluctuations of extensive quantities become negligible in the thermodynamic limit, so that it doesn't matter whether the system is controlled by enforcing a rigid constraint or through the contact with a reservoir [1][2][3]. By the same token lack of equivalence is to be expected when correlations are long ranged. This is a more rare occurrence, but very interesting since new physics obtains by switching from one ensemble to the other within a conjugate pair. Best known and recently much studied is the case of systems with long-range interactions [4].
There is one instance of nonequivalence which stands apart: The one which materializes as an ideal Bose gas (IBG) is driven through the Bose-Einstein condensation (BEC). In the canonical ensemble (CE) fluctuations of the condensate behave normally while in the grand canonical ensemble (GCE) do persist even in the thermodynamic limit [2,5]. Although this is an exact result, it is somewhat puzzling because, dealing with an ideal gas, it is not at all obvious where the long range correlations responsible of the nonequivalence could come from. An unbiased attitude ought to advise to take the facts at face value and to inquiry on possible different mechanisms underlying BEC in the two ensembles. Instead, due to a widespread aversion to macroscopic fluctuations of an extensive quantity, which are not suppressed by lowering the temperature, the GCE result has been variously regarded as unacceptable [6], unphysical [5,[7][8][9] or even wrong [10] and is commonly referred to as the grand canonical catastrophe.
The need to reconsider afresh this matter has been prompted by the recent observation of BEC in the lab [11,12] in a gas of photons under grand canonical conditions, which has changed the outlook by producing evidence for the existence of the macroscopic fluctuations of the condensate. Therefore, after reckoning with the absence of any catastrophe, the challenge is to uncover the mechanism responsible of the nonequivalence. Due to the fundamental character of the question posed, we shall leave the experiment in the background and we shall explore the basic issues in the simplest possible context of the uniform IBG in a box of volume V , aiming primarily to outline the conceptual framework needed to approach this interesting and multifaceted problem.

II. THE PROBLEM
At the phenomenological level the mechanism of BEC appears to be the same in the CE and in the GCE. Denoting by d, d * and d 0 the total density, the density in the excited states and the density in the ground state, respectively, from the obvious identity d = d * + d 0 follows the sum rule which must be satisfied by the average quantities irrespective of the choice of the ensemble where ρ stands for d and the brackets for the average over either ensemble. The condensate density d 0 is called the BEC order parameter. Now, for space dimensionality d > 2 and in the thermodynamic limit, d * is superiorly bounded by a finite critical value ρ c [1,2]. Consequently, keeping T fixed and using ρ as control parameter, from Eq. (1) immediately follows the densitydriven BEC which, we emphasize, holds irrespective of the ensemble. Thus, as far as d 0 is concerned, CE and GCE are equivalent. However, a striking difference between the two emerges when the fluctuations of d 0 are considered, since in the condensed phase, as previously anticipated, one has [5] i.e. normal behavior in the CE and macroscopic fluctuations in the GCE. The crux of the matter is that at this level of observation no insight can be obtained as to the why fluctuations ought to behave so differently in the two ensembles. The point of view that we propose in this paper is that the picture is rationalized by shifting the description to the finer and underlying level of the field-theoretic microscopic degrees of freedom, which however are not directly observable. In order to clarify the interplay of the different levels of description, in the next paragraph we shall exploit the analogy with magnetic systems, where a quite similar and well understood situation arises.

III. SPHERICAL AND MEAN-SPHERICAL MODEL
The IBG in the CE and in the GCE is well known [13][14][15] to be closely related to the spherical and the mean-spherical models of magnetism. Let H(ϕ) = V d r ϕ( r) − 1 2 ∇ 2 ϕ( r) be the energy function of a classical scalar paramagnet [16] in the volume V , where ϕ stands for a configuration of the local unbounded spin variable ϕ( r). Due to its bilinear character the above Hamiltonian can be diagonalized by Fourier transform H = 1 V k k 2 |ϕ k | 2 . In the spherical model (SM) of Berlin and Kac [17] a coupling among the modes is induced by the imposition of an overall constraint on the square magnetization S(ϕ) = V d r ϕ 2 ( r) = 1 V k |ϕ k | 2 . Then, in thermal equilibrium the statistical ensemble reads where Z SM is the partition function, s(ϕ) = 1 V S(ϕ) the square magnetization density and σ a positive number, which usually is set σ = 1, but here will be kept free to vary as a control parameter. In the mean-spherical model (MSM) [13,18] the constraint is imposed in the mean: An S-dependent exponential bias is introduced in place of the δ function and the intensive parameter κ conjugate to S must be adjusted so as to satisfy the requirement Although it is the common usage to refer to these as models, it should be clear from Eqs. (4) and (5) that we are dealing with two conjugate ensembles, distinguished by conserving or letting to fluctuate the density s. Separating the excitations from the ground state contribution s = s * + s 0 , where s * = 1 V 2 k =0 |ϕ k | 2 and s 0 = 1 V 2 ϕ 2 0 , taking the average and using the constraint s = σ, independently from the choice of the model there follows the sum rule analogous to Eq. (1) Therefore, the variables s, s * , s 0 and σ do correspond to the IBG ones d, d * , d 0 and ρ, with the important difference that in the present context these are composite variables, built in terms of the microscopic set of the magnetization components [ϕ k ]. Furthermore, also in this case for d > 2 and in the thermodynamic limit the excitation contribution s * is superiorly bounded by a finite critical value σ c , see Appendices A and B for details. Hence, keeping T fixed and varying σ, from Eq. (7) there follows showing that s 0 behaves like the BEC order parameter and that, as far as s 0 is concerned, the two models are equivalent. However, at the microscopic level a different scenario opens up, since there is no unique way to form a finite expectation s 0 . Let us introduce the probability that s takes the value σ ′ in the MSM, given by K MSM (σ ′ |σ) = dϕ P MSM (ϕ|σ)δ (s(ϕ) − σ ′ ). Then, just as a consequence of definitions, the distributions (4) and (5) are related by The kernel has been worked out by Kac and Thompson [13], obtaining K MSM (σ ′ |σ) = δ(σ ′ − σ) for σ < σ c , which implies that the two distributions coincide and, therefore, that the two models are equivalent below σ c . Conversely, when σ is above σ c the kernel vanishes for σ ′ < σ c , while for σ ′ > σ c is of the spread out form revealing nonequivalence. In the following we shall restrict to the σ > σ c domain, where nontrivial behavior is expected. Integrating out the k = 0 modes from Eq. (9), an identical relation between the marginal probabilities of ψ 0 = 1 V ϕ 0 is obtained. In the left hand side there appears the Gaussian distribution P MSM (ψ 0 |σ) ∝ exp{−βκV ψ 2 0 /2}, as it can be verified by inspection from Eq. (5), since P MSM (ϕ|σ) factorizes in Fourier space. From this follows s 0 = (βκV ) −1 . So, from the second line of Eq. (8) we get which implies Hence, plugging in the explicit expression of K MSM (σ ′ |σ) it is not difficult to verify that Eq. (9) is satisfied by the ansatz where m ± = ± √ σ − σ c is the spontaneous magnetization density which would be obtained, for instance, by switching off an external magnetic field [13,17]. Thus, we have two quite different distributions, as it is clearly illustrated by the plots in Fig. 1.
We are now in the position to draw some conclusions. The BEC-like order parameter s 0 can be computed microscopically as the average composite variable ψ 2 0 . Then, it is straightforward to check from Eqs. (12) and (13) that Eq. (8) is satisfied in both cases. However, it is enough to take a look at Fig. 1 to realize that the numerically identical result ψ 2 0 = (σ − σ c ) for σ > σ c in the two models stands for two different phenomena. The double peaked distribution of the SM case is the familiar one for a ferromagnet in the magnetized phase, each peak being associated to a pure state and with the up-down symmetry of the model spontaneously broken. Namely, the distribution is the even mixture of these two pure states. This means that in the SM the BEC-like transition observed at the level of s 0 is the manifestation of an underlying ordering transition, and that the BEC order parameter is the square of the spontaneous magnetization, i.e. ψ 2 0 = m 2 ± . By contrast, in the MSM case we have the opposite situation, since ψ 2 0 is the variance of a broad Gaussian distribution centered on the origin. Therefore there is no ordering, no breaking of the symmetry. In this case the BEC-like transition undergone by s 0 is the manifestation of the microscopic variable ψ 0 developing finite fluctuations. The reason for this can be grasped intuitively. In the SM, due to the sharp constraint, there is enough nonlinearity to produce ordering. In the MSM framework this cannot be achieved, since the statistics are Gaussian. Then, the only mean to build up the finite value of s 0 needed to saturate the sum rule (7) above σ c is by growing fluctuations in the single degree of freedom ψ 0 . Elsewhere [19][20][21][22][23], this type of transition, characterized by the fluctuations of an extensive quantity condensing into one microscopic component, has been referred to as condensation of fluctuations. The phenomenological picture is completed by the fluctuations of s 0 itself which, as it follows easily from Eqs. (12) and (13), are given by Comparing this with Eq. (3), the analogy is evident. However, now no catastrophical behavior can be envisaged, because the fluctuations of s 0 are trivially a consequence of the different microscopic statistics in the two models.
Having analysed how the nonequivalence unfolds, the remaining task is to clarify where it does to originate from, which ultimately must be in the presence of long range correlations. The explanation is that in the MSM the parameter κ is related to the correlation length ξ by κ = ξ −2 [16] and from Eq. (11) we see that in the thermodynamic limit κ vanishes like 1/V when σ is fixed above σ c . Therefore, in the entire condensed phase the MSM is critical, while the SM is not. Hence, the lack of equivalence. We emphasize that the onset of these critical correlations in the MSM is the unifying thread behind the BEC-like transition accompanied by macroscopic fluctuations of s 0 .

IV. BACK TO THE IBG
We may now go back to the main topic of the IBG with the advantage of hindsight, since we know what to look for: The microscopic variables underlying the phenomenological level, in terms of which we expect to expose both the different mechanisms of BEC in the CE and GCE and the nonequivalence cause. This is accomplished by introducing the creation and destruction operators and by using the representation of the density matrix in the associated coherent state basis. Let us first diagonalize the energy and number operators by where ǫ k is the single particle energy. In the Glauber-Sudarshan P-representation [24,25] the density matrix is given by D(ρ) = d 2 α P (α|ρ)|α α|, where |α are product states k |α k and the k-mode factor |α k is the eigenvector of the annihilation operator a k |α k = α k |α k with complex eigenvalue α k . Then, in the GCE the weight function reads [24] where n k = [e β(ǫ k −µ) − 1] −1 is the usual Bose average occupation number of the state | k [1] and µ stands for the chemical potential. Using the identity |α k | 2 = n k , which easily follows from Eq. (15), the equation fixing µ for the given value of ρ reads 1 V k |α k | 2 = ρ. Since this is nothing but Eq. (1), we may write d = 1 V k |α k | 2 and, consequently, d 0 = |η 0 | 2 , after setting η 0 = 1 √ V α 0 . This allows to identify [α k ] with the sought for set of microscopic variables analogous to [ϕ k ]. Following the magnetic example, we must focus on the statistics of the zero component, keeping in mind however that now is a complex quantity η 0 = |η 0 |e iθ . The starting point is the relation between ensembles analogous to Eq. (9) (see Appendices C and D) where K GCE (ρ ′ |ρ) is the probability in the CGE that the density takes the value ρ ′ . This is known as the Kac function [5], whose form is similar to that of Eq. (10). Namely, K GCE (ρ ′ |ρ) = δ(ρ ′ − ρ) for ρ < ρ c , while when ρ is above ρ c it vanishes for ρ ′ < ρ c and is given by The relation between the η 0 marginal distributions is then obtained by integrating out the k = 0 modes . Inserting in the left hand side the contribution from the first factor of Eq. (15) and substituting for K GCE (ρ ′ |ρ) the above expression, the equation is solved by the ansatz From the plot of the two distributions (18) and (19) in Fig. 2, we see by inspection that we are confronted with a situation qualitatively similar to the one in the magnetic case. In both ensembles |η 0 | develops a nonobservable finite expectation value The observable BEC order parameter exhibits the same numerical value as in Eq. (2) which is achieved through fluctuations in the GCE and without fluctuations in the CE This means that the sum rule (1) in the CE is saturated by fixing the modulus to the precise finite value |η 0 | = √ ρ − ρ c . Because of this freezing of |η 0 |, BEC in the CE fits into the scheme of an ordering transition akin to the ferromagnetic transition in the SM. Conversely, BEC in the GCE does not take place through ordering.
Rather, the saturation of the sum rule is achieved by growing the macroscopic fluctuations of |η 0 |, as Eq. (22) shows. Therefore, in this case BEC fits into the scheme of the condensation transition. Ordering is ruled out because the width of the probability distribution persists in the thermodynamic limit. Moreover, assuming ǫ k ∼ k α , where the power α depends on the dispersion relation, for small k and small µ we may approximate n k −1 ≃ β[k α − µ] and inserting this into Eq. (15) we have that the chemical potential, like κ in the preceding case, is connected to the correlation length by −µ = ξ −α . Since the formation of the condensed phase in the GCE requires µ to vanish in the thermodynamic limit [1], we have that the condensed phase is critical throughout in the GCE but not in the CE. This explains the origin of nonequivalence which, as in the magnetic case, is not revealed by the BEC order parameter but emerges only at the level of the higher cumulant Hence, the phenomenological result of Eq. (3), rather than being pathological, is now accounted for as a byproduct of the critical correlations in the GCE.

V. SUMMARY
In this paper we have investigated the differences arising when BEC in a homogeneous IBG is treated in the CE and in the GCE. The analysis has been carried out by taking advantage of the close analogy with the the spherical and mean spherical models of magnetism. The problem is of particular interest because the ensemble nonequivalence issue encroaches the fundamental question of the nature of BEC. We have shown that ordering takes place in the CE, while condensation takes place in the GCE, whose prominent manifestation are the macroscopic fluctuations of the condensate. Therefore, we suggest that the recent experimental realization of BEC in a gas of photons [11,12,26] ought to be regarded as qualitatively different from other experimental instances of BEC, such as those with cold atoms, precisely because the grand canonical conditions lead to BEC as condensation of fluctuations. Moreover, by retracing the origin of nonequivalence to the onset of critical correlations in the condensed phase of the GCE, we have pointed out that the observable phenomenology follows as a consequence. So, knowledge of the existence of these correlations could possibly serve as a useful guide in the planning of future experiments. As a final remark, notice that the above analysis has involved the modulus, but not the phase of η 0 . This means that the distinction between ordering and condensation is decoupled from the issue of the breaking of the gauge symmetry. This is a separate and important problem which will be the object of future work.

Acknowledgments
Informative and quite useful conversations on photons BEC with Prof. Claudio Conti are gratefully acknowledged. AS acknowledges support from the "Programma VALERE" of University of Campania "L. Vanvitelli" and from MIUR PRIN project "Coarse-grained description for non-equilibrium systems and transport phenomena (CO-NEST)" n. 201798CZLJ.

Appendix A: Spherical and Mean Spherical Model
Let ϕ be a configuration of the magnetization field ϕ( r) ∈ (−∞, +∞) over an hypercube r ∈ V ⊂ R d of side L, whose energy is given by Ensembles, or models, are defined by thermalizing the system and specifying conditions imposed on the overall square magnetization The spherical model (SM) of Berlin and Kac [17], which corresponds to the ensemble canonical with respect to S(ϕ), is obtained by imposing a sharp constraint on the sqaure magnetization density where s(ϕ) = 1 V S(ϕ) and σ is a positive number. This leads to the probability distribution with the partition function The mean spherical model (MSM) of Lewis and Wannier [18], which corresponds to the ensemble grand canonical with respect to S(ϕ), is defined by with the partition function and where the parameter κ is determined self-consistently by imposing the constraint (A3) on average as it is explained in the next section. Notice that in both models β and σ are control parameters.
Imposing periodic boundary conditions and postulating the existence of a microscopic length a 0 , the allowed wave vectors of the Fourier components are given by where N max = L/a 0 supposedly is an integer. The inverse transform reads Using and the energy function (A1) and S(ϕ) take diagonal forms where we have used the reality of ϕ( r) and Eq. (B1), which imply ϕ − k = ϕ * k . The partition function of the MSM can be computed straightforwardly from Eq. (A7) , Re κ > 0, (B7) where N = k 1 is the total number of modes. Using this expression and the mean spherical constraint (A8) reads where we have separated the k = 0 contribution from the rest. Then, for V sufficiently large, the sum over k = 0 can be replaced by an integral and Eq. (B9) can be recast as with the integration measure defined by dµ is the d-dimensional solid angle and Λ ∼ 1/a 0 is a cutoff. The function B(κ) is a monotonic decreasing function of κ, which diverges at κ = 0 for d ≤ 2, while its maximum value at κ = 0 for d > 2 is given by Therefore, for d ≤ 2 the first term in the right hand side of Eq. (B11) can be neglected for any choice of β and σ, and the solution is given by where B −1 is the inverse function of B(κ). Instead, for d > 2, the condition βσ = B(0) defines a critical line on the (β, σ) plane below which the solution is still given by Eq. (B13), while above it is necessary to retain also the first term in the right hand side of Eq. (B11). Thus, keeping β fixed, the critical value of σ is given by and the full solution of Eq. (B11) reads where B 1 is a positive constant. As a consequence of the mode independence, implied by Eq. (B6), the P MSM (ϕ|σ) distribution factorizes. Therefore, introducing the notation ψ 0 = 1 V ϕ 0 for the magnetization density, the k = 0 contribution is given by and inserting the result (B15) for κ, in the V → ∞ limit the result reported in the main text is obtained , for σ > σ c .

(B17)
Appendix C: The SM -MSM connection Using the identity ∞ 0 dσ ′ δ(s(ϕ) − σ ′ ) = 1, the MSM ensemble (A6) can be rewritten as which, using the definition (A5) of the SM partition function, can be further manipulated as Recognizing that in the square bracket there appears P SM (ϕ|σ ′ ), while the last integral (C3) is the probability that s(ϕ) takes the value σ ′ in the MSM parametrized by σ, the probabilities of ϕ in the two models are related by Eliminating the ϕ k =0 components from the above equation by integration, a similar relation is obtained between the ψ 0 marginal distributions The kernel K MSM (σ ′ |σ) has been computed by Kac and Thompson [13] obtaining and for σ > σ c , for σ ′ > σ c .
(C7) Therefore, using the above result together with Eq. (12), one can check that for σ < σ c Eq. (C5) is solved by while for σ > σ c one gets (C10) Appendix D: The CE -GCE connection in the ideal Bose gas The Fock space representation of the Hamiltonian H = k ǫ k a † k a k and number operator N = k a † k a k is given by where n stands for a collection [n k ] of occupation numbers, |n are the product states k |n k and the eigenvalues are given by It is first convenient to lay out the general relation between the Fock-space and the Glauber-Sudarshan representations [24,25] of the density matrix D = n P (n)|n n| = d 2 α P (α)|α α|, where |α = k |α k . The k-mode coherent state |α k is an eigenvector of the annihilation operator a k |α k = α k |α k with complex eigenvalue α k . The weight functions are related by P (n) = d 2 α R(n|α)P (α), with the kernel The canonical and the grand canonical ensembles are obtained by controlling the density of particles ρ either strictly or on average. Then, the corresponding density matrices are given by where the weight functions read In the latter one the chemical potential µ is fixed by the condition By going through the same algebra as in the preceding section it is straightforward to verify that these weights obey the relation analogous to Eq. (C4) where K GCE (ρ ′ |ρ) is the probability that the particle density takes the value ρ ′ in the GCE controlled by the average value ρ, This is given by the Kac function [5], which for ρ < ρ c reads K GCE (ρ ′ |ρ) = δ(ρ ′ − ρ), for ρ < ρ c , while for ρ > ρ c is given by , for ρ ′ > ρ c . (D16) Next, inserting Eq. (D6) into Eq. (D13) and taking into account that the kernel R(n|α) is positive definite, the analogous relation is found to hold between the Prepresentation weight functions Integrating this over all α k =0 , eventually one finds where the left hand side term is given by [24] P GCE (α 0 |ρ) = 1 π n 0 e − |α 0 | 2 n 0 , V (ρ − ρ c ), for ρ > ρ c , (D20) having denoted by µ(ρ) the solution of Eq. (D12) with respect to µ. Therefore, defining η 0 = α 0 / √ V and inserting Eqs. (D15,D16,D18,D19) into Eq. (D17), it is easy to verify that in the large V limit for ρ < ρ c P GCE (|η 0 ||ρ) = P CE (|η 0 ||ρ) = δ(|η 0 |), while for ρ > ρ c and P CE (|η 0 ||ρ) = 1 π δ(|η 0 | 2 − (ρ − ρ c )). (D23)