Thermodynamic instability and first-order phase transition in an ideal Bose gas

We conduct a rigorous investigation into the thermodynamic instability of ideal Bose gas confined in a cubic box, without assuming thermodynamic limit nor continuous approximation. Based on the exact expression of canonical partition function, we perform numerical computations up to the number of particles one million. We report that if the number of particles is equal to or greater than a certain critical value, which turns out to be 7616, the ideal Bose gas subject to Dirichlet boundary condition reveals a thermodynamic instability. Accordingly we demonstrate - for the first time - that, a system consisting of finite number of particles can exhibit a discontinuous phase transition featuring a genuine mathematical singularity, provided we keep not volume but pressure constant. The specific number, 7616 can be regarded as a characteristic number of 'cube' that is the geometric shape of the box.


Introduction
By definition, first-order phase transitions in thermodynamics feature a genuine mathematical singularity. Whether finite systems in Nature can literally exhibit such an infinity is a long standing controversial question [1,2]. When a thermodynamic system is composed of a definite number of particles, say N , and is in contact with a heat reservoir, the key quantity is the canonical partition function: where the sum is over all the quantum states of the N -body system. When the energy eigenvalues depend on the volume of the system, Z N is a function of the volume V , and the temperature through β = 1/(k B T ), where k B denotes Boltzmann constant. Provided the precise canonical partition function, we may compute various physical quantities, which include the pressure, the entropy, the internal energy and the specific heat per particle at constant volume as follows: In particular, ∂ 2 β ln Z N = (E ψ − E) 2 being a standard deviation squared, C V is finite and never negative. Further, the temperature derivative of the probability for the system to occupy a certain quantum state ψ reads As the temperature increases from absolute zero to infinity, the corresponding probability increases if the energy eigenvalue is greater than the average i.e. E ψ > E, and it starts to decrease in the opposite case, E ψ < E.
Since the canonical partition function appears as a positive definite analytical function of its arguments, clearly all the physical quantities listed in (2) and (3) should not feature any singularities. They may do so only in the thermodynamic limit: the limit of N → ∞ and V → ∞ with N/V held fixed [3].
However, strictly speaking, infinite limits are hardly realistic and exist only in theory [1,2,4,5]. The above analysis seems to suggest that Nature does not admit a discontinuous phase transition featuring a genuine mathematical singularity, which is somewhat different from experiments or our daily experiences under the standard pressure 1 atm.
In this paper we pay attention to, among others, the fact that the above finiteness and continuity are for the cases of keeping the volume fixed. Once we switch to an alternative constraint of keeping the pressure constant, we demonstrate, for the first time, that canonical ensembles with finite number of physical degrees may undergo a discontinuous phase transition.
The organization of the present paper is as follows: The final section 4 conveys our discussion. In particular, we comment on the similarity between the permutation symmetry of the identical particle indices and the gauge symmetry in high energy physics.
Appendix contains our numerical verification that ideal Bose or Boltzmann gases under periodic or Neumann boundary conditions exhibit a thermodynamic instability at low temperature near absolute zero.
Although ideal Bose gas has been studied for decades [13][14][15][16][17] and discussed in many textbooks [6,[18][19][20], the implication of the constant pressure constraint to the canonical ensemble of finite N has been rarely explored. 1 To the best of our knowledge the finite N -effect on the canonical ensemble has been addressed only in the case of keeping the volume fixed, rather recently by Kleinert [23] and by Glaum, Kleinert and Pelster [24]. The earlier focus was typically on either grand canonical or micro-canonical ensembles, and the computations often assumed a continuous approximation to convert discrete sums to integrals [25,26], unless an external harmonic potential sets the sum to be taken over a geometric series [27][28][29][30][31][32][33][34][35][36]. 1 The so-called constant pressure ensemble [21,22] fixes the pressure and allows the volume to fluctuate, as its partition function is given by However, just like CV = (kB/N )β 2 ∂ 2 β ln ZN , the specific heat therein is positive definite and finite: As we are interested in the precise change of the volume without allowing any fluctuation of it, in the present paper we focus on the standard canonical ensemble.

First-order phase transitions in finite systems
Here we explain how canonical ensembles with finite number of physical degrees may exhibit a discontinuous phase transition featuring a genuine mathematical singularity when we keep not volume but pressure fixed.
Prior to rigorous analysis, a thought experiment may provide an intuitive clue for this claim: If we fill a rigid box with water to the full capacity and heat it, the temperature will increase but hardly it evaporates. However, once we open the lid, simply it boils.
Explicitly, as pressure being a function of T and V , from (2) we have Hence under constant pressure the temperature derivative acting on any function of β and V can be computed as Note that unless explicitly specified as ∂ ∂T P and ∂

∂T V
, throughout the paper ∂ β always denotes the beta derivative at fixed V and ∂ V is the volume derivative at fixed β, as already taken in (2).
In particular, the specific heat per particle under constant pressure reads We pay attention to the denominator here which is essentially the volume derivative of the pressure: This quantity possesses an indefinite sign. If ∂ V P < 0, the system is stable: it resists against the change of external pressure by adjusting its volume.
On the other hand, the opposite case, ∂ V P ≥ 0, characterizes the firstorder phase transition, as the volume at the phase transition is not single valued [6]. Clearly from (5), when ∂ V P = 0, a singularity develops and every physical quantity should change discontinuously. Moreover, when ∂ V P crosses the vanishing line, C P can be negative. To the best of our knowledge, all the known systems revealing negative specific heat do not keep the volume constant, such as in astrophysics [7,8], melting transitions [9,10], and in real experiments [11,12].
We emphasize that from (5) the only way for a finite canonical ensemble under constant pressure to reveal singularities is through ∂ V P = 0 i.e. the sign of the thermodynamic instability.
We close this section by presenting a sufficient, yet unnecessary, condition for the thermodynamic instability.

• Theorem
At low temperature near absolute zero, a thermodynamic system is un- The proof is simple once we spell the canonical partition function as Z N = n Ω n e −βEn , where E n 's are the possible energy eigenvalues and Ω n is the corresponding degeneracy. By direct manipulation we may express ∂ V P in the following form: where E 0 is the volume independent ground state energy, ∂ V E 0 = 0, with the degeneracy Ω 0 ; E 1 is the first excited state energy satisfying ∂ V E 1 = 0 with the degeneracy Ω 1 . The ellipsis denotes exponentially suppressed terms for large β or low temperature. Clearly at low temperature near absolute zero ∂ V P becomes positive. This completes our proof.
Examples include systems with vanishing ground state energy, such as supersymmetric models, ideal Bose or Boltzmann gases under periodic or Neumann boundary conditions. Appendix contains a numerical verification of the latter.

Ideal Bose gas confined in a cubic box
Based on the general analysis of the previous section, henceforth as a concrete model we focus on ideal Bose gas confined in a cubic box and subject to Dirichlet boundary condition. Since there is no zero mode, the ground state energy depends on the volume and our theorem above is not applicable in this case. Nevertheless, by numerical analysis we show that ∂ V P therein assumes positive values for some interval of temperature if N ≥ 7616.

Algebraic analysis
For non-interacting identical bosonic particle systems, more easily computed than the canonical partition function is the grand canonical partition function: where η denotes a fugacity and n corresponds to a good quantum number valued "vector" which uniquely specifies every quantum state of the single particle system. Taking logarithm and exponentiating back, we acquire an alternative useful expression: From the power series expansion of this, Z = N Z N η N , one can easily read off the canonical partition function, as previously obtained by Matsubara [37] and Feynman [38]: where the sum is over all the partitions of N , given by non-negative integers m a , a = 1, 2, · · · , N satisfying N = N a=1 a m a . In particular, with λ 1 = Z 1 the partition as m 1 = N leads to a conventional approximation [39], or the canonical partition function of ideal Boltzmann gas: This approximation would be only valid if all the particles occupied distinct states, as in the case of high temperature limit. Other partitions then give corrections to such underestimation: Compared to ideal Boltzmann gas, ideal Bose gas has higher probability for the particles to occupy the same quantum state. Yet, according to the Hardy-Ramanujan's estimation, the number of possible partitions grows exponentially like e π √ 2N/3 /(4 √ 3N ), and this would make any numerical computation practically hard for large N . Alternatively, we consider a recurrence relation on Z N , which was first derived by Landsberg [40] and can be easily reproduced here after differentiating Z in (10) by η : Further, if we formally define an ∞ × ∞ triangularized matrix ∆ whose entries are given by the above recurrence relation gets simplified: such that it has a solution given by a particular entry of a certain matrix: Since ∆ k [N +1, 1] vanishes for k ≥ N +1, the triangularized matrix ∆ can be effectively -and happily -truncated to its upper left (N +1)×(N +1) finite block. Then, from the standard recipe [41] to compute the inverse of a matrix, 2 another, novel expression of the canonical partition function follows: where Ω N is an almost triangularized N ×N matrix of which the entries are defined by 2 From det(I − ∆)=1, an intermediate relation between (16) and (17) follows: where ΦN is an N ×N matrix whose entry ΦN [a, b] is given by −λ a−b+1 /a for b ≤ a, unity for b = a+1 and zero otherwise.
In particular, every diagonal entry is unity so that when λ 1 → ∞, we have det(Ω N ) → 1 and hence the reduction: Z N → (Z 1 ) N /N ! as in (12).
In addition to the physical quantities (2) above, by considering −β −1 ∂ E n ln Z = ηe −βE n /(1 − ηe −βE n ) as a trick [42], we can also compute the number of particles occupying the ground state: Each term in the sum above corresponds to the probability for at least k particles to occupy the lowest state. If we denote this probability by p k , the difference p k − p k+1 corresponds to the probability for precisely k particles to occupy the ground state [43,44]. This leads to an alternative derivation of (19) as: Henceforth, exclusively for ideal Bose gas we focus on N particles with mass m, confined in a box of dimension d and length L ≡ V 1/d . Hard, impenetrable walls impose Dirichlet boundary condition [45]. Since we are interested in a finite system, the periodic boundary condition which is somewhat more popular in the literature is not suitable for our purpose. We recall that nevertheless enforcing periodic or Neumann boundary condition leads to a thermodynamic instability at low temperature near absolute zero for arbitrary N (see Figure 3 in Appendix).
In terms of a Jacobi theta function: we get specifically for (10): After all, Z N becomes a function of only one variable q, so that we may put This implies that all the dimensionless physical quantities such as C V /k B , C P /k B , N 0 , etc. are also functions of the single variable q. Consequently, it turns out that the temperature dependence of all these dimensionless quantities can be best analyzed if we introduce the following two dimensionless "temperatures": where Λ de Broglie = 2π/(mk B T ) is the thermal de Broglie wavelength and l eff. = (V /N ) 1/d is the average interparticle distance. While the former in (25) is a monotonically increasing function of q, the latter may be not so as: Any critical value of these quantities will automatically give us the critical temperature at an arbitrarily given volume or pressure.
In a similar fashion, we also define a dimensionless indicator of the thermodynamic instability, φ, as well as dimensionless "volume" and "energy": • Low temperature limit : The variable q lies between zero and one. As q → 0 we have Hence at q = 0, the temperatures vanish τ V = τ P = 0, and In particular, the volume reads at absolute zero, That is to say, despite of the apparent Bose-Einstein condensation i.e. N 0 = N , the volume assumes a finite value which is even not extensive. The finiteness is essentially due to the Heisenberg uncertainty principle: Since the particles are localized in a finite box, the uncertainty principle forbids the ground state energy E 0 to vanish and leads to the nontrivial canonical partition function (28).

Numerical results
Here we present our numerical results of d = 3, i.e. cubic box, at generic temperature. Our analysis is based on a set of recurrence relations which enables us to perform N 2 order computation, enhanced by a parallel computing power -JS20 (PowerPC 970) system. Essentially we utilize (13) 3 and two other relations which can be straightforwardly obtained from the grand canonical partition function expressed in the form (10): where we set ρ n := (β/n)∂ β ln λ n , ζ n := (β 2 /n 2 )∂ 2 β ln λ n .
• Constant volume curves ( Figure 1): As expected, all the physical quantities are smooth single valued functions of the temperature τ V . As N grows, the specific heat C V develops a maximum, at which N 0 drops rapidly. This behaviour is consistent with Ref. [24]. More importantly for us, while φ decreases from infinity at τ V = 0 to one at τ V = ∞, it develops a local minimum which becomes eventually negative if N ≥ 7616. This manifests the thermodynamic instability of the ideal Bose gas confined in a cubic box. 3 Since ZN can be a big number for large N , for the evaluation of the canonical partition function itself, it is convenient to decompose it as and utilize a recurrence relation: which is equivalent to (13).  • Constant pressure curves (Figure 2): Since the thermodynamic instability indicator φ vanishes at two points when N ≥ 7616, there are generically two critical temperatures, τ * P < τ * * P , which we suggest to identify as supercooling and superheating points respectively. As predicted from the general argument (5), all the physical quantities change discontinuously at these points: On each physical quantity versus temperature plane, the constant pressure line zigzags between the two critical points keeping them as two turning points. Accordingly, physical quantities are triple valued between the two points and double valued at the points. This implies the existence of three different phases during the phase transition, say two 'forward' phases and one 'backward' phase.
Specifically, the volume expands abruptly (Figure 2a): below the supercooling point it is almost constant whilst beyond the superheating point it follows the classical ideal gas law (34), which resembles the usual liquid-gas transition. At the same time, the number of particles on the ground state drops abruptly from the full occupancy to the total evacuation (Figure 2c). This clearly realizes a Bose-Einstein condensation taking place both in the momentum and in the coordinate space [47]. The specific heat per particle under constant pressure diverges at the critical points, and when magnified between the supercooling and the superheating points (Figure 2ee), intriguingly it reveals one negative heat capacity in addition to other two positive ones. Because C P /k B = dǫ P /dτ P + dυ P /dτ P , discontinuous changes in both the internal energy and the volume simultaneously contribute to the divergence of the specific heat. 4 Algebraically, while physical quantities are a priori single valued functions of q and hence τ V , the origin of the multi-valuedness and singular behaviours can be all traced back to the zigzagging of the isobar curve on (τ V , τ P ) plane ( Figure  2b). Namely the isobar curves of the ideal Bose gas zigzag on the (V, T ) plane, if N ≥ 7616.
Our numerical results of the supercooling and the superheating points are listed below for selected N : from which, recovering all the dimensionful parameters, we may read off the supercooling and the superheating temperatures at arbitrary constant pressure for each value of N :

Discussion
In summary, finite amount of ideal Bose gas confined in a cubic box reveals a thermodynamic instability, if N ≥ 7616. This implies that, under constant pressure condition, the ideal Bose gas may undergo a first-order phase transition accompanied by genuine mathematical singularities. It is characterized by both the supercooling and the superheating points; Bose-Einstein condensation in both the momentum and the coordinate spaces; and the very fact that all the physical quantities become triple valued between the two points. In particular, one of the three values of the specific heat under constant pressure is negative. While Bose-Einstein condensation appears as a continuous phase transition when the volume is kept constant (Figure 1 b), it may become a discrete phase transition if the pressure is held fixed (Figure 2 c).
Ideal gas consists of featureless point particles. It is such a simple model that the canonical partition function becomes essentially a one-variable function, Z N (q). Assuming an extra internal structure of the particle, e.g. spin or a vibrational mode, will introduce an additional dimensionful parameter into the energy spectrum. This may generate an additional (continuous or discontinuous) phase transition and distinguish the Bose-Einstein condensation to the ground state from the condensation in the coordinate space [49].
The critical number 7616 we report in this paper can be regarded as a characteristic number of 'cube' that is the geometric shape of the box containing the ideal Bose gas. Boxes of different shapes (see e.g. [50]) will have different critical numbers. Thus, our scheme of investigating the thermodynamic instability of ideal Bose gas can provide a novel algorithm to assign a characteristic number to each geometric closed two-dimensional manifold.
Apparently ideal Bose gas has no interaction. However, compared to ideal Boltzmann gas, ideal Bose gas has higher probability for the particles to occupy the same quantum state, as seen from the comparison between (11) and (12). 5 In coordinate space this means that identical bosonic particles have tendency to gather together in comparison to ideal Boltzmann gas. In fact, a path integral representation of the canonical partition function of ideal Bose gas reveals an existence of an attractive effective potential [6]. 6 This provides a physical clue to the condensation in both the momentum and the coordinate spaces: As the temperature decreases, the effective statistical attraction becomes dominant and the system condensates. In this context, it is also worth while to recall the similarity between the permutation symmetry of the identical particle indices and the gauge symmetry of the standard model in high energy physics or matrix models (see e.g. [52]). Although the former is discrete while the latter is continuous, the latter may include the former as a subgroup. The common feature is that they both correspond to nonphysical symmetry [53]. As a matter of fact, up to an overall factor, the canonical partition function of identical harmonic oscillators [27] coincides with that of a massive Yang-Mills quantum mechanics [54], taking the form N n=1 (1 − q n ) −1 . Generically, for a stable matter ∂ V P is negative. What we show by taking ideal 5 As a simple example, consider a two-particle system with quantum states 'up' and 'down' only. The probability for the two identical bosonic particles to occupy the same 'up' state is 1/3 and this is higher than that for the case of two Boltzmann particles, 1/4. Similarly, the probability for the two identical bosonic particles to occupy the two different states is 1/3 which is lower than that of the two Boltzmann particles, 1/2. 6 It is worth while to note that, the nontrivial volume at absolute zero due to the Heisenberg uncertainty principle (30) suggests a statistical repulsive interaction at low temperature. Nevertheless this is valid for ideal Boltzmann gas too. For related subtle issues see e.g. [51].
Bose gas as an exactly solvable model is an explicit demonstration that, if there are sufficiently, yet finitely, many identical bosonic particles, ∂ V P can be positive for at least one interval of temperature. It will be therefore interesting and crucial to see, to what extent interactions can alter this. If not much, as in a weakly interacting system, one first-order phase transition accompanying a discontinuous volume change, or the liquid-gas transition itself, is likely to occur essentially due to the identical nature of particles. In a perturbative analysis of the interaction, ideal Bose gas provides the 'zeroth' order contribution: For each quantum state ψ, as the energy eigenvalue is shifted from E ψ to E ψ + ∆ ψ , the canonical partition function changes from Z N to Z N e −β∆ ψ . Thus, the thermodynamic instability is modified as If the first order correction gives a negative contribution, ∂ 2 V ∆ ψ ≤ 0, the thermodynamic instability persists naturally up to the order. In this case, the average bonding energy between particles will also change discontinuously during the phase transition, but this will be an accompanying side effect rather than a key reason to drive the phase transition.
It will be experimentally challenging to find a corresponding critical number for each molecule to manifest a discontinuous phase transition or its liquid-gas transition under constant pressure. A criterion for the first-order phase transition is to observe the supercooling and the superheating phenomena.

Appendix: Other boundary conditions
Here as an explicit demonstration of the theorem (8), we present numerical results (Figure 3) on the ideal Boltzmann gas (or equivalently ideal Bose gas with N = 1) in a cubic box under boundary condition: Neumann or periodic. This corresponds to taking λ k (q) = [1 + ϑ(q k )] d for the former and λ k (q) = [1 + 2ϑ(q 4k )] d for the latter. Actually we only need λ 1 (q) and set d = 3 as for the cubic box.
As expected, due to the existing zero mode, both boundary conditions lead to a thermodynamic instability i.e. ∂ V P > 0 at low temperature near absolute zero. b: φ versus τ V under periodic boundary condition The description of identical particles appears closely related to a low energy strong coupling limit of Yang-Mills matrix models: The potential therein is generically given by matrix commutator squared, multiplied by a coupling constant. In a strong coupling limit, in order to maintain the energy finite, all the matrices should commute each other and become simultaneously diagonalizable, so that their eigenvalues are effectively only the remaining physical degrees. The unbroken gauge symmetry then corresponds to the permutation of the eigenvalues and can be identified as the permutation symmetry of the identical particle indices. For example,