The unit system

The unit system used in this article is Hartree atomic units: me=kB==4πε0=e=1m_\mrm e=k_\mrm B=\hbar=4\pi\veps_0=e=1, where mem_\mrm e is the electron mass.

In this unit system, the Bohr radius is aB=1a_\mrm B=1, which is of angstrom order. Therefore, I will use 101010^{10} as the order of macroscopic lengths. The Rydberg unit of energy is Ry=1/2\mrm{Ry}=1/2, which is of electronvolt order. Therefore, I will use 10310^3 as the order of inverse room temperature.

One can adjust the units to get results for the cases of other hydrogen-like atoms: use Z2/4πε0=1Z^2/4\pi\veps_0=1 instead of 4πε0=14\pi\veps_0=1, where ZZ is the atomic number.

In this article, I also assume that the mass of the nucleus is infinite. If you want more accuracy, you can use mNme/(mN+me)=1m_\mrm Nm_\mrm e/\p{m_\mrm N+m_\mrm e}=1 instead of me=1m_\mrm e=1, where mNm_\mrm N is the mass of the nucleus.

Terminology about temperatures

I will mainly be working with the inverse temperature β1/kBT\beta\ceq1/k_\mrm BT, where TT is the temperature. However, I will still use “temperature” often to give some physical intuition. To avoid confusion in the context of using β\beta and in appearance of negative temperature, I would avoid using phrases like “high temperature” and “low temperature”. Instead, here are some terminologies that I am going to use:

  • “Cold (positive) temperature” means β+\beta\to+\infty.
  • “Hot positive temperature” means β0+\beta\to0^+.
  • “Cold negative temperature” means β0\beta\to0^-.
  • “Hot negative temperature” means β\beta\to-\infty.

The energy levels of a hydrogen atom are (ignoring fine structures etc.) En=1/2n2E_n=-1/2n^2, with each energy level labeled by nZ+n\in\bZ^+, and each energy level has gnn2g_n\ceq n^2 degeneracy (ignoring spin degeneracy, which merely contributes to an overall factor of the partition function). The partition function is Z ⁣(β)n=1gneβEn=n=1n2eβ/2n2,\fc Z\beta\ceq\sum_{n=1}^\infty g_n\e^{-\beta E_n} =\sum_{n=1}^\infty n^2\e^{\beta/2n^2},(1)\p{1} which diverges for any βC\beta\in\bC (of course, normally we can only have βR\beta\in\bR, but the point of saying that it diverges for any complex β\beta is that there is no way we can analytically continue the function to get a finite result). Does this mean that statistical mechanics breaks down for this system? Not necessarily. Actually, there are multiple ways we can tackle this divergence.

One should notice that, although this article concentrates on regularizing partition functions and that of the hydrogen atom in particular, all the methods are valid for more general divergent sums.

Here is a sentence that is quoted by many literatures on diverging series, so I want to quote it, too:

Divergente Rækker er i det Hele noget Fandenskap, og det er en Skam at man vover at grunde nogen Demonstrasjon derpaa.

—N. H. Abel

It translates to “Divergent series are in general deadly, and it is shameful that anyone dare to base any proof on them.”

The physical answer

A physicist always tell you that one should not be afraid of infinities. Instead, one should look at where the infinity comes out from the seemingly physical model, where there is something sneakily unphysical which ultimately leads to this unphysical divergence. In our case, the divergence comes from high energy levels. It is then a good time to question whether those high energy levels are physical.

There is a radius associated with each energy level in the sense of the Bohr model: rn=n2r_n=n^2. When rnL1010r_n\sim L\ceq10^{10} (which happens at nΛ105n\sim\Lmd\ceq10^5), the orbit is really microscopic now, and the interaction between the electron and the “box” that contains the whole experimental setup is now having significant effects. Or, if there is not a box at all, we can use the size of the universe instead, which is about rnL1036r_n\sim L\ceq10^{36} (Λ1018\Lmd\ceq10^{18}). Use the model of particle in a box for energy levels higher than n=Λn=\Lmd, and we have Z=n=1Λn2eβ/2n2+nx,ny,nz=1exp ⁣(β(nx2+ny2+nz2)π22L2),Z=\sum_{n=1}^\Lmd n^2\e^{\beta/2n^2} +\sum_{n_x,n_y,n_z=1}^\infty\fc\exp{-\beta\fr{\p{n_x^2+n_y^2+n_z^2}\pi^2}{2L^2}}, where LL is the side length of the box (assuming that the box is cubic). If LL is very large, we can approximate the second term as a spherically symmetric integral over the first octant to get L3(2πβ)3/2L^3\p{2\pi\beta}^{-3/2}.

The integral approximation

This is actually the result for Boltzmann ideal gas, so it should be familar, but I still write down the calculation here for completeness.

We can approximate nx,ny,nz=1exp ⁣(β(nx2+ny2+nz2)π22L2)I0d3nexp ⁣(β(nx2+ny2+nz2)π22L2),\sum_{n_x,n_y,n_z=1}^\infty\fc\exp{-\beta\fr{\p{n_x^2+n_y^2+n_z^2}\pi^2}{2L^2}} \approx I\ceq\int_0^\infty\d^3n\,\fc\exp{-\beta\fr{\p{n_x^2+n_y^2+n_z^2}\pi^2}{2L^2}}, where 0d3n\int_0^\infty\d^3n means 000dnxdnydnz\int_0^\infty\int_0^\infty\int_0^\infty\d n_x\,\d n_y\,\d n_z. We can then change the integral to spherical coordinates: I=0184πn2dnexp ⁣(βn2π22L2)=L34π2β3/2dnn2en2/2,I=\int_0^\infty\fr184\pi n^2\,\d n\,\fc\exp{-\beta\fr{n^2\pi^2}{2L^2}} =\fr{L^3}{4\pi^2\beta^{3/2}}\int_{-\infty}^\infty\d n\,n^2\e^{-n^2/2}, where the factor of 1/81/8 is because we only integrate in the first octant, and the second step utilizes the symmetry of the integrand and redefines the integrated variable. This integral is than a familiar Gaussian integral of order unity. The value of it is not important for later discussion because all the arguments that follow only uses orders of magnitude, but I tell you it is 2π\sqrt{2\pi}, which can be evaluated by integrating by parts once and utilizing the famous en2/2dn=2π\int_{-\infty}^{\infty}\e^{-n^2/2}\,\d n=\sqrt{2\pi}. The final result is I=L3(2πβ)3/2I=L^3\p{2\pi\beta}^{-3/2}.

Is this an overestimation or underestimation? It is actually an overestimation. Draw a picture of en2/2\e^{-n^2/2} to convince yourself of this. We do not need to estimate how large the error is, though, because we will see that we only need an upper bound to get the arguments we need.

For the first term, we need to consider how the magnitude of the summand changes with nn. The minimum value of the summand is at n=β/2n=\sqrt{\beta/2}. At room temperature, we have β103\beta\sim10^3, so β/2\sqrt{\beta/2} is well between 11 and Λ\Lmd. Therefore, the largest term is either n=1n=1 or n=Λn=\Lmd. The former is eβ/2\e^{\beta/2}, which is of order 1021710^{217}, while the latter is Λ2\Lmd^2, which is of order 103610^{36} for the case of the size of the universe. We may then be interested in the n=2n=2 term 4eβ/84\e^{\beta/8}, which is of order 105410^{54}. This is much larger than the n=Λn=\Lmd term but much smaller than the n=1n=1 term, so it is second largest term in the sum.

An upper bound of the summation is given by replacing every term except the largest term by the second largest term, which gives Z<eβ/210217+(Λ1)4eβ/81072+L3(2πβ)3/21048eβ/2.Z<\underbrace{\e^{\beta/2}}_{10^{217}} +\underbrace{\p{\Lmd-1}4\e^{\beta/8}}_{10^{72}}+\underbrace{L^3\p{2\pi\beta}^{-3/2}}_{10^{48}}\approx\e^{\beta/2}. Therefore, the n=1n=1 term dominates the entire partition function. This means that the hydrogen atom is extremely likely to be in the ground state (despite the seeming divergence of the partition function). This is intuitive. The probability of the system not being in the ground state is of order 105510^{-55} for the size of the universe and 1015810^{-158} for a typical macroscopic experiment.

More accurate considerations

The usage of the model of particle in a box for energy levels n>Λn>\Lmd gives good enough arguments and results, but one may want to question whether this is appropriate.

What happens if you actually put a hydrogen atom in a box (for simplicity, make the box spherically symmetric)? More accurately, consider the quantum mechanical problem in spherically symmetric potential VV such that Vr1V\sim-r^{-1} for small rr but grows fast and high enough at large rr so that the partition function for bound states is convergent. This is called a confined hydrogen atom. A book chapter The Confined Hydrogen Atom Revisited discusses this problem in detail and cited several papers that did the calculations about the energy levels.

Cutoff regularization

By analyzing the orders of magnitude, we see that we actually do not lose much if we just simply cut off the sum at n=Λn=\Lmd. This corresponds to a regularization method called the simple cutoff: it replaces the infinite sum by a finite partial sum. This can be generalized a little by considering a more general cutoff function χ\chi such that limx0+χ ⁣(x)=1\lim_{x\to0^+}\fc \chi x=1. Then, an infinite sum n=1f ⁣(n)\sum_{n=1}^\infty\fc fn can be written as n=1f ⁣(n)=limλ0+n=1f ⁣(n)χ ⁣(λn).\sum_{n=1}^\infty\fc fn=\lim_{\lmd\to0^+}\sum_{n=1}^\infty\fc fn\fc\chi{\lmd n}. The simple cutoff is then the case where χ ⁣(x)θ ⁣(1x)\fc \chi x\ceq\fc\tht{1-x} and λ1/Λ\lmd\ceq1/\Lmd, where θ\tht is the Heaviside step function. For converging series, this gives the same result as the original sum thanks to the dominated convergence theorem.

For diverging series

For diverging series, this may give a finite result. For example, for f ⁣(n)(1)nnk\fc fn\ceq\p{-1}^nn^k, this method gives η ⁣(k)-\fc\eta{-k} for any complex kk and any smooth enough χ\chi, where η\eta is the Dirichlet eta function. Here is a check for the special case χ ⁣(x)ex\fc\chi x\ceq\e^{-x} (equivalent to the Abel summation). By definition of the polylogarithm, we have n=1(1)nnkeλn=Lik ⁣(eλ).\sum_{n=1}^\infty\p{-1}^nn^k\e^{-\lmd n}=\fc{\mrm{Li}_{-k}}{-\e^{-\lmd}}. Now, substitute λ=0\lmd=0, and utilizing the identity Lis ⁣(1)=η ⁣(s)\fc{\mrm{Li}_s}{-1}=-\fc\eta s, we have the result η ⁣(k)-\fc\eta{-k}.

You may wonder what is the case for f ⁣(n)nk\fc fn\ceq n^k, which is also a diverging series, and it looks much like the case above. However, the limit at λ0+\lmd\to0^+ simply does not exist when Rek1\Re k\ge-1 (i.e., when the series diverges). This is because we have Lis ⁣(1)=ζ ⁣(s)\fc{\mrm{Li}_s}1=\fc\zeta s only for Res>1\Re s>1, where ζ\zeta is the Riemann zeta function, but it is undefined for other values of ss. If you analytically continue the result, you will get the famous Rieman zeta function.

However, although this series may converge for any positive λ\lmd, the limit as λ0+\lmd\to0^+ may not exist. If it diverges because f ⁣(n)\fc fn grows too fast (or decays too slowly) as nn\to\infty, then we should expect that the sum also tends to infinity as λ0+\lmd\to0^+. Assume that we can characterize this divergence by a Laurent series: n=1f ⁣(n)χ ⁣(λn)=k=γkλk.\sum_{n=1}^\infty\fc fn\fc\chi{\lmd n} =\sum_{k=-\infty}^\infty\gma_k\lmd^k.(2)\p{2} If the λ0+\lmd\to0^+ limit converge, we would expect γk<0\gma_{k<0} to be zero, and then the result is simply γ0\gma_0. Therefore, we may also want only γ0\gma_0 when the limit does not exist. To pick out γ0\gma_0, utilize the residue theorem: n=1f ⁣(n)=12πidλλn=1f ⁣(n)χ ⁣(λn),\sum_{n=1}^\infty\fc fn=\fr1{2\pi\i}\oint\fr{\d\lmd}\lmd \sum_{n=1}^\infty\fc fn\fc\chi{\lmd n},(3)\p{3} where the domain of λ\lmd is now analytically continued from R+\bR^+ to a deleted neighborhood of 00. Equation 3 is then a generalized version of Equation 2.

Notice that I have been super slippery in math in the discussion. For example, the Laurent series may not exist at all, and the analytic continuation may not be possible at all; even if they exist, the λ0+\lmd\to0^+ limit may also be different from γ0\gma_0. However, I may claim that we should be able to select smooth enough χ\chi for all of these to work, and the results will be independent of the choice of χ\chi as long as Equation 3 works in this form.

Particularly, one can rigorously prove that for f ⁣(n)nk\fc fn\ceq n^k, the sum obtained by this precedure is ζ ⁣(k)\fc\zeta{-k}, where ζ\zeta is the Riemann zeta function, as long as xkχ ⁣(x)x^k\fc\chi x has bounded (k+2)\p{k+2}th derivative and the sum converges. This is proven in an interesting blog article.

Alternative forms of cutoff regularization

In some cases, one may discover that nf ⁣(n)χ ⁣(λn)\sum_n\fc fn\fc\chi{\lmd n} is not analytic when λ0+\lmd\to0^+ so that the Laurent series expansion is not possible. An example is EnlnlnnE_n\ceq\ln\ln n (for n2n\ge2) with no degeneracies (this system also has a diverging partition function for any complex β\beta). In this case, if you try to use the cutoff function χ ⁣(x)ex\fc\chi x\ceq\e^{-x}, the sum goes like (lnλ)β/λ\p{-\ln\lmd}^{-\beta}/\lmd instead of analytically when λ0+\lmd\to0^+. Proving this is simple. We have Zλ=n=2eλn(lnn)β2eλn(lnn)βdn=1λ(lnλ)β2λexdx(1lnx/lnλ)β,Z_\lmd=\sum_{n=2}^\infty\e^{-\lmd n}\p{\ln n}^{-\beta} \approx\int_2^\infty\e^{-\lmd n}\p{\ln n}^{-\beta}\d n =\fr1{\lmd\p{-\ln\lmd}^\beta}\int_{2\lmd}^\infty \fr{\e^{-x}\,\d x}{\p{1-\ln x/\ln\lmd}^\beta}, where the last step uses the substitution xλnx\ceq\lmd n. Using the binomial theorem, we have Zλ1λ(lnλ)β2λdxexk=0(βk)(lnxlnλ)k,Z_\lmd\approx\fr1{\lmd\p{-\ln\lmd}^\beta}\int_{2\lmd}^\infty\d x\,\e^{-x} \sum_{k=0}^\infty\binom{-\beta}k\p{\fr{\ln x}{-\ln\lmd}}^k, where (βk)\binom{-\beta}k is the binomial coefficient. Note that Γ(k) ⁣(z)=0xk1(lnx)kezxdx\fc{\Gma^{\p k}}z=\int_0^\infty x^{k-1}\p{\ln x}^k\e^{-zx}\d x, where Γ(k)\Gma^{\p k} is the kkth derivative to the Euler Gamma function, so the integral for xx gives a factor Γ(k) ⁣(1)\fc{\Gma^{\p k}}1 in the limit of λ0+\lmd\to0^+. Therefore, Zλ1λ(lnλ)β,Z_\lmd\approx\fr1{\lmd\p{-\ln\lmd}^\beta}, where only the k=0k=0 term in the sum is retained for the leading contribution as λ0+\lmd\to0^+.

However, for any kZ+k\in\bZ^+, one can always choose functions h,χh,\chi so that the sum nf ⁣(n)χ ⁣(λh ⁣(n))\sum_n\fc fn\fc\chi{\lmd\fc hn} goes like λk\lmd^{-k} as λ0+\lmd\to0^+. For example, for χ ⁣(x)ex\fc\chi x\ceq\e^{-x} (equivalent to the Abelian mean or the heat-kernel regularization), we have Zλn0eλh ⁣(n)f ⁣(n)dn=λf ⁣(n0)exf ⁣(h1 ⁣(xλ))dxλh ⁣(h1 ⁣(xλ)).Z_\lmd\approx\int_{n_0}^\infty\e^{-\lmd\fc hn}\fc fn\d n =\int_{\lmd\fc f{n_0}}^\infty\e^{-x}\fc f{\fc{h^{-1}}{\fr x\lmd}}\fr{\d x}{\lmd\fc{h'}{\fc{h^{-1}}{\fr x\lmd}}}. We can choose h ⁣(n)(f ⁣(n)dn)1/k\fc hn\ceq\p{\int\fc fn\d n}^{1/k} so that f ⁣(h1 ⁣(xλ))=k(xλ)k1h ⁣(h1 ⁣(xλ)).\fc f{\fc{h^{-1}}{\fr x\lmd}}=k\p{\fr x\lmd}^{k-1}\fc{h'}{\fc{h^{-1}}{\fr x\lmd}}. Therefore, as λ0+\lmd\to0^+, we have Zλ1λλf ⁣(n0)exk(xλ)k1dxk!λk.Z_\lmd\approx\fr1\lmd\int_{\lmd\fc f{n_0}}^\infty\e^{-x}k\p{\fr x\lmd}^{k-1}\,\d x\approx\fr{k!}{\lmd^k}. However, this does not guarantee that the Laurent series expansion exists. This is a good trial, though. My math capacity does not allow me to confirm whether this is the case for the example of EnlnlnnE_n\ceq\ln\ln n.

Regularizing the hydrogen atom

After saying so much about cutoff regularization in general, what does it say about the partition function of a hydrogen atom? Try multiplying the cutoff function χ ⁣(λn)\fc\chi{\lmd n} to the summand in Equation 1: Zλn=1n2eβ/2n2χ ⁣(λn)=k=0(β/2)kk!n=1n22kχ ⁣(λn)k=0(β/2)kk!ζ ⁣(2k2),Z_\lmd\ceq\sum_{n=1}^\infty n^2\e^{\beta/2n^2}\fc\chi{\lmd n} =\sum_{k=0}^\infty\fr{\p{\beta/2}^k}{k!}\sum_{n=1}^\infty n^{2-2k}\fc\chi{\lmd n} \to\sum_{k=0}^\infty\fr{\p{\beta/2}^k}{k!}\fc\zeta{2k-2},(4)\p{4} where the last step utilizes the result for f ⁣(n)nk\fc fn\ceq n^k, with which we get rid of the dependence on λ\lmd. The last expression is then identified as ZZ.

Now that we get the expression of ZZ, we can get some useful things. However, this time we cannot simply use the summand divided by ZZ to get the probability of each energy level because that will break the normalization of the probability distribution. What we can do, however, is to find the expectation value of the energy using E=dlnZ/dβ\a E=-\d\ln Z/\d\beta. On the other hand, we have Ep1E1+(1p1)E=p1/2\a E\le p_1E_1+\p{1-p_1}E_\infty=-p_1/2, so the probability 1p11-p_1 that the system is not in the ground state is bounded above by 2E+12\a E+1.

The first check to do is to verify that this result is consistent with the known behavior of the system at cold zero temperature, where the system is almost certainly in the ground state; in other words, limβ+E=1/2\lim_{\beta\to+\infty}\a E=-1/2. To get ZZ for large β\beta, we notice that ζ ⁣(+)=1\fc\zeta{+\infty}=1, so Zeβ/2Z\approx\e^{\beta/2}, and this leads to E1/2\a E\approx-1/2 as expected.

Now, we may try to estimate E\a E for finite but large β\beta (e.g., β=103\beta=10^3) and thus give an upper bound for 1p11-p_1. We can study the asymptotic behavior of E\a E for cold positive temperature. It turns out that 1p13e3β/81-p_1\approx3\e^{-3\beta/8}, which is 1016310^{-163} for β=103\beta=10^3. As we can see, without any physical arguments but only with regularization, we get a result that seems sensible and well between the results in the last section for a hydrogen atom confined in a box with a typical macroscopic size or the size of the universe.

Derivation of the asymptotic behavior at cold positive temperature

We have Z=k=0(β/2)kk!ζ ⁣(2k2),dZdβ=12k=0(β/2)kk!ζ ⁣(2k).Z=\sum_{k=0}^\infty\fr{\p{\beta/2}^k}{k!}\fc\zeta{2k-2},\quad \fr{\d Z}{\d\beta}=\fr12\sum_{k=0}^\infty\fr{\p{\beta/2}^k}{k!}\fc\zeta{2k}. Therefore, Z2dZdβ=k=0(β/2)kk!(ζ ⁣(2k2)ζ ⁣(2k)).Z-2\fr{\d Z}{\d\beta}=\sum_{k=0}^\infty\fr{\p{\beta/2}^k}{k!} \p{\fc\zeta{2k-2}-\fc\zeta{2k}}. We can try to find the asymptotic behavior of the coefficient of each term. We have ζ ⁣(2k2)ζ ⁣(2k)=n=1(1n2k21n2k)=n=1n21n2k=322k+O ⁣(132k).\fc\zeta{2k-2}-\fc\zeta{2k}=\sum_{n=1}^\infty\p{\fr1{n^{2k-2}}-\fr1{n^{2k}}} =\sum_{n=1}^\infty\fr{n^2-1}{n^{2k}} =\fr{3}{2^{2k}}+\O{\fr1{3^{2k}}}. We also have ζ ⁣(2k2)=1+O ⁣(22k)\fc\zeta{2k-2}=1+\O{2^{-2k}}, of course. Therefore, 1p12E+1=Z2dZ/dβZ=k(β/2)kk!(322k+O ⁣(132k))k(β/2)kk!(1+O ⁣(122k)).1-p_1\le2\a E+1=\fr{Z-2\d Z/\d\beta}Z =\fr{\sum_k\fr{\p{\beta/2}^k}{k!}\p{\fr3{2^{2k}}+\O{\fr1{3^{2k}}}}}{\sum_k\fr{\p{\beta/2}^k}{k!}\p{1+\O{\fr1{2^{2k}}}}}. These power series are then simply exponential functions. Therefore, 1p13eβ/8+O ⁣(eβ/18)eβ/2+O ⁣(eβ/8)=3e3β/8+O ⁣(e4β/9).1-p_1\le\fr{3\e^{\beta/8}+\O{\e^{\beta/18}}}{\e^{\beta/2}+\O{\e^{\beta/8}}} =3\e^{-3\beta/8}+\O{\e^{-4\beta/9}}.(5)\p{5}

Although the asymptotic behavior at cold temperature (β+\beta\to+\infty) looks good, its behavior is very wrong at some regimes. At some temperature, the monoticity of E\a E reverts, and then it gets even lower than the ground state energy 1/2-1/2 and heads all the way to -\infty at some finite temperature.This is clearly unphysical. This suggests that it is wrong to use the regularized result.


Here is a plot that shows how E\a E starts to decrease with temperature at some point and becomes even lower than the ground state energy:

Plot of  vs. 

Here is a plot that shows how E\a E goes to infinity at different temperatures:

Plot of  vs. 

Here are also plots for ZZ and dZ/dβ\d Z/\d\beta, if you are curious:

Plot of  and  vs. 

The two vertical asymptotes of E\a E corresponds to the two zeros of ZZ, which are β=0\beta=0 and β=1.0721\beta=1.0721. It also has a zero, correponding to the zero of dZ/dβ\d Z/\d\beta at β=0.5530\beta=0.5530. The point where E=1/2\a E=-1/2 is β=11.2486\beta=11.2486, and the point where E\a E has a local maximum is β=13.8021\beta=13.8021.

Another aspect where we can see that this result is wrong is that, if we look at the hot negative temperature limit β\beta\to-\infty, although we have E0=supnEn\a E\to0=\sup_nE_n as expected, it is approaching from the wrong side. In fact, because Z>0Z>0 while dZ/dβ<0\d Z/\d\beta<0 for β<0\beta<0, we have E>0\a E>0 for β<0\beta<0, exceeding the supremum of the energy levels, which is unphysical.

Derivation of the hot negative temperature limit

Here is a non-rigorous derivation. We can rewrite the regularized ZZ in a similar form as Z=β4+n=1n2(eβ/2n21β2n2)=limN((14+N2)β16N(1+N)(1+2N)+n=1Nn2eβ/2n2).\begin{align*} Z&=-\fr\beta4+\sum_{n=1}^\infty n^2\p{\e^{\beta/2n^2}-1-\fr\beta{2n^2}}\\ &=\lim_{N\to\infty}\p{-\p{\fr14+\fr N2}\beta-\fr16N\p{1+N}\p{1+2N}+\sum_{n=1}^Nn^2\e^{\beta/2n^2}}. \end{align*} For finite NN, it has a straight line asymptote as β\beta\to-\infty. The envelope of this family of straight lines (parametrized by NN) is Z=(16β)3/2/363Z=\p{1-6\beta}^{3/2}/36\sqrt3, which means that Z(β)3/2Z\sim\p{-\beta}^{3/2} as β\beta\to-\infty, where “\sim” means that the ratio of the two sides approaches a positive constant. Similarly, we have dZ/dβ(β)1/2\d Z/\d\beta\sim-\p{-\beta}^{1/2}. Therefore, Eβ1\a E\sim-\beta^{-1} as β\beta\to-\infty.

Another regularization special to the hydrogen atom

Here is a special regularization method for the hydrogen atom which is not applicable to general systems. Consider the second derivative d2Z/dβ2\d^2Z/\d\beta^2 by differentiating the summand twice w.r.t. β\beta in Equation 1, and then take twice antiderivative w.r.t. β\beta. This gives Z=A+Bβ+k=0(β/2)kk!ζ ⁣(2k2)=A+Bβ+n=1n2(eβ/2n21β2n2),Z=A+B\beta+\sum_{k=0}^\infty\fr{\p{\beta/2}^k}{k!}\fc\zeta{2k-2} =A+B\beta+\sum_{n=1}^\infty n^2\p{\e^{\beta/2n^2}-1-\fr{\beta}{2n^2}}, where A,BA,B are integration constants. The result from the cutoff regularization and the zeta function regularization is simply A=0A=0, B=1/4B=-1/4. What is interesting about this is that it already determines the asymptotic behavior of 1p11-p_1 at cold temperature, which is 1p13e3β/81-p_1\approx3\e^{-3\beta/8} (see Equation 5), no matter what A,BA,B are.

Zeta function regularization

For a series nf ⁣(n)\sum_n\fc fn, if it diverges, we can instead consider nf ⁣(n)s\sum_n\fc fn^{-s} for some ss whose real part is big enough for the series to converge. Then, we can try to analytically continue to s=1s=-1 to get a finite result for the original series. This is called the zeta function regularization.

When zeta function regularization fails

For the zeta function regularization to work, the asymptotic behavior of f ⁣(n)\fc fn needs to be a non-trivial power law as n+n\to+\infty. Otherwise, the sum may not converge for any ss. For example, consider En=lnlnnE_n=\ln\ln n (with no degeneracies). The partition function with zeta function regularization is Zsn=2(lnn)βs.Z_s\ceq\sum_{n=2}^\infty\p{\ln n}^{\beta s}. This series is divergent for any complex ss.

A famous example is f ⁣(n)n\fc fn\ceq n, which gives nn=ζ ⁣(1)=1/12\sum_nn=\fc\zeta{-1}=-1/12. Generally, for f ⁣(n)nk\fc fn\ceq n^k, we have nnk=ζ ⁣(k)\sum_nn^k=\fc\zeta{-k}. This is the same as the result for the simple cutoff regularization. This raises the question of whether the results obtained from those two methods are necessary the same whenever they both exist. I do not have a rigorous proof, but a strong argument is that both of them are the result of some analytic continuation, so they should be the same by the uniqueness of analytic continuation.

We can check this with the hydrogen atom. We have, for s>1/2s>1/2 and β\beta real, Zsn=1n2sesβ/2n2=k=0(sβ/2)kk!n=1n2s2k=k=0(sβ/2)kk!ζ ⁣(2s+2k).Z_s\ceq\sum_{n=1}^\infty n^{-2s}\e^{-s\beta/2n^2} =\sum_{k=0}^\infty\fr{\p{-s\beta/2}^k}{k!}\sum_{n=1}^\infty n^{-2s-2k} =\sum_{k=0}^\infty\fr{\p{-s\beta/2}^k}{k!}\fc\zeta{2s+2k}. Analytically continue this result to s=1s=-1, and then this gives the same result as Equation 4. The rest will be the same as the last section.

Can we trust this result?

However, can we trust this result, though? Everything is becoming fishy. Probabilities are no longer well-defined because how we normally derive them using ZZ is causing a divergent sum of probabilities and thus invalid. Yet somehow we are trying to estimate the bound of the probability of the system not being in the ground state and getting an expected result. You must have been feeling uncomfortable about this.

The first thing to ask is what we mean by “the expectation value” when the probability distribution is not even well-defined. If it means nothing physical, can we still trust its expression? The simple answer is no.

As we already see, although the result at cold temperature is sensible, the result at some regimes is clearly unphysical. We can also see similar problems with other systems. Consider the system that has energy levels En=lnnE_n=\ln n (with no degeneracies). We can easily get Z=nnβ=ζ ⁣(β)Z=\sum_nn^{-\beta}=\fc\zeta\beta, and thus there is a absi at β=1\beta=1. For β>1\beta>1, ZZ converges, and everything looks good. For β<1\beta<1, the system is so hot that ZZ diverges. Previous arguments suggest that, in this region, the regularized ZZ is still ζ ⁣(β)\fc\zeta\beta. However, we then have E<0\a E<0 in this region, which is lower than the ground state energy. This clearly should not be trusted.

In another aspect, we should note that since the estimation for 1p11-p_1 does not depend on the size of the box confining the hydrogen atom, its rough agreement with the result in the last section should be considered a coincidence.

Another thing to note is that the result of the regularizations depend on whether we “flatten” the energy levels. We can “flatten” all the energy levels: pretend no degeneracies exist. For example, suppose a system with gnng_n\ceq n and EnnE_n\ceq n. However, we can rewrite the same system as En1,2,2,3,3,3,E_n\ceq1,2,2,3,3,3,\ldots (or equivalently En2n+1/2E_n\ceq\floor{\sqrt{2n}+1/2})1, with no degeneracies. This “re-grouping” of the energy levels can affect the result of regularizations and whether a zeta function regularization exists. For an immediate example, if we flatten the energy levels of the hydrogen atom, the zeta function regularization does not exist. Another simple example is that, for a system with En=constE_n=\mrm{const}, we can essentially re-group the all-degenerate states to have any positive integer sequence gng_n to get very arbitrary results for the partition function.

Abscissa of convergence

Forget about the hydrogen atom, and let us consider a general system with (ever-increasing) energy levels EnE_n and degeneracies gng_n. For a given system, there is an abscissa of convergence βc\beta_\mrm c, below (hotter than) which the partition function diverges. In other words, ZngneβEnZ\ceq\sum_ng_n\e^{-\beta E_n} converges for Reβ>βc\Re\beta>\beta_\mrm c and diverges for Reβ<βc\Re\beta<\beta_\mrm c. For most physical systems, we have βc=0\beta_\mrm c=0, meaning that it can have any positive temperature, which sounds sensible. The hydrogen atom has βc=+\beta_\mrm c=+\infty, and a two-level system has βc=\beta_\mrm c=-\infty. A system with En=lnnE_n=\ln n and no degeneracy has βc=1\beta_\mrm c=1.

The term “abscissa of convergence” is borrowed from the study of general Dirichlet series. The form of ZZ is indeed very much like a general Dirichlet series, but a general Dirichlet series requires E=+E_\infty=+\infty, which is not true for the hydrogen atom. However, the existence of an abscissa of convergence is still true for the more general case.

What does it mean physically to have an abscissa of convergence βc\beta_\mrm c? First, if βc=\beta_\mrm c=-\infty, then the system is well behaved at any temperature, which is good and does not need further care.

If βc<\v{\beta_\mrm c}<\infty, normally one should say the system cannot reach a certain temperature: the system can never be in equilibrium with a heat bath hotter than βc\beta_\mrm c. Thermodynamically, one can say that the system needs to absorb an infinite amount of heat to reach this temperature. One can see this easily by considering any sensible system, which has βc=0\beta_\mrm c=0: for β\beta to go below zero means to make the temperature hotter than infinity, which of course needs an infinite amount of heat intuitively. One may want to see whether it is possible to regularize ZZ to get a finite result for Reβ<βc\Re\beta<\beta_\mrm c. A valid claim to make is that, if ZZ can be analytically continued to the half real axis to the left of βc\beta_\mrm c, then any sensible regularization of ZZ there will give the same result as the analytic continuation. Actually, the analytic continuation is exactly the zeta function regularization if there is no degeneracy (or regarding degenerate states as different energy levels). However, it is possible that the analytic continuation does not exist. There may be a branch cut or a natural boundary. For example, if EnlnpnE_n\ceq\ln p_n with no degeneracy, where pnp_n is the nnth prime number, then ZZ is the prime zeta function, which has a natural boundary at Reβ=0\Re\beta=0. Even if such a regularization exists, it should be questioned whether it is physical.

If βc=+\beta_\mrm c=+\infty, then the system is not well behaved at any temperature. This is the case for the hydrogen atom. Physically, this means that the system cannot be in equilibrium with a heat bath at any temperature. The problem with regularization is the same as the case with βc<\v{\beta_\mrm c}<\infty.

In a previous article about statistical ensembles, when I defined the partition function, I briefly mentioned that it is only defined for those intensive variables (β\beta in the context of this article) such that the partition function converges. I did not talk about what to do with the partition function when it diverges, but what that article implied is that it is simply undefined and that no physical meaning should be assigned to it in principle. The existence of an abscissa of convergence tells us that there is a “hottest possible temperature” for any given system. The hydrogen atom is symply the case where the hottest possible temperature coincides with the coldest possible temperature (which is the absolute zero). For most sensible systems, the hottest possible temperature is just the positive hot limit. For systems such as EnlnnE_n\ceq\ln n, the hottest possible temperature is a finite positive temperature, which is at 3.16×105K3.16\times10^5\,\mrm K, resulting from βc=1\beta_\mrm c=1. This can be conterintuitive at first, but one should realize that it is not essentially different from the more common case of βc=0\beta_\mrm c=0.

  1. This is A002024 on OEIS. Coincidentally, the OEIS number of this sequence is the same as the year in which I am writing this article.↩︎