Statistical Mechanics of Ideal Bose Systems

from Statistical Physics using Mathematica
© James J. Kelly, 1996-2002

Introduction

According to the Bose-Einstein distribution function, which is derived within the notebook occupy.nb, the mean occupancy for a single-particle orbital with energy ε is

\(n\&_[ε]\  = \ 
      1\/\(Exp[\(ε - μ\)\/\(\(k\_B\) T\)]\  - \ 1\)\)

where chemical potential μ is a function of density and temperature.  Positivity of the mean occupancy requires the chemical potential μ to be less than or equal to the energy of the lowest single-particle state; hence, since the ground-state is normally placed at (Subscript[ε, 0]) = 0, we require μ ≤  0.  However, if μ →  (Subscript[ε, 0]) for a finite temperature (Subscript[T, c]), the mean occupancy of the ground state appears to diverge.  Under these conditions the ground state can display macroscopic occupancy, in which an appreciable fraction of the particles collect in a single state, and the phenomenon of Bose-Einstein condensation occurs; special care is required to obtain consistent results in the regime T  ≤  (Subscript[T, c]).

The dependence of the chemical potential upon density and temperature is determined by integrating the occupation numbers to obtain the total particle number

\(𝒩\  = \ ∫\_0\%∞ n\&_[ε]\ 𝒟[ε] dε\)

where for large systems the density of states, 𝒟[ε], is approximated by a continuous function of energy.  For an ordinary nonrelativistic gas in three dimensions, the single-particle energy is simply \(TraditionalForm\`ε\  = \ \(\(ℏ\^2\) k\^2\)\/\(2\ m\)\) and the density of momentum states is 𝒟[(Overscript[k, ⇀])] (d^3) k =  
        g V  (((d^3) k)/(((2  π))^3)) where g = 2  s + 1 is the intrinsic degeneracy for a particle with spin s and V is the volume of the system.  Generalizations to one- or two-dimensional systems or to relativistic kinematics are developed in the problems at the end of the notebook.  The chemical potential enters the above integral through  \(TraditionalForm\`n\&_[ε]\) in a rather nontrivial fashion that requires numerical evaluation.  Similarly, other thermodynamic properties can be expressed in terms of integrals of this type.  Therefore, it is useful to define the bose function of order n as

\(bose[n, 
        z]\  = \ \(1\/Γ[
              n]\) \(∫\_0\%∞\( x\^\(n - 1\)\/\(\(z\^\(-1\)\) Exp[x] - 1\)\) dx\)\)

where z  =  Exp[(μ/((Subscript[k, B]) T))] is known as the fugacity.  The gamma function is included to ensure a convenient normalization, namely bose[n, z] -> z for small z.  The properties of bose functions are investigated below in considerable detail.

In this notebook we emphasize the properties of ideal nonrelativistic Bose gases in three dimensions with particular attention to the phenomenon of Bose condensation; generalizations to lower spatial dimensionality or to relativistic systems are explored in the problems at the end.  Thus, we show below that the grand potential for an ideal nonrelativistic Bose gas in three dimensions becomes

\(𝒢 = \(\(-p\)\ V\  = \ 
        g\ k\_B\ T\ \((\ 
            Log[1 - z]\ \  - \ \(\(\ \)\(V\ bose[5\/2, z]\)\)\/λ\^3\ )\)\)\)

where  \(TraditionalForm\`λ = \@\(\(2\ π\ ℏ\^2\)\/\(m\ k\_B\ T\)\)\) is the thermal wavelength.  The fugacity, or equivalently the chemical potential, is determined by the equation

(Subscript[n, Q]) =  ((N (λ^3))/(g V)) = 
        bose[(3/2), z]

where (Subscript[n, Q]) is the quantum concentration.  The phase transition for Bose-Einstein condensation occurs when the fugacity reaches its limiting value of unity for a critical concentration (Subscript[n, c]), such that

(Subscript[n, c])  =  ((N ((Subscript[λ, c])^3))/(g V)) = bose[(3/2), 
            1]  =  ζ[(3/2)] ⟹  (Subscript[T, c]) = 
            3.312502009395621` ((((N/(g  V))))^(2/
                    3)) (( (ℏ^2))/((Subscript[k, B]) m))

For temperatures below the critical temperature, a macroscopic occupation of the ground state develops such that

\(τ < 
        1\ \ \ \ ⟹\ \ \ \(N\_gs\/N\) = \ \(1 - τ\^\(3/2\)\ \ \ \ \ \ \ \ \ \ N\_exc\/N = \ τ\^\(3/2\)\)\)

where τ = T/(Subscript[T, c]) is the reduced temperature expressed in units of the critical temperature.  Since the ground state does not contribute to the entropy, energy, or pressure, most of the thermodynamic functions display discontinuities in either slope or value at the critical temperature.

Most textbook presentations of the thermodynamics of Bose gases employ simple box potentials because the mathematics is relatively simple, but experiments studying Bose-Einstein condensation do not use box potentials.  A much better approximation to experimental conditions of ultracold atomic traps is provided by harmonic confining potentials.  Although the analysis proceeds in a conceptually similar fashion, there are important modifications to the thermodynamic functions.  Nevertheless, the Bose functions developed here remain useful.  The thermodynamics of harmonic traps is developed in a series of problems at the end of this notebook.  The student is strongly encouraged to perform a thorough investigation of these phenomena, using the problems as guidance.

There are also several closely related distributions which describe the behavior of Bose systems which are not subject to a constraint upon the total particle number.  In the absence of such a constraint upon particle number, the chemical potential is identically zero.  For example, electromagnetic radiation in an enclosure consists of photons, which obey Bose-Einstein statistics but are not constrained to a fixed number.  The photon gas obeys the Planck distribution, which is a variation of the Bose-Einstein distribution for μ = 0, and is developed in the notebook planck.nb.  Similarly, phonons also obey Bose-Einstein statistics with μ = 0 and are studied in the notebook debye.nb.  The chemical potential, which governs the sharing of a conserved particle number between a system and a reservoir of particles, vanishes for photons or phonons because quanta of excitation are not conserved.  Instead, excitation quanta represent sharing of energy between a system and a heat bath and their number is already governed by temperature.

Glossary

𝒩	= number of particles
	V	= volume
	v	= volume/particle
	T	= temperature
	(Subscript[k, B])	= Boltzmann constant
	g	= energy-level degeneracy
	ε	= single-particle energy
	k	= single-particle wave number
	𝒟[ε]	= density of states wrt energy
	𝒟[k]	= density of states wrt momentum
	λ	= thermal wavelength
	(Subscript[n, Q])	= quantum concentration
	μ	= chemical potential
	z	= fugacity
	(Subscript[𝒵, 1])	= grand partition function for single orbital
	\(TraditionalForm\`n\&_\)	= mean occupation number
	(Subscript[T, c])	= critical temperature for Bose-Einstein condensation
	τ	= reduced temperature (T/(Subscript[T, c]))
	𝒢	= grand potential
	U	= internal energy
	p	= pressure
	S	= entropy
	F	= free energy
	H	= enthalpy
	(Subscript[C, V])	= isochoric heat capacity
	(Subscript[C, p])	= isochoric heat capacity
	γ  =  ((Subscript[C, p])/(Subscript[C, V]))	= ratio of principal heat capacities
	bose[n, z]	= Bose-Einstein function of order n
	α  =  
        (1/V) (Subscript[((∂(V/∂)T)), p]) 	= isobaric expansivity
	\(TraditionalForm\`κ\_T = \ \ \(-\(1\/V\)\)\ \((∂V\/∂p)\)\_T\)	= isothermal compressibility
	\(TraditionalForm\`κ\_S = \ \ \(-\(1\/V\)\)\ \((∂V\/∂p)\)\_S\)	= adiabatic compressibility

Initialization

Defaults, packages, and symbols

ClearAll["Global`*"];
Off[General::spell, General::spell1];

$DefaultFont = {"Times", 12};
$TextStyle = {FontFamily -> "Times", FontSize -> 12, 
          FontSlant -> "Italic"};

Needs["Utilities`Notation`"];
Needs["Miscellaneous`PhysicalConstants`"];
Needs["Miscellaneous`Units`"];
Needs["Graphics`Master`"]

Symbolize[(Subscript[𝒵, 1]^BE)]

Symbolize[(Subscript[k, B])];Symbolize[(Subscript[n, x])_];Symbolize[(Subscript[λ, x])_];Symbolize[(Subscript[T, x])_];Symbolize[(Subscript[C, x])_];Symbolize[(Subscript[v, x])_];Symbolize[(Subscript[V, x])_]

Symbolize[(Subscript[N, x])_];Symbolize[\(n\&_\)];Symbolize[\(\(n\_x_\)\&_\)]

Symbolize[(Subscript[k, max])];Symbolize[(Subscript[κ, x])_]

Symbolize[\(B\&~\)];Symbolize[\(p\&~\)];Symbolize[\(U\&~\)];Symbolize[\(S\&~\)];Symbolize[\(F\&~\)];Symbolize[\(G\&~\)];Symbolize[\(C\&~\_x_\)];Symbolize[\(γ\&~\)]

SetAttributes[{(Subscript[k, B]),  ℏ,  g,  m},  Constant];

FundamentalConstants  =  {(Subscript[k, B]) -> 
            (BoltzmannConstant/((Joule/Kelvin))), ℏ -> 
            (PlanckConstantReduced/(Joule Second))};

Memory management

This notebook consumes enough memory that it becomes advantageous to run the memory conservation utility.  However, we are not interested in seeing the messages reporting the actions of that process.

Needs["Utilities`MemoryConserve`"];
Off[MemoryConserve::start,  MemoryConserve::end];

Error messages

Several annoying error messages can be suppressed, if desired, but should be enabled when developing or debugging code.

Off[General::ovfl];
Off[Solve::ifun];
Off[NIntegrate::inum];
Off[NIntegrate::ncvb];
Off[NIntegrate::slwcon];
Off[NIntegrate::precw];

Rules for changing variables

TtoLambda = {T -> ((2 π (ℏ^2))/(m (Subscript[k, B]) (λ^2))),  β -> (1/((Subscript[k, B]) T))};
lambdaToT = {λ -> (Sqrt[((2 π (ℏ^2))/(m (Subscript[k, B]) T))])};
lambdaToZ = {λ -> (((((g V bose[(3/2), z])/𝒩)))^(1/3))};
boseToLambda  =  {bose[(3/2), 
              z] -> ((𝒩 (λ^3))/(g V))};

zToTau = z -> fugacity[τ];
muToZ = μ  -> (Subscript[k, B]) T Log[z];
zToT  =  z -> Exp[(μ/((Subscript[k, B]) T))];

Simplification rules

\(\(MyAssumptions = {m > 0, V > 0, g > 0, k\_B > 0, n_ > 0, n\&_ > 0, 
          T > 0, ℏ > 0, λ > 0, β > 0, 𝒩 > 
            0};\)\)

MySimplify = 
        Simplify[#, MyAssumptions] &;
MyFullSimplify = FullSimplify[#, MyAssumptions] &;

Density of states

Evaluate the density of states for nonrelativistic gas in 3 dimensions in terms of both momentum and energy.

𝒟[
          k] = ((4 π g V (k^2))/(((2 π))^3));

(Subscript[ε, k]) = ((((ℏ k))^2)/(2 m));

𝒟[ε] = 𝒟[k]/(∂_k) (Subscript[ε, k]) //. {k -> ((Sqrt[2 m ε])/ℏ)}  // PowerExpand

Bose-Einstein functions

In this section we explore some of the symbolic and numerical properties of Bose-Einstein functions and develop some simplification rules that will be useful in later sections.

Definitions

Bose-Einstein functions are also known as polylogarithm functions and are represented by Mathematica as PolyLog[ν,z]

bose[ν_, z_] = (1/Gamma)[ν] 
        Integrate[((x^(ν - 1))/(Exp[x]/z - 1)), {x, 0, ∞}, 
          Assumptions -> {Re[ν] > 0}]

The Bose-Einstein functions for {ν, (1/2), (3/2)} are plotted as functions of z.  In addition, the limiting case of ν -> ∞, where bose[ν, z] -> z, is shown.  Note that functions with ν ≤ 1 diverge as z -> 1, bose[(3/2), z]  remains finite at z = 1 but with infinite slope, while functions with ν > (3/2)  have finite values and slopes at z = 1.  Hence, for ν ≤ (3/2) special attention is needed to handle the vicinity of z -> 1.

\(\(Plot[
        Evaluate[{Table[bose[ν, z], {ν, 1\/2, 3, 1\/2}], z}], {z, 0, 
          0.999}, PlotRange -> {{0, 1}, {0, 2.5}}, 
        GridLines -> Automatic, Frame -> True, 
        FrameLabel -> {"z", \*""\!\(g\_ν\)[z]""}, 
        PlotLabel -> "Bose-Einstein functions"];\)\)

Small z expansion

Inspection of the Taylor series around z -> 0

Series[bose[ν, z], {z, 0, 3}]

suggests that the general power series takes the form

boseZexpansion[ν_, 
        z_, (Subscript[k, max]) : _] := (∑_(k = 1)^((Subscript[k, max])))((z^k)/(k^ν))

Below we verify that the infinite series does reproduce the desired function.

bose[ν, z] == 
      boseZexpansion[ν, z, ∞]

Thus, we recognize that Bose-Einstein functions are closely related to the Zeta function

(∑_(k = 1)^∞) ((z^k)/(k^ν))

bose[ν, 1]

It will be useful to define a rule for small z as follows.

zSmallRule = 
      bose[ν_, z] -> boseZexpansion[ν, z, 3]

Asymptotic expansion for z near 1

The behavior of Bose-Einstein functions near z -> 1 must be studied with more care because the integrand is singular at x = 0 for ν ≤ 1 and is undefined (0/0) for larger ν.  In the following we develop an expansion for z near unity, following the method of J.E. Robinson, Phys. Rev. 83, 678 (1951).  The method is based upon the Mellin transform of the expansion in powers of z developed in the previous section.  We define α by z = (e^(-α)) such that the Mellin transform takes the following form.

Clear[α]

mellinTransform[ν, 
          t] = (∑_(k = 1)^∞)(((∫_0^∞) (α^(t - 1)) Exp[-k α] \[DifferentialD]α)/(k^ν));

The inverse transform is obtained by integration in the complex t-plane along a line with positive Real[t] that is parallel to the imaginary axis.  The contour is closed with an infinite semicircular arc, to the left, where it can be shown that the integrand vanishes sufficiently rapidly.  The inverse transform is then

\(inverseMellinTransform = \(1\/\(2  π\ i\)\ \(∫\_\(c - i\ ∞\)\%\(c + i\ ∞\)α\^\(-t\)\ mellinTransform[ν, 
                t] dt\)\  = ∑\+poles α\^\(-t\)\ Residue[mellinTransform[ν, pole]]\)\)

Gamma[t] has simple poles at t = -m, where m is a nonnegative integer, with residues ((((-1))^m)/(m!)).  Zeta[t + n] has a simple pole with unit residue at t = 1 - n.  Therefore, when ν ≤ 0 or ν is nonintegral the inverse transform becomes:

inverseMellinTransform[ν, α, 
        (Subscript[k, max])] = (∑_(m = 0)^((Subscript[k, max])))(((α^m) Zeta[ν - m] (((-1))^m))/(m!)) + 
        Gamma[1 - ν] (α^(ν - 1))

Although the present formulation is indeterminate when ν is a positive integer, Robinson used analytic continuation to demonstrate that the result above is, in fact, valid for all ν because the singularities cancel.  Below we evaluate the first few terms numerically for several n.

boseAlphaExpansion[ν_, α_, (Subscript[k, max]) : _] = 
      inverseMellinTransform[ν, α, (Subscript[k, max])]

boseAlphaExpansion[ν, α, 3]

N[boseAlphaExpansion[(1/2), α, 3]]

N[boseAlphaExpansion[(3/2), α, 3]]

N[boseAlphaExpansion[(5/2), α, 3]]

N[boseAlphaExpansion[(3/2), α, 3]/bose[(3/2), z] //. {α -> 0.9, z -> Exp[-α]}]

It will be useful to define an expansion rule for z -> 1 as follows.

zLargeRule[mm_]  =  
      bose[ν_, z_]  -> 
        boseAlphaExpansion[ν, -Log[z], mm]

Recursion relations for symbolic derivatives

It is also useful to note the following downward recursion relation for derivatives.

(∂_z) bose[ν, z]

z (∂_z) bose[ν, z] == bose[ν - 1, z]

Rules for manipulation of Bose-Einstein integrals

Unfortunately, the Mathematica Integrate function still fails to recognize several simple variants of the Bose-Einstein integral.  Furthermore, we would like to employ symbolic integrals without needing to specify the parameters.  Therefore, in order to guide Mathematica toward proper evaluation of these integrals, we need to temporarily unprotect the Integrate routine so that we may install our own rules for the relevant types of integrals.  It is necessary to provide some rather redundant rules to ensure that all equivalent variations are recognized.

Unprotect[Integrate];

\(∫\_0\%∞\(\( A_\ x_\^n_\)\/\(Exp[β_. \ x_] - 
                1\)\) \[DifferentialD]x_ := 
      A\ \(Gamma[n + 1]\ Zeta[n + 1]\)\/β\^\(n + 1\)\)

\(∫\_0\%∞\(\( 
                A_\ Exp[\(-β_. \)\ x_]\ x_\^n_\)\/\(1\  - \ 
                Exp[\(-β_. \)\ x_]\)\) \[DifferentialD]x_ := 
      A\ \(Gamma[n + 1]\ Zeta[n + 1]\)\/β\^\(n + 1\)\)

\(∫\_0\%∞\(\( 
                A_\ x_\^n_\)\/\(Exp[β_. \ x_]\/z - 
                1\)\) \[DifferentialD]x_ := 
      A\ \(Gamma[n + 1]\ bose[n + 1, z]\)\/β\^\(n + 1\)\)

\(∫\_0\%∞\(\( A_\ x_\^n_\)\/\(Exp[β_. \ x_] - 
                z\)\) \[DifferentialD]x_ := 
      A\ \(Gamma[n + 1]\ bose[n + 1, z]\)\/\(z\ β\^\(n + 1\)\)\)

\(∫\_0\%∞\(\( 
                A_\ Exp[\(-β_. \)\ x_]\ x_\^n_\)\/\(1\  - \ 
                z\ Exp[\(-β_. \)\ x_]\)\) \[DifferentialD]x_ := 
      A\ \(Gamma[n + 1]\ bose[n + 1, z]\)\/\(z\ β\^\(n + 1\)\)\)

sumRule = (∫_0^∞)(a_ + 
                  b_) \[DifferentialD]x_ :> (∫_0^∞) a \[DifferentialD]x + (∫_0^∞) 
                b \[DifferentialD]x;

\(\(partsRule1 = ∫\_0\%∞\( x_\^m_\ y_\) \[DifferentialD]x_ :> ∫\_0\%∞\(\(-m\)\ x\^\(m - 1\)\ \(∫y \[DifferentialD]x\)\) \[DifferentialD]x;\)\)

\(\(partsRule2 = ∫\_0\%∞\( x_\^m_\ y_\) \[DifferentialD]x_ :> ∫\_0\%∞\(-\(\(x\^\(m + 1\)\ ∂\_x y\)\/\(m + 1\)\)\) \[DifferentialD]x;\)\)

Protect[Integrate];

gammaRule = n_ Gamma[n_] -> Gamma[n + 1];

Occupation numbers and chemical potential

Mean occupation number

The mean occupation number, as a function of energy and chemical potential,  is determined from the grand partition function.  In addition, the total number of particles in excited states is obtained by integrating the density of states against the occupation number distribution.  However, this method does not determine the ground state occupancy because the density of states is inaccurate for low energies, giving no weight to the ground state.

(Subscript[𝒵, 1]^BE) = (∑_(n = 0)^∞) 
          Exp[-n β (ε - μ)]

\(\(n\_BE\)\&_ = 
      k\_B\ T\ ∂\_\(\(\ \)\(μ\)\)\((Log[𝒵\_1\%BE] /. {β -> 
                    1\/\(k\_B\ T\)})\)\  // FullSimplify\)

\(\(n\_BE\)\&_ /. muToZ // FullSimplify\)

The fugacity, for the Bose-Einstein distribution must be limited to the range 0 ≤ z ≤ 1 so that the mean occupation number is nonnegative.

Evaluate the total number of particles in excited states

Although direct summation of the mean occupancy for all single-particle states should give the total number of particles in the system, for practical reasons we normally replace the sum by an integration over a continuous density of states.  However, the density of states formula for a nonrelativistic three-dimensional gas gives no weight to the ground state because the analytic density of states is proportional to (Sqrt[ε]) and approaches zero.  Therefore, integration of the mean occupancy against the density of states only gives the total number of particles, (Subscript[N, exc]), that occupy excited single-particle states.  The number of particles in the ground state, (Subscript[N, gs]), can then be obtained directly from the mean occupancy as \(TraditionalForm\`N\_gs = \ \(\(\(N\)\(\ \)\(-\)\(\ \)\(N\_exc\)\(\ \)\) = \ \ \ z\/\(1 - z\)\)\).  Therefore, as z → 1 the occupancy of the ground state becomes macroscopic.

\(N\_exc[
        z_]\  = \(Integrate[𝒟[ε]\ \(n\_BE\)\&_ /. 
                muToZ, {ε, 0, ∞}, 
              Assumptions -> {Re[1\/\(k\_B\ T\)] > 0}] /. TtoLambda // 
          PowerExpand\) // Simplify\)

This expression shows that the total number of particles found in excited states depends upon temperature, volume, chemical potential, and the particle type (mass, spin degeneracy).  However, the chemical potential is itself a function of temperature and density that is required to be nonpositive.  Hence, if the chemical potential reaches zero for a finite critical temperature, (Subscript[T, c]), the number of particles in excited states has an explicit value that depends only on temperature and volume (and particle type).  Therefore, any additional particles that are added to the system must enter the ground state because there is no room for them in excited states.  The ground state then develops a macroscopic occupancy that is described as Bose-Einstein condensation.

Critical temperature

The importance of quantum properties of the system is governed by the quantum concentration, (Subscript[n, Q]), defined by the product of the density and the volume corresponding to the thermal de Broglie wave length, λ.  It is also useful to include the intrinsic degeneracy of single-particle states, g.  Bose concentration occurs when the chemical potential reaches zero and it is no longer possible to accommodate the particle number N among the excited states.   At this point, the occupation of the ground state must become macroscopic and hence diverges for N -> ∞.

(Subscript[n, Q]) = ((𝒩 (λ^3))/(g V));

(Subscript[n, c]) = (((Subscript[n, Q]) (Subscript[N, exc])[1])/(𝒩 ))

(Subscript[n, Q]) /. lambdaToT

(Subscript[T, c]) = T/.Solve[(Subscript[n, Q]) == (Subscript[n, c]) /. lambdaToT, 
            T][[1]]

(Subscript[λ, c]) = λ /. Solve[
            (Subscript[n, Q]) == 
              (Subscript[n, c]), λ][[1]]

{(Subscript[n, c]), (Subscript[λ, c]), (Subscript[T, c])}  // N

As expected, the critical concentration for nonrelativistic Bose gases in three dimensions is of order unity.  The critical temperature is then determined by the mass and density.  However, it is also important to note that Bose condensation does not occur for relativistic systems or for smaller spatial dimensionality because with (Subscript[N, exc]) proportional to bose[ν, z] a smaller ν allows (Subscript[N, exc]) to diverge at z = 1 so that the excited states can accommodate essentially infinite particle number.

Ground state occupation

Now it is clear that below the critical temperature, where z is constant, the fraction of the total number of particles that are found in excited states must scale with (τ^(3/2)), where τ  =  T/(Subscript[T, c]), with the remainder in the ground state.

(Subscript[N, gs])[τ_ /; τ < 
            1] := 𝒩 (1 - (τ^(3/2)))
(Subscript[N, gs])[τ_ /; τ ≥ 1] := 0

\(\(Plot[{N\_gs[τ]\/𝒩, 
          1 - N\_gs[τ]\/𝒩}, {τ, 0, 1.5}, 
        PlotStyle -> {{}, Dashing[{0.02, 0.02}]}, 
        PlotLabel -> "Ground-State Population", Frame -> True, 
        FrameLabel -> {\*""T/(Subscript[T, c])"", \*""(Subscript[N, gs])/N, (Subscript[N, exc])/N""}, 
        PlotLegend -> {\*""(Subscript[N, gs])/N"", \ \*""(Subscript[N, exc])/N""}, LegendShadow -> None, 
        LegendPosition -> {0.3, \(-0.1\)}, 
        LegendSize -> {0.5, 0.3}];\)\)

We can now determine the actual value of z for T < (Subscript[T, c]) more precisely.

zsol  =  Solve[(Subscript[N, gs])[τ]  ==  g  (z/(1 - z)) , 
          z][[1]]

Series[ z /. zsol,  {(Subscript[N, gs])[τ], ∞, 1}]

In the immediate vicinity of the phase transition we can expand about τ = 1, but must be careful interpreting this result because 𝒩 is usually extremely large.

(Series[ z /. zsol,  {(Subscript[N, gs])[τ], ∞, 1}] /. 
              (Subscript[N, gs])[τ] -> 𝒩 (1 - (τ^(3/2))) ) // Normal // Series[#, {τ, 1, 2}] &  // 
      Collect[#, {𝒩}] &

Temperature dependence of chemical potential and occupancy

In this section we construct an interpolation function to represent the chemical potential as a function of reduced temperature.  This interpolation function can then be used to display the temperature dependence of thermodynamic functions.  Note that we approximate z by unity for T < (Subscript[T, c]), which is appropriate for large systems.

Here we express energy and chemical potential in units of (Subscript[k, B]) (Subscript[T, c]) and temperature as τ = T/(Subscript[T, c]).

eq1[
        z_, τ_] = 1 == (Subscript[N, exc])[z]/𝒩 /. 
            lambdaToT /. {T -> τ (Subscript[T, c])}  // 
        PowerExpand

zSolution[τ_] := 
      z /. FindRoot[Evaluate[eq1[z, τ]], {z, 0.9, 0.85}]

pointsNearTc = 
        Prepend[Evaluate[
            Table[{τ, τ Log[zSolution[τ]]}, {τ, 1.2, 
                1.5, 0.1}]], {1. , 0. }];

pointsAtLargerT = 
        Table[{τ, τ Log[zSolution[τ]]}, {τ, 2, 10, 
            0.5}];

points = Join[pointsNearTc, pointsAtLargerT];

Clear[chemicalPotential]
chem = 
        Interpolation[points, InterpolationOrder -> 5];
chemicalPotential[τ_ /; τ > 1] := chem[τ];
chemicalPotential[τ_ /; τ ≤ 1] := 
        0;

\(\(DisplayTogether[Plot[chemicalPotential[τ], {τ, 0, 5}], 
        ListPlot[points, PlotStyle -> AbsolutePointSize[4]], 
        PlotRange -> {{0, 5}, {\(-8\), 1}}, 
        PlotLabel -> "Chemical Potential", 
        GridLines -> Automatic, Frame -> True, 
        FrameLabel -> {\*""T/(Subscript[T, c])"", \*""μ/((Subscript[k, B])(Subscript[T, c]))""}];\)\)

fugacity[τ_ /; τ > 1] := 
        Exp[chem[τ]/τ];
fugacity[τ_ /; τ ≤ 1] := 1

\(\(Plot[fugacity[τ], {τ, 0, 5}, 
        PlotLabel -> "Fugacity", GridLines -> Automatic, 
        Frame -> True, 
        FrameLabel -> {\*""T/\!\(T\_c\)"", "z"}];\)\)

\(\(Show[
        GraphicsArray[{Plot[chemicalPotential[τ], {τ, 0, 5}, 
              PlotStyle -> {AbsoluteThickness[1.6]}, 
              PlotRange -> {{0, 5}, {\(-8\), 1}}, 
              PlotLabel -> "Chemical Potential", Frame -> True, 
              FrameLabel -> {\*""T/(Subscript[T, c])"", \*""μ/((Subscript[k, B])(Subscript[T, c]))""}, GridLines -> Automatic, 
              DisplayFunction -> Identity], 
            Plot[fugacity[τ], {τ, 0, 5}, 
              PlotStyle -> {AbsoluteThickness[1.6]}, 
              PlotLabel -> "Fugacity", Frame -> True, 
              FrameLabel -> {\*""T/(Subscript[T, c])"", "z"}, 
              GridLines -> Automatic, 
              DisplayFunction -> Identity]}, 
          GraphicsSpacing -> {\(-0.2\), 0}]];\)\)

Display mean occupation numbers for several temperatures

Clear[temp];
\(\(temp[τ_] := \(n\_BE\)\&_ //. {k\_B -> 1, 
            T -> τ, μ -> 
              chemicalPotential[τ]};\)\)
\(\(nBEplot[τ_] := 
        Plot[temp[τ], {ε, 0, 2}, 
          PlotRange -> {{0, 2}, {0, 5}}, 
          PlotLabel -> "Occupation Number", GridLines -> None, 
          Frame -> True, 
          FrameLabel -> {\*""ε/((Subscript[k, B])(Subscript[T, c]))"", \*""\!\(n\&_\)""}, 
          DisplayFunction -> Identity];\)\)

Show[{nBEplot[1.1], nBEplot[2], nBEplot[5]}, 
        DisplayFunction -> $DisplayFunction];

Rules for derivatives with respect to temperature

dzdT = -((3 z bose[(3/2), z])/(2 T bose[(1/2), 
                  z]));

dLambdaDT = 
        Dt[λ, T, 
            Constants -> {𝒩, 
                V}] -> -(λ/(2 T));
dzDThighT = 
        Dt[z, T, Constants -> {𝒩, V}] -> 
          dzdT;
dzDTlowT = 
        Dt[z, T, Constants -> {𝒩, V}] -> 
          0;

eq2 = ((𝒩 (λ^3))/(g V)) == 
        bose[(3/2), z]

dzHighTrules = 
      Flatten[{Solve[
              Dt[eq2, V, 
                Constants -> {𝒩, T, λ}], 
              Dt[z, V, Constants -> {𝒩, T, λ}]], 
            Solve[Dt[eq2, T, Constants -> {𝒩, V}], 
              Dt[z, T, Constants -> {𝒩, V}]]} /. 
          dLambdaDT]

Thermodynamic functions

Grand potential

The contribution of the ground state is treated separately because the ground state plays an important role at low temperatures, but the use of an approximate density of states would have given it zero weight.  However, because some of the thermodynamic functions become singular at the critical temperature, we insert a parameter, δ, to manipulate the ground-state contributions.

It is convenient to express the grand potential in terms of the so-called q potential defined as \(TraditionalForm\`q\  = \ \(\(p\ V\)\/\(\(k\_B\) T\) = \(ln\ 𝒵\ \ ⟹\ \ 𝒢 = \(-k\_B\) 
              T\ q\)\)\).

q = -(((∫_0^∞) 𝒟[ε] Log[1 - 
                              z Exp[-β ε]] \[DifferentialD]ε + δ g Log[1 - z] /. partsRule2) //. 
                TtoLambda  // PowerExpand) /. lambdaToZ

𝒢 = -(Subscript[k, B]) T q  // Simplify

Recognizing that 1-z  ≈  
        (1/(Subscript[N, gs]))[τ] below the critical temperature, it is clear that the term proportional to δ cannot compete with the term proportional to N for macroscopic systems (N  ∼  (10^24)) and may be omitted safely.  Thus, we can set δ -> 0 when evaluating the thermodynamic functions.  Furthermore, above the critical temperature we can express λ in terms of bose[ (3/2), z].  However, that replacement is not valid below the critical temperature because z reduces to unity whereas λ  ∝  (T^(-3/2)).

𝒢  /. lambdaToZ /. δ -> 0

Mechanical equation of state

Degeneracy pressure

The mechanical equation of state is obtained directly from the grand potential, 𝒢  =  -p V, omitting the ground-state contribution because without any kinetic energy it does not produce pressure.  It is useful to express the mechanical equation of state in the form ((p V)/(N (Subscript[k, B]) T)) = B[τ], where B is sometimes called the compressibility factor.

p = -(𝒢/V)  /. δ -> 0

B = ((p V)/(𝒩 (Subscript[k, B]) T))

Limit[p, z -> 0]

Note that in these expressions the dependence upon density is implicit in z.  Therefore, care must be exercised in using these functions to distinguish between the behaviors above and below the transition temperature for a given density.  Recognizing that only those particles occupying excited states contribute to the pressure and that (Subscript[N, exc])  =   N (τ^(3/2)) below the critical temperature, we find that the pressure of the system is proportional to (T^(5/2)) and is independent of density.

p /. boseToLambda /. lambdaToT /. z -> 1 // 
      PowerExpand

Therefore, the compressibility factor must scale with (τ^(3/2)) for temperatures below the critical temperature.  The function \(TraditionalForm\`B\&~\) is designed to handle temperatures both above and below the critical temperature in a form suitable for numerical evaluation.

Bhigh[τ_] = B /. lambdaToZ /. 
          zToTau;
Blow[τ_] = Bhigh[1] (τ^(3/2));
\(\(B\&~[τ_] = 
        If[τ > 1, Bhigh[τ], Blow[τ]];\)\)

\(\(plot[compressibility] = 
        Plot[B\&~[τ], {τ, 0, 5}, PlotPoints -> 50, 
          GridLines -> Automatic, Frame -> True, 
          FrameLabel -> {"τ", \*""pV/(𝒩(Subscript[k, B])T)""}, PlotLabel -> "Compressibility Factor"];\)\)

The compressibility factor for an ideal Bose gas approaches the classical limit from below at high temperatures.  The tendency for bosons to collect in occupied states reduces the pressure for given density and temperature compared to classical expectations.  This pressure reduction should be compared with the corresponding pressure increase in ideal Fermi gases that arises from the "quantum repulsion" due to the Pauli exclusion principle.   The pressure increase for Fermi gases is often called the degeneracy pressure.  Hence, the analogous degeneracy pressure for Bose gases is negative.  As z -> 1, the degeneracy of the Bose gas reduces the pressure for given density and temperature relative to classical expectations.

B /. z -> 1  // N

Transition line

For temperatures below the condensation temperature, the pressure becomes independent of both volume and particle number, becoming simply a function of temperature.  If any additional particles are added to the system, they enter the ground state without changing the pressure or volume of the system.  In that sense the volume per particle is associated entirely with excited states, with the volume contributed by the condensed phase vanishing.  Thus, the specific volume scales with (τ^(3/2)) in the transition region.   Since the critical temperature itself scales with ((( (N/V)))^(2/3)), the transition line between normal and condensed behavior has the form p ((( (V/N)))^(5/3))  =  constant.  Therefore, the transition line which delineates the onset of Bose condensation satisfies p (Subscript[V, c]^(5/3)) = constant.  Let (Subscript[v, c])[T]  be the specific volume as a function of temperature along the transition line.

(Subscript[v, c])[
        T] = (( (λ^3))/(g bose[(3/2), 1])) /. lambdaToT // 
        PowerExpand

pTransition[T] = p /. boseToLambda /. lambdaToT /. z -> 1 // 
        PowerExpand

PowerExpand[pTransition[T] (Subscript[v, c])[T]^(5/3)]

The latent heat of transition can now be obtained using the Clausius-Clapeyron equation, recognizing that the difference in the specific volumes of the vapor and condensed phases is just the specific volume at the transition temperature.

(∂_T) pTransition[T] T (Subscript[v, c])[T]  // N

Therefore, Bose condensation can be described as a first-order phase transition.

Dimensionless forms

The dependence is displayed in a dimensionless form.  The transition line, p (V^(5/3)) = constant, marks the onset of Bose condensation and defines the scale for a particular system.

pReduced[T_, 
        V_] = p /. boseToLambda /. zToTau /. {g -> 1, 
          (Subscript[k, B]) -> 1, λ -> (1/(Sqrt[T])), τ -> 
            T (V^(2/3))}

transitionLine[V_] = 
      pReduced[1, 1] (V^(-5/3 ))

Plot[
        Evaluate[{transitionLine[V], 
            Table[pReduced[T, V], {T, 0.5, 2, 0.5}]}], {V, 0.1, 5}, 
        PlotStyle -> {{RGBColor[1, 0, 0], 
              Dashing[{0.02, 0.02}]}, {}, {}, {}, {}, {}}, 
        GridLines -> Automatic, Frame -> True, 
        FrameLabel -> {"V", "p"}, 
        PlotLabel -> "Isotherms for ideal Bose systems"];

Compressibility and expansivity

\(κ\_T = \ \(\(\(\(\(-\(1\/\(V\ Dt[p, V, 
                        Constants -> {𝒩, 
                            T, λ}]\)\)\) /. dzHighTrules\) /. 
            dLambdaDT\) /. lambdaToZ // Simplify\)\(
\)
      \)\)

α = (Subscript[κ, T]) Dt[p, T, 
                    Constants -> {𝒩, V}]  /. 
                dzHighTrules /. dLambdaDT /. lambdaToZ  // Simplify  // 
        Apart

Virial expansion

When T ≫ (Subscript[T, c]),  z is small and an expansion about the classical limit shows

Series[B , {z, 0, 3}]

Alternatively, it is useful to develop an expansion for B in terms of the quantum concentration, x = ((N (λ^3))/(g V)).  First we expand x in terms of z and invert that result to determine the dependence of fugacity upon x for small density.  Then we develop the virial expansion itself.  Since it is important to maintain all terms of the same order, we truncate all expressions at (z^2).

zToX = Solve[x == Normal[Series[bose[(3/2), z], {z, 0, 2}]], 
          z][[1]]

virialExpansion = 
      Series[(Normal[Series[B, {z, 0, 2}]] /. zToX), {x, 0, 
          2}]

virialExpansion  // N // Chop

Thermal equation of state

The thermal equation of state for a nonrelativistic three-dimensional gas is obtained directly from the scaling relation U  =  (3/2) p V; other systems display similar scaling relations with different coefficients.

U = (3/2) p V

Recognizing the relationship between B and U, a form suitable for numerical evaluation over a wide range of temperatures is given below.

Uhigh[τ_] = (3/2) τ Bhigh[τ];
Ulow[τ_] = (3/2) τ Blow[τ];
\(\(U\&~[τ_] = 
        If[τ > 1, Uhigh[τ], Ulow[τ]];\)\)

\(\(plot[U] = 
        Plot[{\(3\ τ\)\/2, U\&~[τ]}, {τ, 0, 5}, 
          PlotStyle -> {Dashing[{0.02, 0.02}], {}}, 
          PlotPoints -> 50, GridLines -> Automatic, 
          Frame -> True, 
          FrameLabel -> {"τ", \*""U/((Subscript[Nk, B])(Subscript[T, c]))""}, 
          PlotLabel -> "Internal Energy"];\)\)

The internal energy for an ideal Bose gas is approaches the classical limit from below, and below the transition temperature is substantially lower than the classical prediction.   At the transition temperature the internal energy is only about half the classical prediction.

Isochoric Heat Capacity

It is easy to evaluate the isochoric heat capacity above the critical temperature

(Subscript[C, V]) = Dt[U, T, Constants -> {𝒩, V}] /. 
                dzDThighT /. dLambdaDT /. lambdaToZ  // 
          Collect[#1, {𝒩, (Subscript[k, B])}] &  // Simplify

CVhigh[τ_] = (Subscript[C, V]) /. zToTau

but below the critical temperature we must take care to recognize that only particles in excited states contribute to the heat capacity.  The simplest procedure uses the low-temperature form of the energy per particle as follows

CVlow[τ_] = 𝒩 (Subscript[k, B]) 
        D[Ulow[τ], τ]

CVlow[τ] // N

Hence, we combine these results into a single function and display the temperature dependence of the isochoric heat capacity below.

\(\(C\&~\_V[τ_] = 
        If[τ > 1, CVhigh[τ], CVlow[τ]];\)\)

\(\(plot[C\_V] = 
        Plot[C\&~\_V[τ]\/\(𝒩\ k\_B\), {τ, 0, 5}, 
          PlotRange -> {{0, 5}, {0, 2}}, GridLines -> Automatic, 
          Frame -> True, 
          FrameLabel -> {"τ", \*""(Subscript[C, V])/((Subscript[Nk, B]))""}, PlotLabel -> "Isochoric Heat Capacity"];\)\)

Isobaric Heat Capacity

The isobaric heat capacity is defined as the derivative of enthalpy with respect to temperature for constant pressure.  The enthalpy, U  +  p V, for an ideal nonrelativistic gas is quite simple, reducing to ((5 p V)/2), which allows the isobaric heat capacity to be obtained easily.  Alternatively, it can also be obtained using a more general approach which relates difference between (Subscript[C, p]) and (Subscript[C, V]) to the isobaric expansivity and the isothermal compressibility of the system.  The latter approach is more appropriate here.  However, since the pressure in the transition region depends only on temperature, the isobaric heat capacity becomes meaningless below the transition temperature.  Furthermore, both the isothermal compressibility and the isobaric expansivity diverge in the transition region.

(Subscript[C, p]) = 
      (Subscript[C, V]) + ((V T (α^2))/(Subscript[κ, T]))  // Simplify

\(\(C\&~\_p[τ_] = C\_p /. zToTau;\)\)

\(\(plot[C\_p] = 
        Plot[C\&~\_p[τ]\/\(𝒩\ k\_B\), {τ, 1, 5}, 
          PlotRange -> {{1, 5}, {0, 5}}, GridLines -> Automatic, 
          Frame -> True, 
          FrameLabel -> {\*""T/(Subscript[T, c])"", \*""(Subscript[C, p])/((Subscript[Nk, B]))""}, 
          PlotLabel -> "Isobaric Heat Capacity"];\)\)

Compare (Subscript[C, p]) and (Subscript[C, V])

For a classical ideal gas we would find (Subscript[C, V ])  =  (3/2) N (Subscript[k, B]), (Subscript[C, p])  =  (5/2) N (Subscript[k, B]), and ((Subscript[C, p])/(Subscript[C, V]))  =  (5/3)where p (V^γ) is constant for adiabats with γ = (5/3).  For the Bose gas we find a somewhat more complicated expression for γ in terms of Bose functions, but simple expansions can be developed for the degenerate or the dilute limits.

γ = ((Subscript[C, p])/(Subscript[C, V]))  // Simplify

\(\(plot[γ] = 
        Plot[Evaluate[γ /. zToTau], {τ, 1, 5}, 
          PlotRange -> {{1, 5}, {0, 3}}, GridLines -> Automatic, 
          Frame -> True, 
          FrameLabel -> {\*""T/(Subscript[T, c])"", \*""(Subscript[C, p])/(Subscript[C, V])""}, 
          PlotLabel -> "Ratio of Heat Capacities"];\)\)

γ /. zToTau /. τ -> 10

Entropy, free energy, and free enthalpy

Symbolic forms

S = ((U + p V - 𝒩 μ)/T) /. muToZ // 
          Simplify // Collect[#1, {𝒩, (Subscript[k, B])}] &

F = U  -  T S // Simplify // 
        Collect[#1, {𝒩, (Subscript[k, B])}] &

G = F  +  p V // Simplify // 
        Collect[#1, {𝒩, (Subscript[k, B])}] &

Plot reduced forms

We plot the entropy, free energy, and free enthalpy in dimensionless form.  Note that the behaviors above and below the transition temperature must be handled differently, but the proper scaling is clear from the results above.  When T < (Subscript[T, c]), we must handle the condensation phenomenon explicitly.  Recognizing that only those particles occupying excited states contribute to the entropy and that (Subscript[N, exc])  =   N (τ^(3/2)) below the critical temperature, we find that the entropy of the system is also proportional to N (τ^(3/2)).  Similarly, the Helmholtz free energy is proportional to N (τ^(5/2)) while the Gibbs free enthalpy vanishes (recall that 1 - z is of order (N^(-1))) in the thermodynamic limit.

Shigh[τ_] = 
        (S/(𝒩 (Subscript[k, B]))) /. zToTau // 
          Collect[#1, {𝒩, (Subscript[k, B])}] &;
Slow[τ_] = 
        (S/(𝒩 (Subscript[k, B]))) /. z -> 1;
\(\(S\&~[τ_] = 
        If[τ > 1, Shigh[τ], Slow[τ]];\)\)

Fhigh[τ_] = ((τ F)/(𝒩 (Subscript[k, B]) T)) /. 
            zToTau // Collect[#1, {𝒩, (Subscript[k, B])}] &;
Flow[τ_] = ((τ F)/(𝒩 (Subscript[k, B]) T)) /. 
          z -> 1;
\(\(F\&~[τ_] = 
        If[τ > 1, Fhigh[τ], Flow[τ]];\)\)

Ghigh[τ_] = ((τ G)/(𝒩 (Subscript[k, B]) T)) /. 
            zToTau // Collect[#1, {𝒩, (Subscript[k, B])}] &;
Glow[τ_] = ((τ G)/(𝒩 (Subscript[k, B]) T)) /. 
          z -> 1;
\(\(G\&~[τ_] = 
        If[τ > 1, Ghigh[τ], Glow[τ]];\)\)

\(\(plot[S] = 
        Plot[S\&~[τ], {τ, 0, 5}, 
          PlotLabel -> "Reduced Entropy", 
          GridLines -> Automatic, Frame -> True, 
          FrameLabel -> {"τ", \*""S/(\!\(Nk\_B\))""}, 
          DisplayFunction -> Identity];\)\)

\(\(plot[F] = 
        Plot[F\&~[τ], {τ, 0, 5}, 
          PlotLabel -> "Reduced Free Energy", 
          GridLines -> Automatic, Frame -> True, 
          FrameLabel -> {"τ", \*""F/((Subscript[Nk, B])(Subscript[T, c]))""}, 
          DisplayFunction -> Identity];\)\)

\(\(plot[G] = 
        Plot[G\&~[τ], {τ, 0, 5}, 
          PlotLabel -> "Reduced Free Enthalpy", 
          GridLines -> Automatic, Frame -> True, 
          FrameLabel -> {"τ", \*""G/((Subscript[Nk, B])(Subscript[T, c]))""}, 
          DisplayFunction -> Identity];\)\)

Show[ 
        GraphicsArray[{{plot[U], plot[S]},  {plot[(Subscript[C, V])], plot[F]}}, 
          GraphicsSpacing -> {-0.1, 0.05}], 
        DisplayFunction -> $DisplayFunction ] ;

Latent heat of condensation

Since the entropy of the condensed phase vanishes, the entropy in the transition region is due entirely to the particles in excited states.  Furthermore, since the fraction of the particles in excited states scales with (τ^(3/2)), the reduced entropy for the gas (uncondensed) phase is constant in the transition region.  Hence, the latent heat of transition is simply T times the entropy per particle at the critical temperature.

latentHeat = ((T S)/𝒩) /. lambdaToZ /. 
          z -> 1  // N

Therefore, as expected, we obtain the same result from the entropy as we did previously from the Clausius-Clapeyron equation.

Examples

Liquid helium

Although liquid helium is not an ideal fluid, the superfluid transition occurs at a temperature close to the Bose condensation temperature and the superfluid state exhibits many properties which are similar to the condensed phase of an ideal Bose fluid.  The following numerical calculations employ SI units, with tags removed, unless stated otherwise.  The volume is converted to molar volume.  The experimental value of the liquid helium density is used to estimate the theoretical critical temperature at atmospheric pressure.

values = 
        Join[{g -> 1, 𝒩 -> AvogadroConstant Mole, 
            density -> (146/m), 
            m -> ((4 Convert[AtomicMassUnit, Kilogram])/Kilogram)},  
          FundamentalConstants];
liqHeDensity = V -> (𝒩/density);

{(Subscript[V, He]), (Subscript[T, He])} =  {V, (Subscript[T, c])} /. liqHeDensity //. values

Thus, the Bose-Einstein critical temperature is intriguingly close to the superfluid transition temperature, 2.17 kelvins, for liquid (^4)He.  However, the thermodynamic properties of liquid helium do not conform to those of an ideal Bose fluid.  For example, the isochoric heat capacity displays a marked discontinuity at the superfluid transition whereas (Subscript[C, V])[T] for an ideal Bose fluid would be continuous (despite its kink).  Nevertheless, the similarities suggest that Bose-Einstein condensation plays an important role in the superfluid transition and that deviations from ideality can be attributed to interatomic interactions.

pHe[T_, 
          V_] = p /. boseToLambda /. lambdaToT /. 
                  zToTau /. τ -> (T/(Subscript[T, He])) ((((V/(Subscript[V, He]))))^(2/3)) /. values // PowerExpand // 
          Simplify;

HeTransitionLine[V_] = 
        pHe[(Subscript[T, He]), V] (((((Subscript[V, He])/V)))^(5/3));

\(\(Plot[
        Evaluate[{HeTransitionLine[V]\/10\^6, 
              Table[pHe[T, V]\/10\^6, {T, 2, 6, 2}]} /. {V -> 
                x\/10\^4}], {x, 0.01, 1.0}, 
        PlotRange -> {{0, 1}, {0, 3}}, 
        PlotStyle -> {{RGBColor[1, 0, 0]}, {}, {}, {}, {}, {}}, 
        GridLines -> Automatic, Frame -> True, 
        FrameLabel -> {\*""V  [(10^(-4)) (m^3)/Mole]"", \*""p  [(10^6) Pa]""}, 
        PlotLabel -> "pV diagram for ideal helium"];\)\)

Although this model correctly describes many of the qualitative features of the state diagram for liquid helium, it is not very accurate.  For example, the pressure computed below for the critical temperature at the experimental density for atmospheric pressure is actually much larger than one atmosphere.  This increased pressure is partly due to the difference between the theoretical and experimental transition temperatures.  Furthermore, part of the discrepancy also reflects the repulsive core of the interatomic interaction; the importance of this effect can be gauged by comparing the molar volume to the van der Waals volume parameter.

pHe[(Subscript[T, He]), (Subscript[V, He])] Convert[Pascal, Atmosphere]

pHe[2.17, (Subscript[V, He])] Convert[Pascal, Atmosphere]

vdWaals = {a -> 0.003457, 
          b -> 23.7×(10^(-6))};

((Subscript[V, He])/b) /. vdWaals

Problems

(Subscript[k, B]) (Subscript[T, c]^0) ≪ (Subscript[ε, 1]) ⟹  \(T\_c\/T\_c\%0\)≈1  -  (2/(3  ζ[(3/2)])) 
          Exp[ - ((Subscript[ε, 1])/((Subscript[k, B]) 
                    (Subscript[T, c]^0)))]

(Subscript[k, B]) (Subscript[T, c]^0) ≫ (Subscript[ε, 1]) ⟹  \(T\_c\/T\_c\%0\)≈ \(\((1\/2)\)\^\(2/3\)\ \((1\  + \ \(2\^\(4/3\)\/\(3  ζ[
                      3\/2]\)\) \((\(π\ ε\_1\)\/\(\(k\_B\) T\_c\%0\))\)\^\(1/2\))\)\)

v  =  (Sqrt[(Subscript[((∂(p/∂)ρ)), S])])

\(v\^2 = \ \(\(\(5 \( k\_B\) T\)\/\(3  m\)\) 
          bose[5\/2, z]\/bose[3\/2, z]\  = \ \(5\/9\) \(u\^2\)\&_\)\)