Thermal corrections for positronium

The development of quantum mechanics (QM) and quantum electrodynamics (QED) methods for a detailed description of the processes occurring in atomic systems and the assessment of the corresponding relativistic QED corrections to bound energies [1, 2, 3, 4, 5, 6, 7] play a key role in modern physics. Theoretical calculations in conjunction with the growing experimental accuracy put a way to test our understanding of physics up to the level of $4.2\times 10^{-15}$ [8, 9] or even better $2.1\times 10^{-18}$ [10, 11]. The experiments with such extra-ordinary precision required theoretical calculations of various QED effects at the $\alpha^{6}m^{2}/M$ and $\alpha^{7}m$ levels, see [12] and references therein, where $\alpha$ is the fine structure constant, $m$ and $M$ are the electron and nuclear masses, respectively.

The excellent agreement of theoretical and experimental results up to this level is in some cases violated [13, 14, 15, 16], which directs scientific attention to the verification of fundamental interactions, which, in turn, pursues the goal of eliminating such contradictions [17, 18, 19]. In such cases, a critical analysis of theoretical calculations of energies, fine structure, isotope shifts, etc is required for states in various atomic systems, revealing unresolved discrepancies with experiment, to determine fundamental constants or interactions [20], to set constraints to dark matter [21].

Many theoretical and experimental efforts have gone into the study of the positronium (Ps) atom. Such an atomic system (the bound state of an electron and a positron) is the most attractive for theorists since it is the lightest and the simplest hydrogenlike atom. Representing the purely leptonic system the Ps atom does not depend on hadronic effects and its properties are described by the bound-state QED theory, see, for example, [4]. Thus, a detailed comparison of theory and experimental measurements of the lifetimes [22], transition frequencies [23] or hyperfine energy splitting [24] can be used as a sensitive tool for testing lepton interactions. Probing the positronium atom ipso facto brings periodically exposed disagreements of theory with experimental data [25, 23]. The ’purification’ [25] combined with increasing the accuracy of the experiment stimulates theoretical development and computing of corresponding magnitudes, creating a basis for searching and testing various physical effects and hypotheses.

Special attention is paid to thermal-induced effects aimed at the accurate calculation of the binding energies in the atom. For an atomic system placed in a heat bath (blackbody radiation) the known thermal effects consist in presence of a Stark shift of energy levels and an induced line broadening. In nowadays, evaluation of the Stark shift induced by the blackbody radiation (BBR) is of fundamental importance in atomic clocks, see [26, 27, 28, 29, 30], establishing the most significant restraint on the accuracy of operating transitions [31, 32, 33].

Calculations of the BBB-induced Stark shift based on the QM approach [34] are generally extended to a many-electron atomic system, the theory of which includes evaluation of the static polarizability and dynamical corrections to it [35]. Theoretical predictions made within the framework of the multipolar decomposition method [35] was later verified experimentally for ytterbium clocks [33]. However, the theory of thermal action on atomic systems can be defined within the framework of the QED approach allowing the rigorous consideration of various effects (for example, relativistic, radiative, or finite lifetimes). By this means, the QED derivation of the Stark shift induced by the BBR field was recently performed in [36], where, in particular, it was found that the real and imaginary parts of the thermal one-loop self-energy correction represent BBR-Stark shift and line broadening, respectively.

More recently, the QED approach has been used to study the interaction induced by the BBR for two charges [37]. The lowest-order thermal correction obtained in the nonrelativistic limit turned out to be cubic in temperature, while the Stark shift is proportional to the fourth power. Such a difference shows a fundamental difference in the effects, and the results for the shift of atomic levels are more significant, see [37] for one-electron atomic systems and [38] for helium. In particular, the thermal interaction of a bound electron with a nucleus can lead to an energy shift exceeding the corresponding Stark shift. The hypothesis put forward in [37, 38] is mainly based on the thermal quantum electrodynamics (TQED) [39, 40, 41], and the existence of a corresponding correction can be indirectly proven by a theoretical prediction made in [41], as well as the experimentally observed thermal effect called as $T^{2.7}$ in [42]. Finally, the relativistic thermal corrections to the thermal photon exchange between two charges were obtained in [43]. Revealing these effects is still a matter of experimental examination.

In view of the close attention to the verification of fundamental physical interactions in experiments with simple atomic systems, the derivation of thermal effects leading to the correction of the lowest order, fine and hyperfine level splitting is of considerable interest for positronium. This problem can be solved using the formalism presented in [37, 43]. In this paper, the thermal corrections arising from the scalar and transversal parts of the thermal photon propagator are evaluated for the Ps atom. All the derivations are performed within the framework of rigorous quantum electrodynamics at finite temperatures. In addition, thermal corrections to the probabilities of two- and three-photon annihilation of the positronium atom are briefly discussed in an attempt to provide an exhaustive description of the thermal impact in the lowest order.

Starting with the description of the interaction of two charges, one can use the relation from textbooks (see, for example, [2]) connecting the nuclear current, $j^{\nu}(x^{\prime})$, with the field, $A_{\mu}(x)$, it creates:

$\displaystyle A_{\mu}(x)=\int d^{4}xD_{\mu\nu}(x,x^{\prime})j^{\nu}(x^{\prime}),$

(1)

where $x=(t,\vec{r})$ represents the four-dimensional coordinate vector ($t$ represents time and $\vec{r}$ denotes a space vector), $D_{\mu\nu}(x,x^{\prime})$ is the Green’s function of the photon, and $\mu$, $\nu$ are the indices running the values $0,1,2,3$. Then, the zero component of $A_{\mu}(x)$ corresponds to the Coulomb interaction, and the components $1,2,3$ are the transversal part, which gives the interaction of retardation and advance. According to [39, 40, 41], the photon Green’s function (photon propagator) is represented by the sum of two contributions, which are the result of expectation value with the states of zero and heated vacuum, $D_{\mu\nu}(x,x^{\prime})=D_{\mu\nu}^{0}(x,x^{\prime})+D^{\beta}_{\mu\nu}(x,x^{\prime})$, respectively.

Thermal interaction can be introduced by analogy, see [37], when the ’ordinary’ photon Green’s function is replaced by the thermal one, $D_{\mu\nu}^{\beta}(x_{1},x_{2})$, [39, 40, 41]. In [37] it was found that the thermal part of photon propagator $D_{\mu\nu}^{\beta}(x,x^{\prime})$ admits a different (equivalent) form:

$\displaystyle D_{\mu\nu}^{\beta}(x,x^{\prime})=-4\pi g_{\mu\nu}\int\limits_{C_{1}}\frac{d^{4}k}{(2\pi)^{4}}\frac{e^{ik(x-x^{\prime})}}{k^{2}}n_{\beta}(|\vec{k}|),$

(2)

where $g_{\mu\,\nu}$ is the metric tensor, $k^{2}=k_{0}^{2}-\vec{k}^{2}$ and $n_{\beta}$ is the Planck’s distribution function. The contour of integration in $k_{0}$-plane for Eq. (2) is given in Fig. 1.

Thermal photon propagator in the form Eq. (2) has the advantage of allowing the simple introduction of gauges, see [37]. In the Coulomb gauge the function $D_{\mu\nu}^{\beta}$ recasts into

$\displaystyle D_{00}^{\beta}(x,x^{\prime})$

$\displaystyle=$

$\displaystyle 4\pi i\int\limits_{C_{1}}\frac{d^{4}k}{(2\pi)^{4}}\frac{e^{ik(x-x^{\prime})}}{\vec{k}^{2}}n_{\beta}(\omega)$

(3)

for the Coulomb part, and the transversal part is

$\displaystyle D_{ij}^{\beta}(x,x^{\prime})$

$\displaystyle=$

$\displaystyle 4\pi i\int\limits_{C_{1}}\frac{d^{4}k}{(2\pi)^{4}}\frac{e^{ik(x-x^{\prime})}}{k^{2}}n_{\beta}(\omega)\left(\delta_{ij}-\frac{k_{i}k_{j}}{\vec{k}^{2}}\right).\,\,\,$

(4)

Concentrating first on the thermal Coulomb interaction, i.e. the zero component of the thermal photon propagator, the substitution of Eq. (3) into the expression (1) produces the thermal potential for the point-like nucleus in the static limit. The appropriate analytical calculations [37] of the integral yield

$\displaystyle V^{\beta}(r)=-\frac{4e^{2}}{\pi}\left(-\frac{\gamma}{\beta}+\frac{i}{2r}\ln\left[\frac{\Gamma\left(1+\frac{ir}{\beta}\right)}{\Gamma\left(1-\frac{ir}{\beta}\right)}\right]\right),$

(5)

where $\beta\equiv 1/(k_{B}T)$ ($k_{B}$ is the Boltzmann constant and $T$ is the temperature in kelvin), $r$ is the modulus (length) of the corresponding radius vector for the interparticle distance, $\Gamma$ is the gamma function and $\gamma$ is the Euler-Mascheroni constant, $\gamma\simeq 0.577216$.

Expanding the thermal potential Eq. (5) in terms of $\beta\rightarrow\infty$ (low-temperature regime) or in the small distance limit (the nonrelativistic limit for bound particles), $r\rightarrow 0$, the thermal corrections of lowest order are

$\displaystyle V^{\beta}(r)\approx\frac{4e^{2}\zeta(3)}{3\pi\beta^{3}}r^{2}-\frac{4e^{2}\zeta(5)}{5\pi\beta^{5}}r^{4}+\dots,$

(6)

where $\zeta(s)$ is the Riemann zeta function. The parameter estimation of expression (6) is as follows. In relativistic units $e^{2}\sim\alpha$, $r\sim 1/(m\alpha Z)$ and $1/\beta=k_{B}(r.u.)T$, where $\alpha$ is the fine structure constant, $m$ is the particle mass and $T$ is the temperature in kelvin. Then, combining these parametrizations, one can find for the first term $\frac{1}{\alpha m^{2}Z^{2}}(k_{B}T)^{3}$ r.u.$=\frac{\alpha^{3}}{Z^{2}}(k_{B}T)^{3}$ in atomic units, and $\frac{1}{\alpha^{3}m^{4}Z^{4}}(k_{B}T)^{5}$ r.u.$=\frac{\alpha^{5}}{Z^{4}}(k_{B}T)^{5}$ in atomic units for the second contribution. Here we used the ratio $k_{B}(r.u.)=m\alpha^{2}k_{B}(a.u.)$ and conversion pre-factor $m\alpha^{2}$ between relativistic and atomic units.

The potential (5) was used to find thermal corrections to the energy of a bound electron in the hydrogen atom [37] and helium [38], where the thermal correction turned out to be of the order of the experimental accuracy [44]. As it should be, in the lowest order the heat bath environment removes the orbital momentum degeneracy in the hydrogen atom.

Relativistic corrections arising through the transversal part of Eq. (4) (thermal Breit interaction) to interaction (5) were recently discussed in [43]. To find them (corrections to the fine and hyperfine structure), one should turn to the Pauli approximation or determine the relativistic corrections proportional to $1/c^{2}$, where $c$ is the speed of light, see [2, 7, 6, 4]. The particle interaction operator (in the case of zero vacuum) in momentum representation and Coulomb gauge [4] reads

	$\displaystyle U(\vec{p}_{1},\vec{p}_{2},\vec{k})=4\pi e^{2}\left[\frac{1}{\vec{k}^{2}}-\frac{1}{8m_{1}^{2}c^{2}}-\frac{1}{8m_{2}^{2}c^{2}}\qquad\right.$		(7)
	$\displaystyle+\left.\frac{(\vec{k}\vec{p}_{1})(\vec{k}\vec{p}_{2})}{m_{1}m_{2}c^{2}\vec{k}^{4}}-\frac{\vec{p}_{1}\vec{p}_{2}}{m_{1}m_{2}c^{2}\vec{k}^{2}}+\frac{i\vec{\sigma}_{1}[\vec{k}\times\vec{p}_{1}]}{4m_{1}^{2}c^{2}\vec{k}^{2}}-\frac{i\vec{\sigma}_{2}[\vec{k}\times\vec{p}_{2}]}{4m_{2}^{2}c^{2}\vec{k}^{2}}\right.$
	$\displaystyle\left.-\frac{i\vec{\sigma}_{1}[\vec{k}\times\vec{p}_{2}]}{2m_{1}m_{2}c^{2}\vec{k}^{2}}+\frac{i\vec{\sigma}_{2}[\vec{k}\times\vec{p}_{1}]}{2m_{1}m_{2}c^{2}\vec{k}^{2}}+\frac{(\vec{\sigma}_{1}\vec{k})(\vec{\sigma}_{2}\vec{k})}{4m_{1}m_{2}c^{2}\vec{k}^{2}}-\frac{\vec{\sigma}_{1}\vec{\sigma}_{2}}{4m_{1}m_{2}c^{2}}\right].$

Here $\vec{p}$ is the electron momentum operator, $\vec{\sigma}$ is the Pauli matrix, $m$ is the particle mass and the index 1 or 2 refers to the corresponding particle.

The subsequent evaluation repeats calculations performed in [43], i.e. Fourier transform should be applied to Eq. (7). Integrating with $\int\frac{d^{3}k}{(2\pi)^{3}}e^{i\vec{k}\vec{r}}$ the scattering amplitude $U(\vec{p}_{1},\vec{p}_{2},\vec{k})$, one can be find a coordinate representation that will contain relativistic corrections from the scalar and transverse parts of the photon propagator $D_{00}(x,x^{\prime})$ and $D_{ij}(x,x^{\prime})$. Then the operators $\vec{p}_{1}$ and $\vec{p}_{2}$ should be replaced by the $\vec{p}_{1}=-i\nabla_{1}$ and $\vec{p}_{2}=-i\nabla_{2}$. In the thermal case, however, the Fourier transform to coordinate representation is given as $2\int\frac{d^{3}k}{(2\pi)^{3}}n_{\beta}(|\vec{k}|)e^{i\vec{k}\vec{r}}$, where the factor 2 occurs in a result of integration along the contour $C_{1}$.

The Fourier transform of the first term in Eq. (7) results in the expression (5), where the regularization of divergent contribution at $|\vec{k}|\rightarrow 0$ was performed by introducing a coincidence limit [37] (see also Appendix A). For a positronium atom, the masses $m_{1}$ and $m_{2}$ should be equated to each other. Then the second and third terms in Eq. (7), after integration over angles, convert to

$\displaystyle(2)+(3)\rightarrow-\frac{e^{2}}{\pi}\frac{1}{m^{2}c^{2}}\int\limits_{0}^{\infty}d\kappa\,n_{\beta}(\kappa)\frac{\kappa\sin\kappa r_{12}}{r_{12}}.$

(8)

Truncation of the Taylor series by the first two terms yields a constant contribution and a correction proportional to the square $r_{12}$. The constant contribution does not depend on the state and, therefore, vanishes for the energy difference of atomic states. It can be also found that the coincidence limit [37] regularizes such contributions along with divergences. In other words, subtracting the $r_{12}\rightarrow 0$ limit, the regular expression is

$\displaystyle(2)+(3)\rightarrow\frac{4e^{2}}{\pi}\frac{\zeta(5)}{\beta^{5}m^{2}c^{2}}r_{12}^{2},$

(9)

where we have integrated over $\kappa\equiv|\vec{k}|$.

The fourth term can be integrated using the equality:

	$\displaystyle\frac{4\pi e^{2}}{m_{1}m_{2}c^{2}}2\int\frac{d^{3}k}{(2\pi)^{3}}n_{\beta}(\kappa)e^{i\vec{k}\vec{r}_{12}}\frac{(\vec{k}\vec{p}_{1})(\vec{k}\vec{p}_{2})}{\vec{k}^{4}}=$		(10)
	$\displaystyle\frac{4\pi e^{2}}{m^{2}c^{2}}2\int\frac{d^{3}k}{(2\pi)^{3}}n_{\beta}(\kappa)\frac{(\vec{\nabla}_{1}\vec{p}_{1})(\vec{\nabla}_{2}\vec{p}_{2})}{\vec{k}^{4}}e^{i\vec{k}\vec{r}_{12}}.$

Then, integrating over angles and acting by the gradient operators, the formula (10) reduces to

	$\displaystyle\frac{4e^{2}}{\pi m^{2}c^{2}}\int\limits_{0}^{\infty}d\kappa\,\frac{n_{\beta}(\kappa)}{\kappa^{2}}\left[\frac{\cos\kappa r_{12}}{r_{12}^{2}}-\frac{\sin\kappa r_{12}}{\kappa r_{12}^{3}}\right](\vec{p}_{1}\vec{p}_{2})$		(11)
	$\displaystyle-\frac{4e^{2}}{\pi m^{2}c^{2}}\int\limits_{0}^{\infty}d\kappa\,\frac{n_{\beta}(\kappa)}{\kappa^{2}}\left[\frac{3\cos\kappa r_{12}}{r_{12}^{4}}-\frac{3\sin\kappa r_{12}}{\kappa r_{12}^{5}}\right.$
	$\displaystyle\left.+\frac{\kappa\sin\kappa r_{12}}{r_{12}^{3}}\right](\vec{r}_{12}\vec{p}_{1})(\vec{r}_{12}\vec{p}_{2}).$

The coincidence limit $r_{12}\rightarrow 0$ in this case is

$\displaystyle-\frac{4e^{2}}{3\pi m^{2}c^{2}}\int\limits_{0}^{\infty}d\kappa\,n_{\beta}(\kappa),$

(12)

which cancels the divergence in Eq. (11). Finally, the Fourier transform of fourth term in Eq. (7) reduces to

$\displaystyle\frac{4\zeta(3)e^{2}}{15\pi\beta^{3}m^{2}c^{2}}\left[r_{12}^{2}(\vec{p}_{1}\vec{p}_{2})+2(\vec{r}_{12}\vec{p}_{1})(\vec{r}_{12}\vec{p}_{2})\right].$

(13)

Evaluation of the fifth contribution in Eq. (7) results in

$\displaystyle-\frac{\vec{p}_{1}\vec{p}_{2}}{m^{2}c^{2}\vec{k}^{2}}\rightarrow\frac{4\zeta(3)e^{2}}{3\pi\beta^{3}m^{2}c^{2}}r_{12}^{2}(\vec{p}_{1}\vec{p}_{2}).$

(14)

The Fourier transform for the next four terms can be performed using the substitution $\vec{k}\rightarrow-i\vec{\nabla}_{1}$. Acting by the gradient operator on the expression arising after angular integration, we find

	$\displaystyle\frac{i\vec{\sigma}_{1}[\vec{k}\times\vec{p}_{1}]}{4m^{2}c^{2}\vec{k}^{2}}\rightarrow-\frac{2\zeta(3)e^{2}}{3\pi\beta^{3}m^{2}c^{2}}\left(\vec{\sigma}_{1}[\vec{r}_{12}\times\vec{p}_{1}]\right),$
	$\displaystyle-\frac{i\vec{\sigma}_{2}[\vec{k}\times\vec{p}_{2}]}{4m^{2}c^{2}\vec{k}^{2}}\rightarrow\frac{2\zeta(3)e^{2}}{3\pi\beta^{3}m^{2}c^{2}}\left(\vec{\sigma}_{2}[\vec{r}_{12}\times\vec{p}_{2}]\right),$		(15)
	$\displaystyle-\frac{i\vec{\sigma}_{1}[\vec{k}\times\vec{p}_{2}]}{2m_{1}m_{2}c^{2}\vec{k}^{2}}\rightarrow\frac{4\zeta(3)e^{2}}{3\pi\beta^{3}m^{2}c^{2}}\left(\vec{\sigma}_{1}[\vec{r}_{12}\times\vec{p}_{2}]\right),$
	$\displaystyle\frac{i\vec{\sigma}_{2}[\vec{k}\times\vec{p}_{1}]}{2m_{1}m_{2}c^{2}\vec{k}^{2}}\rightarrow-\frac{4\zeta(3)e^{2}}{3\pi\beta^{3}m^{2}c^{2}}\left(\vec{\sigma}_{2}[\vec{r}_{12}\times\vec{p}_{1}]\right).$

Similar calculations for the last two contributions lead to

	$\displaystyle\frac{(\vec{\sigma}_{1}\vec{k})(\vec{\sigma}_{2}\vec{k})}{4m_{1}m_{2}c^{2}\vec{k}^{2}}\rightarrow-\frac{4\zeta(5)e^{2}}{5\pi\beta^{5}m^{2}c^{2}}\left[r_{12}^{2}(\vec{\sigma}_{1}\vec{\sigma}_{2})+2(\vec{r}_{12}\vec{\sigma}_{1})(\vec{r}_{12}\vec{\sigma}_{2}\right]$
	$\displaystyle-\frac{\vec{\sigma}_{1}\vec{\sigma}_{2}}{4m_{1}m_{2}c^{2}}\rightarrow\frac{4\zeta(5)e^{2}}{\pi\beta^{5}m^{2}c^{2}}r_{12}^{2}(\vec{\sigma}_{1}\vec{\sigma}_{2}).\qquad$		(16)

We emphasize that the replacement $\vec{k}\rightarrow-i\vec{\nabla}_{1}$ assumes first the action of the gradient operator and then going to the coincidence limit.

These calculations, however, represent the transformation of the scattering (’direct’) diagram. There is also a second independent contribution corresponding to the ’exchange’ or ’annihilation’ diagram since the wave function of the electron-positron system need not be asymmetric [4]. In the latter, it is convenient to use the Feynman gauge for the thermal part of the photon propagator, Eq. (2). In the approximation of ’almost nonrelativistic’ particles, it can be reduced to

$\displaystyle D_{\mu\nu}^{\beta}(x,x^{\prime})=-\frac{\pi g_{\mu\nu}}{m^{2}c^{2}}\int\limits_{C_{1}}\frac{d^{4}k}{(2\pi)^{4}}e^{ik(x-x^{\prime})}n_{\beta}(|\vec{k}|).$

(17)

Repeating the calculations given in [4], we find the annihilation amplitude in a coordinate space as

$\displaystyle U^{\mathrm{(ann)}}=\frac{\pi e^{2}}{2m^{2}c^{2}}\int\limits_{C_{1}}\frac{d^{4}k}{(2\pi)^{4}}e^{ik(x-x^{\prime})}n_{\beta}(|\vec{k}|)\left[3+\vec{\sigma}_{1}\vec{\sigma}_{2}\right].$

(18)

Performing the remaining integrations, we obtain

	$\displaystyle U^{\mathrm{(ann)}}=\frac{ie^{2}\left[3+\vec{\sigma}_{1}\vec{\sigma}_{2}\right]}{4\pi m^{2}c^{2}\beta^{2}\,r_{12}}\left[\psi^{(1)}\left(1+\frac{ir}{\beta}\right)-\psi^{(1)}\left(1-\frac{ir}{\beta}\right)\right]$
	$\displaystyle\approx\left[\frac{e^{2}\zeta(3)}{\pi\beta^{3}c^{2}m^{2}}-\frac{2e^{2}\zeta(5)}{\pi\beta^{5}c^{2}m^{2}}r_{12}^{2}\right]\left(3+\vec{\sigma}_{1}\vec{\sigma}_{2}\right).\qquad$		(19)

Here $\psi^{(1)}(x)$ represents the first derivative of the Digamma function [45]. The first term in square brackets is state-independent (in the sense of $r_{12}$) and is, therefore, canceled by the coincidence limit. The second contribution is $\beta^{2}$ times smaller and, according to the estimations given after Eq. (6), insignificant.

The total contribution to the binding energy in the lowest order in temperature for a positronium atom can be written as

	$\displaystyle U(\vec{p}_{1},\vec{p}_{2},\vec{r}_{12})=-\frac{4\zeta(3)e^{2}}{3\pi\beta^{3}}r_{12}^{2}+\frac{8\zeta(3)e^{2}}{5\pi\beta^{3}m^{2}c^{2}}r_{12}^{2}(\vec{p}_{1}\vec{p}_{2})$
	$\displaystyle+\frac{8\zeta(3)e^{2}}{15\pi\beta^{3}m^{2}c^{2}}\vec{r}_{12}(\vec{r}_{12}\vec{p}_{2})\vec{p}_{1}-\frac{2\zeta(3)e^{2}}{3\pi\beta^{3}m^{2}c^{2}}\left(\vec{\sigma}_{1}[\vec{r}_{12}\times\vec{p}_{1}]\right)$
	$\displaystyle+\frac{2\zeta(3)e^{2}}{3\pi\beta^{3}m^{2}c^{2}}\left(\vec{\sigma}_{2}[\vec{r}_{12}\times\vec{p}_{2}]\right)+\frac{4\zeta(3)e^{2}}{3\pi\beta^{3}m^{2}c^{2}}\left(\vec{\sigma}_{1}[\vec{r}_{12}\times\vec{p}_{2}]\right)$
	$\displaystyle-\frac{4\zeta(3)e^{2}}{3\pi\beta^{3}m^{2}c^{2}}\left(\vec{\sigma}_{2}[\vec{r}_{12}\times\vec{p}_{1}]\right).\qquad$		(20)

In the center-of-mass system, the electron and positron momentum operators in positronium are $\vec{p}_{1}=-\vec{p}_{2}\equiv\vec{p}$, where $\vec{p}=-i\vec{\nabla}$ is the operator of the momentum of relative motion corresponding to relative position vector $\vec{r}\equiv\vec{r}_{12}=\vec{r}_{1}-\vec{r}_{2}$. Then the thermal contribution Eq. (II) reduces to

	$\displaystyle U=-\frac{4\zeta(3)e^{2}}{3\pi\beta^{3}}r^{2}-\frac{8\zeta(3)e^{2}}{5\pi\beta^{3}m^{2}c^{2}}r^{2}p^{2}$
	$\displaystyle-\frac{8\zeta(3)e^{2}}{15\pi\beta^{3}m^{2}c^{2}}\vec{r}(\vec{r}\vec{p})\vec{p}-\frac{4\zeta(3)e^{2}}{\pi\beta^{3}m^{2}c^{2}}\left(\vec{S}\vec{l}\right),$		(21)

where we have introduced the operators of the total spin $\vec{S}=\frac{1}{2}(\vec{\sigma}_{1}+\vec{\sigma}_{2})$ and the orbital angular momentum $\vec{l}=\left[\vec{r}\times\vec{p}\right]$ (details of the corresponding simplifications can be found in [4, 2]). In principle, the operator (II) is sufficient for calculations, but it can be further simplified using the relation $\vec{l}^{2}\equiv\left[\vec{r}\times\vec{p}\right]^{2}=r^{2}p^{2}-(\vec{r}\vec{p})^{2}+i(\vec{r}\vec{p})$. Then, employing the commutator $\left[r_{i},p_{j}\right]=i\delta_{ij}$ in the third term of Eq. (II), we find

	$\displaystyle U=-\frac{4\zeta(3)e^{2}}{3\pi\beta^{3}}r^{2}-\frac{16\zeta(3)e^{2}}{15\pi\beta^{3}m^{2}c^{2}}r^{2}p^{2}$
	$\displaystyle-\frac{8\zeta(3)e^{2}}{15\pi\beta^{3}m^{2}c^{2}}\vec{l}^{2}-\frac{4\zeta(3)e^{2}}{\pi\beta^{3}m^{2}c^{2}}\left(\vec{S}\vec{l}\right).$		(22)

Subsequent calculations correspond to the evaluation of the average values of the operators given by Eq. (II). Since the Bohr’s radius in the positronium atom is estimated as $\langle r\rangle\sim 1/(\frac{1}{2}m\alpha)$, we use the hydrogen ratio $\langle a|r^{2}|a\rangle=\frac{n_{a}^{2}}{2}(5n_{a}^{2}+1-3l_{a}(l_{a}+1))$, which converts to $\langle a|r^{2}|a\rangle=2n_{a}^{2}(5n_{a}^{2}+1-3l_{a}(l_{a}+1))$ for an arbitrary $a$-state in positronium. The average value of the third term in Eq. (II) gives $l_{a}(l_{a}+1)$ and the fourth can be found as $\langle a|\left(\vec{S}\vec{l}\right)|a\rangle=\frac{1}{2}\left[j_{a}(j_{a}+1)-l_{a}(l_{a}+1)-S_{a}(S_{a}+1)\right]$.

Finally, to evaluate the second term, we take into account that $p^{2}\psi_{a}=(E+1/r)\psi_{a}$, which follows from the Schrödinger equation for positronium, where $E$ represents the energy levels of positronium: $E_{n_{a}}=-1/4n_{a}^{2}$. Then, the average value of $(r^{2}p^{2})\psi_{a}=\left(-r^{2}/(4n_{a}^{2})+r\right)\psi_{a}$ can be easily calculated using $\langle a|r|a\rangle=3n_{a}^{2}-l_{a}(l_{a}+1)$, and, therefore, $\langle a|r^{2}p^{2}|a\rangle=\frac{1}{2}\left[n^{2}_{a}-1+l_{a}(l_{a}+1)\right]$. In total we have

	$\displaystyle\langle a|U|a\rangle=-\frac{8\zeta(3)\alpha^{3}}{3\pi\beta^{3}}n_{a}^{2}\left[5n_{a}^{2}+1-3l_{a}(l_{a}+1)\right]$
	$\displaystyle-\frac{8\zeta(3)\alpha^{5}}{15\pi\beta^{3}}\left[n_{a}^{2}-1+2l_{a}(l_{a}+1)\right]\qquad$		(23)
	$\displaystyle-\frac{2\zeta(3)\alpha^{5}}{\pi\beta^{3}}\left[j_{a}(j_{a}+1)-l_{a}(l_{a}+1)-S_{a}(S_{a}+1)\right],$

This expression is written in atomic units in conjunction with $1/\beta=k_{B}T=3.16681\times 10^{-6}T$.

The numerical results for some low-lying states in the positronium atom are given in Table 1. It should be emphasized here that the $4\pi$ factor, corresponding to the Heaviside definition of charge, was lost in [37, 38, 43], i.e. in these works, one should additionally divide by $4\pi$.

The values listed in Table 1 demonstrate that the effects described above are beyond the precision of modern laboratory experiments. The scale factor $T^{3}$ can be prolonged to astrophysical conditions, taking into account the lowest order thermal correction, the result is $60$ kHz for the Ly_α line at a recombination temperature of the Universe, $3000$ K.

In this part of the work, we briefly describe the thermal Stark shift for Ps. The corresponding derivations can be attributed to the earlier work [34], where a quantum mechanical description of the ac-Stark shift and the transition rate induced by blackbody radiation was given. However, here we apply the QED formalism discussed in [36] for the appropriate derivations and further calculations in the positronium atom. According to [36], in this case, within the framework of the QED approach at finite temperatures, it is sufficient to evaluate the one-loop self-energy correction (see also [37]). Then, replacing the ’ordinary’ photon line by a thermal one, the real part of this correction gives the ac-Stark effect, while the imaginary part is the level width (the sum of all partial transition to the lower and upper states) induced by the BBR.

After several successive conversions of the thermal photon propagator Eq. (2) (see [37] for details), the representation used in [36] can be found:

$\displaystyle D^{\beta}_{\mu\,\nu}=-\frac{g_{\mu\,\nu}}{\pi r_{12}}\int\limits_{-\infty}^{+\infty}d\omega n_{\beta}(|\omega|)\sin{|\omega|r_{12}}e^{-i\omega(t_{1}-t_{2})}.\qquad$

(24)

Then the energy shift for an arbitrary state $a$ is

$\displaystyle\Delta E_{a}^{\beta}=\frac{e^{2}}{\pi}\sum\limits_{n}\left(\frac{1-\vec{\alpha}_{1}\vec{\alpha}_{2}}{r_{12}}I^{\beta}_{na}(r_{12})\right)_{anna},$

(25)

where

$\displaystyle I^{\beta}_{na}(r_{12})=\int\limits_{-\infty}^{+\infty}d\omega n_{\beta}(|\omega|)\frac{\sin{|\omega|r_{12}}}{E_{n}(1-i0)-E_{a}+\omega}.$

(26)

Here the sum runs over the entire spectrum $n$, including the continuum, and the matrix element is to be understood as $\left(\hat{A}(12)\right)_{abcd}\equiv\langle a(1)b(2)|\hat{A}|c(1)d(2)\rangle$ [6].

Omitting the description of the effect associated with the finite lifetime of states (see [36] for details), the result can be obtained using the Sokhotski-Plemelj theorem:

$\displaystyle\lim\limits_{\epsilon\rightarrow 0}\frac{1}{x\pm i\epsilon}=\mathrm{P.V.}\left(\frac{1}{x}\right)\mp i\pi\delta(x),$

(27)

where $\mathrm{P.V.}$ means the principal value. Then, one can find

	$\displaystyle I^{\beta}_{na}(r_{12})=\sum\limits_{\pm}\mathrm{P.V.}\int\limits_{0}^{\infty}d\omega\,n_{\beta}(\omega)\frac{\sin\omega r_{12}}{E_{n}-E_{a}\pm\omega}$		(28)
	$\displaystyle+i\,\pi n_{\beta}(|E_{an}|)\sin|E_{an}|r_{12},$

where the notation $E_{na}=E_{n}-E_{a}$ was introduced, and $\sum\limits_{\pm}$ means the sum of two contributions with $-$ and $+$ before $\omega$ in the energy denominator.

Equation (28) already demonstrates the existence of real and imaginary contributions for the energy shift (25). To give them a physical interpretation, it is useful to consider the nonrelativistic limit, which we arrive at employing Taylor’s series expansion of the function $\sin\omega r_{12}\approx\omega r_{12}-\frac{1}{6}(\omega r_{12})^{3}$. Then the imaginary part (after successive but ordinary calculations [36, 37]) is

$\displaystyle\Gamma_{a}^{\beta}\equiv-2{\rm Im}\Delta E_{a}^{\beta}=\frac{4}{3}e^{2}\sum\limits_{n}\left|\langle a|\vec{r}|n\rangle\right|^{2}n_{\beta}(|\omega_{an}|)\omega_{an}^{3}.\qquad$

(29)

The real part should be considered more carefully. By truncating the $\sin$ Taylor series with two terms, the real part is reduced to

	$\displaystyle{\rm Re}\Delta E_{a}^{\beta}\approx\frac{e^{2}}{\pi}\sum\limits_{n}\mathrm{P.V.}\int\limits_{0}^{\infty}d\omega n_{\beta}(\omega)\left[\frac{1}{E_{n}-E_{a}-\omega}\qquad\right.$		(30)
	$\displaystyle\left.+\frac{1}{E_{n}-E_{a}+\omega}\right]\left(\omega-\omega\vec{\alpha}_{1}\vec{\alpha}_{2}-\frac{1}{6}\omega^{3}r_{12}^{2}\right)_{an\,na}.\qquad$

Hereinafter, we use that the sum in square brackets is $2E_{na}/(E_{na}^{2}-\omega^{2})$. Then the first term is equal to zero due to the orthogonality property of wave functions and the presence of $E_{na}=0$ in numerator for $n=a$.

Applying the nonrelativistic limit to matrix elements with $\vec{\alpha}$-matrix ($(\vec{\alpha}_{1}\vec{\alpha}_{2})_{an\,na}=E_{an}^{2}(\vec{r}_{1}\vec{r}_{2})$ and ratio $r_{12}^{2}=r_{1}^{2}+r_{2}^{2}-2(\vec{r}_{1}\vec{r}_{2})$, see [6]), we obtain

	$\displaystyle{\rm Re}\Delta E_{a}^{\beta}=\frac{e^{2}}{\pi}\sum\limits_{n}\mathrm{P.V.}\int\limits_{0}^{\infty}d\omega\frac{2E_{na}n_{\beta}(\omega)}{E_{na}^{2}-\omega^{2}}\times$		(31)
	$\displaystyle\left(-\omega E_{na}^{2}+\frac{1}{3}\omega^{3}\right)(\vec{r}_{1}\vec{r}_{2})_{an\,na}.$

The expression (31) can be simplified by substituting the zero contribution $\pm\omega^{3}$. Then,

	$\displaystyle{\rm Re}\Delta E_{a}^{\beta}=\frac{4e^{2}}{3\pi}\sum\limits_{n}\mathrm{P.V.}\int\limits_{0}^{\infty}d\omega\frac{E_{an}n_{\beta}(\omega)}{E_{an}^{2}-\omega^{2}}\left|\langle a|\vec{r}|n\rangle\right|^{2}\qquad$		(32)
	$\displaystyle+\frac{2e^{2}}{\pi}\sum\limits_{n}\mathrm{P.V.}\int\limits_{0}^{\infty}d\omega\frac{E_{na}n_{\beta}(\omega)}{E_{na}^{2}-\omega^{2}}\left(\omega^{3}-\omega E_{na}^{2}\right)\left|\langle a|\vec{r}|n\rangle\right|^{2}.$

The first contribution represents well-known ac-Stark shift induced by the blackbody radiation field:

$\displaystyle\Delta E_{a}^{\rm Stark}=\frac{4e^{2}}{3\pi}\sum\limits_{n}\mathrm{P.V.}\int\limits_{0}^{\infty}d\omega\frac{E_{an}n_{\beta}(\omega)\omega^{3}}{E_{an}^{2}-\omega^{2}}\left|\langle a|\vec{r}|n\rangle\right|^{2}.\qquad$

(33)

Parametric estimation of Eqs. (29) and (33) arises as indicated above and is given as $m\alpha^{5}(k_{B}T)^{4}/Z^{4}$ in relativistic units, where the Boltzmann should be taken in atomic units. To get result in completely atomic units it is necessary to divide this estimate by $m\alpha^{2}$. Thus, for evaluating ac-Stark, ${\rm Re}\Delta E_{a}^{\beta}$, and level width, $\Gamma_{a}^{\beta}$, in positronium atom it is sufficient to take into account the coefficient $1/2$, arising from the reduced mass. We do not provide the corresponding values here, assuming that they can be easily found using the results of [34] and concluding, that the ac-Stark shift does not exceed a few Hz at $T=300$ K for low-lying states, while the line broadening remains insignificant at room temperature.

The second term in Eq. (32) (indicated by a cross below) is more delicate. First of all, the energy denominator is canceled by the numerator, which gives

$\displaystyle\Delta E_{a}^{\beta,\times}=\frac{2e^{2}}{\pi}\sum\limits_{n}\mathrm{P.V.}\int\limits_{0}^{\infty}d\omega\,\omega n_{\beta}(\omega)E_{an}\left|\langle a|\vec{r}|n\rangle\right|^{2}.\qquad$

(34)

Then, using the sum rule for the oscillator strength, one can find that this contribution is constant and independent of states. Thus, it represents an immeasurable contribution to the atomic energy of the bound electron. Basically, we can just throw it away [46]. However, evaluation of the coincidence limit [37] shows the same result and, therefore, cancels this contribution and results in a gauge invariant expression (33).