Research at PIDE: Key Messages

arXiv:cond-mat/0505316v1 [cond-mat.stat-mech] 12 May 2005 Ideal Fermi gases in harmonic oscillator potential traps David J. Toms School of Mathematics and Statistics, University of Newcastle Upon Tyne, Newcastle Upon Tyne, United Kingdom NE1 7RU Abstract We study the thermodynamic properties of an ideal gas of fermions in a harmonic oscillator confining potential. The analogy between this problem and the de Haasvan Alphen effect is discussed and used to obtain analytical results for the chemical potential and specific heat in the case of both isotropic and anisotropic potentials. Step-like behaviour in the chemical potential, first noted in numerical studies, is obtained analytically and shown to result in an oscillatory behaviour of the specific heat when the particle number is varied. The origin of these oscillations is that part of the thermodynamic potential responsible for the de Haas-van Alphen-type effect. At low temperatures we show analytically that there are significant deviations in the specific heat from the expected linear temperature dependence, again as a consequence of the de Haas-van Alphen part of the thermodynamic potential. Results are given for one, two, and three spatial dimensions. In the anisotropic case we show how the specific heat jumps as the ratio of oscillator frequencies varies. Key words: Trapped Fermi gas, de Haas-van Alphen, specific heat, chemical potential, oscillations PACS: 03.75.Ss, 05.30.Fk, 71.10.Ca Preprint submitted to Elsevier Science 31 October 2018 1 Introduction The experimental realization [1] of Bose-Einstein condensation in confined gases of atoms at low temperatures has been the stimulus for a wide range of experimental and theoretical investigations. (See Ref. [2] for a review.) By ingenious experimental techniques it is possible to prepare the atoms as either bosons or fermions. There have now been a number of experiments on trapped gases of fermions, as well as mixtures of bosons and fermions [3]. For a gas of fermions, it is a good approximation to model the system as a collection of non-interacting particles obeying Fermi-Dirac statistics confined by a simple harmonic oscillator potential. Unlike the case of bosons, the dominant s-wave scattering channel is suppressed making the effects of interactions less important [4]. With this simple theoretical model for the system, it is easy to obtain exact and simple expressions for the single-particle energy levels using the standard quantum mechanical result for the simple harmonic oscillator. This can serve as a starting point for a theoretical analysis of the thermodynamic behaviour of the system. Because the energy levels for a set of particles (bosons or fermions) confined by a simple harmonic oscillator potential form a discrete set, in order to calculate the thermodynamic properties we must be able to perform sums over the quantum numbers labelling the states; this can lead to difficulties in evaluation. By restricting the parameters of the system (temperature and oscillator frequencies for example) it may be possible to argue that the sums can be approximated with integrals, an approximation that renders the computation Email address: [email protected] (David J. Toms). URL: http://www.staff.ncl.ac.uk/d.j.toms/ (David J. Toms). 2 of analytical results easier. We will refer to the approximation of sums with integrals as the continuum approximation, since it is equivalent to regarding the discrete energy spectrum as continuous. This approximation was first applied to the Fermi gas in the now classic paper of Butts and Rokhsar [4]. Even when one goes beyond this approximation, as we will show below, the continuum approximation can still be very accurate in some temperature regimes, and can capture the leading order behaviour if the temperature is not too low. The first study that showed there might be more features present beyond what is obtained by the continuum approximation was given by Schneider and Wallis [5]. These authors computed the thermodynamic expressions numerically, and their results showed that there were a number of step-like features present, not predicted by the continuum approximation. In particular, the chemical potential, which is a smooth function in the continuum approximation, has a series of step-like jumps when a numerical analysis of the exact result is performed. This in turn can lead to step-like, or oscillatory, behaviour of other thermodynamic properties, such as the specific heat. Another important consequence of the work presented in Ref. [5] is that at sufficiently low temperatures there are significant deviations in the specific heat from the linear temperature dependence predicted [4] by the continuum approximation. The main purpose of the present paper is to analyze the situation analytically to obtain the step-like behaviour of the chemical potential, and to apply this method to the evaluation of the specific heat. In addition to the two studies already mentioned [4,5], there have been a number of other related theoretical works. Ref. [6] generalized the calculations of Ref. [5] to a gas with two spin states with an interaction between the two states. (We will make some brief comments on the effects of interactions in 3 Sec. 6.) A number of authors [7,8,9,10,11,12,13,14] have obtained analytical approximations in a number of cases using path integral, Green’s function, or density matrix methods. The approach described in the present paper was used to study partially confined gases in Ref. [15] In the present paper we will concentrate on the behaviour of the chemical potential and specific heat, and will use the evaluation of the thermodynamic potential as the basic starting point. The key observation is that for a simple harmonic oscillator confining potential, the evaluation of the thermodynamic potential is mathematically very similar to that of a gas of charged fermions in a constant magnetic field. Landau [16] was the first to show that the thermodynamic properties of electrons in a constant magnetic field undergo oscillations whose period is determined by the inverse of the magnetic field strength. This gave a theoretical basis for the observations of de Haas and van Alphen, now referred to as the de Haas-van Alphen effect. (See Ref. [17] for a good discussion of the de Haas-van Alphen effect and the important role it plays in condensed matter physics.) What we will show here is that for the trapped Fermi gas the step-like behaviour found in Ref. [5] is completely analogous to the de Haas-van Alphen effect. Just as in the de Haas-van Alphen effect, where use of the continuum approximation would miss the oscillations, here too we must go beyond the continuum approximation. The outline of our paper is as follows. In Sec. 2 we describe the general method following the classic analysis of Sondheimer and Wilson [18]. There are some differences with the de Haas-van Alphen case considered in Ref. [18] that we will describe. We will show how the continuum approximation comes about in the method, and even find the next order correction to it, beyond the leading order result given in Ref. [4]. In Sec. 3 we analyze the de Haas-van Alphen 4 contribution to the thermodynamic potential for a 3-dimensional isotropic harmonic oscillator potential, and obtain approximate analytical results for the chemical potential. This is then used to examine the specific heat. Results for the 1- and 2-dimensional gases are given in Sec. 4 and the results compared with the 3-dimensional case. In Sec. 5 we analyze the situation when the trapping potential is anisotropic. Finally, Sec. 6 contains a brief summary and a short discussion of the main results. 2 General method We will begin with the grand canonical ensemble and the thermodynamic potential Ω for a system of fermions with energy levels En defined by Ω = −T X f (En ) , (2.1) n where f (E) = ln [1 + exp(−β(E − µ))] . (2.2) We use the usual notation β = 1/T , with T the temperature, and choose units such that the Boltzmann constant is equal to one. Take the Laplace transform of f (E), and call it ϕ(β), so that ϕ(β) = Z∞ dE e−βE f (E) . (2.3) 0 The inverse Laplace transform of Eq. (2.3) reads c+i∞ 1 Z f (E) = dβ eβE ϕ(β) , 2πi (2.4) c−i∞ 5 where c is an arbitrary constant chosen so that the integration path lies to the right of any singularities of ϕ(β). We can now use Eq. (2.4) in Eq. (2.1). The sum over the energy levels that results can be related to the partition function for the µ = 0 system, defined by Z(β) = X e−βEn . (2.5) n We can then write Ω defined by Eq. (2.1) as T Ω=− 2πi c+i∞ Z dβ Z(−β)ϕ(β) . (2.6) c−i∞ It is necessary here to regard Z(β) as a function of β with β viewed as a complex variable. Z(−β) is the result obtained from this complex function by analytic continuation of the definition Eq. (2.5) to ℜ(β) < 0. We now define Z(E) to be the Laplace transform of β −2 Z(β). (The factor of β −2 is for later convenience [18].) This results in the definition Z(β) = β 2 Z∞ dE e−βE Z(E) . (2.7) 0 The inverse Laplace transform gives 1 Z(E) = 2πi c+i∞ Z dβ eβE β −2 Z(β) . (2.8) c−i∞ Substitution of Eq. (2.7) into Eq. (2.6) followed by the use of Eq. (2.4) results in Ω = −T Z∞ 0 dE Z(E) ∂ 2 f (E) . ∂E 2 (2.9) The Fermi-Dirac distribution function F (E) is 6 h F (E) = eβ(E−µ) + 1 = −T i−1 (2.10) ∂ f (E) , ∂E (2.11) if we use Eq. (2.2). This leads to Ω= Z∞ dE Z(E) 0 ∂ F (E) , ∂E (2.12) giving the key starting point in the Sondheimer-Wilson [18] analysis of the de Haas-van Alphen effect. The method rests on an evaluation of Z(E) defined by Eq. (2.7) or Eq. (2.8) and its use in Eq. (2.12). To evaluate Z(E) in Eq. (2.8) we require the partition function, or at least some information about the singularities of β −2Z(β) in the complex β-plane. For the case of a simple harmonic oscillator potential this information is very easy to obtain. We will consider the D-dimensional oscillator potential D 1 X ωj2x2j , V (x) = m 2 j=1 (2.13) with x = (x1 , . . . , xD ) the D spatial coordinates. The energy levels are simply given by (in ~ = 1 units) En = D  X nj + j=1 1 ωj , 2  (2.14) where n = (n1 , . . . , nD ) and nj = 0, 1, 2, . . . for j = 1, . . . , D. We then have Z(β) in Eq. (2.5) expressed as a product of geometric series that are easily summed to give Z(β) = D Y e−βωj /2 . −βωj ) j=1 (1 − e (2.15) 7 Because of the neglect of interactions, and the simple form of the potential, it has been possible to evaluate the partition function in closed form. To obtain Z(E) we need to know the singularity structure of β −2 Z(β). It is now obvious from Eq. (2.15) that β −2 Z(β) is a meromorphic function. There is a pole at β = 0 of order D + 2 as well as a series of poles along the positive and negative imaginary β-axis at β = βkj = 2πikj /ωj with kj = ±1, ±2, . . .. The order of these poles depends on the relative ratios of the harmonic oscillator frequencies. If none of the oscillator frequencies are rational multiples of any of the others, all of the poles, apart from the one at β = 0, will be simple. However, if some of the oscillator frequencies are rational multiples of the others, some of the poles away from the origin will be higher order. This makes a determination of the residues of poles along the positive and negative imaginary β-axis for a general harmonic oscillator potential an inelegant treatment of special cases. We will return to this problem later. We can write Z(E) in Eq. (2.8) as Z(E) = Z0 (E) + Zr (E) , (2.16) where we close the contour in the left hand side of the complex plane and use Z0 (E) to denote the contribution to the integral from the pole at β = 0, and Zr (E) the contribution coming from the rest of the poles along the positive and negative imaginary β-axis. When Eq. (2.16) is used in Eq. (2.12) we obtain the thermodynamic potential as Ω = Ω0 + Ωr (2.17) 8 in an obvious way. As in the de Haas-van Alphen effect, the oscillations in thermodynamic quantities will come from Ωr , but we will first examine Ω0 . Ωr will be evaluated in Sec. 3 for the isotropic potential, and in Sec. 5.1 for the anisotopic potential. To obtain Z0 (E) we need the residue of β −2 eβE Z(β) at β = 0. For general values of D and arbitrary frequencies ωj this is difficult to write down in any simple way. We will concentrate on just the cases D = 1, 2, 3 here, although there is no intrinsic reason why the method cannot be extended to other values. A straightforward calculation yields E2 ω − for D = 1, (2.18) 2ω 24  E  2 4E − ω12 − ω22 for D = 2, (2.19) Z0 (E) = 24ω1ω2 h 1 Z0 (E) = 240E 4 − 120E 2 (ω12 + ω22 + ω32 ) + 7(ω14 + ω24 + ω34 ) 5760ω1ω2 ω3 i +10(ω12ω22 + ω22ω32 + ω32 ω12 ) for D = 3. (2.20) Z0 (E) = These results can now be used in Eq. (2.12) to obtain Ω0 = Z∞ dE Z0 (E) 0 ∂ F (E) . ∂E (2.21) As T → 0, meaning that βµ → ∞, the Fermi-Dirac distribution function approaches a step-function, F (E) → θ(µ − E), so that ∂ F (E) ∂E → −δ(E − µ). This crude approximation results in Ω0 → −Z0 (µ) , (2.22) and is really the first term in a systematic expansion due to Sommerfeld. (See Refs. [17,19,20] for three derivations of the Sommerfeld expansion.) In our case 9 we find Ω0 ≃ −Z0 (µ) − π 2 2 ′′ 7π 4 4 ′′′′ T Z0 (µ) − T Z0 (µ) . 6 360 (2.23) Although the general Sommerfeld expansion contains higher order terms, these vanish here because Z0 (E) is no more than quartic in E due to our restriction that D ≤ 3. For arbitrary D, the expansion Eq. (2.23) will contain more terms with increasing powers of T . For low values of T these should be less important than those indicated in Eq. (2.23) in any case. Using Eqs. (2.182.20) in Eq. (2.23) results in 1 2 π2 2 ω (µ + T ) + , for D = 1 , 2ω 3 24 µ µ3 + (ω 2 + ω22 − 4π 2 T 2 ) , for D = 2 , Ω0 ≃ − 6ω1 ω2 24ω1 ω2 1 µ4 µ2 Ω0 ≃ − + (ω 2 + ω22 + ω32 − 4π 2 T 2 ) 24ω1 ω2 ω3 48ω1 ω2 ω3 1 h 1 7(ω14 + ω24 + ω34) + 10(ω12 ω22 + ω22 ω32 + ω32 ω12 ) − 5760ω1ω2 ω3 i −40π 2 T 2 (ω12 + ω22 + ω32 ) + 112π 4 T 4 , for D = 3 . Ω0 ≃ − (2.24) (2.25) (2.26) The thermodynamic properties of the system follow from a knowledge of the thermodynamic potential. In particular, the average particle number N is given by N =− ∂Ω ∂µ ! , (2.27) T,ω and the internal energy U is given by U= ∂ . (βΩ) ∂β βµ,ω (2.28) 10 In the most physically interesting case we take D = 3. (The cases of D = 1, 2 will be given later in Sec. 4.) If we temporarily ignore the contribution Ωr to Ω, we have N ≃ N0 where N0 = − ∂Ω0 ∂µ ! (2.29) T,ω µ3 µ ≃ − (ω12 + ω22 + ω32 − 4π 2 T 2 ) . 6ω1 ω2 ω3 24ω1 ω2 ω3 (2.30) For large values of µ, if we keep only the term in Eq. (2.30) of order µ3 , we find (with N ≃ N0 ), µ3 ≃ 6ω1 ω2 ω3 N , (2.31) showing that µ ∝ N 1/3 . For large values of N we will have µ much larger than the average oscillator frequency, consistent with the assumptions made in deriving these results. It is easy to obtain the next order correction to Eq. (2.31), (ω 2 + ω22 + ω32 − 4π 2 T 2 ) µ3 ≃ 6N 1 + 1 ω1 ω2 ω3 4(6Nω1 ω2 ω3 )2/3 ( ) , (2.32) showing that the correction to Eq. (2.31) becomes increasingly unimportant for large values of N. If we use Eq. (2.32) and rewrite the result in terms of the Fermi energy defined by EF = µ(T = 0) , (2.33) it is easy to see that µ ≃ EF π2T 2 1− 3EF2 ! . (2.34) 11 This agrees with the result in Ref. [4] who used the continuum approximation and only the T -dependent part of the second term in Eq. (2.30). We have demonstrated here that Eq. (2.34) holds true even when the next order approximation for the density of states is included in the continuum approximation. However, as we will see in the next section there are important corrections to this from terms that come from beyond the use of the continuum approximation and support the numerical investigations of Ref. [5]. We can also calculate the contribution to the internal energy, that we call U0 , coming from Ω0 . Using Eq. (2.26) in Eq. (2.28) and then eliminating µ in favour of N with Eq. (2.32) results in the following approximation : 3 (ω 2 + ω22 + ω32 + 4π 2 T 2 ) 2/3 U0 ≃ (6ω1 ω2 ω3 )1/3 N 4/3 + 1 N , 4 8(6ω1ω2 ω3 )1/3 (2.35) assuming that N is large. (The intermediate result for U0 in terms of µ can be found in Eq. (3.25) below for the isotropic potential.) The specific heat can be found using C= ∂U ∂T ! . (2.36) N,ω The contribution coming from the approximate result in Eq. (2.35) is then seen to be C ≃ π 2 (6ω1 ω2 ω3 )−1/3 N 2/3 T . (2.37) As before this agrees with the result quoted in Ref. [4] demonstrating a linear dependence on the temperature. Again we will find that at low temperatures the part of Ω that has not yet been considered, namely Ωr , alters this ex12 pected behaviour, confirming the numerical results of Ref. [5]. We consider the evaluation of Ωr and its effect in the next section. 3 Isotropic 3-dimensional potential 3.1 Thermodynamic potential We will examine the case of the isotropic harmonic oscillator potential with ω1 = ω2 = ω3 = ω , (3.1) in which case the results of the previous section simplify. Our aim here is to calculate Ωr that has been neglected up to now from the analysis. To do this we need to evaluate the contribution to Ω coming from the poles of β −2 eβE Z(β) along the imaginary β-axis away from β = 0. From Eq. (2.15), because Z(β) = e−3βω/2 , (1 − e−βω )3 (3.2) there will be poles of order 3 at β = βk defined by βk = 2πi k ω (3.3) for k = ±1, ±2, . . .. It is a straightforward matter to evaluate the residues of β −2 eβE Z(β) at β = βk and show that their contribution to Eq. (2.8) gives (−1)k ω 4π 2 k 2 E 2 2 2 Zr (E) = 6+π k − cos(2πkE/ω) 4 4 ω2 k=1 16π k ! ∞ X + (−1)k E sin(2πkE/ω) . 3 3 k=1 2π k ∞ X 13 (3.4) It now remains to use this in Eq. (2.12) and to try to extract something useful from the result. We have ∞ ∞ βω X (−1)k (−1)k β X A − Bk , k 64π 4 k=1 k 4 8π 3 k=1 k 3 (3.5) Ak = Z∞ 4π 2 k 2 E 2 dE 6 + π k − ω2 (3.6) Bk = Z∞ dE E Ωr = − where 0 ! 2 2 0 cos(2πkE/ω) , cosh2 [ 12 β(E − µ)] sin(2πkE/ω) . cosh2 [ 12 β(E − µ)] (3.7) In these last two integrals we can make the change of variable E =µ+ 2 θ, β (3.8) defining a new integration variable θ. The lower limits on the integrals in Eqs. (3.6) and (3.7) then become −βµ/2. If we look at low enough temperatures, specifically T << µ as we have already assumed, then to a good approximation (up to exponentially small terms) we can replace the lower limits on the integrals defining Ak and Bk with −∞. The approximate results then become 2 Ak ≃ β Z∞ −∞ 4π 2 k 2 µ2 16π 2 k 2 µθ 16π 2 k 2 θ2 dθ 2 2 6 + π k − − − ω2 βω 2 β 2ω2 cosh2 θ " 2πkµ 4πkθ + × cos ω βω 2 Bk ≃ β Z∞ −∞ ! ! , 2θ 2πkµ 4πkθ dθ µ+ sin + 2 β ω βω cosh θ 14 # (3.9) ! . (3.10) The integrals in Eqs. (3.9) and (3.10) may be evaluated exactly using residues. All of the expressions required may be related to the basic integral Z∞ −∞ πa cos b cos(aθ + b) dθ = , 2 sinh( π2 a) cosh θ (3.11) and derivatives with respect to the parameters a and b. After some straightforward calculation, using the results for Ak and Bk found from Eqs. (3.9) and (3.10) as described, it can be shown that ∞  2  nh X 1 (−1)k 2 2π k 4 2 Ωr ≃ − 2 3 2 8π k csch 2 8π β ω k=1 k 3 sinh ( 2πβωk ) βω +4π 2 kβω coth  2π 2 k  βω + 2β 2 ω 2 i +π 2 k 2 (4π 2 + β 2 ω 2 − 4β 2 µ2 ) cos(2πkµ/ω) h +4πkβµ βω + 2π 2 k coth  2π 2 k i βω o sin(2πkµ/ω) , (3.12) The presence of the trigonometric functions in Eq. (3.12) is responsible for the oscillations that occur in thermodynamic quantities if we go beyond the leading order continuum approximation. The presence of the hyperbolic functions in Eq. (3.12) with argument 2π 2 k βω means that unless βω > 1 these oscillations will be suppressed. The oscillatory behaviour should therefore show up for T < ω and become more prominent as T is reduced. The inclusion of Ωr in the expression used for the thermodynamic potential will lead to corrections to the thermodynamic behaviour of the system beyond what is found using the continuum approximation of Sec. 2. We will look first at how the chemical potential is affected. 15 3.2 Chemical potential We can define (from Eq. (2.27) ) Nr = − ∂Ωr ∂µ ! , (3.13) T,ω and then use Eq. (3.12) to obtain ∞ 2 i (−1)k nh 2 2 π X 2 2π k 2 2 2 2 ( 4β µ − 4π − β ω − 8π csch ) Nr ≃ 3 3 2 4β ω k=1 sinh ( 2πβωk ) βω × sin(2πkµ/ω) + 8πβµ coth ( o 2π 2 k ) cos(2πkµ/ω) . βω (3.14) This will make an oscillatory correction to the contribution N0 for the average number of particles found using the continuum approximation in Sec. 2. The chemical potential may be found by solving N = N0 + Nr (3.15) for µ, where N0 is given by (see Eq. (2.30) with ω1 = ω2 = ω3 = ω) µ3 4π 2 µ N0 ≃ 3 − − 6ω 3 24ω β 2ω2 ! . (3.16) Obviously the complicated dependence on µ in Nr as given in Eq. (3.14) renders an analytical evaluation of µ difficult. We can simplify by using the assumption that µ is large ( since we have already assumed that βµ ≫ 1 and µ ≫ ω), and keep only the leading term in µ. This gives Nr ≃ ∞ (−1)k sin(2πkµ/ω) πµ2 X . 2 βω 3 k=1 sinh ( 2πβωk ) 16 (3.17) The continuum approximation in Eq. (3.16) still gives the leading contribution to N for large µ; however, Nr in Eq. (3.17) can be more important than the sub-leading term (proportional to µ) in Eq. (3.16) if the temperature is low enough. For βω of the order of 1 or less, the sum in Eq. (3.17) is well approximated by simply the first term and it is clear that the oscillations, although present, will have a very small amplitude. We would therefore expect that for the temperature range T ≥ ω, Eq. (2.31) or Eq. (2.32) would provide a good approximation for the chemical potential. However as the temperature is reduced, so that βω ≫ 1, the amplitude of the oscillations coming from Eq. (3.17) become more pronounced and must be taken into effect. It is possible to find an asymptotic expansion for Nr in Eq. (3.17) valid for βω ≫ 1 (or T ≪ ω). After some calculation, it can be shown that µ2 βω Nr ≃ tanh (2µ̄ − 1) + 1 − 2µ̄ 2 4ω 4 ( ! ) , (3.18) where µ̄ = µ µ − ω ω   , (3.19) with [x] denoting the largest integer whose value is less than or equal to x. (Thus 0 ≤ µ̄ < 1 can be assumed in Eq. (3.18).) Similar asymptotic expansions can be found for the sub-leading terms in Nr given in Eq. (3.14); however we will not give them here. It can be shown by numerically evaluating the sum in Eq. (3.17) and comparison with the analytical approximation in Eq. (3.18) that Eq. (3.18) does give an accurate result for large values of βω. 17 If we take the limit βω → ∞ in Eq. (3.18), the result simplifies further to give Nr ≃ µ2 2ω 2  µ µ 1 − + ω 2 ω   . (3.20) Strictly speaking, this further approximation is only valid if h µ ω + 1 2 i is not equal to an integer. In this special case, Nr in Eq. (3.17) or Eq. (3.18) can be seen to vanish, and we must look at the sub-leading contributions to Nr that follow from Eq. (3.14). If we use Eq. (3.20) for Nr in Eq. (3.15) along with Eq. (3.16) for N0 , it is possible to show that the solution for µ is given by the simple expression µ 1 ≃ [(6N)1/3 ] + . ω 2 (3.21) The presence of the greatest integer function in Eq. (3.21) leads to the steplike behaviour first found in the numerical studies of [5]. Our results provide a confirmation of this behaviour by analytical means. For large values of N, our result in Eq. (3.21) shows that these steps will occur roughly for N ≃ ℓ3 /6 where ℓ is an integer. This agrees with the “magic numbers” found by [5] that occurred for N = ℓ(ℓ + 1)(ℓ + 2)/6 if we make N, and hence ℓ, large enough. We have therefore seen how the step-like behaviour of the chemical potential comes about in an analytical way, and traced its origin back to the same type of terms that are responsible for the de Haas-van Alphen effect. In Fig. 1 we show the result for µ/ω plotted as a function of N over a range of N. We have taken N large, but the same type of behaviour can be found for smaller values as well. The continuum approximation for µ is shown as the smooth solid, almost straight, line, and the simple analytical approximation in Eq. (3.21) is superimposed on it. The steps occur at the magic numbers as 18 39 βω = 10 38 µ ω βω = 20 37 βω = 100 36 35 8000 8200 8400 8600 8800 9000 9200 9400 9600 9800 10000 N Fig. 1. (color online) This plot shows µ/ω plotted over a range of N for three sample temperatures. The solid, almost straight, line shows the result found using the continuum approximation for the particle number resulting in µ/ω = (6N )1/3 . The step-function superimposed shows the result found from the very low temperature approximation of Eq. (3.21). The smooth sinuous curves give the results for T = 0.1ω, T = 0.05ω and T = 0.01ω. The last two curves have been displaced from the first one for clarity. predicted. We have displaced the curves for different values of βω away from each other for clarity to show the trend towards the step-like behaviour as the temperature is reduced. As the temperature is increased, the amplitudes of the oscillations decreases as mentioned above. 19 3.3 Specific heat We will first calculate that part of the energy that arises from Ωr . Define Ur = ∂ , (βΩr ) ∂β βµ,ω (3.22) and use Eq. (3.12) for Ωr . After some calculation it can be shown that ∞ πµ3 X 2πkµ (−1)k Ur ≃ sin 2k 2π 3 βω k=1 sinh ( βω ) ω ! 2 ∞ (−1)k cosh ( 2π k ) 2πkµ 3π 2 µ2 X βω cos + 2 3 2 2π 2 k β ω k=1 sinh ( βω ) ω ! (3.23) where we only include the two leading order terms in µ since it is these two terms we will need to calculate the leading contribution to the specific heat. The specific heat was defined in Eq. (2.36). It is straightforward to show that this expression is equivalent to C= ∂U ∂T ! µ,ω −  ∂U ∂µ   T,ω  ∂N ∂µ ∂N ∂T  µ,ω , (3.24) T,ω which is more useful for explicitly calculating C. The presence of the second term on the right hand side complicates the evaluation of C, but we will obtain an expansion for C in powers of µ. If we look at the first term on the right hand side of Eq. (3.24) and use U = U0 + Ur with Ur given by Eq. (3.23) and U0 = ∂ (βΩ0 ) ∂β βµ,ω µ2 µ4 (ω 2 − 4π 2 T 2 ) = 3− 3 8ω 16ω 20 + 1 (17ω 4 + 40π 2 ω 2 T 2 − 112π 4T 4 ) 3 1920ω (3.25) we obtain ∂U ∂T ! µ,ω π 2 µ2 T ≃ + 2ω 3 ∂Ur ∂T ! . (3.26) µ,ω From Eq. (3.23), counting powers of µ, it can be observed that the second term on the right hand side of Eq. (3.26) will contain expressions in µ3 , µ2, . . .; thus, it appears as if the leading order behaviour of the specific heat will be µ3 , rather than µ2 as predicted by the continuum approximation. However, we will show that the µ3 part of Eq. (3.26) cancels with a similar expression coming from the second term in Eq. (3.24) leaving the overall leading behaviour of the specific heat as µ2 rather than µ3 . In any case, we must work consistently to order µ2 , and to this order only the two contributions to Ur that we have written down in Eq. (3.23) are necessary. We will drop all terms that result in explicit factors of µ in C that are of order µ and lower. Using Eq. (3.25) we have ∂U ∂µ ! T,ω µ3 ≃ + 2ω 3 ∂Ur ∂µ ! , (3.27) T,ω and it is noted that the second term on the right hand side contains terms involving µ3 , µ2 , . . .. From Eq. (3.15) and Eq. (3.16) we have ∂N ∂T ! µ,ω π 2 µT ≃ + 3ω 3 ∂Nr ∂T ! , (3.28) µ,ω with the second term on the right hand side involving µ2 , µ, . . ., as well as ∂N ∂µ ! T,ω µ2 ≃ + 2ω 3 ∂Nr ∂µ ! , T,ω 21 (3.29) with the second term on the right hand side involving µ2 , µ, . . .. A simple counting of powers of µ, using the expressions for Ur and Nr , shows that the second term on the right hand side of Eq. (3.24) does begin at order µ3 . It now remains to use the explicit results for Ur given in Eq. (3.23) and Nr given in Eq. (3.14) and evaluate C for large µ keeping terms of order µ3 and µ2 . After some calculation it can be shown that the order µ3 terms cancel leaving C≃ π 2 T µ2 n Σ22 o 1 + Σ − 12 , 1 6ω 3 Σ3 (3.30) with π 2 kT 2 1 cosh θk − + Σ1 = 12 (−1) 2 ω sinh θk sinh3 θk sinh θk k=1   µ , × cos 2πk ω ( )   ∞ X 1 µ 2π 2 kT cosh θk k Σ2 = (−1) sin 2πk − , sinh θk ω sinh2 θk ω k=1 ∞ X k !) ( µ (−1)k 4π 2 k cos 2πk Σ3 = 1 + βω sinh θk ω k=1 ∞ X   , (3.31) (3.32) (3.33) where θk = 2π 2 k βω (3.34) has been defined to save a bit of writing. The continuum approximation of Sec. 2 can be regained by dropping all of the terms in Σ1,2,3 that have arisen from the de Haas-van Alphen part of the thermodynamic potential. This gives the familiar linear dependence on temperature [4]. The numerical calculations of [5] showed that as the temperature got sufficiently small there was a significant departure from the linear temperature 22 dependence in the specific heat. Once again we will show that this follows from our analytical method, and that the origin of this behaviour is in the de Haasvan Alphen part of the thermodynamic potential. To do this we will evaluate the asymptotic expansion of the three sums defined in Eqs. (3.31–3.33) for large values of βω. Leaving out the technical details of this for brevity, we find the approximate forms Σ1 ≃ 3β 3ω 3 (2µ̄ − 1)2 16π 2 cosh2 Σ2 ≃ − Σ3 ≃  βω (2µ̄ 4 − 1) β 2 ω 2 (2µ̄ − 1) 16π cosh2  βω (2µ̄ 4 βω 2 4 cosh  βω (2µ̄ 4  − 1) − 1)  , −1,  , (3.35) (3.36) (3.37) with µ̄ defined as in Eq. (3.19). The results in Eqs. (3.35–3.37) can be checked against a numerical evaluation of the sums defined by Eqs. (3.31–3.32) and found to be accurate for βω ≃ 10 and µ̄ not too close to 0 or 1. Once βω ≃ 100 the results become very accurate even for µ̄ close to 0 and 1. Thus for T ≤ ω/100, the simple expressions in Eqs. (3.35–3.37) become reliable approximations for Σ1,2,3 . If we use Eqs. (3.35–3.37) in the expression for the specific heat in Eq. (3.30) the result can be shown to vanish. The de Haas-van Alphen contribution to the specific heat cancels the continuum approximation to the leading order we are working to. The specific heat therefore vanishes as T → 0 faster than T . This is completely consistent with the numerical results found in [5]. Because the de Haas-van Alphen approximation, as well as the asymptotic evaluation of the sums leading to Eqs. (3.35–3.37) neglect terms that are exponentially suppressed, we suspect that the specific heat vanishes like e−αω/T for some constant α as T → 0, but we have not been able to establish this in any 23 simple way. Further support for this belief follows from the result we are able to establish for the 1-dimensional gas in Sec. 4.1 below. A more refined estimate of the sums, as well as the de Haas-van Alphen contribution would reveal the exact nature of the T → 0 limit. 1.4 1.2 (6N )−2 / 3C 1 T = 0.5ω T = 0 .3 ω 0.8 0.6 T = 0 .2 ω 0.4 0.2 0 T = 0.1ω 100 200 300 400 500 NN Fig. 2. (color online) This plot shows C/(6N )2/3 plotted over a range of N for four sample temperatures T = 0.1ω, T = 0.2ω, T = 0.3ω and T = 0.5ω. As the temperature increases the oscillation amplitude decreases and the curves approach the continuum limit of 1 in the scaled specific heat. As the temperature decreases there are significant deviations from the result found from using the continuum limit, and for very small temperatures the specific heat starts to become vanishingly small. For T ≥ ω/100, but still small, the results in Eqs. (3.35–3.37) start to become less reliable. As a check on our results against the numerical ones of [5] we plot the specific heat as found from Eq. (3.30) to demonstrate the de Haasvan Alphen oscillations. This is shown in Fig. 2. As the temperature increases we do find, as expected, that the contribution from Σ1,2,3 becomes smaller exponentially, and the linear behaviour with temperature is regained. 24 4 One and two dimensions We will examine the thermodynamics of trapped Fermi gases in one and two spatial dimensions, since these cases may be of relevance as limiting cases of the 3-dimensional gas in some situations. Because the methods used are similar to those described above in the 3-dimensional case we will be brief here and only exhibit key results that show a difference with results already obtained. 4.1 One dimension With D = 1, Z(β) in Eq. (2.15) has simple poles at β = βk with βk still defined by Eq. (3.3). It is easy to show that the contribution from the poles with k 6= 0 to Z(E) is   ∞ ω X E (−1)k Zr (E) = − 2 . cos 2πk 2π k=1 k 2 ω (4.1) The contribution from the pole at β = 0 was given in Eq. (2.18). After some calculation, making the approximation µ ≫ T , it can be shown that Ωr ≃ ∞ 1X (−1)k cos(2πkµ/ω) . 2 β k=1 k sinh ( 2πβωk ) (4.2) The calculation is very similar to the D = 3 case, so details will not be given here. The continuum approximation for Ω, called Ω0 , was given in Eq. (2.24). 25 Using Eq. (2.24) and Eq. (4.2) in the general expression for the average particle number N in Eq. (2.27) results in N≃ ∞ µ 2π X (−1)k sin(2πkµ/ω) . + 2 ω βω k=1 sinh ( 2πβωk ) (4.3) The first term on the right hand side comes from Ω0 , and the second term from Ωr , so that the continuum approximation would yield simply µ ≃ ωN (4.4) in this case. Again for large particle numbers we expect µ ≫ ω to be a valid approximation. For T ≪ ω we can again find the asymptotic form of the sum in Eq. (4.3) and obtain ( 1 1 βω µ + + tanh N≃ ω 2 2 4    2µ µ −1 −2 ω ω   ) , (4.5) and this proves to be very accurate for low values of T . The internal energy can be found to be U≃ µ2 π2 2 ω + T + + Ur , 2ω 6ω 24 (4.6) with ( ∞ 2π X π cosh θk µ sin(2πkµ/ω) sin(2πkµ/ω) + (−1)k Ur ≃ β k=1 ω sinh θk βω sinh2 θk ) (4.7) and θk defined by Eq. (3.34). We can calculate the specific heat using Eq. (3.24) with the result C≃ Σ2 o π2T n 1 + Σ1 − 12 2 , 3ω Σ3 (4.8) 26 The sums Σ1,2,3 are the same as those given in Eqs. (3.31–3.34). The result in Eq. (4.8) is very similar to that found in Eq. (3.30) in the 3-dimensional case apart from the overall factor that can be recognized as the continuum approximation for the specific heat of the 1-dimensional gas. We can therefore conclude immediately that as T → 0, the specific heat vanishes faster than the linear temperature dependence deduced from the continuum approximation. 1.6 1.4 1.2 1 C π2T 3ω 0.8 0.6 0.4 0.2 0 0.1 0.2 0.3 0.4 0.5 T /b ω Fig. 3. (color online) This plot shows the specific heat plotted over a range of temperature for the 1-dimensional gas. As the temperature increases the specific heat approach the continuum limit that behaves linearly on the temperature. As the temperature decreases the result deviates substantially from the linear result and the curve decays exponentially fast as found in Eq. (4.9). Because of the simplicity of the D = 1 case, we can shed some light on the behaviour of the specific heat as T → 0. We note that the continuum approximation for the specific heat in Eq. (4.4) is also a solution when the de Haasvan Alphen contributions to N are included, since the sin term vanishes for integral N. The sum Σ2 defined by Eq. (3.32) vanishes as well if µ = ωN, and 27 we can use the approximation for Σ1 given by Eq. (3.35), valid for T ≪ ω to obtain from Eq. (4.8) the simple result C≃ ω2 ω exp − 2 16T 2T   . (4.9) This demonstrates that the specific heat does vanish exponentially fast as T → 0, not linearly. The role of the de Haas-van Alphen part of the thermodynamic potential is again responsible for this behaviour. As a check on this conclusion we have plotted the specific heat as a function of temperature in Fig. 3. The results can be seen to be consistent with our analytical result in Eq. (4.9) as T → 0. For larger values of the temperature the specific heat approaches the linear temperature dependence predicted by the continuum approximation. This can also be seen from the expression obtained in Eq. (4.8) since as βω becomes small, the sums Σ1,2,3 start to vanish as a consequence of the hyperbolic functions present in the expressions. 4.2 Two dimensions Using D = 2 in Eq. (2.15) we find Zr (E) = − ∞ o ω 1n 1 X E cos(2πkE/ω) − sin(2πkE/ω) . 2π 2 k=1 k 2 πk (4.10) This results in ( ∞ 1 X 1 µ cos(2πkE/ω) Ωr ≃ 2πk 2 2πβ k=1 k ω sinh θk " # cosh θk 1 + θk sin(2πkE/ω) − sinh θk sinh2 θk 28 ) , (4.11) with θk defined by Eq. (3.34) if we make the same approximations as in the D = 3 case. The (−1)k factor in the summand, present for D = 1, 3, is absent here, a result that is true for any even dimension. The average particle number can be shown to be µ2 1 π2 − + + Nr , 2ω 2 12 6β 2 ω 2 N≃ (4.12) with Nr ≃ ∞ n o µ sin(2πkµ/ω) 2π X π cosh θk cos(2πkµ/ω) + βω k=1 ω sinh θk βω sinh2 θk (4.13) the contribution coming from the de Haas-van Alphen part of Ω. In the continuum approximation (dropping the term in Nr ) we find µ ≃ (2N)1/2 , ω (4.14) to leading order for large N. However the de Haas-van Alphen part of the particle number can be shown to lead to the step-like behaviour we saw previously in the 3-dimensional case. In the low temperature limit (T ≪ ω) an asymptotic analysis of Eq. (4.13) can be used to show that 1 µ N≃ 2 ω  2 ( µ βω 3 µ +2 tanh − 2 ω ω 2      µ µ − ω ω  ) . (4.15) Solving this expression for µ gives rise to the approximate analytical solution 1 µ . ≃ (2N)1/2 + ω 2   (4.16) (Again we remind the reader that the square brackets in Eqs. (4.15) and (4.16) denote the greatest integer function.) We therefore conclude that the 29 step-like behaviour found earlier for the chemical potential when D = 3 is also found for the 2-dimensional gas. We have checked this result by solving Eq. (4.12) numerically and found that the approximation in Eq. (4.16) becomes increasingly accurate as T is reduced. Because the results resemble the similar behaviour exhibited in Fig. 1 we will not show them here. The jumps in the chemical potential occur, for large N, when N ≃ 21 ℓ2 for integral ℓ. (The exact result is 12 ℓ(ℓ + 1) analogously to the 3-dimensional case studied in Ref. [5].) It is straightforward to obtain expressions for the internal energy and specific heat. Once again the results are similar to those found in the 3-dimensional case, this time shown in Fig. 2, so need not be exhibited in detail. As T → 0 the result for the specific heat vanishes much faster than the linear approximation C ≃ π 2 T (2N)1/2 /(3ω) predicted by the continuum approximation. 5 Anisotropic 3-dimensional potential 5.1 Thermodynamic potential As in the isotropic case, we will apply the analysis used by Sondheimer and Wilson [18] for the de Haas-van Alphen effect to obtain the thermodynamic potential. This allows a clear separation of the terms responsible for the continuum approximation and those that lead to oscillatory behaviour in thermodynamic expressions. The calculation of the thermodynamic potential will be used in the next two sections to obtain the chemical potential, particle number, and the specific heat. Because we have already outlined the details of the Sondheimer-Wilson [18] 30 method above, we will begin with the main result that expresses the thermodynamic potential as (see (2.12)) ∞ 1 Z Z(E) Ω = − β dE 2 1 4 cosh 2 β(E − µ) (5.1) 0 where Z(E) is related to the inverse Laplace transform of the partition function Z(β) for the gas with µ = 0. (See (2.5) and (2.8).) To evaluate Z(E) in Eq. (2.8) we require the partition function. For the case of a free Fermi gas confined by a simple harmonic oscillator potential we have the result in Eq. (2.15) with D = 3 taken there. The heart of the calculation consists of using Eq. (2.15) in Eq. (2.8), closing the contour in the left-hand side of the complex plane, and evaluating the result by residues as in the isotropic case considered in Sec. 3. As explained above, and as should be evident from Eq. (2.15), the poles of the integrand of Eq. (2.8) depend on the ratios of the three oscillator frequencies. There will be a pole of order 5 at β = 0 as well as poles along the imaginary axis. It is easy to see from Eq. (2.15) that there are poles away from the origin at β = βj where βj = 2πikj ωj (5.2) with j = 1, 2, 3 and kj = 0, ±1, ±2, . . .. The ratio of oscillator frequencies determines the order of these poles, and the order affects the residue. If none of the oscillator frequencies are rational multiples of each other, then the poles for β1,2,3 will all be distinct and therefore simple. However if any of the trap frequencies are rational multiples of any of the others, then some of the poles along the imaginary axis will no longer be simple. Therefore a 31 complete analysis of the problem requires investigating a variety of different cases depending on the relative ratios of the three oscillator frequencies. In order to avoid too many different cases, and to stop the paper being overly lengthy, we will only look at the case of an axisymmetric trap characterized by ω1 = ω2 = ω , ω3 = λω , (5.3) where λ is some real number that characterizes the anisotropy. There are then only two distinct cases to consider: λ rational, and λ irrational. The isotropic results can be recovered in the special case λ = 1. Regardless of whether λ is rational or irrational, there is a pole of order 5 at β = 0 for the integrand of Eq. (2.8). If we call Ω0 the contribution to the thermodynamic potential that arises from this pole, then i µ4 µ2 h 2 2 2 2 + (2 + λ )ω − 4π T 24λω 3 48λω 3 i h 1 4 4 2 2 2 2 4 4 , 3(8 + 3λ )ω − 80π T ω (2 + λ ) + 112π T − 5760λω 3 Ω0 ≃ − (5.4) as follows from the result in (2.26) with (5.3) used. The first term on the right-hand side of Eq. (5.4) can be recognized as the continuum (or semiclassical) approximation for the thermodynamic potential that arises from approximating the sums over the discrete energy levels with integrals, and is equivalent to results given in Ref. [4]. The other two terms are the next order corrections to the leading order behaviour. As we will discuss below, even in situations where the continuum approximation is valid, the leading order term may not be accurate if the trap becomes too anisotropic. If λ becomes large, it is possible for the next to leading order terms, as counted in terms of µ, to 32 become non-negligible. 5.1.1 λ irrational When λ is irrational it is easy to see from Eq. (2.15) that the integrand of Eq. (2.8) has double poles at β = βk = 2πik ω (5.5) and simple poles at β = β̃k = 2πik λω (5.6) where k = ±1, ±2, . . . in both cases. (Because λ is assumed irrational here, the β̃k poles are all distinct from the βk poles.) It is a straightforward calculation to show that these poles make a contribution to Z(E), that we will call Zr (E), of ∞ n 1 E 1 X E sin 2πk Zr (E) = − 2 2 4π k=1 k sin(πkλ) ω o  h ω i λω E + + cot(πkλ) cos 2πk πk 2  ω  2πkE ∞ cos − πk λω X λω + 2 . 2 2 8π k=1 k sin (πk/λ)   (5.7) This can be used in Eq. (5.1), and again making the approximation T ≪ µ we find Ωr ≃ ∞ nh i 1 1 X 1 + πλ cot(πkλ) + θ coth θ k k 4πβ k=1 k 2 sin(πkλ) sinh θk     µ µ µ o × cos 2πk + 2πk sin 2πk ω ω ω   2πkµ ∞ cos λω − πk 1 X , − 4β k=1 k sin2 (πk/λ) sinh(θk /λ) 33 (5.8) where 2π 2 k θk = . βω (5.9) At this stage we have not made any assumptions about the magnitude of βω in relation to βµ. A simple observation that can be made here is that for small values of βω the presence of the hyperbolic functions in Eq. (5.8) tends to suppress the contribution of Ωr to Ω. Thus the oscillatory part of Ω would be expected to become less important as the temperature is increased (say T ∼ ω), and more important at lower temperatures (T < ω). For values of T larger than ω the leading part of Ω would be expected to be well described by the continuum approximation. given by Ω0 . This will be verified in the next two sections. 5.1.2 λ rational For rational values of λ we will define, without any loss of generality, λ= p , q (5.10) where p and q are two integers that are relatively prime. (That is, p and q have no non-trivial common factors.) The isotropic case considered above is a special case of this, corresponding to p = q = 1. The poles of the integrand of Eq. (2.8) still occur at values of β equal to βk and β̃k given in Eq. (5.5) and Eq. (5.6), but this time some of the values of k result in some of the β̃k poles coinciding with those of βk . It is easy to show that Z(β) has triple poles at 2πiqk/ω for k = ±1, ±2, . . .; double poles at 2πik/ω for k = ±1, ±2, . . ., but k 6= ±q, ±2q, . . .; simple poles at 2πqki/(pω) for k = ±1, ±2, . . . but 34 k 6= ±p, ±2p, . . .. After some calculation, assuming that µ ≫ T as before, it can be shown that Ωr ≃ ∞ nh 1X µ2 p2 + 2q 2 − β k=1 2ω 2 pk sinh(qθk ) 24pq 2 k sinh(qθk )  π2  1 2 + 2β 2 ω 2pk sinh(qθk ) sinh3 (qθk ) i 1 cosh(qθk ) cos(2πqkµ/ω − πkp) − − 2βpqk 2 sinh2 (qθk ) 4π 2 pq 2 k 3 sinh(qθk ) o h µ qθk cosh(qθk ) i 1 − sin(2πqkµ/ω − πkp) + 2πpqωk 2 sinh(qθk ) sinh2 (qθk ) ∞ ′′ n 1 X 1 [1 + πλ cot(πλk) + θk coth θk ] + 4πβ k=1 k 2 sin(πλk) sinh θk o 2πkµ × cos(2πkµ/ω) + sin(2πkµ/ω) ω   2πkµ ′ ∞ − πk cos X 1 λω . (5.11) − 2 4β k=1 k sin (πk/λ) sinh(θk /λ) − Here ∞ ′′ X denotes that terms with k = q, 2q, 3q, . . . are omitted from the sum, k=1 ∞ ′ and X denotes that the terms with k = p, 2p, 3p, . . . are omitted. In the k=1 isotropic case, where p = q = 1, these last two restricted sums both vanish, and we recover the result for the isotropic case above. 5.2 Particle number and chemical potential Having obtained the thermodynamic potential we can use it to calculate all of the thermodynamic properties of interest. In this section we will study the chemical potential and its relation to the particle number. This can be obtained by solving (2.27) for µ. If we split Ω = Ω0 + Ωr as in (2.17) with Ω0 given by Eq. (2.7) and Ωr given by either Eq. (2.11) or Eq. (2.14), we can write N = N0 + Nr in an obvious 35 way to find N0 = i µ3 µ h 2 2 2 2 . (2 + λ )ω − 4π T − 6λω 3 24λω 3 (5.12) By temporarily neglecting the contribution Nr to N coming from Ωr , we obtain µ 4π 2 T 2 1 2 + λ2 − (6λN)−1/3 , ≃ (6λN)1/3 + ω 12 ω2 ! (5.13) valid for large values of λN. For λN ≫ 1 we have µ ≫ ω. The result found in Eq. (5.13) contains the continuum approximation of Ref. [4] and the next to leading order correction. We will see that the next to leading order term can have a significant effect for some traps. The approximation in Eq. (5.13) has neglected any effects that arise from the de Haas-van Alphen part of the thermodynamic potential. As in the isotropic case discussed in Sec. 3, the de Haas-van Alphen terms can make a significant alteration to the continuum result for the chemical potential for certain temperatures, and lead to the step-like behaviour discovered in the numerical results of Ref. [5]. Because of the different expressions for Ω found in the last section depending on the rational or irrational nature of λ, we will look at two separate cases. 5.2.1 λ rational With λ defined as in Eq. (2.13) we find, Nr ≃ ∞ nh qµ2 p2 + 2q 2 π X − βω k=1 pω 2 sinh(qθk ) 12pq sinh(qθk ) + i π2q  2 1 sin(2πqkµ/ω − πkp) + pβ 2 ω 2 sinh(qθk ) sinh3 (qθk ) 36 o 2πqµ cosh(qθk ) cos(2πqkµ/ω − πkp) pβω 2 sinh2 (qθk ) ∞ ′′ n 1 X 1 + [πλ cot(πλk) + θk coth θk ] sin(2πkµ/ω) 2βω k=1 k sin(πλk) sinh θk o 2πkµ − cos(2πkµ/ω) ω   ′ ∞ − πk sin 2πkµ 1 X λω . (5.14) − 2βλω k=1 sin2 (πk/λ) sinh(θk /λ) + For µ ≫ ω the leading contribution would be expected to be given simply by the first term in the first sum on the right hand side of Eq. (5.14). It is the presence of the trigonometric functions in Eq. (5.14) that result in the step-like behaviour of the chemical potential. Because of the presence of sinh(qθk ) in the summand, where θk = 2π 2 kT /ω, it can be seen that the sums converge very rapidly if T ∼ ω. For T close to ω very good analytical approximations can be achieved by merely taking the first term of the sums. Because the temperature dependence of the sums is dictated by the magnitude of exp(−2π 2 qT /ω), the oscillations, although present, would be expected to be negligible and unobservable for T close to ω. It is only once T is reduced below ω that the oscillatory, or step-like, behaviour will become observable. For values of T < ω, it is possible to rewrite the sums that occur in Eq. (5.14) in a way that is much more rapidly convergent. (As T becomes smaller than ω, the sums in Eq. (5.14) converge much more slowly rendering their evaluation more difficult. In addition, it is not sufficient to retain the first term, or even few terms, and achieve a reliable analytical approximation.) One way to obtain a reliable analytical approximation for T small in comparison with ω is to make use of the Poisson summation formula, familiar from the well-known Jacobi inversion formula for the θ-function [21]. It can be shown that 37 ∞ X 1 1 1 sin(2πkm − πkp) = tanh (m̄ − p̄/2) − (m̄ − p̄/2) 2 sinh(2π bk) 4πb 2b 2πb k=1 +  sinh h 1 (m̄ b i ∞  − p̄/2) X 2πb k=1  1 cosh(k/b) + cosh (m̄ − p̄/2) b  −1 , (5.15) where m̄ = m − [m] , p̄ = p mod 2 . (5.16) (5.17) Here [x] denotes the greatest integer that does not exceed x. The result in Eq. (5.15) is sufficient to evaluate the leading term on the right hand side of Eq. (5.14). Other sums that occur in Eq. (5.14) may be evaluated in a similar manner. The advantage of this transformation is that for small values of b (where for Eq. (5.14) we have b ∝ T /ω) the slowly convergent sum on the left hand side of Eq. (5.15) is converted into a rapidly convergent sum on the right hand side. In fact, for b ≪ 1, we can safely neglect the sum that occurs on the right hand side of Eq. (5.15) and obtain a reliable analytical approximation from the first two terms on the right hand side, giving a very simple approximation. It is worth noting that the appearance of m̄ and p̄ as defined in Eqs. (5.16) and (5.17) on the right hand side of Eq. (5.15) is necessary to preserve the periodic structure in m and p found in the original sum on the left hand side; it is incorrect to simply take m and p on the right hand side of Eq. (5.15). An alternate way of obtaining the analytical terms on the right hand side of Eq. (5.15) makes use of the Euler-Maclaurin summation formula. This will not give the last term on the right hand side of Eq. (5.15) as it neglects exponentially small terms. If we take b ≪ 1, and neglect the sum on the right hand side of Eq. (5.15), it is possible to approximate the result in an even simpler way. It is easy to 38 show that ∞ X sin(2πkm − πkp) 1 ≃ {1 − p̄ + 2 ([m + p̄/2] − m)} . 2 sinh(2π bk) 4πb k=1 (5.18) ([x] is the greatest integer function mentioned above.) If we keep only the term of order µ2 on the right hand side of Eq. (5.14), that would be expected to dominate the result for large values of µ ≫ ω, the following simple approximation can be found: Nr ≃ µ2 1 − p̄ + 2 4pω 2   qµ qµ p̄ − + ω 2 ω   . (5.19) For the isotropic oscillator (p = q = 1) this reduces to the result in Eq. (3.20). When T < ω, but not necessarily such that T ≪ ω, it is more reliable to use the result of Eq. (5.15). When the temperature is low enough so that Eq. (5.19) is a reliable approximation, it is possible to solve for µ in terms of N and obtain the following analytical result: µ 1 ≃ ω q  6pq 2 N 1/3 + p̄ 1 − p̄ + 2 2   . (5.20) Because of the greatest integer function in Eq. (5.20) the result shows the step-like behaviour found in the numerical studies of Ref. [5]. In obtaining this result we have approximated Nr with Eq. (5.19) and only included the leading order continuum approximation, the order µ3 term, from Eq. (5.12). This means that the result would be expected to break down in situations where the neglected terms can become comparable in magnitude to those retained. One case of direct physical significance occurs when we allow the trap anisotropy parameter λ to become large. In this case the subleading 39 T = 0.1ω 106 (6λN )1 / 3  1   1 (6λN )1 / 3 1 + (2 + λ2 ) (6λN )−2 / 3  +    2  12  T = 0.05ω 104 µ ω T = 0.01ω 102 100 98 96 8000 8200 8400 8600 8800 9000 9200 9400 9600 9800 10000 N Fig. 4. (color online) This figure shows the chemical potential in units of ω plotted as a function of the particle number N for a range of temperatures T = 0.1ω, 0.05ω, 0.01ω with the anisotropy parameter λ = 20. The curves with T = 0.05ω and 0.01ω have been displaced vertically by 1.5 and 3 units respectively for clarity. The almost straight line shows the continuum approximation (6λN )1/3 , and the step function over the T = 0.1ω curve shows the analytical approximation   1 (6λN )1/3 + 12 (2 + λ2 )(6λN )−1/3 + 12 found from Eq. (5.21) with p = 20, q = 1. (The temperature correction is negligible here and has been omitted.) part of Eq. (5.12) can lead to a non-trivial correction to the standard result of µ/ω ≃ (6λN)1/3 as obtained in Eq. (5.13). It is not difficult to see that a better approximation for the chemical potential, that takes the de Haas-van Alphen part of Ω into account, is µ 1 ≃ ω q 1 4π 2 T 2 1+ q(6λN) 2 + λ2 − (6λN)−2/3 12 ω2 # ) p̄ 1 − p̄ + , + 2 2 (" ! 1/3 ! (5.21) in place of Eq. (5.20). (Recall that λ = p/q here.) For small λ and T ≪ ω 40 the result in Eq. (5.20) is a good approximation. For larger values of λ it is necessary to use the more accurate estimation Eq. (5.21). In Fig. 4 we show the transition from smooth oscillations to step-like behaviour of the chemical potential for a range of N and decreasing temperatures. (Other ranges of N show a similar picture.) The trend towards the step-like behaviour as the temperature is reduced, predicted by the result in Eq. (5.21) is evident. We have chosen a large anisotropy parameter λ = 20 to illustrate how the leading order continuum approximation becomes more inaccurate as the anisotropy increases. Had we plotted the improved expression of Eq. (5.13) it would have been centered on the chemical potential curves; it therefore provides a better estimate for the chemical potential, although still missing the step-like features that are present. For smaller values of the anisotropy parameter, the result (6λN)1/3 becomes more accurate, and the plots look similar to those presented in the isotropic case described earlier. 5.2.2 λ irrational For irrational values of λ we find ∞ n 1 1 X [πλ cot(πλk) + θk coth θk ] sin(2πkµ/ω) 2βω k=1 k sin(πλk) sinh θk o 2πkµ − cos(2πkµ/ω) ω   2πkµ ∞ 1 X sin λω − πk − . (5.22) 2βλω k=1 sin2 (πk/λ) sinh(θk /λ) Nr ≃ This expression is simply the last two terms of Eq. (5.14) with the restrictions on the sums removed. For irrational values of λ, the confluence of poles in Z(β) that gave rise to the triple pole cannot occur here, and this is why the first term of Eq. (5.14) is not present. The leading order contribution to Nr is 41 proportional to µ rather than µ2 as we found in the case of rational λ above. This leads to quantitatively different behaviour for the chemical potential. It has proven very difficult to obtain reliable asymptotic approximations for the sums in Eq. (5.22), even for the leading order part of Nr because of the presence of the sin(πλk) factor in the summand. The techniques that prove so useful in the rational case do not work here, as they all involve the conversion of the sum into an ordinary integral at some stage, and this requires the integral to converge. The sin(πλk) factor prevents this here. As a result, our results are less complete than for the rational case examined earlier. Because of the lack of analytical expressions when λ is irrational we will be brief here and only mention the main conclusions. For values of λ ∼ 1 (we √ used λ = 3) we found that although the oscillations in µ were present, they only became noticeable at very low values of the temperature, typically T = 10−2 ω. This is expected as a consequence of the de Haas-van Alphen contribution to the particle number involving µ, rather than µ2 as in the rational case. In addition the amplitude of the oscillations is much smaller than in the rational case. A conclusion is that the continuum approximation for µ as given by the first term on the right hand side of Eq. 5.13, is more reliable, even for temperatures as low as T = 0.1ω, than in the rational case. √ For larger values of λ (we used λ = 401 to compare with λ = 20 in the rational case) similar conclusions were found provided that the extra term in the continuum approximation as given in Eq. 5.13 was included. 42 5.3 Specific heat The internal energy U is defined in terms of the thermodynamic potential Ω by (2.28). The specific heat C is then defined in terms of U by (2.36). With Ω = Ω0 + Ωr it can be readily seen that U0 = ≃ ∂ (βΩ0 ) ∂β βµ,ω i µ4 µ2 h 2 2 2 2 − (2 + λ )ω − 12π T 8λω 3 48λω 3 (5.23) in agreement with what was found in (3.25). As in the previous section, the result for Ur , coming from taking Ω = Ωr in Eq. (2.28), will be different depending upon whether λ is rational or irrational. The results for Ur will be given explicitly in Secs. 5.3.1 and 5.3.2 below. Although the internal energy U has a simple representation as U = U0 + Ur , the same is not true for the specific heat because of N being held fixed in Eq. (2.36). It is easy to see that Eq. (2.36) can be expressed as in (3.24). The second term on the right hand side of (3.24) frustrates any attempt to express C simply as C = C0 + Cr with Cr the correction to the continuum approximation coming from the de Haas-van Alphen part of Ω. In order to see the effects that the de Haas-van Alphen contributions make to the specific heat we will initially ignore them, and work out just the continuum approximation. We will define C0 to be the continuum approximation for the specific heat where we take U = U0 and N = N0 in Eq. (3.24). Using Eqs. (5.23) and (3.3) we find (for T ≪ ω) π 2 µ2 T . C0 ≃ 6λω 3 (5.24) 43 The conclusion based on the continuum approximation is that at low temperatures the specific heat varies linearly with the temperature. However, as in the isotropic case, and the numerical studies of Ref. [5], this behaviour turns out to be significantly altered once the de Haas-van Alphen part of Ω is included. 5.3.1 λ rational If we concentrate on the leading order contribution to the specific heat for large µ, it would be expected, based on Eq. (5.24), that we would find Cr ∝ µ2 . We will therefore work consistently to order µ2 when we evaluate Cr . If we use Ωr as given in Eq. (2.14) for Ω in the definition Eq. (2.28) we obtain Ur ≃ ∞ n X k=1 πqµ3 sin(2πqkµ/ω − πkp) pβω 3 sinh(qθk ) o 3π 2 qµ2 cosh(qθk ) cos(2πqkµ/ω − πkp) pβ 2 ω 3 sinh2 (qθk ) ∞ ′′ X πµ2 cos(2πkµ/ω) − +··· . βω 2 sin(πλk) sinh θk k=1 + (5.25) Terms of order µ or less have not been kept in this expression. (Recall that ′′ on the second sum means that the terms with k = q, 2q, . . . are omitted and that λ = p/q here. θk was defined in Eq. (2.12).) ! ∂Ur , that is needed to find the complete expression When we compute ∂T µ,ω for C, it can be seen that there will be a term of order µ3 present rather than µ2 as expected; however, it turns out that the term of order µ3 coming from the first term on the right hand side of Eq. (3.24) cancels with an identical contribution coming from the second term, leaving C ∝ µ2 as predicted. After some calculation, assuming that µ is large, it can be shown that C≃ π 2 T µ2 n Σ22 o , 1 + Σ − 12 1 6λω 3 Σ3 (5.26) 44 with Σ1 = 12 ∞ X k=1 ( cosh(qθk ) π 2 qk − βω sinh2 (qθk ) 2 1 + 3 sinh(qθk ) sinh (qθk ) ! !) 2πqkµ − πkp , × cos ω ( ) ! ∞ X 1 2πqkµ cosh(qθk ) Σ2 = sin − qθk − πkp , ω sinh2 (qθk ) k=1 sinh(qθk ) ∞ 2πqkµ 4π 2 q X k cos − πkp Σ3 = 1 + βω k=1 sinh(qθk ) ω ! . (5.27) (5.28) (5.29) It is noteworthy that the dependence on the restricted summation in Eq. (5.25) that could have occurred in the end result has cancelled. The continuum approximation of Eq. (5.24) is regained if we set Σ1 , Σ2 and Σ3 to zero; this is expected because this is equivalent to dropping the contribution from the de Haas-van Alphen part of Ω. For very low temperatures it is possible to obtain reliable analytic approximations for Σ1,2,3 . One way to do this is to use the Poisson summation formula as we did in Sec. 3.1 to approximate Nr . If we assume that qT ≪ ω, it is possible to derive the following approximations: 1 + Σ1 ≃ 4π 2 q cosh2 Σ2 ≃ − Σ3 ≃  3β 3 ω 3 m − 1q [qm] −  βω 2    p̄ 2 2q m − 1q [qm] − β 2 ω 2 m − 1q [qm] − 8πq cosh2 4q cosh 2   βω 2 βω 2   p̄ 2q p̄ 2q  m − 1q [qm] −  p̄ 2q βω m − 1q [qm] − p̄ 2q  ,  ., (5.30) , (5.31) (5.32) As before, [x] is the greatest integer function and p̄ = p mod 2. We have not written out the terms in Eqs. (5.30–5.32) that behave like exp(−α/T ) for constant α as T → 0. It is now possible to see that as T becomes much smaller than ω/q the specific heat behaves like T exp(−α/T ), demonstrating 45 analytically that there is a deviation from the linear dependence predicted by the continuum approximation. As T increases, the terms in Σ1,2,3 in Eq. (5.26) begin to decrease in importance and the linear behaviour starts to become a good approximation. T = 0.5ω 400 T = 0.4ω T = 0.3ω 300 C T = 0.2ω 200 100 T = 0.1ω 0 8000 8200 8400 8600 8800 9000 9200 9400 9600 9800 10000 N Fig. 5. (color online) This figure shows the specific heat plotted as a function of the particle number N over a range of N . The different curves show the results for five different temperatures. The anisotropy parameter λ = 20. The behaviour of the specific heat with particle number and with temperature are shown in Figs. 5 and 6 respectively. In Fig. 5 the oscillations in the specific heat are clearly visible. The amplitude of the oscillations starts to decrease very rapidly as the temperature is increased. Once T becomes greater than about 0.5ω they all but disappear. As the temperature is reduced there is the competing effect that the specific heat starts to vanish exponentially fast as discussed above and illustrated in Fig. 5. Thus there is really a narrow range of temperature when both the specific heat and the amplitude of the oscillations are large enough to be potentially observable. In Fig. 5 it can be 46 500 400 300 C 200 100 0 0.1 0.2 0.3 0.4 0.5 T /ω Fig. 6. (color online) This figure shows the specific heat plotted as a function of the temperature for λ = 20 and N = 105 both fixed. The straight line indicates the linear π 2 N 2/3 T temperature dependence given by the continuum approximation C = . The (6λ)1/3 ω other curve shows the effect if the de Haas-van Alphen contribution is included. noted that once T is about 0.1ω there is a dramatic fall-off in the specific heat from the expected linear behaviour that is fully consistent with the analytical expressions obtained above. 5.3.2 λ irrational As in the calculation of the chemical potential, we were less successful at obtaining simple analytical expressions here than in the case of rational λ. Using Eq. (5.8) for Ωr in Eq. (2.28) we find Ur ≃ nh 1 2π 3 k π3k + β 3 ω 2 sinh3 θk β 3 ω 2 sinh θk k=1 k sin(πkλ) ∞ X + π 2 λ cot(πkλ) cosh θk πkµ2 i − cos(2πkµ/ω) βω 2 sinh θk 2β 2 ω sinh2 θk 47 250 200 150 C 100 T = 0.5ω 50 T = 0.3ω T = 0.15ω 0 20 40 60 80 100 p Fig. 7. (color online) This figure shows the specific heat plotted as a function of the anisotropy parameter with N = 103 held fixed. λ = p here with p an integer. The result is shown for three different temperatures T = 0.15ω, 0.3ω and 0.5ω.  2π 2 k cosh θk i µ πλ cot(πkλ) + βω sinh θk β 2 ω 2 sinh2 θk 2βω sinh θk × sin(2πkµ/ω) ∞ n   2πkµ X πµ − − πk sin 2 λω k=1 2λβω sin (πk/λ) sinh(θk /λ) + h π 2 kµ cosh θ k + +  2πkµ o π 2 cosh(θk /λ) cos − πk . (5.33) λω 2λβ 2 ω sin2 (πk/λ) sinh2 (θk /λ) This expression can be contrasted with that found in the rational case in Eq. (5.25). The de Haas-van Alphen contribution to the internal energy begins at order µ2 rather than µ3 . This has a direct consequence for the specific heat. If we use Eq. (3.24) to calculate the specific heat to leading order we find that to leading order in µ the specific heat is given by the continuum expression in Eq. (5.24). For irrational values of λ the de Haas-van Alphen contribution cancel out of the specific heat to leading order. There is a contribution from the de Haas-van Alphen part of the internal energy at sub-leading order µ; 48 however we will not give the expression here as it has only a negligible effect. 6 Discussion and conclusions The main conclusion of this paper is that the origin of the step-like features found in the numerical calculations of Schneider and Wallis [5] have the same origin as the periodicity found in the de Haas-van Alphen effect and that they can be approximated by analytical means. The general approach of Sondheimer and Wilson [18], so useful in the analysis of the de Haas-van Alphen effect, can be used to great effect here to obtain analytical results for various thermodynamic quantities. In addition, it is possible to see that within this approach the continuum approximation used by Butts and Rokhsar [4] corresponds to the neglect of an infinite set of poles in the Laplace transform of the partition function. It is also possible to recover the continuum approximation using the methods described in the general framework of Ref. [22]. In Sec. 3 we obtained approximate analytical results for the 3-dimensional, isotropic harmonic oscillator potential. In the case T ≪ ω, a very simple result (Eq. (3.21)) was found for the chemical potential that clearly exhibited the step-like features found in Ref. [5]. This was used to evaluate the specific heat, and it was shown (in agreement with Ref. [5]) that there were significant deviations from the linear temperature dependence predicted by the continuum approximation. The specific heat was seen to exhibit a periodic structure in the particle number. In Sec. 4 we studied the trapped Fermi gas in 1 and 2 spatial dimensions. For the 1-dimensional gas we were able to show that as T → 0 the specific 49 heat vanished exponentially fast in the inverse temperature (Eq. (4.9)). The 2-dimensional gas was similar in many ways to the 3-dimensional case of Sec. 3. The same methods that were applied in the isotropic case can be used to examine anisotropic potentials. The technical details are more involved due to the structure of the poles in the inverse Laplace transform of the partition function as described above. In addition to a periodicity in the particle number, there can also be a periodicity when the trapping potential is altered. This latter effect was not found for the isotropic case and it is worth commenting on why this occurs, in contrast with what might be expected from the de Haasvan Alphen effect where the thermodynamics shows a periodicity in the inverse magnetic field strength. The difference between the two cases is related to the relationship between the chemical potential and the particle number. In both cases (de Haas-van Alphen and trapped Fermi gas) the periodic structure results from a trigonometric dependence on µ/ω. (For the de Haas-van Alphen effect ω is the cyclotron frequency associated with the magnetic field strength.) In the de Haas-van Alphen effect, when µ is solved for in terms of the particle number the leading order contribution turns out to be independent of ω. Thus the trigonometric functions that involve µ/ω exhibit the familiar de Haasvan Alphen oscillations as the magnetic field is varied. For the trapped Fermi gas in an isotropic potential we found µ ∝ ω, so that µ/ω is independent of ω, although there is a dependence on the particle number. We claim that this is an artifact of the simplicity of the isotropic potential, and that the situation for anisotropic potentials is more interesting. The treatment presented here has only been performed for the free Fermi gas. As already mentioned in the introduction, the neglect of interactions for a single component gas is a good approximation [4,5]. In Ref. [6] a two component 50 model was studied that allowed for an interaction between the two spin components. A numerical study showed that as the interaction strength is increased the step-like behaviour, like that we have been discussing, is increasingly suppressed. For the single component gas, we expect that an analysis similar to that given by Luttinger [23] in the de Haas-van Alphen effect could be used to show that both the amplitude and period of the oscillations are not affected by interactions to leading order. (This conclusion does not hold for the non-oscillatory part of Ω, that we have called Ω0 .) It would be of interest to study the Luttinger analysis for a multi-component Fermi gas in more detail to make contact with the results of Ref. [6]. References [1] M. H. Anderson, J. R. Ensher, M. R. Matthews, C. E. Wieman, and E. A. Cornell, Science 269, 198 (1995); C. C. Bradley, C. A. Sackett, J. J. Tollett, and R. G. Hulet, Phys. Rev. Lett. 75, 1687 (1995); K. B. Davis, M.-O. Mewes, M. R. Andrews, N. J. van Druten, D. S. Durfee, D. M. Kurn, and W. Ketterle, Phys. Rev. Lett. 75, 3969 (1995). [2] C. J. Pethick and H. Smith, Bose-Einstein Condensation in Dilute Gases, (Cambridge University Press, Cambridge, 2002). [3] B. DeMarco and D. S. Jin, Science 285, 1703 (1999); A. G. Truscott, K. E. Strecker, W. I. McAlexander, G. B. Partridge, and R. G. Hulet, Science 291, 2570 (2001); F. Schreck, L. Khaykovich, K. L. Corwin, G. Ferrari, T. Bourdel, J. Cubizolles, and C. Salomon, Phys. Rev. Lett. 87, 080403 (2001); S. R. Granade, M. E. Gehm, K. M. O’Hara, and J. E. Thomas, Phys. Rev. Lett. 88, 120405 (2002); Z. Hadzibabic, C. A. Stan, K. Dieckmann, S. Gupta, M.W. Zwierlein, 51 A. Gorlitz, and W. Ketterle, Phys. Rev. Lett. 88, 160401 (2002); G. Roati, F. Riboli, G. Modugno, and M. Inguscio, Phys. Rev. Lett. 89, 150403 (2002). [4] D. A. Butts and D. S. Rokhsar, Phys. Rev. A 55, 4346 (1997). [5] J. Schneider and H. Wallis, Phys. Rev. A 57, 1253 (1998). [6] G. M. Bruun and K. Burnett, Phys. Rev. A 58, 2427 (1998). [7] F. Brosens, J. T. Devreese, and L. F. Lemmens, Phys. Rev. E 57, 3871 (1998). [8] P. Vignolo, A. Minguzzi, and M. P. Tosi, Phys. Rev. Lett. 85, 2850 (2000). [9] F. Gleisberg, W. Wonneberger, U. Schlöder, and C. Zimmermann, Phys. Rev. A 62, 063602 (2000). [10] M. Brack and B. P. van Zyl, Phys. Rev. Lett. 86, 1574 (2001). [11] X. Z. Wang, J. Phys. A 35, 9601 (2002). [12] B. P. van Zyl, R. K. Bhaduri, A. Suzuki, and M. Brack, Phys. Rev. A 67, 023609 (2003). [13] M. Brack and M. V. N. Murthy, J. Phys. A 36, 1111 (2003). [14] B. P. van Zyl, Phys. Rev. A 68, 033601 (2003). [15] D. J. Toms, J. Phys. A 37, 3111 (2004). [16] L. Landau, Z. Phys. 64, 629 (1930). [17] N. W. Ashcroft and N. D. Mermin, Solid State Physics, (Holt, Rinehart, and Winston, Philadelphia, 1976). [18] E. H. Sondheimer and A. H. Wilson, Proc. Roy. Soc. (London) A210, 173 (1951). [19] K. Huang, Statistical Mechanics, (J. Wiley, New York, 1987). 52 [20] R. K. Pathria, Statistical Mechanics, (Pergamon, London, 1972). [21] E. T. Whittaker and G. N. Watson, A Course of Modern Analysis, (Cambridge University Press, Cambridge, 1928). [22] K. Kirsten and D. J. Toms, Phys. Lett. A 222, 148 (1996); K. Kirsten and D. J. Toms, Phys. Rev. E 59, 158 (1999). [23] J. M. Luttinger, Phys. Rev. 121, 1251 (1961). 53