On the derivation of the entropy of ideal quantum gases confined in a three-dimensional harmonic pot
本站小编 Free考研考试/2022-01-02
Le Tri Dat,, Vinh N T Pham,41Laboratory of Applied Physics, Advanced Institute of Materials Science, Ton Duc Thang University, Ho Chi Minh City 756636, Vietnam 2Faculty of Applied Sciences, Ton Duc Thang University, Ho Chi Minh City 756636, Vietnam 3Ho Chi Minh City University of Education, Ho Chi Minh City 748342, Vietnam
First author contact:4Author to whom any correspondence should be addressed. Received:2019-10-2Revised:2020-02-3Accepted:2020-02-12Online:2020-03-27
Abstract In this work, the entropy functions of ideal quantum gases in a three-dimensional harmonic trap are analytically calculated using temperature as an explicit variable. Afterward, the applicability of the analytical formulas is validated by comparison with the numerical calculation. The results illustrate that the obtained functions could be applied for the whole temperature regime with a maximum relative deviation of less than 7.5% in the vicinity of the critical temperature Tc in the case of Bose gases. Meanwhile, for Fermi gases, although the analytical formula fits well at very low- and high-temperature regimes, it cannot be applied at temperature in the range $[0.3-0.5]{T}_{{\rm{F}}}$, where TF is the Fermi temperature. In addition, the consistency between our formulas and classical ones at significantly high temperatures is also discussed. Keywords:thermodynamic properties;entropy;3D harmonic trap;ideal quantum gases
PDF (455KB)MetadataMetricsRelated articlesExportEndNote|Ris|BibtexFavorite Cite this article Le Tri Dat, Vinh N T Pham. On the derivation of the entropy of ideal quantum gases confined in a three-dimensional harmonic potential. Communications in Theoretical Physics, 2020, 72(4): 045701- doi:10.1088/1572-9494/ab76fd
1. Introduction
The observation of the Bose–Einstein condensation (BEC) in dilute 85Rb gas has opened a new era in the research of ultra-cold atoms [1] and motivated subsequent experiments on BEC in other bosonic atoms [2–8], molecules [9], and exciton-polariton [10]. By using a similar cooling technique as that of bosonic atoms, physicists also achieved interesting phenomena with fermions such as the onset of observing 40K at $0.5{T}_{{\rm{F}}}$ [11] (TF is the Fermi temperature), the observation of Fermi pressure [12], and of strong interacting of degenerate Fermi gas [13]. These ultra-cold atoms play a crucial role in modern physics since they can be used to interpret many interesting phenomena.
Thermodynamic quantities such as the chemical potential, the total energy, and the heat capacity of trapped quantum gases have been extensively studied by several authors [14–23]. However, the entropy, which plays a crucial role in the interpretation of phase transition in physics, has not been mentioned in detail, especially for Fermi gas where other functions are usually considered [15, 19, 21]. In comparison to the entropy of harmonically trapped Fermi gases, that of Bose gases receives less concern [14, 16, 18, 22, 23]. Grether et al have recently examined the properties of entropy in such a system; however, they expressed the results in terms of polylogarithm functions of chemical potential [17]. In fact, in experiments, such property is controlled using temperature as variable. Therefore, it is vital to conduct a comprehensive study to better understand the entropy of ideal quantum gases confined harmonically.
In this work, we propose a mathematical scheme to obtain the entropy of a harmonically trapped ideal quantum gas as a function of temperature rather than the chemical potential since this expression is more natural and closed to experiments. We establish the grand potential (Ω) of N ideal bosons and fermions trapped in a three-dimensional (3D) harmonic trap. Then, using the relation between the entropy and Ω in the grand canonical ensemble, the entropy functions are analytically derived. Finally, a numerical calculation is conducted in order to validate the correctness and the applicable range of the obtained analytical formulas.
The article is organized as follows. In section 2, we introduce the mathematical scheme to derive analytical formulas that describe the entropy of N non-interacting Bosons or Fermions trapped in a 3D harmonic potential. In addition, the numerical results are also presented to assess the exactness of the analytical formula as well as the applicable range. Section 3 summarizes and concludes the work. Throughout this paper, the term NkB where N is the number of particles in the system and kB denotes the Boltzmann constant is chosen as the unit of the entropy, the critical temperature Tc and the Fermi temperature TF are chosen as the units of the temperature associated with the bosonic and fermionic systems, respectively.
2. Results and discussion
2.1. Ideal Bose gases
For an ideal Bose gas confined harmonically, the grand potential is given by (see details in appendix)$\begin{eqnarray}{{\rm{\Omega }}}_{{\rm{Bose}}}=-\displaystyle \frac{1}{6{\left({\hslash }\omega \right)}^{3}}{\int }_{0}^{\infty }\displaystyle \frac{{\varepsilon }^{3}{\rm{d}}\varepsilon }{\exp [\beta (\varepsilon -\mu )]-1},\end{eqnarray}$where ℏ is the reduced Planck constant, ϵ and ω denote the energy and the frequency of the trap, respectively, and $\beta =1/({k}_{{\rm{B}}}T)$. According to the grand canonical ensemble, the entropy S can be obtained by taking derivative the grand potential with respect to the absolute temperature$\begin{eqnarray}\begin{array}{rcl}{S}_{{\rm{Bose}}} & = & -{\left(\displaystyle \frac{\partial {{\rm{\Omega }}}_{{\rm{Bose}}}}{\partial T}\right)}_{\mu }\\ & = & \displaystyle \frac{{k}_{{\rm{B}}}{\beta }^{2}}{6{\left({\hslash }\omega \right)}^{3}}{\int }_{0}^{\infty }\displaystyle \frac{{\varepsilon }^{3}\exp [\beta (\varepsilon -\mu )](\varepsilon -\mu )}{{\left(\exp [\beta (\varepsilon -\mu )]-1\right)}^{2}}{\rm{d}}\varepsilon .\end{array}\end{eqnarray}$
Owing to the phase transition at a critical temperature Tc of bosons (the BEC), the thermodynamic properties are treated with two separate phases, condensed and normal phases [14, 16, 23]. The entropy S, therefore, is also proceeded. Such phases are separated by the critical temperature ${T}_{c}$ where all bosons are in the excited states with the maximum chemical potential μ=0 [24–26]. The critical temperature has been derived in [14, 16, 22–26]$\begin{eqnarray}{T}_{c}=\displaystyle \frac{{\hslash }\omega }{{k}_{{\rm{B}}}}\sqrt[3]{N\zeta (3)},\end{eqnarray}$where N, which is fixed, is the total number of particles and ζ(3) denotes the Riemann zeta function. In the condensed phase, T<Tc, the chemical potential remains zero, thus we have$\begin{eqnarray}{S}_{{\rm{Bose}}}=\displaystyle \frac{{k}_{{\rm{B}}}}{6{\left(\beta {\hslash }\omega \right)}^{3}}{\int }_{0}^{\infty }\displaystyle \frac{{\left(\beta \varepsilon \right)}^{4}\exp (\beta \varepsilon )}{{[\exp (\beta \varepsilon )-1]}^{2}}{\rm{d}}(\beta \varepsilon ).\end{eqnarray}$
Integrating equation (4), the entropy in condensed phase is straightforwardly obtained$\begin{eqnarray}\displaystyle \frac{{S}_{{\rm{Bose}}}}{{{Nk}}_{{\rm{B}}}}=\displaystyle \frac{4\zeta (4)}{\zeta (3)}{\left(\displaystyle \frac{T}{{T}_{c}}\right)}^{3}.\end{eqnarray}$
In the normal phase (high-temperature regime) the Taylor expansion is used to approximately compute equation (2)$\begin{eqnarray}\displaystyle \frac{1}{\exp (x)-1}\approx \exp (-x)+\exp (-2x)+\ldots .\end{eqnarray}$
Although the detailed explanation of the term $\exp (\beta \mu )$ in normal phase has been given in the [23], we still briefly present here the procedure to derive this form from the grand potential. According to the grand canonical ensemble, the total number of particles is derived by partially differentiating the grand potential Ω with respect to the chemical potential μ, hence$\begin{eqnarray}N=-\displaystyle \frac{\partial {{\rm{\Omega }}}_{{\rm{Bose}}}}{\partial \mu }=\displaystyle \frac{1}{2{\left({\hslash }\omega \right)}^{3}}{\int }_{0}^{\infty }\displaystyle \frac{{\varepsilon }^{2}}{\exp \left[\beta (\varepsilon -\mu )\right]-1}{\rm{d}}\varepsilon .\end{eqnarray}$
Again, we approximately expand the above equation by equation (6) up to the second order due to the complexity of solving third- and higher-order terms$\begin{eqnarray}\exp (\beta \mu )=-4+\sqrt{16+8\zeta (3){\left(\displaystyle \frac{{T}_{c}}{T}\right)}^{3}}.\end{eqnarray}$
Finally, we obtain the formula describing the dependence of entropy on the absolute temperature in the normal phase$\begin{eqnarray}\displaystyle \frac{{S}_{{\rm{Bose}}}}{{{Nk}}_{{\rm{B}}}}=4-\displaystyle \frac{\zeta (3)}{2}{\left(\displaystyle \frac{{T}_{c}}{T}\right)}^{3}-\mathrm{log}\left[-4+\sqrt{16+8\zeta (3){\left(\displaystyle \frac{{T}_{c}}{T}\right)}^{3}}\right].\end{eqnarray}$
We compare the analytical results to the numerical results in order to validate the correctness and the applicable range of the approximated formulas. In our work, the numerical results are obtained by firstly integrating the equation (A8) for whole temperature regime, then numerically computing the entropy using the relation $S=-{(\partial {\rm{\Omega }}/\partial T)}_{\mu }$ for both Fermi and Bose systems. The detailed results are depicted in figure 1. It is obvious that in the condensed phase, the result computed by the analytical method is well-matched to the numerical calculation. The relative deviation is nearly zero, due to the fact that in this temperature regime, equation (4) can be integrated exactly. In the case of normal phase, in high-temperature regime, the approximated formula describes well the dependence of the entropy on the temperature. While the relative deviation between the approximated formula and the numerical calculation in the vicinity of the critical temperature Tc is quite large, around 7.5%; however, it is still an acceptable deviation. The relative error between numerical and analytical approaches in the normal phase is natural and mainly due to the Taylor expansion. This expansion is considered to work well as $\beta \mu \ll 1$ associated with sufficiently high temperature. Thus the deviation is large around Tc where the temperature is not adequately high and rapidly decreases as the temperature increases as expected. We will correct the formula to achieve a better match at Tc in our next work.
Figure 1.
New window|Download| PPT slide Figure 1.(a) Entropy of an ideal Bose gas in a 3D harmonic trap as a function of the temperature and (b) the relative deviation between the numerical and approximated results. The black solid curve, red squares, and blue circles indicate the results obtained from the numerical results and analytical results in condensed phase and normal phase, respectively.
Obviously, once the temperature is significantly high, quantum gases behave classically. Therefore, we assume that the temperature is high enough to compute the limit of equation (10), then compare to the classical formula derived by Ligare [27] in order to validate our formula, hence$\begin{eqnarray}{\left(\displaystyle \frac{{S}_{{\rm{Bose}}}}{{{Nk}}_{{\rm{B}}}}\right)}_{{\rm{classical}}}=4+\mathrm{log}\left[\displaystyle \frac{1}{\zeta (3)}{\left(\displaystyle \frac{T}{{T}_{c}}\right)}^{3}\right].\end{eqnarray}$
Substituting the definition of the critical temperature Tc in equation (3) back into equation (11), we can rederive the formula for the entropy of a classical gas trapped in a 3D harmonic potential, see equation (24) in [27] which is$\begin{eqnarray}{\left(\displaystyle \frac{S}{{{Nk}}_{{\rm{B}}}}\right)}_{{\rm{Ligare}}}=3+\mathrm{log}\left[\displaystyle \frac{e}{N}{\left(\displaystyle \frac{U}{3{\hslash }\omega }\right)}^{3}\right],\end{eqnarray}$where $U=3{{Nk}}_{{\rm{B}}}T$ is the total energy of a classical gas in a 3D harmonic trap. For significantly high temperature, the Bose gas behaves as classical one, thus the deviation between the analytical formula and numerical result has to approach to zero as expected.
2.2. Ideal Fermi gas
In this section, the entropy function of a system of N ideal fermions confined in a 3D harmonic trap is taken into account. Traditionally, the temperature of a Fermi gas is considered high if $T\gg {T}_{{\rm{F}}}$ and low if $T\ll {T}_{{\rm{F}}}$. It is important to note that the Fermi temperature is the temperature scale of the Fermi energy that is defined mathematically as the chemical potential at T=0 or the energy level of the last-filled state [24, 25]. In case of fixed N ideal fermions trapped in a 3D harmonic potential, the Fermi temperature is given by$\begin{eqnarray}{T}_{{\rm{F}}}=\displaystyle \frac{{\hslash }\omega }{{k}_{{\rm{B}}}}\sqrt[3]{6N}.\end{eqnarray}$
Different from the case of Bose gas where the chemical potential μ equals to zero in low-temperature regime corresponding to condensed phase, for Fermi gas, $\mu \ne 0$ for whole range of temperature since there is no existence of condensation. Hence it is challenging to correctly derive the formula for the entropy at low temperature in this case. Instead, we use the relation$\begin{eqnarray}S={\int }_{0}^{T}\displaystyle \frac{C({T}^{{\prime} })}{{T}^{{\prime} }}{{\rm{d}}{T}}^{{\prime} },\end{eqnarray}$to compute the entropy in lieu of using the Sommerfeld expansion since we have already known the heat capacity [15, 17, 19, 21] giving$\begin{eqnarray}\displaystyle \frac{C}{{{Nk}}_{{\rm{B}}}}={\pi }^{2}\displaystyle \frac{T}{{T}_{{\rm{F}}}}.\end{eqnarray}$
Using equations (14) and (15), the fermion entropy at low-temperature regime is obtained$\begin{eqnarray}\displaystyle \frac{{S}_{{\rm{Fermi}}}}{{{Nk}}_{{\rm{B}}}}={\pi }^{2}\displaystyle \frac{T}{{T}_{{\rm{F}}}}.\end{eqnarray}$
We proceed to derive the fermion entropy in high-temperature regime. The grand potential of a Fermi gas trapped harmonically reads$\begin{eqnarray}\begin{array}{rcl}{{\rm{\Omega }}}_{{\rm{Fermi}}} & = & -\displaystyle \frac{1}{{\left(\beta {\hslash }\omega \right)}^{3}}{{\rm{Li}}}_{4}\left[-\exp (\beta \mu )\right]\\ & = & \displaystyle \frac{1}{{\left(\beta {\hslash }\omega \right)}^{3}}{f}_{4}\left[\exp (\beta \mu )\right].\end{array}\end{eqnarray}$where fν(z) denotes the Fermi–Dirac integrals, giving$\begin{eqnarray}{f}_{v}[z]=\displaystyle \frac{1}{{\rm{\Gamma }}(\nu )}{\int }_{0}^{\infty }\displaystyle \frac{{x}^{v-1}{\rm{d}}{x}}{{z}^{-1}\exp (x)+1}=-\sum _{k=1}^{\infty }\displaystyle \frac{{\left(-z\right)}^{k}}{{k}^{v}}.\end{eqnarray}$
Differentiating the equation (17) with respect to the chemical potential μ, and combining with the relation between the Fermi–Dirac integrals$\begin{eqnarray}\displaystyle \frac{\partial {f}_{\sigma }(z)}{\partial z}=\displaystyle \frac{1}{z}{f}_{\sigma -1}(z),\end{eqnarray}$the fermion entropy is given by$\begin{eqnarray}{S}_{{\rm{Fermi}}}=4\displaystyle \frac{{f}_{4}[\exp (\beta \mu )]}{{f}_{3}[\exp (\beta \mu )]}-\mathrm{log}[\exp (\beta \mu )].\end{eqnarray}$
The Fermi–Dirac integrals in equation (18) are expanded up to second order by equation (18), hence$\begin{eqnarray}{S}_{{\rm{Fermi}}}=4\displaystyle \frac{\exp (\beta \mu )-\tfrac{{[\exp (\beta \mu )]}^{2}}{16}}{\exp (\beta \mu )-\tfrac{{[\exp (\beta \mu )]}^{2}}{8}}-\mathrm{log}[\exp (\beta \mu )].\end{eqnarray}$
Once again, we need the expression of the term $\exp (\beta \mu )$ on the temperature in high-temperature regime. It is a similar procedure as that for Bose gases, thus, the final result is presented straightforwardly$\begin{eqnarray}\exp (\beta \mu )=\displaystyle \frac{1}{6}{\left(\displaystyle \frac{{T}_{{\rm{F}}}}{T}\right)}^{3}+\displaystyle \frac{1}{288}{\left(\displaystyle \frac{{T}_{{\rm{F}}}}{T}\right)}^{6}.\end{eqnarray}$
Note that the second term in equation (22) is our contribution compared to other studies [15, 19, 21] in order to avoid the accumulated error and then improve the accuracy of the final result. Substituting equation (22) into equation (21), we find the approximated formula describing the fermion entropy in high-temperature regime$\begin{eqnarray}\displaystyle \frac{{S}_{{\rm{Fermi}}}}{{{Nk}}_{{\rm{B}}}}=4+\displaystyle \frac{1}{24}{\left(\displaystyle \frac{{T}_{{\rm{F}}}}{T}\right)}^{3}-\mathrm{log}\left[\displaystyle \frac{1}{6}{\left(\displaystyle \frac{{T}_{{\rm{F}}}}{T}\right)}^{3}+\displaystyle \frac{1}{288}{\left(\displaystyle \frac{{T}_{{\rm{F}}}}{T}\right)}^{6}\right].\end{eqnarray}$
The results of numerical and analytical works from equations (16) and (23) are shown in figure 2. As we could see in figure 2, for $T\lt 0.3{T}_{{\rm{F}}}$, the analytical formula for low-temperature regime has a relative deviation below 10% compared to the numerical result. In this regime, the analytical formula departs from the numerical calculation as the temperature increases. The root of this behavior is the implication of Sommerfeld expansion in considering the total energy of the system before deriving the heat capacity based on the relation $C=\partial E/\partial T$ [28]. This is well-known that Sommerfeld expansion works well for sufficiently large β (low temperature). Thus the relative error of entropy for $T\gt 0.3{T}_{{\rm{F}}}$ in this regime is accumulated from the derivative procedure. Meanwhile, the approximated formula for the fermion entropy in high-temperature regime has the maximum relative error below 10% for $T\gt 0.5{T}_{{\rm{F}}}$, and the deviation also falls rapidly as the temperature rises. It is straightforward to understand the high relative error for $T\lt 0.5{T}_{{\rm{F}}}$. The expansion of $\exp (\beta \mu )$ in equation (22) works well for adequately high temperature. We note that the approximated formulas in this present work are not able to be applied in the temperature interval of $[0.3-0.5]{T}_{{\rm{F}}}$, since the relative deviation is large, over 10%, which means that the correction has to be studied in our next work.
Figure 2.
New window|Download| PPT slide Figure 2.(a) Entropy of an ideal Fermi gas in the harmonic trap as a function of the reduced temperature and (b) the relative error between the exact and approximated results. The black solid curve, red squares, and blue circles indicate the results obtained from the numerical and analytical approach in low- and high-temperature regime, respectively.
Finally, we verify the correctness of the high-temperature fermion entropy function again by computing the limit at higher temperature then compare to classical results in [27]. The fermion entropy in the significantly high-temperature limit is given by$\begin{eqnarray}{\left(\displaystyle \frac{{S}_{{\rm{Fermi}}}}{{{Nk}}_{{\rm{B}}}}\right)}_{{\rm{classic}}}=4+\mathrm{log}\left[6{\left(\displaystyle \frac{T}{{T}_{{\rm{F}}}}\right)}^{3}\right].\end{eqnarray}$
Once again, we also rederive the classical entropy by substituting the definition of the Fermi temperature into the equation (24). Similarly to the case of Bose gas, for significant high temperature, the Fermi gas again acts classically, therefore the deviation between analytical formula and numerical results has to approach to zero as expected.
3. Conclusion
In this paper, we provide a mathematical scheme to derive analytically the grand potential of ideal quantum gases confined harmonically, then the entropy is derived by differentiating the grand potential with respect to the absolute temperature. For ideal Bose gas, since the integrals are exactly calculable in condensed phase, the entropy is derived without any approximation. Meanwhile, the analytical formula describing the entropy in normal phase is approximated by using Taylor expansion thus working perfectly for sufficiently high temperature. Regarding the ideal Fermi gas, we use the relation between the entropy and the heat capacity instead of directly using Sommerfeld approximation to derive the entropy at low temperature. The deviation between analytical formula and numerical result in this regime is due to the nature of Sommerfeld expansion implied in the derivation of total energy of the system. In high-temperature regime, the entropy is treated by using expansion of polylogarithm function, then again has good performance at adequately high temperature. Finally, the analytical entropy formulas are validated by numerical approach. The results indicate that for ideal Bose gas, the analytical formulas could describe well the entropy in both condensed and normal phases with an acceptable deviation (< 7.5%) at close to the critical temperature. Although the approximated formulas of ideal Fermi gas could be applied for very low- and high- temperature regime, they are not able to be used in the temperature interval of $[0.3-0.5]{T}_{{\rm{F}}}$ since the deviation is relatively large. The correction has to be studied in our next work.
Acknowledgments
The authors highly appreciate Mr Tran Duong Anh-Tai in Okinawa Institute of Science and Technology for his fruitful discussion.
Appendix. Derivation of the integral form of the grand canonical potential of ideal quantum gases
We introduce a mathematical scheme for establishing the grand potential of ideal quantum gases in a 3D harmonic potential. The grand potential in the Grand Canonical ensemble can be written in form of$\begin{eqnarray}{\rm{\Omega }}=-\displaystyle \frac{1}{\beta a}\sum _{i=0}^{\infty }\mathrm{log}\left(1+a\exp \left[-\beta ({\varepsilon }_{i}-\mu )\right]\right),\end{eqnarray}$where $a=-1$ for bosons, a=1 for fermions, $\beta =1/{k}_{{\rm{B}}}T,{k}_{{\rm{B}}}$ is the Boltzmann constant, and T is the absolute temperature. We assume that the number of particles in the system is large enough so that the semiclassical approximation can be applied. In the equation (A1), the summation is replaced by ${\int }_{0}^{\infty }g(\varepsilon ){\rm{d}}\varepsilon $ , where g(ϵ) is called the density of state function because we suppose that the system behaves as the semiclassical approximation for the adequately large number of N particles. In the case of a 3D harmonic oscillator, the density of state is given by [26]$\begin{eqnarray}g(\varepsilon )=\displaystyle \frac{{\varepsilon }^{2}}{2{\left({\hslash }\omega \right)}^{3}},\end{eqnarray}$where ℏ is the reduced Planck constant, ϵ and ω denote the energy and the frequency of the trap, respectively. Hence, the equation (A1) is rewritten as follow$\begin{eqnarray}{\rm{\Omega }}=-\displaystyle \frac{1}{\beta a}\displaystyle \frac{1}{2{\left({\hslash }\omega \right)}^{3}}{\int }_{0}^{\infty }{\varepsilon }^{2}\mathrm{log}\left(1+a\exp \left[-\beta (\varepsilon -\mu )\right]\right){\rm{d}}\varepsilon .\end{eqnarray}$
The Taylor–Maclaurin formula of logarithm function that valids for $| x| \lt 1$ is used to expand the equation (A3)$\begin{eqnarray}\mathrm{log}(1+x)=-\sum _{k=1}^{\infty }\displaystyle \frac{{\left(-x\right)}^{{\ell }}}{{\ell }}.\end{eqnarray}$
We consider the integral in the equation (A5) and then the sum is calculated sequentially. Hence, the equation (A5) becomes$\begin{eqnarray}{\rm{\Omega }}=\displaystyle \frac{1}{{\beta }^{4}{\left({\hslash }\omega \right)}^{3}}\sum _{{\ell }=1}^{\infty }\displaystyle \frac{1}{a}\displaystyle \frac{{[-a\exp (\beta \mu )]}^{{\ell }}}{{{\ell }}^{4}}.\end{eqnarray}$
According to the generalized definition of polylogarithm function$\begin{eqnarray}\begin{array}{rcl}-{{a}{\rm{L}}{\rm{i}}}_{\sigma }(-{az}) & = & -\displaystyle \frac{1}{a}\sum _{{\ell }=1}^{\infty }\displaystyle \frac{{\left(-{az}\right)}^{{\ell }}}{{{\ell }}^{\sigma }}\\ & = & \displaystyle \frac{1}{{\rm{\Gamma }}(\sigma )}{\displaystyle \int }_{0}^{\infty }\displaystyle \frac{{x}^{\sigma -1}{\rm{d}}{x}}{{z}^{-1}\exp (x)+a},\end{array}\end{eqnarray}$the infinite summation in equation (A6) can be treated as polylogarithm function $-a{\mathrm{Li}}_{4}[-a\exp (\beta \mu )]$, since integrals do not depend on the variable of integration, we will choose x=βϵfor simplicity. The integral form of the grand potential of ideal quantum gases confined harmonically is obtained$\begin{eqnarray}{\rm{\Omega }}=-\displaystyle \frac{1}{6{\left({\hslash }\omega \right)}^{3}}{\int }_{0}^{\infty }\displaystyle \frac{{\varepsilon }^{3}{\rm{d}}\varepsilon }{\exp [\beta (\varepsilon -\mu )]+a}.\end{eqnarray}$