1. Introduction
One of the remarkable consequences of the quantization of electromagnetic fields is the existence of vacuum fluctuations as a result of the Heisenberg uncertainty principle, which has profoundly changed our usual understanding of a vacuum. The vacuum fluctuations are deemed as the origin of various physical phenomena, such as the spontaneous emission [1] and the Lamb shift [2], the attractive force between two perfectly conducting neutral plates in a vacuum, which is known as the Casimir effect [3], as well as a force felt by a neutral atom near a perfectly conducting plate, due to the coupling between the atom and the vacuum electromagnetic fluctuations modified by the presence of the boundary [4].
When two atoms are involved, a vacuum-fluctuation-induced interatomic interaction then appears, which can be physically understood as follows. One atom interacts with fluctuating vacuum electromagnetic fields, and a radiative field is subsequently produced, which then acts on the other atom, and vice versa. As a result, an interaction potential, which is also generally referred to as the van der Waals or Casimir–Polder interaction potential, is induced [4–11]. Based on the standard perturbation theory, London found that, for two ground-state atoms with an interatomic distance L much smaller than the characteristic transition wavelength of the atoms, the interatomic interaction potential behaves as L−6 [5], and Casimir and Polder further showed that in the far region where the interatomic distance L is much larger than the atomic transition wavelength, the potential decays more rapidly with a behavior of L−7 rather than L−6, due to the effect of retardation [4]. A few years later, Lifshitz et al developed a macroscopic theory for the interaction of dielectric bodies, and by going to the limiting case of rarefied media, they recovered both the nonretarded and retarded interatomic interaction potential [12, 13].
The early works mentioned above are about the interatomic interaction energy of two ground-state inertial atoms coupled to the fluctuating electromagnetic field in a trivial flat space in a vacuum. Over the past decades, endeavors have also been devoted to this fluctuation-induced interaction in more complex circumstances, such as atoms in a bounded space [14, 15], in noninertial motions [16, 17], in a thermal bath [14, 18–24], and even in a curved spacetime [25, 26]. For the effects of a thermal bath, let us note that after the work of Lifshitz and his collaborators [12, 13], McLachlan revisited the long-range dispersion forces in a non-uniform dielectric at a finite temperature [18]. He derived most of the results of Lifshitz et al in a more elementary way and discovered that the interaction energy of two atoms in a thermal bath at a low temperature T, which is much smaller than the energy gap of the atoms ω0, scales as L−6 at small interatomic distances and as L−7 at large interatomic distances when the retardation effect is important; while at sufficiently long distances, L−1 ≪ T ≪ ω0, the temperature effects mask the retardation effect and the interatomic interaction energy scales as TL−6 [18]. Later, this temperature and distance dependence was also shown by Boyer in the theory of classical electrodynamics with classical electromagnetic zero-point radiation [19], by Milonni and others via considering the energy of an induced dipole moment of a molecule in an electric field correlated with the second molecule [20, 21], and by Wennerstr$\ddot{{\rm{o}}}{\rm{m}}$ and his collaborators with two different methods (the Lifshitz theory and the generalized fourth-order quantum-electrodynamical perturbation theory) [22]. In addition, by taking the dilute medium limit of the Lifshitz theory [12, 13] for two parallel dielectric slabs in a thermal environment, Ninham and Daicic also recovered the London potential in the nonretarded region and revealed that as the retarded limit is approached with increasing interatomic separation, temperature effects become increasingly significant and ultimately the London description must break down, even at a fixed temperature [23]. Recently, the study of the thermal effects on the dispersion interaction energy of two ground-state atoms has also been extended to the case of atoms out of thermal equilibrium [24].
When a two-atom system in a thermal bath is concerned, there are three characteristic physical parameters, i.e. the interatomic separation, L, the transition wavelength of the atoms, $\lambda \sim {\omega }_{0}^{-1}$, and the wavelength of thermal photons, β = T−1, and so there are in general six typical regions for the problem under consideration. However, the previous studies on this topic are only concerned with one or a few regions, and to the best of our knowledge, a complete analysis of the interatomic interaction at a finite temperature is still lacking. In this paper, we aim to fill this gap, i.e. we will give a complete study for the interaction potential of two ground-state atoms in a thermal bath at an arbitrary temperature.
We exploit the formalism proposed by Dalibard, Dupont-Roc, and Cohen-Tannoudji (the DDC formalism) [27, 28], which separates the contributions of field fluctuations and the contributions of radiation reaction of atoms, to deal with the interaction potential. It is worth pointing out that the DDC formalism has been widely used, during the past years, to study various second-order perturbation quantum effects, such as the energy shifts of a single atom [29–36] and the resonance interaction energy of two entangled atoms [37–42], and very recently, it was generalized from the second order to the fourth order to study the interatomic interaction between two ground-state atoms in vacuum [16, 26, 43]. In this paper, as an attempt to utilize the DDC formalism to deal with the effects of a thermal bath on the interatomic interaction between two ground-state atoms, we first generalize the fourth-order DDC formalism in vacuum [43] to a thermal bath. For simplicity, we start with the model of atoms linearly coupled to a massless scalar field and derive basic formulas for the contribution of thermal fluctuations (tf-contribution) and that of the radiation reaction of atoms (rr-contribution), with which we then calculate the interatomic interaction energy from both contributions and analyze its behaviors in all typical regions in detail. Here let us note that the scalar field model has been widely used in the study of radiative properties of atoms in various circumstances in recent years [25, 26, 29, 30, 33, 34, 36–40, 43–49].
2. General formalism of the contributions of thermal fluctuations and radiation reaction of atoms to interatomic interaction
We consider that, two spatially separated atoms are interacting with massless scalar fields in a thermal bath at an arbitrary temperature T. Each atom, labeled by ξ (ξ = A or B), is assumed to be of two energy levels with an energy gap ωξ, and ∣gξ⟩ and ∣eξ⟩ represent atomic ground state and excited state, with energies $-\tfrac{{\omega }_{\xi }}{2}$ and $+\tfrac{{\omega }_{\xi }}{2}$ respectively. The total Hamiltonian describing the ‘atoms + field' system is given by1 )-(4 ), we shall first establish the general fourth-order DDC formalism for the contributions of thermal fluctuations and radiation reaction of atoms to the interatomic interaction energy.
$\begin{eqnarray}H(\tau )={H}_{S}(\tau )+{H}_{F}(\tau )+{H}_{I}(\tau ),\end{eqnarray}$
with respect to the proper time of the atoms τ. Here, HS(τ) and HF(τ) are the free Hamiltonian of the two-atom system and that of the scalar fields, which can be written as $\begin{eqnarray}{H}_{S}(\tau )={\omega }_{A}{R}_{3}^{A}(\tau )+{\omega }_{B}{R}_{3}^{B}(\tau ),\end{eqnarray}$
$\begin{eqnarray}{H}_{F}(\tau )=\int {{\rm{d}}}^{3}{\boldsymbol{k}}\,{\omega }_{{\boldsymbol{k}}}{a}_{{\boldsymbol{k}}}^{\dagger }(t(\tau )){a}_{{\boldsymbol{k}}}(t(\tau ))\displaystyle \frac{{\rm{d}}t}{{\rm{d}}\tau },\end{eqnarray}$
where ${R}_{3}^{\xi }=\tfrac{1}{2}(| {e}_{\xi }\rangle \langle {e}_{\xi }| -| {g}_{\xi }\rangle \langle {g}_{\xi }| )$, and ${a}_{{\boldsymbol{k}}}^{\dagger }$ and ak are respectively the creation and annihilation operators of the field modes with momentum k. The four-dimensional coordinates are denoted by x = (t, x), and the atom-field interaction Hamiltonian HI(τ) is expressed as $\begin{eqnarray}{H}_{I}(\tau )=\mu {R}_{2}^{A}(\tau )\phi ({x}_{A}(\tau ))+\mu {R}_{2}^{B}(\tau )\phi ({x}_{B}(\tau )).\end{eqnarray}$
Here μ is a very small coupling constant, ${R}_{2}^{\xi }=\tfrac{{\rm{i}}}{2}({R}_{-}^{\xi }-{R}_{+}^{\xi })$, with ${R}_{+}^{\xi }=| {e}_{\xi }\rangle \langle {g}_{\xi }| $ and ${R}_{-}^{\xi }=| {g}_{\xi }\rangle \langle {e}_{\xi }| $ being atomic raising and lowering operators, and φ(x) is the operator of scalar fields, $\begin{eqnarray}\phi (x)=\int {{\rm{d}}}^{3}{\boldsymbol{k}}\,{g}_{{\boldsymbol{k}}}[{a}_{{\boldsymbol{k}}}(t){f}_{{\boldsymbol{k}}}({\boldsymbol{x}})+{a}_{{\boldsymbol{k}}}^{\dagger }(t){f}_{{\boldsymbol{k}}}^{* }({\boldsymbol{x}})],\end{eqnarray}$
where ${g}_{{\boldsymbol{k}}}={\left[2{\omega }_{{\boldsymbol{k}}}{\left(2\pi \right)}^{3}\right]}^{-1/2}$, and fk(x) is the scalar-field mode with momentum k. With the Hamiltonians of the ‘atoms+field' system in equations (Following similar procedures we have adopted in [43] for the interatomic interaction energy between two ground-state atoms in a vacuum, we obtain and solve the Heisenberg equations of motion of the dynamical variables of the atoms and the scalar fields, and separate each solution into a free part and a source part. We then perturbatively expand each source part with respect to the small coupling constant μ. By exploiting the symmetric ordering between the operators of the atoms and the fields [27, 28], we evaluate the contributions of the free field and the source field to the variation rate of the Hamiltonian of the two-atom system. After taking the average of them over the thermal state of the scalar fields ∣β⟩, we can pick out the following effective Hamiltonian of the thermal fluctuations contribution (tf-contribution) caused by the free field, and that of the radiation reaction contribution (rr-contribution) caused by the source field,6 ) and (7 ), over the two-atom state, ∣gAgB⟩, and doing some simplifications, we finally arrive at the general formula for the tf-contribution to the interatomic interaction energy for two ground-state atoms,
$\begin{eqnarray}\begin{array}{l}{\left({H}_{S}(\tau )\right)}_{\mathrm{tf}}^{\mathrm{eff}}\\ =\displaystyle \frac{1}{8}{\rm{i}}{\mu }^{4}{\displaystyle \int }_{{\tau }_{0}}^{\tau }{\rm{d}}{\tau }_{1}{\displaystyle \int }_{{\tau }_{0}}^{{\tau }_{1}}{\rm{d}}{\tau }_{2}{\displaystyle \int }_{{\tau }_{0}}^{{\tau }_{2}}{\rm{d}}{\tau }_{3}\\ \times \,\left\langle \beta \left|\left\{{\phi }^{f}({x}_{A}(\tau )),{\phi }^{f}({x}_{B}({\tau }_{3}))\right\}\left[{\phi }^{f}({x}_{B}({\tau }_{2})),{\phi }^{f}({x}_{A}({\tau }_{1}))\right]\right|\beta \right\rangle \\ \times \,\left[{R}_{2}^{A,f}(\tau ),{R}_{2}^{A,f}({\tau }_{1})\right]\left[{R}_{2}^{B,f}({\tau }_{3}),{R}_{2}^{B,f}({\tau }_{2})\right]\\ +\,A\rightleftharpoons B\,\mathrm{term},\end{array}\end{eqnarray}$
and $\begin{eqnarray}\begin{array}{l}{\left({H}_{S}(\tau )\right)}_{\mathrm{rr}}^{\mathrm{eff}}\\ =\displaystyle \frac{1}{8}{\rm{i}}{\mu }^{4}{\displaystyle \int }_{{\tau }_{0}}^{\tau }{\rm{d}}{\tau }_{1}{\displaystyle \int }_{{\tau }_{0}}^{{\tau }_{1}}{\rm{d}}{\tau }_{2}{\displaystyle \int }_{{\tau }_{0}}^{{\tau }_{2}}{\rm{d}}{\tau }_{3}\\ \times \,\left\langle \beta \left|\left[{\phi }^{f}({x}_{B}({\tau }_{2})),{\phi }^{f}({x}_{A}({\tau }_{1}))\right]\left[{\phi }^{f}({x}_{B}({\tau }_{3})),{\phi }^{f}({x}_{A}(\tau ))\right]\right|\beta \right\rangle \\ \times \left\{{R}_{2}^{A,f}({\tau }_{1}),{R}_{2}^{A,f}(\tau )\right\}\left[{R}_{2}^{B,f}({\tau }_{2}),{R}_{2}^{B,f}({\tau }_{3})\right]\\ +\displaystyle \frac{1}{8}{\rm{i}}{\mu }^{4}{\displaystyle \int }_{{\tau }_{0}}^{\tau }{\rm{d}}{\tau }_{1}{\displaystyle \int }_{{\tau }_{0}}^{{\tau }_{1}}{\rm{d}}{\tau }_{2}{\displaystyle \int }_{{\tau }_{0}}^{{\tau }_{2}}{\rm{d}}{\tau }_{3}\\ \times \,\left\langle \beta \left|\left[{\phi }^{f}({x}_{A}({\tau }_{1})),{\phi }^{f}({x}_{B}({\tau }_{3}))\right]\left[{\phi }^{f}({x}_{B}({\tau }_{2})),{\phi }^{f}({x}_{A}(\tau ))\right]\right|\beta \right\rangle \\ \times \left\{{R}_{2}^{A,f}({\tau }_{1}),{R}_{2}^{A,f}(\tau )\right\}\left[{R}_{2}^{B,f}({\tau }_{2}),{R}_{2}^{B,f}({\tau }_{3})\right]\\ +\displaystyle \frac{1}{8}{\rm{i}}{\mu }^{4}{\displaystyle \int }_{{\tau }_{0}}^{\tau }{\rm{d}}{\tau }_{1}{\displaystyle \int }_{{\tau }_{0}}^{{\tau }_{1}}{\rm{d}}{\tau }_{2}{\displaystyle \int }_{{\tau }_{0}}^{{\tau }_{2}}{\rm{d}}{\tau }_{3}\\ \times \,\left\langle \beta \left|\left[{\phi }^{f}({x}_{A}({\tau }_{3})),{\phi }^{f}({x}_{B}({\tau }_{2}))\right]\left[{\phi }^{f}({x}_{B}({\tau }_{1})),{\phi }^{f}({x}_{A}(\tau ))\right]\right|\beta \right\rangle \\ \times \left\{{R}_{2}^{A,f}({\tau }_{3}),{R}_{2}^{A,f}(\tau )\right\}\left[{R}_{2}^{B,f}({\tau }_{1}),{R}_{2}^{B,f}({\tau }_{2})\right]\\ +\displaystyle \frac{1}{8}{\rm{i}}{\mu }^{4}{\displaystyle \int }_{{\tau }_{0}}^{\tau }{\rm{d}}{\tau }_{1}{\displaystyle \int }_{{\tau }_{0}}^{\tau }{\rm{d}}{\tau }_{2}{\displaystyle \int }_{{\tau }_{0}}^{{\tau }_{2}}{\rm{d}}{\tau }_{3}\\ \times \,\left\langle \beta \left|\left[{\phi }^{f}({x}_{B}({\tau }_{1})),{\phi }^{f}({x}_{A}(\tau ))\right]\left[{\phi }^{f}({x}_{B}({\tau }_{3})),{\phi }^{f}({x}_{A}({\tau }_{2}))\right]\right|\beta \right\rangle \\ \times \left[{R}_{2}^{A,f}(\tau ),{R}_{2}^{A,f}({\tau }_{2})\right]\left\{{R}_{2}^{B,f}({\tau }_{1}),{R}_{2}^{B,f}({\tau }_{3})\right\}\\ +\displaystyle \frac{1}{8}{\rm{i}}{\mu }^{4}{\displaystyle \int }_{{\tau }_{0}}^{\tau }{\rm{d}}{\tau }_{1}{\displaystyle \int }_{{\tau }_{0}}^{{\tau }_{1}}{\rm{d}}{\tau }_{2}{\displaystyle \int }_{{\tau }_{0}}^{\tau }{\rm{d}}{\tau }_{3}\\ \times \,\left\langle \beta \left|\left\{{\phi }^{f}({x}_{B}({\tau }_{2})),{\phi }^{f}({x}_{A}({\tau }_{3}))\right\}\left[{\phi }^{f}({x}_{B}({\tau }_{1})),{\phi }^{f}({x}_{A}(\tau ))\right]\right|\beta \right\rangle \\ \times \left[{R}_{2}^{A,f}(\tau ),{R}_{2}^{A,f}({\tau }_{3})\right]\left[{R}_{2}^{B,f}({\tau }_{2}),{R}_{2}^{B,f}({\tau }_{1})\right]\\ +\displaystyle \frac{1}{8}{\rm{i}}{\mu }^{4}{\displaystyle \int }_{{\tau }_{0}}^{\tau }{\rm{d}}{\tau }_{1}{\displaystyle \int }_{{\tau }_{0}}^{{\tau }_{1}}{\rm{d}}{\tau }_{2}{\displaystyle \int }_{{\tau }_{0}}^{{\tau }_{1}}{\rm{d}}{\tau }_{3}\\ \times \,\left\langle \beta \left|\left[{\phi }^{f}({x}_{B}({\tau }_{2})),{\phi }^{f}({x}_{A}({\tau }_{1}))\right]\left[{\phi }^{f}({x}_{B}({\tau }_{3})),{\phi }^{f}({x}_{A}(\tau ))\right]\right|\beta \right\rangle \\ \times \left\{{R}_{2}^{A,f}({\tau }_{1}),{R}_{2}^{A,f}(\tau )\right\}\left[{R}_{2}^{B,f}({\tau }_{3}),{R}_{2}^{B,f}({\tau }_{2})\right]\\ +\,A\rightleftharpoons B\,\mathrm{terms},\end{array}\end{eqnarray}$
where $\begin{eqnarray}\begin{array}{rcl}{\phi }^{f}(x) & = & \displaystyle \int {{\rm{d}}}^{3}{\boldsymbol{k}}\ {g}_{{\boldsymbol{k}}}\left[{a}_{{\boldsymbol{k}}}({t}_{0}){{\rm{e}}}^{-{\rm{i}}{\omega }_{{\boldsymbol{k}}}(t-{t}_{0})}{f}_{{\boldsymbol{k}}}({\boldsymbol{x}})\right.\\ & & \left.+{a}_{{\boldsymbol{k}}}^{\dagger }({t}_{0}){{\rm{e}}}^{{\rm{i}}{\omega }_{{\boldsymbol{k}}}(t-{t}_{0})}{f}_{{\boldsymbol{k}}}^{* }({\boldsymbol{x}})\right],\end{array}\end{eqnarray}$
is the free scalar-field operator, and $\begin{eqnarray}\begin{array}{rcl}{R}_{2}^{\xi ,f}(\tau ) & = & \displaystyle \frac{{\rm{i}}}{2}\left[{R}_{-}^{\xi }({\tau }_{0}){{\rm{e}}}^{-{\rm{i}}{\omega }_{\xi }(\tau -{\tau }_{0})}\right.\\ & & \left.-{R}_{+}^{\xi }({\tau }_{0}){{\rm{e}}}^{{\rm{i}}{\omega }_{\xi }(\tau -{\tau }_{0})}\right],\end{array}\end{eqnarray}$
is the free part of the atomic operator ${R}_{2}^{\xi }(\tau )$. After taking the expectation values of the two effective Hamiltonians, ${\left({H}_{S}(\tau )\right)}_{\mathrm{tf}}^{\mathrm{eff}}$ and ${\left({H}_{S}(\tau )\right)}_{\mathrm{rr}}^{\mathrm{eff}}$ in equations ( $\begin{eqnarray}\begin{array}{rcl}{\left(\delta E\right)}_{\mathrm{tf}} & = & \langle {g}_{A}{g}_{B}| {\left({H}_{S}(\tau )\right)}_{\mathrm{tf}}^{\mathrm{eff}}| {g}_{A}{g}_{B}\rangle \\ & = & 2{\rm{i}}{\mu }^{4}{\displaystyle \int }_{{\tau }_{0}}^{\tau }{\rm{d}}{\tau }_{1}{\displaystyle \int }_{{\tau }_{0}}^{{\tau }_{1}}{\rm{d}}{\tau }_{2}{\displaystyle \int }_{{\tau }_{0}}^{{\tau }_{2}}{\rm{d}}{\tau }_{3}{C}_{\beta }^{F}\\ & & \times \,({x}_{A}(\tau ),{x}_{B}({\tau }_{3})){\chi }_{\beta }^{F}({x}_{A}({\tau }_{1}),{x}_{B}({\tau }_{2}))\\ & & \times \,{\chi }^{A}(\tau ,{\tau }_{1}){\chi }^{B}({\tau }_{2},{\tau }_{3})\\ & & +A\rightleftharpoons B\,\mathrm{term},\end{array}\end{eqnarray}$
and that of the rr-contribution, $\begin{eqnarray}\begin{array}{l}{\left(\delta E\right)}_{\mathrm{rr}}\\ =\langle {g}_{A}{g}_{B}| {\left({H}_{S}(\tau )\right)}_{\mathrm{rr}}^{\mathrm{eff}}| {g}_{A}{g}_{B}\rangle \\ =2{\rm{i}}{\mu }^{4}{\displaystyle \int }_{{\tau }_{0}}^{\tau }{\rm{d}}{\tau }_{1}{\displaystyle \int }_{{\tau }_{0}}^{{\tau }_{1}}{\rm{d}}{\tau }_{2}{\displaystyle \int }_{{\tau }_{0}}^{{\tau }_{2}}{\rm{d}}{\tau }_{3}\\ \times \,{\chi }_{\beta }^{F}({x}_{A}(\tau ),{x}_{B}({\tau }_{3})){\chi }_{\beta }^{F}\left({x}_{A}({\tau }_{1}),{x}_{B}({\tau }_{2})\right)\\ \times \,{C}^{A}(\tau ,{\tau }_{1}){\chi }^{B}({\tau }_{2},{\tau }_{3})\\ +2{\rm{i}}{\mu }^{4}{\displaystyle \int }_{{\tau }_{0}}^{\tau }{\rm{d}}{\tau }_{1}{\displaystyle \int }_{{\tau }_{0}}^{{\tau }_{1}}{\rm{d}}{\tau }_{2}{\displaystyle \int }_{{\tau }_{0}}^{{\tau }_{2}}{\rm{d}}{\tau }_{3}\\ \times \,{\chi }_{\beta }^{F}({x}_{A}({\tau }_{1}),{x}_{B}({\tau }_{3})){\chi }_{\beta }^{F}\left({x}_{B}({\tau }_{2}),{x}_{A}(\tau )\right)\\ \times \,{C}^{A}(\tau ,{\tau }_{1}){\chi }^{B}({\tau }_{2},{\tau }_{3})\\ +2{\rm{i}}{\mu }^{4}{\displaystyle \int }_{{\tau }_{0}}^{\tau }{\rm{d}}{\tau }_{1}{\displaystyle \int }_{{\tau }_{0}}^{{\tau }_{1}}{\rm{d}}{\tau }_{2}{\displaystyle \int }_{{\tau }_{0}}^{{\tau }_{2}}{\rm{d}}{\tau }_{3}\\ \times \,{\chi }_{\beta }^{F}({x}_{A}({\tau }_{3}),{x}_{B}({\tau }_{2})){\chi }_{\beta }^{F}({x}_{B}({\tau }_{1}),{x}_{A}(\tau ))\\ \times \,{C}^{A}(\tau ,{\tau }_{3}){\chi }^{B}({\tau }_{1},{\tau }_{2})\\ +2{\rm{i}}{\mu }^{4}{\displaystyle \int }_{{\tau }_{0}}^{\tau }{\rm{d}}{\tau }_{1}{\displaystyle \int }_{{\tau }_{0}}^{\tau }{\rm{d}}{\tau }_{2}{\displaystyle \int }_{{\tau }_{0}}^{{\tau }_{2}}{\rm{d}}{\tau }_{3}\\ \times \,{\chi }_{\beta }^{F}({x}_{A}({\tau }_{2}),{x}_{B}({\tau }_{3})){\chi }_{\beta }^{F}({x}_{A}(\tau ),{x}_{B}({\tau }_{1})){\chi }^{A}(\tau ,{\tau }_{2})\\ \times \,{C}^{B}({\tau }_{1},{\tau }_{3})\\ +2{\rm{i}}{\mu }^{4}{\displaystyle \int }_{{\tau }_{0}}^{\tau }{\rm{d}}{\tau }_{1}{\displaystyle \int }_{{\tau }_{0}}^{{\tau }_{1}}{\rm{d}}{\tau }_{2}{\displaystyle \int }_{{\tau }_{0}}^{\tau }{\rm{d}}{\tau }_{3}\\ \times \,{C}_{\beta }^{F}({x}_{B}({\tau }_{2}),{x}_{A}({\tau }_{3})){\chi }_{\beta }^{F}({x}_{B}({\tau }_{1}),{x}_{A}(\tau ))\\ \times \,{\chi }^{A}({\tau }_{3},\tau ){\chi }^{B}({\tau }_{1},{\tau }_{2})\\ +2{\rm{i}}{\mu }^{4}{\displaystyle \int }_{{\tau }_{0}}^{\tau }{\rm{d}}{\tau }_{1}{\displaystyle \int }_{{\tau }_{0}}^{{\tau }_{1}}{\rm{d}}{\tau }_{2}{\displaystyle \int }_{{\tau }_{0}}^{{\tau }_{1}}{\rm{d}}{\tau }_{3}\\ \times \,{\chi }_{\beta }^{F}({x}_{A}(\tau ),{x}_{B}({\tau }_{3})){\chi }_{\beta }^{F}({x}_{A}({\tau }_{1}),{x}_{B}({\tau }_{2}))\\ \times \,{C}^{A}(\tau ,{\tau }_{1}){\chi }^{B}({\tau }_{3},{\tau }_{2})\\ +A\rightleftharpoons B\,\mathrm{terms},\end{array}\end{eqnarray}$
where ${C}^{\xi }(\tau ,\tau ^{\prime} )$ and ${\chi }^{\xi }(\tau ,\tau ^{\prime} )$ are the statistical functions of the atoms in the ground state, defined by $\begin{eqnarray}{C}^{\xi }(\tau ,\tau ^{\prime} )\equiv \displaystyle \frac{1}{2}\langle {g}_{\xi }| \{{R}_{2}^{\xi ,f}(\tau ),{R}_{2}^{\xi ,f}(\tau ^{\prime} )\}| {g}_{\xi }\rangle ,\end{eqnarray}$
$\begin{eqnarray}{\chi }^{\xi }(\tau ,\tau ^{\prime} )\equiv \displaystyle \frac{1}{2}\langle {g}_{\xi }| [{R}_{2}^{\xi ,f}(\tau ),{R}_{2}^{\xi ,f}(\tau ^{\prime} )]| {g}_{\xi }\rangle ,\end{eqnarray}$
and ${C}_{\beta }^{F}({x}_{A}(\tau ),{x}_{B}(\tau ^{\prime} ))$ and ${\chi }_{\beta }^{F}({x}_{A}(\tau ),{x}_{B}(\tau ^{\prime} ))$ are the correlation functions of the scalar fields in the thermal state ∣β⟩, defined by $\begin{eqnarray}{C}_{\beta }^{F}({x}_{A}(\tau ),{x}_{B}(\tau ^{\prime} ))\equiv \displaystyle \frac{1}{2}\langle \beta | \{{\phi }^{f}({x}_{A}(\tau )),{\phi }^{f}({x}_{B}(\tau ^{\prime} ))\}| \beta \rangle ,\end{eqnarray}$
$\begin{eqnarray}{\chi }_{\beta }^{F}({x}_{A}(\tau ),{x}_{B}(\tau ^{\prime} ))\equiv \displaystyle \frac{1}{2}\langle \beta | \left[{\phi }^{f}({x}_{A}(\tau )),{\phi }^{f}({x}_{B}(\tau ^{\prime} ))\right]| \beta \rangle .\end{eqnarray}$
Note here that, in order to obtain equations (10 ) and (11 ), we have simplified every scalar-field four-point correlation function of the effective Hamiltonians in equations (6 ) and (7 ), into products of two two-point correlation functions, and the details on the simplifications are given in the appendix .
3. Contributions of thermal fluctuations and radiation reaction to the interaction of two static ground-state atoms
In this section, we shall use the general formalism established in the proceeding section to calculate the interatomic interaction energy between two static ground-state atoms in a thermal bath. The trajectories of the two atoms at rest are depicted by14 ) and (15 ) are respectively
$\begin{eqnarray}{x}_{A}=(\tau ,0,0,0),\quad {x}_{B}=(\tau ,0,0,L),\end{eqnarray}$
where L (> 0) denotes the interatomic separation, and along these trajectories of the atoms, the correlation functions of the scalar fields defined in equations ( $\begin{eqnarray}\begin{array}{rcl}{C}_{\beta }^{F}({x}_{A}(\tau ),{x}_{B}(\tau ^{\prime} )) & = & \displaystyle \frac{1}{4{\pi }^{2}L}{\displaystyle \int }_{0}^{\infty }{\rm{d}}\omega \left(1+\displaystyle \frac{2}{{{\rm{e}}}^{\omega /T}-1}\right)\\ & & \times \,\sin (\omega L)\cos [\omega (\tau -\tau ^{\prime} )],\end{array}\end{eqnarray}$
$\begin{eqnarray}\begin{array}{rcl}{\chi }_{\beta }^{F}({x}_{A}(\tau ),{x}_{B}(\tau ^{\prime} )) & = & -\displaystyle \frac{{\rm{i}}}{4{\pi }^{2}L}{\displaystyle \int }_{0}^{\infty }{\rm{d}}\omega \,\sin (\omega L)\\ & & \times \,\sin [\omega (\tau -\tau ^{\prime} )].\end{array}\end{eqnarray}$
For atomic statistical functions defined in equations (12 ) and (13 ), we obtain, with equation (9 ),17 )–(20 ) into equations (10 ) and (11 ), and performing the triple integrations with respect to τ1, τ2 and τ3, we have respectively, for an infinitely long time interval τ − τ0, the tf-contribution and the rr-contribution23 ) and (24 ) has a prefactor ∼L−2, and it arises from the prefactor L−1 in the field correlation functions (equations (17 ) and (18 )) as the tf- and rr-contribution depend on the product of two field correlation functions (see equations (10 ) and (11 )). Generally, the tf- and rr-contribution are distinctive due to that the first term in equation (23 ) and equation (24 ) which is oscillatory is of the same magnitude but an opposite sign. When added up, the first terms in these two equations cancel out perfectly, and it means that this part of the effects of zero-point fluctuations in the tf- and rr-contribution is cancelled. As a result, the total interaction energy follows,
$\begin{eqnarray}{C}^{\xi }(\tau ,\tau ^{\prime} )=\displaystyle \frac{1}{4}\cos [{\omega }_{\xi }(\tau -\tau ^{\prime} )],\end{eqnarray}$
$\begin{eqnarray}{\chi }^{\xi }(\tau ,\tau ^{\prime} )=-\displaystyle \frac{{\rm{i}}}{4}\sin [{\omega }_{\xi }(\tau -\tau ^{\prime} )],\end{eqnarray}$
which are valid when the two-level atom is in its ground state ∣gξ⟩. Substituting equations ( $\begin{eqnarray}\begin{array}{rcl}{\left(\delta E\right)}_{\mathrm{tf}} & = & -\displaystyle \frac{{\mu }^{4}}{64{\pi }^{4}{L}^{2}}{\displaystyle \int }_{0}^{\infty }{\rm{d}}{\omega }_{1}{\displaystyle \int }_{0}^{\infty }{\rm{d}}{\omega }_{2}\\ & & \times \,\displaystyle \frac{{\omega }_{A}{\omega }_{B}{\omega }_{2}\sin ({\omega }_{1}L)\sin ({\omega }_{2}L)}{({\omega }_{1}^{2}-{\omega }_{A}^{2})({\omega }_{1}^{2}-{\omega }_{B}^{2})({\omega }_{2}^{2}-{\omega }_{1}^{2})}\\ & & \times \,\left(1+\displaystyle \frac{2}{{{\rm{e}}}^{{\omega }_{1}/T}-1}\right),\end{array}\end{eqnarray}$
and $\begin{eqnarray}\begin{array}{rcl}{\left(\delta E\right)}_{\mathrm{rr}} & = & {\left(\delta E\right)}_{\mathrm{tf}}-\displaystyle \frac{{\mu }^{4}}{32{\pi }^{4}{L}^{2}}{\displaystyle \int }_{0}^{\infty }{\rm{d}}{\omega }_{1}{\displaystyle \int }_{0}^{\infty }{\rm{d}}{\omega }_{2}\\ & & \times \,\displaystyle \frac{{\omega }_{1}{\omega }_{2}({\omega }_{1}^{2}-{\omega }_{A}^{2}-{\omega }_{A}{\omega }_{B}-{\omega }_{B}^{2})}{({\omega }_{A}+{\omega }_{B})({\omega }_{1}^{2}-{\omega }_{A}^{2})({\omega }_{1}^{2}-{\omega }_{B}^{2})({\omega }_{2}^{2}-{\omega }_{1}^{2})}\\ & & \times \,\sin ({\omega }_{1}L)\sin ({\omega }_{2}L).\end{array}\end{eqnarray}$
The above results can be further simplified by performing the ω2-integration with the residue theorem to yield the final expressions of the tf-contribution and the rr-contribution $\begin{eqnarray}\begin{array}{rcl}{\left(\delta E\right)}_{\mathrm{tf}} & = & -\displaystyle \frac{{\mu }^{4}}{512{\pi }^{2}{L}^{2}}\displaystyle \frac{{\omega }_{B}\cos (2{\omega }_{A}L)-{\omega }_{A}\cos (2{\omega }_{B}L)}{{\omega }_{A}^{2}-{\omega }_{B}^{2}}\\ & & -\displaystyle \frac{{\mu }^{4}}{256{\pi }^{3}{L}^{2}}{\displaystyle \int }_{0}^{\infty }{\rm{d}}u{\alpha }_{A}({\rm{i}}u){\alpha }_{B}({\rm{i}}u){{\rm{e}}}^{-2{uL}}\\ & & -\displaystyle \frac{{\mu }^{4}}{128{\pi }^{3}{L}^{2}}{\displaystyle \int }_{0}^{\infty }{\rm{d}}u{\alpha }_{A}(u){\alpha }_{B}(u)\displaystyle \frac{\sin (2{uL})}{{{\rm{e}}}^{u/T}-1}\end{array}\end{eqnarray}$
and $\begin{eqnarray}\begin{array}{rcl}{\left(\delta E\right)}_{\mathrm{rr}} & = & \displaystyle \frac{{\mu }^{4}}{512{\pi }^{2}{L}^{2}}\displaystyle \frac{{\omega }_{B}\cos (2{\omega }_{A}L)-{\omega }_{A}\cos (2{\omega }_{B}L)}{{\omega }_{A}^{2}-{\omega }_{B}^{2}}\\ & & -\displaystyle \frac{{\mu }^{4}}{256{\pi }^{3}{L}^{2}}{\displaystyle \int }_{0}^{\infty }{\rm{d}}u{\alpha }_{A}({\rm{i}}u){\alpha }_{B}({\rm{i}}u){{\rm{e}}}^{-2{uL}}\\ & & -\displaystyle \frac{{\mu }^{4}}{128{\pi }^{3}{L}^{2}}{\displaystyle \int }_{0}^{\infty }{\rm{d}}u{\alpha }_{A}(u){\alpha }_{B}(u)\displaystyle \frac{\sin (2{uL})}{{{\rm{e}}}^{u/T}-1},\end{array}\end{eqnarray}$
with ${\alpha }_{\xi }({\rm{i}}u)={\omega }_{\xi }/({\omega }_{\xi }^{2}+{u}^{2})$, both of which are composed of three terms with the first two temperature-independent terms reflecting effects of zero-point fluctuations and the third temperature-dependent term representing thermal revisions. It is also worth noting that every term of equations ( $\begin{eqnarray}\begin{array}{l}{\left(\delta E\right)}_{\mathrm{tot}}=-\displaystyle \frac{{\mu }^{4}}{128{\pi }^{3}{L}^{2}}{\displaystyle \int }_{0}^{\infty }{\rm{d}}u\displaystyle \frac{{\omega }_{A}{\omega }_{B}\ }{({\omega }_{A}^{2}+{u}^{2})({\omega }_{B}^{2}+{u}^{2})}{{\rm{e}}}^{-2{uL}}\\ -\displaystyle \frac{{\mu }^{4}}{64{\pi }^{3}{L}^{2}}{\displaystyle \int }_{0}^{\infty }{\rm{d}}u\displaystyle \frac{{\omega }_{A}{\omega }_{B}\sin (2{uL})}{({\omega }_{A}^{2}-{u}^{2})({\omega }_{B}^{2}-{u}^{2})({{\rm{e}}}^{u/T}-1)},\quad \end{array}\end{eqnarray}$
which contains two terms with the first temperature-independent term being exactly the interaction energy of two static ground-state atoms in a vacuum and the other temperature-dependent term giving the thermal corrections.It is difficult to analytically evaluate the integrals in equations (23 )–(25 ). Fortunately, however, analytical approximate results are obtainable in some special cases. For simplicity, we assume in the following discussions that the two atoms are identical with the transition frequency denoted by ω0, and divide our analysis into the van der Waals regime where the interatomic separation L is much smaller than the atomic transition wavelength λ, defined by $\lambda =2\pi {\omega }_{0}^{-1}$, and the Casimir–Polder regime where L is much larger than the atomic transition wavelength λ. As indicated in the general expressions of the tf- and rr-contribution (23 ) and (24 ), for the interaction energy of two atoms in a thermal bath, besides the two length scales L and λ, another length scale, i.e. the characteristic wavelength of photons in the thermal bath β singles out. As a result, the van der Waals region L ≪ λ and the Casimir–Polder region L ≫ λ are respectively further divided into three typical regions, and the tf- and rr-contribution (equations (23 ) and (24 )), which can be expressed in terms of some complicated special functions after doing the integrations, can be further approximated in polynomials of L. We next discuss those results one by one.
3.1. The van der Waals interaction in a thermal bath
This regime can be further divided into three sub-regimes depending on the relationship among three characteristic length scales, i.e. the interatomic separation L, the thermal wavelength of the thermal bath, β = T−1, and atomic transition wavelength λ. We now start with L ≪ λ ≪ β.
3.1.1. The van der Waals interaction in the region L ≪ λ ≪ β
When L ≪ λ ≪ β, the tf-contribution and the rr-contribution to the interatomic van der Waals interaction, according to equations (23 ) and (24 ), are respectively approximated by
$\begin{eqnarray}\begin{array}{rcl}{\left(\delta E\right)}_{\mathrm{tf}} & \simeq & \displaystyle \frac{{\mu }^{4}}{256{\pi }^{3}L}-\displaystyle \frac{{\mu }^{4}{T}^{2}}{384\pi {\omega }_{0}^{2}L}\\ & & -\displaystyle \frac{{\mu }^{4}{\omega }_{0}^{2}L}{192{\pi }^{3}}\mathrm{ln}({\omega }_{0}L),\end{array}\end{eqnarray}$
$\begin{eqnarray}\begin{array}{rcl}{\left(\delta E\right)}_{\mathrm{rr}} & \simeq & -\displaystyle \frac{{\mu }^{4}}{512{\pi }^{2}{\omega }_{0}{L}^{2}}+\displaystyle \frac{{\mu }^{4}}{256{\pi }^{3}L}\\ & & -\displaystyle \frac{{\mu }^{4}{T}^{2}}{384\pi {\omega }_{0}^{2}L}-\displaystyle \frac{{\mu }^{4}{\omega }_{0}^{2}L}{192{\pi }^{3}}\mathrm{ln}({\omega }_{0}L).\end{array}\end{eqnarray}$
In each of these two equations, the first term on the right plays a dominant role, and the thermal corrections dependent on the temperature appear in higher orders, which are negligible. Taking the derivative of the above leading terms with respect to the interatomic separation, we find that the thermal fluctuations lead to a repulsive van der Waals force between two atoms, which is proportional to L−2, while the radiation reaction of atoms yields an attractive force which is proportional to L−3. Let us note that the different scaling of L in these two leading terms is caused by the difference in the sign of oscillatory terms in equations (23 ) and (24 ). To be specific, the oscillatory term in the tf-contribution (equation (23 )) cancels the other temperature-independent term, while that in the rr-contribution (equation (24 )) coincides it and they jointly lead to the first term in equation (27 ). It is easy to see that ∣δE∣tf ≪ ∣δE∣rr, which means that the rr-contribution dominates over the tf-contribution in this case. Adding up equations (26 ) and (27 ), we get
$\begin{eqnarray}\begin{array}{rcl}{\left(\delta E\right)}_{\mathrm{tot}} & \simeq & -\displaystyle \frac{{\mu }^{4}}{512{\pi }^{2}{\omega }_{0}{L}^{2}}+\displaystyle \frac{{\mu }^{4}}{128{\pi }^{3}L}\\ & & -\displaystyle \frac{{\mu }^{4}{T}^{2}}{192\pi {\omega }_{0}^{2}L}-\displaystyle \frac{{\mu }^{4}{\omega }_{0}^{2}L}{96{\pi }^{3}}\mathrm{ln}({\omega }_{0}L),\end{array}\end{eqnarray}$
for the total van der Waals interaction energy of two ground-state atoms in a thermal bath, which is attractive.3.1.2. The van der Waals interaction in the region L ≪ β ≪ λ
When L ≪ β ≪ λ, we simplify the tf-contribution and the rr-contribution in equations (23 ) and (24 ) into the following two expressions,26 ) and (27 ), we find that the first (leading) term of the tf-contribution in equation (29 ), as well as the second (subleading) term of the rr-contribution in equation (30 ), which are independent of the temperature, remarkably change sign. This sign-changing phenomenon results from that the sole temperature-dependent term in equations (23 ) and (24 ), in the limit of L ≪ β ≪ λ also generates the temperature-independent behaviors of L−1, and thus suggests that, in sharp contrast to the low-temperature case of λ ≪ β, the thermal radiation can also induce the temperature-independent corrections to the van der Waals interaction in the high-temperature case, i.e. when the thermal wavelength β is much smaller than the transition wavelength of atoms λ. Consequently, the tf-contribution in the L ≪ β ≪ λ region leads to an attractive van der Waals force, as opposed to a repulsive force in the region L ≪ λ ≪ β, and the rr-contribution leads to a slight enhancement for the van der Waals interaction energy, rather than a slight reduction in the region L ≪ λ ≪ β.
$\begin{eqnarray}{\left(\delta E\right)}_{\mathrm{tf}}\simeq -\displaystyle \frac{{\mu }^{4}}{256{\pi }^{3}L}-\displaystyle \frac{{\mu }^{4}{\omega }_{0}^{2}}{384{\pi }^{3}{T}^{2}L}\mathrm{ln}\left(\displaystyle \frac{{\omega }_{0}}{T}\right),\end{eqnarray}$
$\begin{eqnarray}\begin{array}{rcl}{\left(\delta E\right)}_{\mathrm{rr}} & \simeq & -\displaystyle \frac{{\mu }^{4}}{512{\pi }^{2}{\omega }_{0}{L}^{2}}-\displaystyle \frac{{\mu }^{4}}{256{\pi }^{3}L}\\ & & -\displaystyle \frac{{\mu }^{4}{\omega }_{0}^{2}}{384{\pi }^{3}{T}^{2}L}\mathrm{ln}\left(\displaystyle \frac{{\omega }_{0}}{T}\right).\end{array}\end{eqnarray}$
Comparing the above results with their counterparts in the case of L ≪ λ ≪ β, i.e. equations (The summation of the above two equations, equations (29 ) and (30 ), gives rise to the total van der Waals interaction energy in the region L ≪ β ≪ λ,30 ), which however is unaffected by the temperature. As a result, the total interatomic van der Waals force in this case is, to the leading order, still an attractive force proportional to L−3. In the present region and the region of L ≪ λ ≪ β, the contribution of radiation reaction of atoms to the interaction energy is always larger than that of field fluctuations, and this is physically understandable as the propagation distance of the photons of the radiative field emitted from one atom to the other atom is very short when two atoms are very close.
$\begin{eqnarray}\begin{array}{rcl}{\left(\delta E\right)}_{\mathrm{tot}} & \simeq & -\displaystyle \frac{{\mu }^{4}}{512{\pi }^{2}{\omega }_{0}{L}^{2}}-\displaystyle \frac{{\mu }^{4}}{128{\pi }^{3}L}\\ & & -\displaystyle \frac{{\mu }^{4}{\omega }_{0}^{2}}{192{\pi }^{3}{T}^{2}L}\mathrm{ln}\left(\displaystyle \frac{{\omega }_{0}}{T}\right).\end{array}\end{eqnarray}$
Obviously, though the dominant behavior of the tf-contribution is altered by the temperature, the leading term of the total interaction energy comes from the rr-contribution in equation (3.1.3. The van der Waals interaction in the region β ≪ L ≪ λ
The third region for the van der Waals interaction we are to examine is β ≪ L ≪ λ, within which the tf-contribution and the rr-contribution in equations (23 ) and (24 ) are approximated as
$\begin{eqnarray}{\left(\delta E\right)}_{\mathrm{tf}}\simeq -\displaystyle \frac{{\mu }^{4}}{256{\pi }^{3}L}-\displaystyle \frac{{\mu }^{4}{\omega }_{0}^{2}{L}^{2}T}{384{\pi }^{2}},\end{eqnarray}$
$\begin{eqnarray}\begin{array}{rcl}{\left(\delta E\right)}_{\mathrm{rr}} & \simeq & -\displaystyle \frac{{\mu }^{4}}{512{\pi }^{2}{\omega }_{0}{L}^{2}}-\displaystyle \frac{{\mu }^{4}}{256{\pi }^{3}L}\\ & & -\displaystyle \frac{{\mu }^{4}{\omega }_{0}^{2}{L}^{2}T}{384{\pi }^{2}},\end{array}\end{eqnarray}$
and they jointly result in the following total van der Waals interaction energy, $\begin{eqnarray}\begin{array}{rcl}{\left(\delta E\right)}_{\mathrm{tot}} & \simeq & -\displaystyle \frac{{\mu }^{4}}{512{\pi }^{2}{\omega }_{0}{L}^{2}}-\displaystyle \frac{{\mu }^{4}}{128{\pi }^{3}L}\\ & & -\displaystyle \frac{{\mu }^{4}{\omega }_{0}^{2}{L}^{2}T}{192{\pi }^{2}}.\end{array}\end{eqnarray}$
In order to elucidate the behaviors of the tf-contribution and the rr-contribution to the interatomic van der Waals interaction energy in the β ≪ L ≪ λ region, we further divide this region into the following three different subregions, i.e. the subregion of $\beta \ll L\ll \sqrt[3]{{\lambda }^{2}\beta }$, $\sqrt[3]{{\lambda }^{2}\beta }\ll L\ll \sqrt[4]{{\lambda }^{3}\beta }$, and $\sqrt[4]{{\lambda }^{3}\beta }\ll L\ll \lambda $.In the first subregion $\beta \ll L\ll \sqrt[3]{{\lambda }^{2}\beta }$, we obtain29 ) and (30 )]. And the tf-contribution and the rr-contribution, in the second subregion $\sqrt[3]{{\lambda }^{2}\beta }\ll L\ll \sqrt[4]{{\lambda }^{3}\beta }$, can be expressed as
$\begin{eqnarray}{\left(\delta E\right)}_{\mathrm{tf}}\simeq -\displaystyle \frac{{\mu }^{4}}{256{\pi }^{3}L},\end{eqnarray}$
$\begin{eqnarray}{\left(\delta E\right)}_{\mathrm{rr}}\simeq -\displaystyle \frac{{\mu }^{4}}{512{\pi }^{2}{\omega }_{0}{L}^{2}},\end{eqnarray}$
which are exactly the same in the leading order as those in the region L ≪ β ≪ λ [see equations ( $\begin{eqnarray}{\left(\delta E\right)}_{\mathrm{tf}}\simeq -\displaystyle \frac{{\mu }^{4}{\omega }_{0}^{2}{L}^{2}T}{384{\pi }^{2}},\end{eqnarray}$
$\begin{eqnarray}{\left(\delta E\right)}_{\mathrm{rr}}\simeq -\displaystyle \frac{{\mu }^{4}}{512{\pi }^{2}{\omega }_{0}{L}^{2}},\end{eqnarray}$
respectively. Here the van der Waals force attributed to the tf-contribution is a repulsive force, with a new behavior of TL. It is easy to verify that $| {\left(\delta E\right)}_{\mathrm{rr}}| \gg | {\left(\delta E\right)}_{\mathrm{tf}}| $ in the first two subregions, i.e. the rr-contribution weighs much more than the tf-contribution. The third subregion comes with $\sqrt[4]{{\lambda }^{3}\beta }\ll L\ll \lambda $, and now the tf-contribution plays an equally important role as the rr-contribution does, i.e. $\begin{eqnarray}{\left(\delta E\right)}_{\mathrm{tf}}\simeq {\left(\delta E\right)}_{\mathrm{rr}}\simeq -\displaystyle \frac{{\mu }^{4}{\omega }_{0}^{2}{L}^{2}T}{384{\pi }^{2}}.\end{eqnarray}$
It is interesting to note that the interatomic van der Waals energy leads to an attractive force in the first two subregions which scales as L−3, and remarkably a repulsive force in the third subregion, i.e. the $\sqrt[4]{{\lambda }^{3}\beta }\ll L\ll \lambda $ region, which scales as TL. The crossover of the van der Waals force from attractive to repulsive when the interatomic separation enters the third subregion is a result of the fact that both the tf- and rr-contribution are dominated by the temperature-dependent terms which weigh equally, whereas, in the first two regions L ≪ λ ≪β, L ≪ β ≪ λ, or the first two subregions of β ≪ L ≪ λ, although the tf-contributions may be dominated by temperature-dependent terms, the rr-contribution, which is often dominated by temperature-independent terms, is much greater than the tf-contribution, and consequently the total interaction energy is still dominated by temperature-independent terms, resulting in an overall attractive force.Some comments are in order on whether or not this repulsive van der Waals force is experimentally detectable. For the remarkable crossover of the van der Waals force from attractive to repulsive to occur, β ≪ λ has to be satisfied, i.e. the thermal-bath temperature T has to be much higher than the critical temperature Tc ≡ ℏω0/kB (in the International System of Units). For natural atoms with transition frequencies ω0 ∈ [1014, 1017]s−1, Tc ∼ 764 K at least, which is much higher than the room temperature T ∼ 300 K. Thus this repulsive van der Waals force of two natural atoms does not show up near the room temperature. One may use the atoms with relatively lower transition frequencies to reduce the temperature required for a repulsive van der Waals force to happen. However, a smaller ω0 and a smaller T means a larger interatomic separation since $L\gg \sqrt[4]{{\lambda }^{3}\beta }$ is required. For instance, it can be estimated that at room temperature, if a repulsive van der Waals force is to be generated between two atoms with a transition frequency ω0 ∼ 1012 s−1, the interatomic separation has to be approximately hundreds of micrometers, which makes the force too tiny to detect at present. So, experimental detection of the repulsive van der Waals force relies on technological advances in precision of measurement in the future.
3.2. The Casimir–Polder interaction in a thermal bath
In this subsection, we turn to explore how the interatomic Casimir–Polder interaction energy is affected by the presence of a thermal bath. Similar to the discussions on the van der Waals interaction energy above, our analytical results for the interatomic Casimir–Polder interaction energy are approximately obtained in the regions of λ ≪ L ≪ β, λ ≪ β ≪ L, and β ≪ λ ≪ L.
3.2.1. The Casimir–Polder interaction in the region λ ≪ L ≪ β
In the region λ ≪ L ≪ β, we approximate the tf-contribution and the rr-contribution to the interatomic Casimir–Polder interaction energy, equations (23 ) and (24 ), as40 ) and that of equation (41 ), which are temperature-independent, oscillate with the interatomic separation and they respectively correspond to the first term of the general expression of the tf-contribution equation (23 ) and that of the rr-contribution equation (24 ) in the limit of ωξ → ω0 and L ≫ λ. The second term is also temperature-independent, and it can be derived by replacing αξ(iu) in the second term of equations (23 ) and (24 ) with αξ(0) which is reasonable as the exponential e−2uL in the integrand and the condition λ ≪ L in the present case ensure the rapid decrease of the integrand with the increase of u. The oscillatory term together with the second term in the above results are manifestations of the effects of zero-point fluctuations. The third terms which are proportional to T2L−1, are thermal revisions and they result from the third terms of equations (23 ) and (24 ). Comparing the amplitude of the oscillatory term in each of the above two equations with non-oscillatory terms, we find that the former is much greater, indicating that the tf- and rr-contribution to the interaction energy in the region λ ≪ L ≪ β can be either positive or negative and even be null, signaling an either attractive or repulsive and even vanishing force.
$\begin{eqnarray}\begin{array}{rcl}{\left(\delta E\right)}_{\mathrm{tf}} & \simeq & \displaystyle \frac{{\mu }^{4}}{512{\pi }^{2}L}\sin (2{\omega }_{0}L+{\theta }_{1})\\ & & -\displaystyle \frac{{\mu }^{4}}{512{\pi }^{3}{\omega }_{0}^{2}{L}^{3}}-\displaystyle \frac{{\mu }^{4}{T}^{2}}{384\pi {\omega }_{0}^{2}L},\end{array}\end{eqnarray}$
$\begin{eqnarray}\begin{array}{rcl}{\left(\delta E\right)}_{\mathrm{rr}} & \simeq & -\displaystyle \frac{{\mu }^{4}}{512{\pi }^{2}L}\sin (2{\omega }_{0}L+{\theta }_{1})\\ & & -\displaystyle \frac{{\mu }^{4}}{512{\pi }^{3}{\omega }_{0}^{2}{L}^{3}}-\displaystyle \frac{{\mu }^{4}{T}^{2}}{384\pi {\omega }_{0}^{2}L},\end{array}\end{eqnarray}$
with ${\theta }_{1}=\arcsin {\left(1+4{\omega }_{0}^{2}{L}^{2}\right)}^{-1/2}$. In sharp contrast to those results in the three van der Waals regions previously discussed, here the first term of equation (However, this oscillatory separation-dependent behavior does not remain in the total Casimir–Polder interaction energy, as the signs of the oscillatory terms in ${\left(\delta E\right)}_{\mathrm{tf}}$ and ${\left(\delta E\right)}_{\mathrm{rr}}$ are opposite. As a result, we obtain the following monotonic behavior of the interatomic Casimir–Polder interaction energy in the region λ ≪ L ≪ β,40 ) and (41 ). Here, the change from the van der Waals behavior of L−2 in equation (28 ) to the Casimir–Polder behavior of L−3 in equation (42 ) shows the retardation effect as a result of finite light speed.
$\begin{eqnarray}{\left(\delta E\right)}_{\mathrm{tot}}\simeq -\displaystyle \frac{{\mu }^{4}}{256{\pi }^{3}{\omega }_{0}^{2}{L}^{3}}-\displaystyle \frac{{\mu }^{4}{T}^{2}}{192\pi {\omega }_{0}^{2}L},\end{eqnarray}$
which are equally contributed by the tf-contribution and the rr-contribution, equations (3.2.2. The Casimir–Polder interaction in the region λ ≪ β ≪ L
When λ ≪ β ≪ L, we obtain, with equations (23 ) and (24 ), the tf-contribution and rr-contribution to the Casimir–Polder interaction energy,23 ) and (24 ) since an exponential factor ${\left({{\rm{e}}}^{u/T}-1\right)}^{-1}$ is included in the integrand and meanwhile ω0 ≫ T for the present case. As compared with the temperature-independent behavior of L−3 in the vacuum case or the thermal case with λ ≪ L ≪ β, this Casimir–Polder interaction energy is dominated by thermal corrections and it displays the same separation-dependence, L−2, just as the van der Waals interaction energy in the vacuum case, suggesting that temperature effects in this case suppress the retardation effect.
$\begin{eqnarray}\begin{array}{rcl}{\left(\delta E\right)}_{\mathrm{tf}} & \simeq & \displaystyle \frac{{\mu }^{4}}{512{\pi }^{2}L}\sin (2{\omega }_{0}L+{\theta }_{1})\\ & & -\displaystyle \frac{{\mu }^{4}T}{256{\pi }^{2}{\omega }_{0}^{2}{L}^{2}},\end{array}\end{eqnarray}$
$\begin{eqnarray}\begin{array}{rcl}{\left(\delta E\right)}_{\mathrm{rr}} & \simeq & -\displaystyle \frac{{\mu }^{4}}{512{\pi }^{2}L}\sin (2{\omega }_{0}L+{\theta }_{1})\\ & & -\displaystyle \frac{{\mu }^{4}T}{256{\pi }^{2}{\omega }_{0}^{2}{L}^{2}}.\end{array}\end{eqnarray}$
We find that, the same oscillatory terms as those in the region λ ≪ L ≪ β, also appear in the present case, and their amplitudes are also much larger than the monotonic temperature-dependent terms. The summation of the above tf-contribution and rr-contribution, as a result of the perfect cancellation of the oscillatory terms, yields $\begin{eqnarray}{\left(\delta E\right)}_{\mathrm{tot}}\simeq -\displaystyle \frac{{\mu }^{4}T}{128{\pi }^{2}{\omega }_{0}^{2}{L}^{2}},\end{eqnarray}$
which gives the total Casimir–Polder interaction energy of two ground-state atoms in the λ ≪ β ≪ L region. Obviously, this interaction energy equally comes from the tf- and rr-contribution, and it can be easily derived by replacing αξ(u) with αξ(0) in the integrand of equations (3.2.3. The Casimir–Polder interaction in the region β ≪ λ ≪ L
Finally, in the β ≪ λ ≪ L region, the tf-contribution and the rr-contribution to the interatomic interaction energy are respectively43 ) and (44 )], lies in that the leading oscillatory terms of the tf-contribution and the rr-contribution are temperature-dependent and they do not cancel out with each other, leading to the total interatomic Casimir–Polder interaction which is oscillatory,
$\begin{eqnarray}\begin{array}{rcl}{\left(\delta E\right)}_{\mathrm{tf}} & \simeq & {\left(\delta E\right)}_{\mathrm{rr}}\simeq \displaystyle \frac{{\mu }^{4}T}{256{\pi }^{2}{\omega }_{0}L}\sin (2{\omega }_{0}L+{\theta }_{2})\\ & & -\displaystyle \frac{{\mu }^{4}T}{256{\pi }^{2}{\omega }_{0}^{2}{L}^{2}},\end{array}\end{eqnarray}$
with ${\theta }_{2}=\arcsin {\left(1+{\omega }_{0}^{2}{L}^{2}\right)}^{-1/2}$. A remarkable difference, as compared to the case for λ ≪ β [refer to equations ( $\begin{eqnarray}\begin{array}{rcl}{\left(\delta E\right)}_{\mathrm{tot}} & \simeq & \displaystyle \frac{{\mu }^{4}T}{128{\pi }^{2}{\omega }_{0}L}\sin (2{\omega }_{0}L+{\theta }_{2})\\ & & -\displaystyle \frac{{\mu }^{4}T}{128{\pi }^{2}{\omega }_{0}^{2}{L}^{2}}.\end{array}\end{eqnarray}$
So, the total interatomic Casimir–Polder interaction oscillates with a decreasing amplitude which is proportional to TL−1, and as a result, the Casimir–Polder force can be either attractive or repulsive, and even vanish, depending on the concrete value of the separation. For some special interatomic separations, in which the leading oscillatory term vanishes, the Casimir–Polder interaction energy is proportional to TL−2 and corresponds to an attractive force, as in the λ ≪ β ≪ L region [see equation (45 )].
So far, we have derived the interaction energy of two ground-state atoms in a thermal bath and analysed in detail the approximate analytical results in three van der Waals regions and three Casimir–Polder regions. It is worth pointing out that for two atoms with an arbitrary separation in a thermal bath at an arbitrary temperature, though the interaction energy equations (23 )–(25 ) are generally difficult to simplify, they can be numerically computed after the removal of the removable singularity u = ωξ. We have done such numerical computations, compared the numerical results with all the approximate analytical results in six typical regions we have reported in the present work and found that they agree very well.
4. Summary
In this paper, we first make a generalization of the fourth-order DDC formalism from vacuum to a thermal bath, and with this general formalism, we then study the interatomic interaction energy between two static ground-state atoms coupled to massless scalar fields at finite temperature T. We analytically calculate and analyze, in six regions depending on the relative magnitudes among three characteristic physical parameters of the system, i.e. the interatomic separation L, the atomic transition wavelength λ, and the thermal wavelength β, the separate contribution of thermal fluctuations (tf-contribution) and that of the radiation reaction of atoms (rr-contribution) to the interatomic van der Waals and Casimir–Polder interaction energies, and some important results are summed up as follows.
The analytical results for the van der Waals interaction energy are approximately obtained in the following three regions: L ≪ λ ≪ β, L ≪ β ≪ λ and β ≪ L ≪ λ, and the third region β ≪ L ≪ λ is further divided into three subregions: $\beta \ll L\ll \sqrt[3]{{\lambda }^{2}\beta }$, $\sqrt[3]{{\lambda }^{2}\beta }\ll L\ll \sqrt[4]{{\lambda }^{3}\beta }$, and $\sqrt[4]{{\lambda }^{3}\beta }\ll L\ll \lambda $. We find that, in the first two regions as well as the first two subregions of the third region, the thermal corrections to the van der Waals interaction are of a high-order effect, and the interatomic interaction in the leading order is the same as that in a vacuum, where the rr-contribution plays a dominant role and leads to an attractive force with the behavior of L−3. Remarkably, a crossover of the van der Waals force from attractive to repulsive occurs, when the temperature of the thermal bath is sufficiently high such that $\sqrt[4]{{\lambda }^{3}\beta }\ll L\ll \lambda $. In this subregion within the β ≪ L ≪ λ region, the tf-contribution and the rr-contribution play an equally important role, and they jointly result in a repulsive van der Waals force which scales as TL.
Similarly, analytical results for the Casimir–Polder interaction energy are obtained in regions λ ≪ L ≪ β, λ ≪ β ≪ L and β ≪ λ ≪ L, and the tf-contribution and the rr-contribution in every region are found to be oscillatory and are equally important. However, these oscillatory separation-dependent behaviors do not carry onto the total interaction energy in the first two regions because of a perfect cancellation of the two oscillatory contributions, resulting in an attractive interatomic interaction force with a monotonic separation-dependence. Concretely, the force scales as L−4 in the first region, which is the same as its counterpart in a vacuum, and scales as TL−3 in the second region, which is dominated by the thermal corrections. In the third region, the total interaction energy oscillates with a decreasing amplitude proportional to TL−1 and gives rise to a force that can be either attractive or repulsive and even vanishing, depending on the concrete value of the interatomic separation.
Appendix A. Details on the simplifications of the four-point correlation functions of the scalar fields in equations (6 ) and (7 )
As we have mentioned below in equation (15 ), in order to obtain equations (10 ) and (11 ), we have simplified the four-point correlation functions of the free fields in the effective Hamiltonians of equations (6 ) and (7 ), into the products of two two-point correlation functions of the scalar fields. We now give the details on the simplifications, i.e. how ⟨β∣φf(x1)φf(x2)φf(x3)φf(x4)∣β⟩ can be expanded in terms of the product of two two-point correlation functions.
The four-point correlation function of the scalar fields in a thermal state, ⟨β∣φf(x1)φf(x2)φf(x3)φf(x4)∣β⟩, is generally given by8 ). After further simplifications, we find that there are only six nonvanishing terms left on the right of this equation, the coefficients of which are respectively given by
$\begin{eqnarray}\begin{array}{l}\langle \beta | {\phi }^{f}({x}_{1}){\phi }^{f}({x}_{2}){\phi }^{f}({x}_{3}){\phi }^{f}({x}_{4})| \beta \rangle \\ =\displaystyle \frac{1}{Z}\sum _{{{\boldsymbol{k}}}_{1},{{\boldsymbol{k}}}_{2},{{\boldsymbol{k}}}_{3},{{\boldsymbol{k}}}_{4}}\mathrm{Tr}\left\{{{\rm{e}}}^{-\beta {H}_{F}}\left[{a}_{{{\boldsymbol{k}}}_{1}}^{f}{f}_{{{\boldsymbol{k}}}_{1}}({x}_{1})+{\rm{H}}.\ {\rm{c}}.\right]\right.\\ \times \,\left[{a}_{{{\boldsymbol{k}}}_{2}}^{f}{f}_{{{\boldsymbol{k}}}_{2}}({x}_{2})+{\rm{H}}.\ {\rm{c}}.\right]\\ \left.\times \,\left[{a}_{{{\bf{k}}}_{3}}^{f}{f}_{{{\bf{k}}}_{3}}({x}_{3})+{\rm{H}}.\ {\rm{c}}.\right]\left[{a}_{{{\boldsymbol{k}}}_{4}}^{f}{f}_{{{\boldsymbol{k}}}_{4}}({x}_{4})+{\rm{H}}.{\rm{c}}.\right]\right\},\end{array}\end{eqnarray}$
where ${H}_{F}={\sum }_{{\boldsymbol{k}}}{\omega }_{{\boldsymbol{k}}}{a}_{{\boldsymbol{k}}}^{f\dagger }{a}_{{\boldsymbol{k}}}^{f}$ is the free Hamiltonian of the scalar fields, $Z\equiv \mathrm{Tr}({{\rm{e}}}^{-\beta {H}_{F}})={\prod }_{{\boldsymbol{k}}}{\left(1-{{\rm{e}}}^{-\beta {\omega }_{{\boldsymbol{k}}}}\right)}^{-1}$, and on the right of this equation, we have used the expansion of the free fields, equation ( $\begin{eqnarray}\left\{\begin{array}{l}\displaystyle \frac{1}{Z}\mathrm{Tr}({{\rm{e}}}^{-\beta {H}_{F}}{a}_{{{\boldsymbol{k}}}_{1}}^{\dagger }{a}_{{{\boldsymbol{k}}}_{2}}^{\dagger }{a}_{{{\boldsymbol{k}}}_{3}}{a}_{{{\boldsymbol{k}}}_{4}})\\ ={N}_{{{\boldsymbol{k}}}_{1}}{N}_{{{\boldsymbol{k}}}_{2}}({\delta }_{{{\boldsymbol{k}}}_{1}{{\boldsymbol{k}}}_{3}}{\delta }_{{{\boldsymbol{k}}}_{2}{{\boldsymbol{k}}}_{4}}+{\delta }_{{{\boldsymbol{k}}}_{1}{{\boldsymbol{k}}}_{4}}{\delta }_{{{\boldsymbol{k}}}_{2}{{\boldsymbol{k}}}_{3}}),\\ \displaystyle \frac{1}{Z}\mathrm{Tr}({{\rm{e}}}^{-\beta {H}_{F}}{a}_{{{\boldsymbol{k}}}_{1}}^{\dagger }{a}_{{{\boldsymbol{k}}}_{2}}{a}_{{{\boldsymbol{k}}}_{3}}^{\dagger }{a}_{{{\boldsymbol{k}}}_{4}})\\ ={N}_{{{\boldsymbol{k}}}_{1}}{N}_{{{\boldsymbol{k}}}_{3}}{\delta }_{{{\boldsymbol{k}}}_{1}{{\boldsymbol{k}}}_{2}}{\delta }_{{{\boldsymbol{k}}}_{3}{{\boldsymbol{k}}}_{4}}+{N}_{{{\boldsymbol{k}}}_{1}}({N}_{{{\boldsymbol{k}}}_{3}}+1){\delta }_{{{\boldsymbol{k}}}_{1}{{\boldsymbol{k}}}_{4}}{\delta }_{{{\boldsymbol{k}}}_{2}{{\boldsymbol{k}}}_{3}},\\ \displaystyle \frac{1}{Z}\mathrm{Tr}({{\rm{e}}}^{-\beta {H}_{F}}{a}_{{{\boldsymbol{k}}}_{1}}^{\dagger }{a}_{{{\boldsymbol{k}}}_{2}}{a}_{{{\boldsymbol{k}}}_{3}}{a}_{{{\boldsymbol{k}}}_{4}}^{\dagger })\\ ={N}_{{{\boldsymbol{k}}}_{1}}({N}_{{{\boldsymbol{k}}}_{4}}+1)({\delta }_{{{\boldsymbol{k}}}_{1}{{\boldsymbol{k}}}_{2}}{\delta }_{{{\boldsymbol{k}}}_{3}{{\boldsymbol{k}}}_{4}}+{\delta }_{{{\boldsymbol{k}}}_{1}{{\boldsymbol{k}}}_{3}}{\delta }_{{{\boldsymbol{k}}}_{2}{{\boldsymbol{k}}}_{4}}),\\ \displaystyle \frac{1}{Z}\mathrm{Tr}({{\rm{e}}}^{-\beta {H}_{F}}{a}_{{{\boldsymbol{k}}}_{1}}{a}_{{{\boldsymbol{k}}}_{2}}^{\dagger }{a}_{{{\boldsymbol{k}}}_{3}}^{\dagger }{a}_{{{\boldsymbol{k}}}_{4}})\\ ={N}_{{{\boldsymbol{k}}}_{2}}({N}_{{{\boldsymbol{k}}}_{3}}+1){\delta }_{{{\boldsymbol{k}}}_{1}{{\boldsymbol{k}}}_{3}}{\delta }_{{{\boldsymbol{k}}}_{2}{{\boldsymbol{k}}}_{4}}+({N}_{{{\boldsymbol{k}}}_{2}}+1){N}_{{{\boldsymbol{k}}}_{3}}{\delta }_{{{\boldsymbol{k}}}_{1}{{\boldsymbol{k}}}_{2}}{\delta }_{{{\boldsymbol{k}}}_{3}{{\boldsymbol{k}}}_{4}},\\ \displaystyle \frac{1}{Z}\mathrm{Tr}({{\rm{e}}}^{-\beta {H}_{F}}{a}_{{{\boldsymbol{k}}}_{1}}{a}_{{{\boldsymbol{k}}}_{2}}^{\dagger }{a}_{{{\boldsymbol{k}}}_{3}}{a}_{{{\boldsymbol{k}}}_{4}}^{\dagger })\\ =({N}_{{{\boldsymbol{k}}}_{2}}{N}_{{{\boldsymbol{k}}}_{4}}+{N}_{{{\boldsymbol{k}}}_{2}}+{N}_{{{\boldsymbol{k}}}_{4}}+1){\delta }_{{{\boldsymbol{k}}}_{1}{{\boldsymbol{k}}}_{2}}{\delta }_{{{\boldsymbol{k}}}_{3}{{\boldsymbol{k}}}_{4}}+{N}_{{{\boldsymbol{k}}}_{2}}({N}_{{{\boldsymbol{k}}}_{4}}+1){\delta }_{{{\boldsymbol{k}}}_{1}{{\boldsymbol{k}}}_{4}}{\delta }_{{{\boldsymbol{k}}}_{2}{{\boldsymbol{k}}}_{3}},\\ \displaystyle \frac{1}{Z}\mathrm{Tr}({{\rm{e}}}^{-\beta {H}_{F}}{a}_{{{\boldsymbol{k}}}_{1}}{a}_{{{\boldsymbol{k}}}_{2}}{a}_{{{\boldsymbol{k}}}_{3}}^{\dagger }{a}_{{{\boldsymbol{k}}}_{4}}^{\dagger })\\ =({N}_{{{\boldsymbol{k}}}_{3}}{N}_{{{\boldsymbol{k}}}_{4}}+{N}_{{{\boldsymbol{k}}}_{3}}+{N}_{{{\boldsymbol{k}}}_{4}}+1)({\delta }_{{{\boldsymbol{k}}}_{1}{{\boldsymbol{k}}}_{3}}{\delta }_{{{\boldsymbol{k}}}_{2}{{\boldsymbol{k}}}_{4}}+{\delta }_{{{\boldsymbol{k}}}_{1}{{\boldsymbol{k}}}_{4}}{\delta }_{{{\boldsymbol{k}}}_{2}{{\boldsymbol{k}}}_{3}}),\end{array}\right.\end{eqnarray}$
with ${N}_{{{\boldsymbol{k}}}_{i}}={\left({{\rm{e}}}^{\beta {\omega }_{{{\boldsymbol{k}}}_{i}}}-1\right)}^{-1}$. Hereafter, we omit the subscripts ‘f' of the creation and annihilation operators for brevity.Notice that for the two-point correlation function of the scalar fields in the same thermal state,
$\begin{eqnarray}\begin{array}{rcl}\langle \beta | \phi ({x}_{1})\phi ({x}_{2})| \beta \rangle & \equiv & \displaystyle \frac{1}{Z}\mathrm{Tr}\left[{{\rm{e}}}^{-\beta {H}_{F}}\phi ({x}_{1})\phi ({x}_{2})\right]\\ & = & \displaystyle \frac{1}{Z}\sum _{{{\boldsymbol{k}}}_{1},{{\boldsymbol{k}}}_{2}}\mathrm{Tr}\left\{{{\rm{e}}}^{-\beta {H}_{F}}\left[{a}_{{{\boldsymbol{k}}}_{1}}{f}_{{{\boldsymbol{k}}}_{1}}({x}_{1})+{\rm{H}}.\ {\rm{c}}.\right]\right.\\ & & \left.\times \,\left[{a}_{{{\boldsymbol{k}}}_{2}}{f}_{{{\boldsymbol{k}}}_{2}}({x}_{2})+{\rm{H}}.\ {\rm{c}}.\right]\right\}\\ & = & \displaystyle \sum _{{\boldsymbol{k}}}\left[{N}_{{\boldsymbol{k}}}{f}_{{\boldsymbol{k}}}^{* }({x}_{1}){f}_{{\boldsymbol{k}}}({x}_{2})\right.\\ & & \left.+({N}_{{\boldsymbol{k}}}+1){f}_{{\boldsymbol{k}}}({x}_{1}){f}_{{\boldsymbol{k}}}^{* }({x}_{2})\right].\end{array}\end{eqnarray}$
Here to obtain the third line of the above equation, we have used the following relations: $\begin{eqnarray}\left\{\begin{array}{l}\mathrm{Tr}({{\rm{e}}}^{-\beta {H}_{F}}{a}_{{{\boldsymbol{k}}}_{1}}{a}_{{{\boldsymbol{k}}}_{2}})=\mathrm{Tr}({{\rm{e}}}^{-\beta {H}_{F}}{a}_{{{\boldsymbol{k}}}_{1}}^{\dagger }{a}_{{{\boldsymbol{k}}}_{2}}^{\dagger })=0,\\ \mathrm{Tr}({{\rm{e}}}^{-\beta {H}_{F}}{a}_{{{\boldsymbol{k}}}_{1}}^{\dagger }{a}_{{{\boldsymbol{k}}}_{2}})={\delta }_{{{\boldsymbol{k}}}_{1}{{\boldsymbol{k}}}_{2}}\displaystyle \frac{1}{{{\rm{e}}}^{\beta {\omega }_{{{\boldsymbol{k}}}_{1}}}-1}\prod _{{\boldsymbol{k}}}\displaystyle \frac{1}{1-{{\rm{e}}}^{-\beta {\omega }_{{\boldsymbol{k}}}}},\\ \mathrm{Tr}({{\rm{e}}}^{-\beta {H}_{F}}{a}_{{{\boldsymbol{k}}}_{1}}{a}_{{{\boldsymbol{k}}}_{2}}^{\dagger })={\delta }_{{{\boldsymbol{k}}}_{1}{{\boldsymbol{k}}}_{2}}\left(1+\displaystyle \frac{1}{{{\rm{e}}}^{\beta {\omega }_{{{\boldsymbol{k}}}_{1}}}-1}\right)\prod _{{\boldsymbol{k}}}\displaystyle \frac{1}{1-{{\rm{e}}}^{-\beta {\omega }_{{\boldsymbol{k}}}}}.\end{array}\right.\end{eqnarray}$
Combining equations (A2 )–(A4 ) with equation (A1 ), we can obtain6 ) and (7 ), and take the expectation values of them over the ground state ∣gAgB⟩, and we finally get equations (10 ) and (11 ).
$\begin{eqnarray}\begin{array}{l}\langle \beta \left|\phi ({x}_{1})\phi ({x}_{2})\phi ({x}_{3})\phi ({x}_{4})\right|\beta \rangle \\ =\langle \beta \left|\phi ({x}_{1})\phi ({x}_{2})\right|\beta \rangle \langle \beta \left|\phi ({x}_{3})\phi ({x}_{4})\right|\beta \rangle \\ +\langle \beta \left|\phi ({x}_{1})\phi ({x}_{3})\right|\beta \rangle \langle \beta \left|\phi ({x}_{2})\phi ({x}_{4})\right|\beta \rangle \\ +\langle \beta \left|\phi ({x}_{1})\phi ({x}_{4})\right|\beta \rangle \langle \beta \left|\phi ({x}_{2})\phi ({x}_{3})\right|\beta \rangle .\end{array}\end{eqnarray}$
Then use this relation in equations (