text:
# Vacuum Polarization of Massive Scalar Fields on the Black Hole Horizon

## I Introduction

Since Hawking’s discovery of black holes radiation , it has been an important subject to investigate quantum effects in black hole spacetimes and a great deal of effort has been done. The expectation values of the regularized vacuum polarization $`\varphi ^2`$ and the stress energy tensor $`T_{\mu \nu }`$ are important quantities to understand quantum effects in black hole spacetimes. They represent the effect of quantum fluctuations and back reaction on black hole spacetimes and play an important role in the context of semi-classical theory of quantum gravity.

For a scalar field in the large mass limit, of which Compton length is much smaller than the characteristic curvature scale of background geometry, we can evaluate regularized $`\varphi ^2`$ and $`T_{\mu \nu }`$ by using the DeWitt-Schwinger approximation. This is an adiabatic expansion in terms of inverse power series of the mass of scalar field. One can evaluate divergent and finite parts of $`\varphi ^2`$ and $`T_{\mu \nu }`$ up to arbitrary order. Frolov and Zelnikov  suggest that the finite parts of $`\varphi ^2`$ and $`T_{\mu \nu }`$ calculated by the DeWitt-Schwinger expansion coincides with the value calculated by assuming the Hartle-Hawking state for the scalar field, although the DeWitt-Schwinger expansion does not impose any condition on the state of the scalar field. The justification is confirmed numerically near the horizon.

The Hartle-Hawking state is defined by requiring that the Euclidean Green’s function is regular on the horizon. Its physical meaning is that a black hole is placed in a cavity and is equilibrium with thermal black-body radiation. The thermal state is defined by the Euclidean Green’s function by imposing periodicity in the imaginary time coordinate with a period of the inverse of the temperature. On the order hand, in the DeWitt-Schwinger expansion, the condition which characterizes the thermal state is not imposed. So it is not obvious whether the DeWitt-Schwinger expansion yields the same finite part of the Hartle-Hawking state near the horizon.

In this paper, we consider scalar fields in static spherically symmetric background spacetimes. By assuming that the mass of scalar fields is large, we calculate $`\varphi ^2`$ in the thermal state with arbitrary temperature using point splitting method near the horizon. A similar calculation is already done by Anderson et al., but their expression of the vacuum polarization diverges on the horizon which contradicts with their numerical calculation. Our purpose is to derive the vacuum polarization which is regular on the horizon, and resolve why the finite part of $`\varphi ^2`$ in the Hartle-Hawking state near the horizon can be approximated by the DeWitt-Schwinger expansion which does not assume any condition on the thermal state.

The plan of this paper is as follows. In Sec.II, we review the method to calculate $`\varphi ^2`$ in a thermal state using point splitting regularization. In Sec.III, we present our method to calculate $`\varphi ^2`$ in the large mass limit near the black hole horizon. Sec.IV is devoted to conclusion and discussion. Throughout the paper, we use units such that $`\mathrm{}=c=G=k_B=1`$. Our sign conventions are those of Misner, Thorne and Wheeler .

## II derivation of vacuum polarization $`\varphi ^2`$ in a thermal state

We consider a minimally coupled scalar field with mass $`m`$ in a thermal state at arbitrary temperature $`T`$ in a Reissner-Nordström spacetime. The metric in Euclidean section $`\tau =it`$ is given by

$$ds^2=fd\tau ^2+\frac{1}{f}dr^2+r^2d\mathrm{\Omega }^2,$$
(1)
where $`f=(rr_+)(rr_{})/r^2`$, $`r_+`$ and $`r_{}`$ are the location of inner and outer horizon, respectively. The surface gravity of a black hole is given by $`\kappa =(r_+r_{})/(2r_+^2)`$. The vacuum polarization $`\varphi ^2`$ can be computed from the Euclidean Green’s function by noting the identity

$$\varphi ^2_{\text{unreg}}=\frac{1}{2}G^{(1)}(x,x)=iG_F(x,x)=G_E(x,x),$$
(2)
where $`G^{(1)}`$ is the Hadamard Green’s function, $`G_F`$ is the Feynman Green’s function, and $`G_E`$ is the Euclidean Green’s function, respectively. $`G_E`$ is divergent at coincident limit and must be regularized. In this paper, the covariant point-splitting method is employed for the regularization of ultraviolet divergences. We start from an expression for $`G_E(x,x^{})`$ when the points $`x`$ and $`x^{}`$ are separated. Next we prepare an appropriate counterterm to subtract the divergence, and then take coincident limit. This procedure is shown schematically as

$$\varphi ^2_{\text{reg}}=\underset{x^{}x}{lim}\left[G_E(x,x^{})G_{DS}(x,x^{})\right],$$
(3)
where $`G_{DS}(x,x^{})`$ is a point splitting counterterm needed to regularize $`G_E(x,x^{})`$. The Euclidean Green’s function obeys the equation

$$\left[\mathrm{}m^2\right]G_E(x,x^{})=\frac{\delta ^4(x,x^{})}{g^{1/2}(x)},$$
(4)
where $`m`$ is the mass of the scalar field and d’Alembertian $`\mathrm{}`$ evaluated using the Euclidean metric (1). The point splitting counter term needed to renormalize $`G_E(x,x^{})`$ in arbitrary spacetime is

$$\begin{array}{cc}\hfill G_{DS}(x,x^{})& =\frac{1}{8\pi ^2\sigma }+\frac{m^2+(\xi \frac{1}{6})R}{8\pi ^2}\left[\gamma +\frac{1}{2}\mathrm{ln}\left[\frac{m^2\sigma }{2}\right]\right]\hfill \\ & \frac{m^2}{16\pi ^2}+\frac{1}{96\pi ^2}R_{\alpha \beta }\frac{\sigma ^{,\alpha }\sigma ^{,\beta }}{\sigma },\hfill \end{array}$$
(5)
where $`\xi `$ is coupling constant to the scalar curvature $`R`$, $`\sigma `$ is the one half the square of geodesic distance between the points $`x`$ and $`x^{}`$ and $`\gamma `$ is Euler’s constant.

For simplicity of calculation, the quantity $`\varphi ^2`$ is evaluated by timelike point separation $`\tau ^{}\tau `$. By assuming that the Green’s function is periodic in $`ϵ=\tau \tau ^{}`$ with period $`1/T`$, the expression for the vacuum polarization in the thermal state at temperature $`T`$ is given by

$`\varphi ^2_T`$ $`=\underset{\tau ^{}\tau }{lim}\left[G_E(x,\tau ;x,\tau ^{})G_{DS}(x,\tau ;x,\tau ^{})\right]`$
$`=\underset{\tau ^{}\tau }{lim}[{\displaystyle \frac{T}{4\pi }}{\displaystyle \underset{l=0}{\overset{\mathrm{}}{}}}[(2l+1)p_{0l}q_{0l}{\displaystyle \frac{1}{r\sqrt{f}}}]`$
$`+{\displaystyle \frac{T}{2\pi }}{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}\mathrm{cos}[n2\pi T(\tau \tau ^{})]{\displaystyle \underset{l=0}{\overset{\mathrm{}}{}}}[(2l+1)p_{nl}q_{nl}{\displaystyle \frac{1}{r\sqrt{f}}}]G_{DS}(x,\tau ;x,\tau ^{})]`$
$`\varphi ^2_{n=0}+\varphi ^2_{n1}\varphi ^2_{DS},`$ (6)
where

$$\varphi ^2_{DS}=\frac{m^2}{8\pi ^2}\left(\mathrm{ln}\frac{m\sqrt{f}ϵ}{2}+\gamma \frac{1}{2}\right)\frac{1}{4\pi ^2fϵ^2}+\left(\frac{f^{}^2}{192\pi ^2f}\frac{f^{\prime \prime }}{96\pi ^2}\frac{f^{}}{48\pi ^2r}\right).$$
(7)
The term $`1/(r\sqrt{f})`$ is subtracted to cancel the superficial divergence which comes from the choice of timelike point-splitting. $`p_{nl}`$ and $`q_{nl}`$ are independent solutions of the mode equation

$`{\displaystyle \frac{d^2\chi _{nl}}{dr^2}}+\left[{\displaystyle \frac{2}{r}}+{\displaystyle \frac{1}{f}}{\displaystyle \frac{df}{dr}}\right]{\displaystyle \frac{d\chi _{nl}}{dr}}\left[{\displaystyle \frac{(2\pi nT)^2}{f^2}}+{\displaystyle \frac{l(l+1)}{r^2f}}+{\displaystyle \frac{m^2}{f}}\right]\chi _{nl}=0`$ (8)
and satisfy the Wronskian condition

$$p_{nl}\frac{dq_{nl}}{dr}q_{nl}\frac{dp_{nl}}{dr}=\frac{1}{r^2f}.$$
(9)
Anderson et al. calculated the second and the third term of (II) analytically using the second order WKB approximation . Their result is

$$\begin{array}{cc}\hfill \varphi ^2_{n1}\varphi ^2_{DS}& =\frac{m^2}{16\pi ^2}\mathrm{ln}\left(\frac{m^2f}{16\pi ^2T^2}\right)+\frac{m^2}{16\pi ^2}\hfill \\ & +\frac{T^2}{12f}\frac{f_{}^{}{}_{}{}^{2}}{192\pi ^2f}\frac{f^{\prime \prime }}{96\pi ^2}\frac{f^{}}{48\pi ^2r}.\hfill \end{array}$$
(10)
The expression (10) was obtained by taking into account only the contribution from $`n1`$ mode. In the case of massless scalar fields, excluding the lowest frequency mode $`n=0`$ corresponds to imposing an infrared cutoff. The expression (10) gives good approximation for massless fields both near and far from the horizon in Schwarzschild and Reissner-Nordström spacetimes. This was confirmed by their numerical calculation.

When the result (10) is applied to the massive scalar field, however, two unpleasant features appear. Firstly, the regularized expression always has logarithmic divergence on the horizon $`f=0`$ even if the scalar field is in the Hartle-Hawking state. It is not consistent with their numerical work which shows regular behavior of $`\varphi ^2_T`$ on the horizon. Secondly, the finite terms of $`O(m^2)`$ and $`O(m^0)`$ are left in the regularized expression. This disagrees with the result of the DeWitt-Schwinger approximation by Frolov which does not contain such finite terms. Therefore the expression (10) cannot yield good approximation for the massive scalar field. In the next section, we present the method to improve their approximations and resolve these problems.

## III vacuum polarization near the horizon in the large mass limit

In this section, we present the method to resolve the unwanted behavior arose in the expression of $`\varphi ^2_T`$ near the horizon $`f0`$. We assume that the mass $`m`$ of the scalar field is large:

$$m\frac{1}{r_+}.$$
(11)
We evaluate $`n1`$ and $`n=0`$ contribution to $`\varphi ^2_T`$ separately.

### A $`n1`$ contribution

We can use the WKB approximation for $`n1`$ mode. The zeroth order WKB solution is given by

$`p_n={\displaystyle \frac{1}{\sqrt{2r^2W}}}\mathrm{exp}\left[{\displaystyle 𝑑r\frac{W}{f}}\right],q_n={\displaystyle \frac{1}{\sqrt{2r^2W}}}\mathrm{exp}\left[{\displaystyle 𝑑r\frac{W}{f}}\right],`$
$`;W^2=\kappa ^2n^2+m^2f+{\displaystyle \frac{\left(l+\frac{1}{2}\right)^2}{r^2}}f=(\alpha _n^2+{\displaystyle \frac{\left(l+\frac{1}{2}\right)^2}{r^2}})f,`$ (12)
$`\alpha _n^2=m^2+{\displaystyle \frac{\kappa ^2n^2}{f}},\kappa =2\pi T.`$
The WKB approximation is correct near the horizon. Using this solution, we have

$$\varphi ^2_{n1}=\frac{T}{2\pi }\underset{n=1}{\overset{\mathrm{}}{}}\mathrm{cos}(n\kappa ϵ)I_n,$$
(13)
where

$`I_n`$ $`{\displaystyle \underset{l=0}{\overset{\mathrm{}}{}}}\left[(2l+1)p_nq_n{\displaystyle \frac{1}{r\sqrt{f}}}\right]`$
$`={\displaystyle \frac{1}{r\sqrt{f}}}{\displaystyle \underset{l=0}{\overset{\mathrm{}}{}}}\left\{{\displaystyle \frac{l+\frac{1}{2}}{\left[r^2\alpha _n^2+\left(l+\frac{1}{2}\right)^2\right]^{1/2}}}1\right\}.`$ (14)
$`I_n`$ can be evaluated using the Plana sum formula under the assumption $`f0,\alpha _n\kappa n/f`$:

$`{\displaystyle \underset{l=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{l+\frac{1}{2}}{\left[r^2\alpha _n^2+\left(l+\frac{1}{2}\right)\right]^{1/2}}}`$ $`={\displaystyle \frac{1}{2}}{\displaystyle \frac{\frac{1}{2}}{\left[r^2\alpha _n^2+\frac{1}{4}\right]^{1/2}}}+{\displaystyle _0^{\mathrm{}}}𝑑l{\displaystyle \frac{l+\frac{1}{2}}{\left[r^2\alpha _n^2+\left(l+\frac{1}{2}\right)^2\right]^{1/2}}}`$
$`+i{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{dl}{e^{2\pi l}1}}\left\{{\displaystyle \frac{il+\frac{1}{2}}{\left[r^2\alpha _n^2+\left(il+\frac{1}{2}\right)^2\right]^{1/2}}}(ll)\right\}`$
$`=\left[\left(l+{\displaystyle \frac{1}{2}}\right)^2+r^2\alpha _n^2\right]^{1/2}|_{l=0}^{\mathrm{}}+O\left({\displaystyle \frac{f}{r\kappa n}}\right).`$ (15)
Therefore, for $`n1`$,

$`I_n`$ $`={\displaystyle \frac{1}{r\sqrt{f}}}\left\{\left[{\displaystyle \frac{1}{4}}+r^2\left(m^2+{\displaystyle \frac{r^2\kappa ^2n^2}{f^2}}\right)\right]^{1/2}+O\left({\displaystyle \frac{f}{r\kappa n}}\right)\right\}`$
$`={\displaystyle \frac{\kappa n}{f}}{\displaystyle \frac{m^2}{2\kappa n}}+O({\displaystyle \frac{f}{rn^3}},{\displaystyle \frac{\sqrt{f}}{r^2\kappa n}}).`$ (16)
By taking $`n`$ sum, $`O(f/(rn^3),\sqrt{f}/(r^2\kappa n))`$ terms give finite values which becomes zero on the horizon and we can omit these terms. Using the formula

$$\begin{array}{cc}\hfill \kappa \underset{n=1}{\overset{\mathrm{}}{}}\frac{\mathrm{cos}(n\kappa ϵ)}{n\kappa }& =\frac{1}{2}\mathrm{ln}(\kappa ^2ϵ^2)+O(ϵ^2),\hfill \\ \hfill \kappa \underset{n=1}{\overset{\mathrm{}}{}}n\kappa \mathrm{cos}(n\kappa ϵ)& =\frac{1}{ϵ^2}\frac{\kappa ^2}{12}+O(ϵ^2),\hfill \end{array}$$
(17)
we have

$`\varphi ^2_{n1}`$ $`{\displaystyle \frac{\kappa ^2}{4\pi ^2f}}{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}n\mathrm{cos}(n\kappa ϵ){\displaystyle \frac{m^2}{8\pi ^2}}{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{\mathrm{cos}(n\kappa ϵ)}{n}}`$
$`={\displaystyle \frac{1}{4\pi ^2fϵ^2}}+{\displaystyle \frac{\kappa ^2}{48\pi ^2f}}+{\displaystyle \frac{m^2}{16\pi ^2}}\mathrm{ln}(kappa^2ϵ^2).`$ (18)
Using the DeWitt-Schwinger counter term (7), the regularized expression of the expectation value becomes

$$\varphi ^2_{n1}\varphi ^2_{DS}=\frac{m^2}{16\pi ^2}\mathrm{ln}\left(\frac{m^2f}{16\pi ^2T^2}\right)+\frac{m^2}{16\pi ^2}+\frac{T^2}{12f}\frac{f_{}^{}{}_{}{}^{2}}{192\pi ^2f}\frac{f^{\prime \prime }}{96\pi ^2}\frac{f^{}}{48\pi ^2r}.$$
(19)
This gives the same expression of Anderson et al., who used the second order WKB solution. Near the black hole horizon, it is sufficient to use the zeroth order WKB solution to reproduce the result of the second order WKB approximation. On the horizon $`f=0`$, this expression diverge and we must examine the contribution from $`n=0`$ mode.

### B $`n=0`$ contribution

We must recall that the expression (10) does not include $`n=0`$ mode but there is no reason to exclude this contribution in the case of the massive scalar field. So we must investigate the contribution of $`n=0`$ mode. Since the WKB method breaks down for $`n=0`$ mode near the horizon, we cannot apply the approximation used by Anderson et al.. We solve $`n=0`$ mode function by the following method.

As we assume that the mass of the scalar field is large, a dimensionless constant

$$ϵ\frac{r_+r_{}}{4r_+[m^2r_+^2+l(l+1)]}$$
(20)
becomes smaller than unity for all value of $`l`$:

$$ϵ1(\text{for all }l)$$
(21)
We rescale the radial coordinate $`r`$ as follows:

$$x\sqrt{\frac{rr_+}{ϵr_+}}.$$
(22)
Assuming that $`x`$ is $`O(1)`$ quantity, this rescaling of the radial coordinate means we are concentrate on the region $`(rr_+)O(ϵ)`$, which is the region near the black hole horizon. Using $`x`$ and $`ϵ`$, the radial equation (8) for $`n=0`$ mode can be written as

$`{\displaystyle \frac{d^2\chi _{0l}}{dx^2}}+\left[{\displaystyle \frac{1}{x}}+F(x,ϵ)\right]{\displaystyle \frac{d\chi _{0l}}{dx}}\left[1+G(x,ϵ)\right]\chi _{0l}=0,`$ (23)
where

$`F(x,ϵ)`$ $`={\displaystyle \frac{2ϵx}{2\kappa r_++ϵx^2}},`$ (24)
$`G(x,ϵ)`$ $`={\displaystyle \frac{16ϵ}{2\kappa r_++ϵx^2}}\left[\kappa r_++2m^2r_+^2(2ϵ^2x^2+3ϵ^3x^4)\right].`$ (25)
We can express the approximate solution of the mode equation (23) by the power series expansion with respect to a small parameter $`ϵ`$:

$$\chi _{0l}(x,\eta )=\chi _{0l}^{(0)}(x)+ϵ\chi _{0l}^{(1)}(x)+ϵ^2\chi _{0l}^{(2)}(x)+\mathrm{}.$$
(26)
$`\chi _{0l}^{(n)}`$ obeys

$`{\displaystyle \frac{d\chi _{0l}^{(0)}}{dx^2}}+{\displaystyle \frac{1}{x}}{\displaystyle \frac{d\chi _{0l}^{(0)}}{dx}}\chi _{0l}^{(0)}=0,`$ (27)
$`{\displaystyle \frac{d\chi _0^{(n)}}{dx^2}}+{\displaystyle \frac{1}{x}}{\displaystyle \frac{d\chi _{0l}^{(n)}}{dx}}\chi _{0l}^{(n)}=`$
$`{\displaystyle \underset{k=1}{\overset{n}{}}}{\displaystyle \frac{1}{k!}}\left[{\displaystyle \frac{^kF}{ϵ^k}}|_{ϵ=0}{\displaystyle \frac{d\chi _{0l}^{(nk)}}{dx}}+{\displaystyle \frac{^kG}{ϵ^k}}|_{ϵ=0}\chi _{0l}^{(nk)}\right](\text{for }n1)`$ (28)
We can obtain solutions of these equations by using modified Bessel functions. For example, the lowest order solution is given by

$$p_{l0}^{(0)}=I_0(x),q_{l0}^{(0)}=\frac{2}{r_+r_{}}K_0(x).$$
(29)
Substituting the expansion (26) to

$$\varphi ^2_{n=0}=\frac{T}{4\pi }\underset{l=0}{\overset{\mathrm{}}{}}\left[(2l+1)p_{l0}q_{l0}\frac{1}{r\sqrt{f}}\right]$$
(30)
and taking summation about $`l`$ using the Plana sum formula, $`\varphi ^2_{n=0}`$ is written as the power series expansion with respect to $`(mr_+)^1`$. We calculated $`\varphi ^2_{n=0}`$ up to $`O\left((mr_+)^4\right)`$, which is accomplished by the $`ϵ`$ expansion of the mode function up to $`O(ϵ^3)`$. By taking into the account of the contribution of $`n=0`$ mode, we finally get the form of the vacuum polarization near the horizon as follows:

$`\varphi ^2_T`$ $`=\left(T{\displaystyle \frac{\kappa }{2\pi }}\right)\left\{{\displaystyle \frac{m^2}{16\pi \kappa }}+{\displaystyle \frac{m^2}{8\pi \kappa }}\left[2\gamma +\mathrm{ln}{\displaystyle \frac{m^2(rr_+)}{2\pi T}}\right]{\displaystyle \frac{1}{12\pi r_+}}\right\}`$
$`+\left\{T^2\left({\displaystyle \frac{\kappa }{2\pi }}\right)^2\right\}\left[{\displaystyle \frac{1}{24\kappa (rr_+)}}{\displaystyle \frac{r_{}}{48r_+^3\kappa ^2}}\right]`$ (31)
$`+{\displaystyle \frac{m^2}{16\pi ^2\kappa }}(\kappa +2\pi T)\mathrm{ln}\left({\displaystyle \frac{2\pi T}{\kappa }}\right)`$
$`+{\displaystyle \frac{T}{4\pi }}{\displaystyle \frac{1}{45m^2}}{\displaystyle \frac{1}{(r_+r_{})r_+^2}}\left[3{\displaystyle \frac{6r_{}}{r_+}}+{\displaystyle \frac{4r_{}^2}{r_+^2}}\right]`$
$`+{\displaystyle \frac{T}{4\pi }}{\displaystyle \frac{2}{315m^4}}{\displaystyle \frac{1}{(r_+r_{})r_+^4}}\left[10+39{\displaystyle \frac{r_{}}{r_+}}52{\displaystyle \frac{r_{}^2}{r_+^2}}+23{\displaystyle \frac{r_{}^3}{r_+^3}}\right]+O\left((mr_+)^6\right).`$ (32)
The divergences in $`\varphi ^2_T`$ appears only in the terms of $`O(m^2)`$ and $`O(m^0)`$. Taking $`f0`$ limit, $`O(m^2)`$ term diverge as $`\mathrm{ln}(rr_+)`$ and $`O(m^0)`$ as $`(rr_+)^1`$. These divergence disappear completely if the temperature of the scalar field equals $`T=\kappa /2\pi T_H`$, where $`T_H`$ is the Hawking temperature. Specifying $`T=T_H`$ is equivalent to require the regurality of $`\varphi ^2`$ on the horizon. In this case, the black hole is equilibrium with the thermal scalar field with temperature $`T=T_H`$. This is nothing but a Hartle-Hawking state.

We can evaluate a finite part of $`\varphi ^2_T`$ in the Hartle-Hawking state on the horizon. The leading order of $`\varphi ^2_T`$ at $`T=T_H`$ is $`O(m^2)`$ because not only divergent but also finite terms contained in $`O(m^2)`$ and $`O(m^0)`$ terms vanish. The result is given by

$$\begin{array}{cc}\hfill \varphi ^2_{T=T_H}& =\frac{1}{720\pi ^2m^2}\frac{1}{r_+^4}\left[3\frac{6r_{}}{r_+}+\frac{4r_{}^2}{r_+^2}\right]\hfill \\ & +\frac{1}{2520\pi ^2m^4}\frac{1}{r_+^6}\left[10+39\frac{r_{}}{r_+}52\frac{r_{}^2}{r_+^2}+23\frac{r_{}^3}{r_+^3}\right]+O\left(m^6\right).\hfill \end{array}$$
(33)
This expression agrees with the result of the DeWitt-Schwinger expansion .

## IV Conclusion and discussion

We calculated $`\varphi ^2_T`$ of a massive scalar field near the horizon of a Reissner-Nordström black hole using point splitting method. Our results are as follows: (i)all divergent and finite terms of $`O(m^2)`$ and $`O(m^0)`$ in $`\varphi ^2_T`$ are canceled to be zero if the temperature of the scalar field equals to the Hawking temperature $`T_H`$ of a black hole. (ii)at this temperature, the leading order of the regularized $`\varphi ^2_T`$ becomes $`O(m^2)`$ and its value on the horizon agrees with that of the DeWitt-Schwinger expansion up to $`O(m^4)`$. (i) and (ii) indicate that regularized $`\varphi ^2_T`$ in the Hartle-Hawking state near the black hole horizon is well approximated by the DeWitt-Schwinger expansion. This confirms the results of .

It is not trivial that the finite part of $`\varphi ^2_T`$ in the thermal state at $`T=T_H`$ is reproduced by the DeWitt-Schwinger expansion which does not specify any thermal state of the scalar field. It was shown that the leading order of the DeWitt-Schwinger expansion reproduces the value $`\varphi ^2_{T=0}`$ of scalar field with large mass in general static spherical symmetric spacetime . This result is natural because the scalar field with $`T=0`$ is in a vacuum state and the DeWitt-Schwinger expansion does not incorporate information on a thermal state.

Discrepancy between $`\varphi ^2_T`$ in the Hartle-Hawking state and that in he DeWitt-Schwinger expansion will appear far from the horizon, because the thermal effect dominates the gravitational effect in the asymptotic flat region. If the thermal effect is dominant or comparable to the contribution from curvature effect near the horizon, there would also be discrepancy between $`\varphi ^2_T`$ in the Hartle-Hawking state and that of the DeWitt-Schwinger expansion. However this is not seen from our results. This implies that the finite part of $`\varphi ^2_T`$ in the Hartle-Hawking state near the black hole horizon is dominated by the contribution from curvature effect. Why the contribution from thermal effects to finite parts can be negligible? The characteristic scale for thermal state is the inverse of the temperature $`1/T_HM`$ which is a scale length of thermal fluctuation. On the other hand, the Compton length of the scalar field is $`1/m`$. The contribution of the Hawking radiation to $`\varphi ^2`$ is $`T_H^4/m^2`$ and that of thermal excitation is $`\mathrm{exp}(m/T_H)`$. For $`m1/M`$, the effect of thermal excitation is exponentially small. Furthermore, near the black hole horizon, almost all part of thermally excited particle is absorbed by the black hole. This is the reason of the DeWitt-Schwinger approximation gives the same result of Hartle-Hawking state on the horizon.


format: text

text:
# Causality, delocalization and positivity of energy.

## 1 Introduction

Are there deviations from Einstein’s causality? G. C. Hegerfeldt has written : “Positivity of the Hamiltonian alone is used to show that particles, if initially localized in a finite region, immediately develop infinite tails.” This seems to imply superluminality. One of his examples is the Fermi problem  of two atoms coupled by a radiation field. Consider the initial condition when one of the atoms is in an excited state, the other in the ground state, and no photons are present. The probability to find the second atom in an excited state is non-vanishing immediately after the initial moment, independently of the distance between the atoms ,  \- . Hegerfeldt’s arguments are based on the analyticity of the expectation values of the operator $`N(V)`$, which gives the probability to find a particle inside a finite volume $`V`$. He showed that a state in a quantum system with positive energy localized in a finite volume $`V`$ at the instant $`t=0`$, will develop infinite tails immediately afterwards. Positivity of energy plays an essential role in his proof. In this paper we present an illustration of Hegerfeldt’s theorem, without any appeal to superluminality. We apply Hegerfeldt’s consideration to wave packets. Moreover, we show that Hegerfeldt’s effect appears even for classical fields, if wave packets are constructed from positive frequencies (corresponding positive energy quantum fields).

We first study the positive-frequency solutions of the classical wave equation (section 2). We consider wave packets $`\mathrm{\Phi }(x,t)`$ localized at $`t=0`$. The localization is due to interference of the two complex solutions, each propagating causally

$$\mathrm{\Phi }(x,t)=\psi (xt)+\psi ^{}(x+t)$$
(1)
where “” denotes complex conjugation and we take $`c=1`$. We show that both these wave packets are delocalized. They present long tails, extending to arbitrary distances and decaying according to a power law. As we shall show, the “nonlocal effect” can also be understood from the point of view of the initial conditions. Indeed, in our construction of the solution (1) we shall use two conditions; one is the initial condition of the local shape of the field $`\mathrm{\Phi }(x,t)`$ and the other is the condition of the positivity of frequencies. The frequency positivity replaces the usual initial condition on the time derivative of the field $`\mathrm{\Phi }/t`$. We shall show that our initial condition with positive frequencies leads to the nonlocality of $`\mathrm{\Phi }/t`$ at $`t=0`$.

In section 3 we show that similar conclusions are obtained for the wave packet of a free field in relativistic quantum field theory. We construct an operator reminding of the position operator of Newton-Wigner . The expectation value of this operator with the state corresponding to our wave packet is local at $`t=0`$. However, it has infinite tails which are “hidden” at time $`t=0`$, but emerge immediately afterwards. We may call this effect a “curtain” effect. No superluminal propagation is involved. We note that at the same time other quantities such as the energy density have a nonlocal expectation value in the same state even for $`t=0`$. It should be also pointed out that for the Dirac equation there are no positive energy solutions which can be localized in a finite region (see ). This demonstrates that localization in relativistic quantum field theory cannot be “complete”.

## 2 Classical wave packets

Consider classical wave packets constructed by the solutions of the wave equation with positive frequency and localized at time $`t=0`$. We show that these wave packets will spread immediately over the whole space. Curiously we have not found any reference to this effect in the literature. We start from the wave equation on the real line ($`c=1`$).

$$\left(\frac{^2}{t^2}\frac{^2}{x^2}\right)\mathrm{\Phi }(x,t)=0$$
(2)
The general complex  solution of Eq. (2) is, by the Fourier transform, of the form

$$\mathrm{\Phi }(x,t)=\frac{1}{2\pi }_{\mathrm{}}^{\mathrm{}}𝑑k\left\{\varphi _+(k)e^{i\omega _kt}+\varphi _{}(k)e^{i\omega _kt}\right\}e^{ikx}$$
(3)
where $`\omega _k=|k|`$ and where $`\varphi _\pm (k)`$ are arbitrary functions. To determine $`\varphi _+`$ and $`\varphi _{}`$ one can use the two initial conditions $`\mathrm{\Phi }(x,0)`$ and $`\dot{\mathrm{\Phi }}(x,0)`$. However, one can also consider the special class of positive-frequency solutions to Eq. (2), i.e. $`\varphi _{}(k)0`$ and

$$\mathrm{\Phi }_+(x,t)=\frac{1}{2\pi }_{\mathrm{}}^{\mathrm{}}𝑑k\varphi _+(k)e^{i\omega _kt}e^{ikx}$$
(4)
These positive-frequency solutions are determined by the initial condition $`\mathrm{\Phi }(x,0)`$. Note that relation (4) leads to a complex field for $`t0`$, even if $`\mathrm{\Phi }_+(x,0)`$ or $`\varphi _+(k)`$ are real. Consider as an example a localized (rectangular) wave packet with center $`x_0`$, and width $`2b`$ at time $`t=0`$:

$$\mathrm{\Phi }_{x_0,b}(0)=\frac{1}{2b}\mathrm{\Theta }(b|xx_0|)$$
(5)
The normalization has been chosen so that the integral of this function over $`x`$ is equal to one. Then the function $`\varphi (k)`$ is:

$$\varphi _+(k)=\frac{1}{2b}_{\mathrm{}}^+\mathrm{}𝑑xe^{ikx}\mathrm{\Theta }(b|xx_0|)$$
(6)
where $`\mathrm{\Theta }(x)`$ is the step function, which is $`0`$ for $`x`$ negative, and $`1`$ for $`x`$ positive. Then the function $`\mathrm{\Phi }(x,t)`$ in (4) is given by

$$\mathrm{\Phi }(x,t)=\frac{1}{4\pi b}_{\mathrm{}}^+\mathrm{}𝑑k_{x_0b}^{x_0+b}𝑑x^{}e^{i|k|t+ik(xx^{})}$$
(7)
This is a sum of two functions corresponding to two wave packets moving in opposite directions,

$$\mathrm{\Phi }(x,t)=\psi (xt)+\psi ^{}(x+t)$$
(8)
where

$$\psi (x)=\frac{1}{4\pi b}_{x_0b}^{x_0+b}𝑑x^{}_0^+\mathrm{}𝑑ke^{ik(xx^{})}$$
(9)
To evaluate the integral over $`k`$ we introduce the usual regularisation by adding a positive infinitesimal to $`x`$, which leads to

$$\psi (x)=\frac{1}{4\pi bi}_{x_0b}^{x_0+b}\frac{dx^{}}{xx^{}+i0}$$
(10)
After integration over $`x^{}`$ we obtain:

$$\psi (x)=\frac{i}{4\pi b}[\mathrm{ln}(xx_0+b+i0)\mathrm{ln}(xx_0b+i0)]$$
(11)
The logarithm of a complex number is given by

$$\mathrm{ln}(z)=\mathrm{ln}|z|+i(\mathrm{arg}(z)+2\pi n)$$
(12)
where $`n`$ is an integer. In order to have both terms in (11) on the same branch of the logarithm we take $`n=0`$ for both of them (due to the difference of the two terms in (11) the result does not depend on the particular value of $`n`$). The argument of $`x+i0`$ can be expressed as

$$\mathrm{arg}(x+i0)=\frac{\pi }{2}(1\mathrm{sign}(x))$$
(13)
where $`\mathrm{sign}(x)=x/|x|`$ is the sign of $`x`$. Then, inserting (12) and (13) into (11) we obtain

$$\psi (x)=\frac{1}{8b}\left(\mathrm{sign}(xx_0+b)\mathrm{sign}(xx_0b)\right)+\frac{i}{4\pi b}\mathrm{ln}\left|\frac{xx_0+b}{xx_0b}\right|$$
(14)
We see that the function $`\psi (x)`$ in (14) consists of a local real part ($`\mathrm{sign}`$) and a nonlocal imaginary part ($`\mathrm{log}`$). For $`t0`$ it is sufficient to replace $`x`$ by $`xt`$ in (8). Similar result is obtained for $`\psi ^{}(x+t)`$. As a result, the function $`\mathrm{\Phi }(x,t)`$ is also nonlocal because it is the superposition of the two complex functions $`\psi (xt)`$ and $`\psi ^{}(x+t)`$ in (8), which describe nonlocal objects moving with the speed of light in opposite directions. However, at $`t=0`$ the imaginary parts cancel each other (see Fig. 1), and we recover our localized initial condition (5), because only the real parts of these functions, which are local, remain. In all our figures time $`t`$ is measured in seconds ($`s`$), the coordinate $`x`$ is measured in “light seconds” ($`ls`$) and wave packet amplitudes are dimensionless. Fig. 2 corresponds to $`t=0.25s`$. At this moment the real local parts of $`\psi (xt)`$ and $`\psi ^{}(x+t)`$ have moved in opposite directions. The nonlocal imaginary parts of $`\psi (xt)`$ and $`\psi ^{}(x+t)`$ have also shifted in opposite directions and no more cancel each other. At this time, the two waves overlap and we have

$$|\mathrm{\Phi }(x,t)||\psi (xt)+\psi ^{}(x+t)||\psi (xt)|+|\psi ^{}(x+t)|$$
(15)
At $`t=1s`$ in Fig. 3 the overlapping is small and we have

$$|\mathrm{\Phi }(x,t)||\psi (xt)|+|\psi ^{}(x+t)|$$
(16)
We see that the initial condition $`\mathrm{\Phi }(x,0)`$ is local (Fig. 1) only because at $`t=0`$ the nonlocal parts cancel each other completely by destructive interference. We may describe the appearance of nonlocality as a sort of “curtain effect”. The nonlocal nature of each wave packet $`\psi (xt)`$ and $`\psi ^{}(x+t)`$ is hidden behind a “curtain” at the initial time and emerges immediately afterwards. Each of the nonlocal wave packets is complex and propagates at the speed of light.

In conclusion, we have illustrated Hegerfeldt’s theorem for classical wave packets. We see that the localization of wave packets corresponding to positive frequency is unstable and involves complex space structures.

Note that the localization of a positive frequency wave packet is not “complete”, because the time derivative of the function $`\mathrm{\Phi }(x,t)`$ is nonlocal even at $`t=0`$:

$$\left[\frac{\mathrm{\Phi }(x,t)}{t}\right]_{t=0}=\frac{i}{2\pi b}\left(\frac{1}{xx_0b}\frac{1}{xx_0+b}\right)$$
(17)
The wave equation being of second order demands two initial conditions: for the function itself and for its time derivative. The additional requirement of positivity of frequency replaces the second condition. There are no wave packets containing only positive frequency modes, which are localized together with their time derivative .

## 3 Relativistic quantum field

We turn now to relativistic quantum field theory. We show that the previous discussion is applicable to relativistic quantum particles. In this case the condition $`\omega _k>0`$ appears naturally since the energy $`E=\mathrm{}\omega _k`$ must be positive (we take $`\mathrm{}=1`$). We consider massless particles with no spin (“photons”). To simplify our consideration we use again a 1+1-dimensional spacetime. In terms of second quantization we have the scalar field operator

$$\widehat{\psi }(x,t)=_{\mathrm{}}^+\mathrm{}𝑑\stackrel{~}{k}\left(a_k^{}e^{i(\omega _ktkx)}+a_ke^{i(\omega _ktkx)}\right)$$
(18)
where $`d\stackrel{~}{k}=dk/(4\pi \omega _k)`$ is a relativistic invariant measure and

$$\omega _k=|k|$$
(19)
with $`c=1`$. The creation and annihilation operators $`a_k^{}`$ and $`a_k`$ of the photon with wave vector $`k`$ obey the commutation relation

$$[a_k^{},a_k^{}]=4\pi \omega _k\delta (kk^{})$$
(20)
We construct a wave packet from a linear combination of normal modes.

$$|\mathrm{\Phi }_{x_0,b}(t)=_{\mathrm{}}^+\mathrm{}𝑑\stackrel{~}{k}\varphi _{x_0,b}(k)e^{i\omega _kt}a_k^{}|0$$
(21)
Here $`|0`$ is the vacuum state for the field (18). The fact that our wave packet is obtained by the action of creation operators on the vacuum state implies that the state consists of normal modes with positive energy. As before, we chose the function $`\varphi _{x_0,b}(k)`$ so that the wave packet is localized at the time $`t=0`$ in a domain with center $`x_0`$ and width $`2b`$,

$$\varphi _{x_0,b}(k)=(2\omega _k)^{1/2}_{\mathrm{}}^+\mathrm{}𝑑xe^{ikx}\mathrm{\Phi }_{x_0,b}(x,0)$$
(22)
where

$$\mathrm{\Phi }_{x_0,b}(x,0)=\frac{1}{(2b)^{1/2}}\mathrm{\Theta }(b|xx_0|)$$
(23)
The function $`\mathrm{\Phi }_{x_0,b}(x,0)`$ is normalized to ensure the normalization of the state $`|\mathrm{\Phi }_{x_0,b}(0)`$ in (21).

Let us introduce the operator $`\rho (x)`$:

$$\rho (x)=a^{}(x)a(x)$$
(24)
where $`a^{}(x)`$ and $`a(x)`$ are defined by

$`a^{}(x)`$ $`=`$ $`{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑\stackrel{~}{k}(2\omega _k)^{1/2}e^{ikx}a_k^{}`$
$`a(x)`$ $`=`$ $`{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑\stackrel{~}{k}(2\omega _k)^{1/2}e^{ikx}a_k`$ (25)
These operators satisfy the commutation relation

$$[a(x),a^{}(x^{})]=\delta (xx^{})$$
(26)
This construction follows the ideas of the construction of positions operators by Newton and Wigner . We shall express the localization of our state $`|\mathrm{\Phi }_{x_0,b}(t)`$ in terms of the expectation value of the operator $`\rho (x)`$. We call a state localized if the expectation value of $`\rho (x)`$ in this state vanishes when $`x`$ is outside a finite region. Our choice of $`\varphi _{x_0,b}(k)`$ in (22) and $`\mathrm{\Phi }_{x_0,b}(x)`$ in (23) guarantees that the state $`|\mathrm{\Phi }_{x_0,b}(t)`$ is localized at $`t=0`$ in the domain $`[x_0b,x_0+b]`$, i.e.,

$$\mathrm{\Phi }_{x_0,b}(0)|\rho (x)|\mathrm{\Phi }_{x_0,b}(0)=\frac{1}{2b}\mathrm{\Theta }(b|xx_0|)$$
(27)
Using (25) we obtain the time evolution of this quantity

$$\mathrm{\Phi }_{x_0,b}(t)|\rho (x)|\mathrm{\Phi }_{x_0,b}(t)=_{\mathrm{}}^+\mathrm{}𝑑\stackrel{~}{k}𝑑\stackrel{~}{k}^{}(4\omega _k\omega _k^{})^{1/2}e^{i(kk^{})x}\mathrm{\Phi }_{x_0,b}(t)|a_k^{}a_k^{}|\mathrm{\Phi }_{x_0,b}(t)$$
(28)
where $`\mathrm{\Phi }_{x_0,b}(t)|a_k^{}a_k^{}|\mathrm{\Phi }_{x_0,b}(t)`$ is expressed using our form of the wave packet (21) as follows:

$$\mathrm{\Phi }_{x_0,b}(t)|a_k^{}a_k^{}|\mathrm{\Phi }_{x_0,b}(t)=_{\mathrm{}}^+\mathrm{}𝑑\stackrel{~}{l}𝑑\stackrel{~}{l}^{}\varphi _{x_0,b}^{}(l)\varphi _{x_0,b}(l^{})e^{i(\omega _l\omega _l^{})t}0|a_la_k^{}a_k^{}a_l^{}^{}|0$$
(29)
Using the commutation relation (20) we integrate (29) over $`l`$ and $`l^{}`$ and then insert the result into (28). Taking into account the positivity of energy (19) and the form of the function $`\varphi _{x_0,b}(k)`$ in (22) we obtain

$$\mathrm{\Phi }_{x_0,b}(t)|\rho (x)|\mathrm{\Phi }_{x_0,b}(t)=\left|\frac{1}{2\pi (2b)^{1/2}}_{x_0+b}^{x_0b}𝑑x^{}_{\mathrm{}}^+\mathrm{}𝑑ke^{i|k|t+ik(xx^{})}\right|^2$$
(30)
By comparison with (7) we see that this quantity is equal to the absolute value squared of the classical function $`\mathrm{\Phi }(x,t)`$ up to the normalization constant. Our discussion of nonlocality remains, therefore, also valid in the quantum case, and the expression inside the absolute value in (30) is a superposition of two nonlocal wave packets that move in opposite directions at the speed of light. In the appendix give a second example using an analogy with Fermi’s problem .

Let us note that in quantum field theory localization depends on the observable. If a state is local from the point of view of one observable, it can be nonlocal from the point of view of another. In our example the expectation value of the operator $`\rho (x)`$ is local in the state $`|\mathrm{\Phi }_{x_0,b}(0)`$. At the same time, the energy density of the field in the same state is nonlocal. Indeed, the energy density $`T_{00}(x)`$ of the free massless field (18) is

$$T_{00}(x)=\frac{1}{2}\left(\left(\frac{\widehat{\psi }}{t}\right)^2+\left(\frac{\widehat{\psi }}{x}\right)^2\right)$$
(31)
It contains the derivatives of the field operator. (As we have seen in classical case, the time derivative of the function $`\mathrm{\Phi }(x,t)`$ is nonlocal even at $`t=0`$). To determine the expectation value of $`T_{00}(x)`$ in the state $`|\mathrm{\Phi }_{x_0,b}(0)`$ taking into account the positivity of energy (19) we first calculate this expectation value for a finite $`t`$ and then take the limit $`t0`$. Using (18)-(23) obtain:

$$\mathrm{\Phi }_{x_0,b}(0)|T_{00}(x)|\mathrm{\Phi }_{x_0,b}(0)=\frac{1}{4\pi b}\left(\frac{1}{|xb|}+\frac{1}{|x+b|}\frac{1\mathrm{sign}(xb)\mathrm{sign}(x+b)}{\sqrt{|xb|}\sqrt{|x+b|}}\right)$$
(32)
where we put $`x_0=0`$ to simplify the expression. This quantity is obviously nonlocal.

## 4 Conclusion

Positivity of energy for quantum fields (or frequency for classical fields) leads to a decomposition of localized wave packets in terms of nonlocal wave packets with long tails. The long tails, which cancel each other initially, appear immediately afterwards as the nonlocal wavepackets move in opposite directions. In our examples the long tails decay with the distance $`x`$ according to $`b/x`$ for $`b/x1`$, where $`b`$ is the size of the localized wave packet. They are precursors to the usual wave propagation. Although we may have instant interactions, these are not the result of superluminal propagation, but of “preformed” structures.

We shall study the interaction between nonlocal structures in a separate paper. We have then “contact interactions”, due to the overlapping of the long tails. We shall also show that the photon clouds around atoms and molecules are nonlocal, which leads to the precursor effect and eliminates the apparent deviation from causality in Fermi’s two-atom problem. However, it is true that the two atoms “feel” each other instantaneously. Even inside a relativistic theory (the wave equation is Lorenz invariant) there is place for instantaneous interactions due to nonlocality.

Einstein’s relativistic events are associated to four dimensional points. Here we see nonlocal but still relativistic events that are due to the instability of localization as shown in the examples presented in this paper.

ACKNOWLEDGEMENTS

The authors would like to thank Prof. B. Pavlov and Prof. I. Antoniou for comments and fruitful discussions. This work was carried out with financial support of the International Solvay Institutes, the European Commission ESPRIT Project 28890 NTGONGS. This work was partially supported by the Engineering Research Program of the Office of Basic Energy Sciences at the Department of Energy Grant No. DE-FG03-94ER14465, the Robert A Welch Foundation Grant No. F-0365.

APPENDIX

Here we show an analogy of our problem with Fermi’s two-atom problem. We prepare our wave packet $`|\mathrm{\Phi }_{x_0,b}(t)`$ “localized” at $`t=0`$. At time $`t`$ we project this state on an second wave packet $`|\mathrm{\Phi }_{x_1,b}(0)`$, which is localized in a domain with center $`x_1x_0`$ and width $`2b`$, and plays the role of a measurement device. We choose $`x_1`$ so that $`x_1x_0>2b`$. Then, at $`t=0`$ these two states do not overlap and the scalar product $`\mathrm{\Phi }_{x_1,b}(0)|\mathrm{\Phi }_{x_0,b}(t)`$ vanishes. We consider this scalar product as a function of $`t`$, i.e., at each moment $`t`$ we project the evolving wave packet $`|\mathrm{\Phi }_{x_0,b}(x,t)`$ on the localized state $`|\mathrm{\Phi }_{x_1,b}(0)`$. As we have shown, the initially localized packet, which evolves in time, is nonlocal immediately after $`t=0`$ and our scalar product has a non-vanishing value. This can be interpreted as the second (localized) wave packet “feeling” the existence of the first one even at $`t<x_1x_02b`$ when the causal component of the first wave packet still did not reach the domain of localization of the second wave packet. The scalar product $`\mathrm{\Phi }_{x_1,b}(0)|\mathrm{\Phi }_{x_0,b}(t)`$ grows as the overlapping of the two wave packets increases (see Fig. 4). We expect some essential change of this growth when the main part of the first wave packet corresponding to the position of its local (“causal”) component reaches the domain of localization of the second wave packet. Then, as the moving component goes away, the scalar product decreases. To show this, we write this scalar product using (21) in the following form

$$\mathrm{\Phi }_{x_1,b}(0)|\mathrm{\Phi }_{x_0,b}(t)=_{\mathrm{}}^+\mathrm{}𝑑k𝑑k^{}\varphi _{x_1,b}^{}(k)\varphi _{x_0,b}(k^{})e^{i\omega _kt}0|a_ka_k^{}^{}|0$$
$`(A1)`$
We perform the integration over $`k^{}`$ with the help of the commutation relation (20). Then, inserting (22) and (23) we obtain

$$\mathrm{\Phi }_{x_1,b}(x,0)|\mathrm{\Phi }_{x_0,b}(x,t)=\frac{1}{4\pi b}_{x_1b}^{x_1+b}𝑑x^{}_{x_0b}^{x_0+b}𝑑x^{\prime \prime }_{\mathrm{}}^+\mathrm{}𝑑ke^{i|k|t+ik(x^{}x^{\prime \prime })}$$
$`(A2)`$
Using a regularisation similar to (7), we integrate over $`k`$ and obtain

$$\mathrm{\Phi }_{x_1,b}(0)|\mathrm{\Phi }_{x_0,b}(t)=\frac{1}{4\pi ib}_{x_1b}^{x_1+b}𝑑x^{}_{x_0b}^{x_0+b}𝑑x^{\prime \prime }\left(\frac{1}{tx^{}+x^{\prime \prime }i0}\frac{1}{t+x^{}x^{\prime \prime }i0}\right)$$
$`(A3)`$
The integration over $`x^{}`$ and $`x^{\prime \prime }`$ and the rearrangement of terms give us

$$\mathrm{\Phi }_{x_1,b}(0)|\mathrm{\Phi }_{x_0,b}(t)=\psi _1(x_1x_0+ti0)\psi _1(x_1x_0t+i0)$$
$`(A4)`$
where

$$\psi _1(x)=\frac{1}{4\pi ib}\left((x2b)\mathrm{ln}(x2b)+(x+2b)\mathrm{ln}(x+2b)2x\mathrm{ln}(x)\right)$$
$`(A5)`$
Then, using (12), (13) we come to

$$\mathrm{\Phi }_{x_1,b}(0)|\mathrm{\Phi }_{x_0,b}(t)=\psi _2(x_1x_0t)+\psi _2^{}(x_1x_0+t)$$
$`(A6)`$
where

$`\psi _2(x)`$ $`=`$ $`{\displaystyle \frac{1}{8b}}\left(|x2b|+|x+2b|2|x|\right)`$
$`+`$ $`{\displaystyle \frac{i}{4\pi b}}\left(|x2b|\mathrm{ln}|x2b|+|x+2b|\mathrm{ln}|x+2b|2|x|\mathrm{ln}|x|\right)`$
Fig. 5 shows the real part, the imaginary part and the absolute value of the scalar product $`\mathrm{\Phi }_{x_1,b}(0)|\mathrm{\Phi }_{x_0,b}(t)`$ as functions of $`t`$ for two different values of the wave packet’s width $`b`$. The real component is non-vanishing only in the time interval corresponding to the overlapping of the localized components. In contrast, the imaginary part is non-vanishing immediately after $`t=0`$, because it reflects the overlapping of nonlocal components. Fig. 5 also shows that the “causal” part of the effect is much bigger than the contribution of the long tails, if the size of the wave packet is much less than the distance between the domains of localization.

| 1 | The real part (dashed lines) and the imaginary part (solid lines) |  |
| --- | --- | --- |
|  | of $`\psi (xt)`$ (a), $`\psi ^{}(x+t)`$ (c) and $`\mathrm{\Phi }(x,t)=\psi (xt)+\psi ^{}(x+t)`$ (e) |  |
|  | as functions of $`x`$ at $`t=0`$; the absolute values $`|\psi (xt)|`$ (b), $`|\psi ^{}(x+t)|`$ (d) |  |
|  | and $`|\mathrm{\Phi }(x,t)||\psi (xt)|+|\psi ^{}(x+t)|`$ (f) as functions of $`x`$ at $`t=0`$. |  |
| 2 | The real part (dashed lines) and the imaginary part (solid lines) |  |
|  | of $`\psi (xt)`$ (a), $`\psi ^{}(x+t)`$ (c) and $`\mathrm{\Phi }(x,t)=\psi (xt)+\psi ^{}(x+t)`$ (e) |  |
|  | as functions of $`x`$ at $`t=0.25s`$; the absolute values $`|\psi (xt)|`$ (b), $`|\psi ^{}(x+t)|`$ (d) |  |
|  | and $`|\mathrm{\Phi }(x,t)||\psi (xt)|+|\psi ^{}(x+t)|`$ (f) as functions of $`x`$ at $`t=0.25s`$. |  |
| 3 | The real part (dashed lines) and the imaginary part (solid lines) of $`\psi (xt)`$ (a), |  |
|  | $`\psi ^{}(x+t)`$ (c) and $`\mathrm{\Phi }(x,t)=\psi (xt)+\psi ^{}(x+t)`$ (e) as functions of $`x`$ |  |
|  | at $`t=1s`$; the absolute values $`|\psi (xt)|`$ (b), $`|\psi ^{}(x+t)|`$ (d) and |  |
|  | $`|\mathrm{\Phi }(x,t)||\psi (xt)|+|\psi ^{}(x+t)|`$ (f) as functions of $`x`$ at $`t=1s`$. |  |
| 4 | $`\mathrm{\Phi }_{0,b}(t)|\rho (x)|\mathrm{\Phi }_{0,b}(t)`$ (evolving object) and $`\mathrm{\Phi }_{2,b}(0)|\rho (x)|\mathrm{\Phi }_{2,b}(0)`$ |  |
|  | (right rectangle) with no overlap at $`t=0`$ (a), overlapping only by the |  |
|  | nonlocal tail at $`t=0.75s`$ (b), overlapping also by the local (causal) |  |
|  | component at $`t=1.25s`$ (c). |  |
| 5 | The real part (a) and (b), the imaginary part (c) and (d), the absolute value |  |
|  | (e) and (f) of $`\mathrm{\Phi }_{2,b}(0)|\mathrm{\Phi }_{0,b}(t)`$ as functions of $`t`$ for $`b=0.5ls`$ ((a), (c), (e)) and |  |
|  | $`b=0.01ls`$ ((b), (d), (f)). |  |


text:
# Acknowledgments

## Acknowledgments

This project is supported by the Ministry of Science and Technology of the Republic of Slovenia and by the Rector’s Fund of the University of Maribor. VR acknowledges the support of the work by the grant of the Ministry of Science and Technology of the Republic of Slovenia and the Abdus Salam ICTP (Trieste) Joint Programme.


text:
# Generating moment equations in the Doi model of liquid–crystalline polymers

## I Introduction

Kinetic theory is a powerful analytical tool for describing the dynamics of dilute and semi–dilute solutions of polymers in terms of a diffusion equation for the particle distribution function or, equivalently, by the full system of moment equations. In general, the moment system has to be truncated at some level. The problem of the “closure approximation” is very well–known in the literature, and an enormous amount of suggestions have been analyzed in the case where each moment couples only to a few higher–order moments. However, for some kinetic equations, the time evolution of each moment couples to an infinite set of higher moments, so that further analytical work is often precluded since closure approximations are less studied for this case. In this paper we present a simple method that allows to derive moment equations with a finite coupling valid for a wide class of kinetic equations.

In order to be specific, we consider a particularly important example: the Doi theory of liquid–crystalline polymers (LCP), subject to the Onsager excluded–volume potential . As it is well–known, in this model each moment equation depends on an infinite set of higher–order moments. In the original work , this problem was treated in two steps: First, the Onsager potential was replaced by a different, phenomenological potential of the Maier–Saupe type , which gives raise to a coupling to the next higher moment only. In the second step, the “decoupling” approximation was used to solve the resulting closure problem for the second moment. Subsequent extensive studies were focused on improvements of the second step . At the same time, we are not aware of improvements on the first step and closure approximations are limited to the Maier–Saupe potential up to now. However, it would be desirable to deal with the true Onsager potential, not only because it becomes exact in the limit of low concentrations of perfectly rigid rod–like molecules, but also because it contains no phenomenological parameters and therefore gives more quantitative predictions. In addition, the Onsager potential is preferred in the study of the influence of flow on the isotropic–nematic transition, since it gives a clear–cut prediction of the range of coexistence of the equilibrium isotropic and nematic phase, both in stationary  and non–stationary  flows. The method, which we propose in this work, leads to an approximation of the Onsager potential, which, to the lowest order, is at the same time as simple as the Maier–Saupe potential but also closer to the true Onsager potential. Moreover, corrections to this approximation can be obtained in a systematic manner.

## II The Doi model

Let $`\psi (𝒖;t)`$ be the probability distribution function for a rigid rod–like polymer molecule to be oriented parallel to the unit vector $`𝒖`$. The time evolution of $`\psi `$ in the presence of flow and the Onsager excluded–volume potential was given by Doi  and may be written as:

$$_t\psi =\left[𝒖\times \left(𝜿𝒖\psi \right)\right]+\widehat{D}_\mathrm{r}\psi \left(\frac{\delta A}{\delta \psi (𝒖)}\right).$$
(1)
Here $`=𝒖\times /𝒖`$ is the rotational operator, $`/𝒖`$ the gradient on the unit sphere, $`𝜿`$ the gradient of the velocity, $`\widehat{D}_\mathrm{r}`$ the rotational diffusivity, $`\delta /\delta \psi `$ the functional derivative and $`A=A_0+A_1`$ is the free energy functional per molecule divided by $`k_\mathrm{B}T`$,

$`A_0`$ $`=`$ $`\mathrm{ln}\nu 1+\mathrm{ln}\psi (𝒖)`$ (3)
$`A_1`$ $`=`$ $`{\displaystyle \frac{U}{2}}\sqrt{1(𝒖𝒘)^2}.`$ (4)
$`U=2bL^2\nu `$ is the reduced excluded–volume, $`2b`$ and $`L`$ are the diameter and the length of the rod–like polymeric molecules, respectively, and $`\nu `$ is the number of molecules per unit volume. Here and below we use the following notations for averages: $`f(𝒖)=f(𝒖)\psi (𝒖)𝑑𝒖`$, and $`f(𝒖,𝒘)=f(𝒖,𝒘)\psi (𝒖)\psi (𝒘)𝑑𝒖𝑑𝒘`$. $`A_0`$ describes the loss of entropy with molecular alignment, while $`A_1`$ expresses the Onsager free energy of steric interaction in the second virial approximation . Following Doi and Edwards , the rotational diffusivity is approximated by

$$\widehat{D}_\mathrm{r}\overline{D}_\mathrm{r}=D_\mathrm{r}\left[\frac{4}{\pi }\sqrt{1(𝒖𝒘)^2}\right]^2,$$
(5)
where $`D_\mathrm{r}`$, the rotational diffusion coefficient for a rod in an isotropic, semi–dilute solution of like rods, is related to the rotational diffusion constant for a dilute solution, $`D_{\mathrm{r0}}`$, by $`D_\mathrm{r}=cD_{\mathrm{r0}}(\nu L)^2`$ with an empirical coefficient $`c`$. Nonlinearity of Eq. (1) in $`\psi `$ brought about by the potential (4) reflects the mean–field nature of the Onsager theory of the excluded–volume effect. The self–consistent potential, identified by Doi , is related to the free energy of interaction, $`V(𝒖)=k_\mathrm{B}T\delta A_1/\delta \psi (𝒖)`$. Various phases of the LCP are conveniently described by the order parameter $`\mathrm{SS}=𝒖𝒖(1/3)\mathrm{𝟏}`$, where $`\mathrm{𝟏}`$ is the unit tensor. It is reasonable therefore to look for approximate formulations of the dynamics in terms of the order parameter alone. However, as mentioned above, the time evolution equation for $`\mathrm{SS}`$ couples to an infinite number of moments of $`\psi `$. In the derivation given by Doi, this difficulty was circumvented by replacing the Onsager potential (4) by a different, phenomenological expression of the Maier–Saupe type :

$$A_1^{\mathrm{MS}}=a_0\frac{a_1}{2}U\mathrm{SS}:\mathrm{SS},$$
(6)
where $`a_0`$ and $`a_1`$ are parameters independent of $`\psi `$. A further separate treatment of the diffusivity (5) is also necessary. A compact presentation of the entire development is given by Doi and Edwards . The Doi model with the Maier–Saupe potential (6) constitutes the basic kinetic model of LCP used by many authors for analytical studies to derive equations for the order parameter. As is well–known, the kinetic equation (1) with the potential (6) does not give a closed equation for the order parameter but contains also the higher–order moment $`𝒖𝒖𝒖𝒖`$, and therefore constitutes a further problem of closure. The original Doi approach was based on the decoupling approximation for the fourth–order moments of $`\psi `$ in terms of $`\mathrm{SS}`$. Improvements on the decoupling approximation are currently under active research .

## III Generating moment equations

In this communication, we demonstrate that a different self–consistent treatment of the kinetic equation (1) is possible. Modifications concern only the relaxational part of the Eq. (1), specifically, the excluded–volume potential (4) and the diffusivity (5), and therefore we consider the case $`𝜿=0`$ in the sequel to simplify notations. Specifically, we employ the cumulant expansion of the potential (4) and the diffusivity (5). The leading term of this expansion results in an effective potential that differs from the Maier–Saupe potential (6), and which contains a non-polynomial dependence on the order parameter $`\mathrm{SS}`$.

In the second virial approximation, the free energy of interaction, $`A_1`$ can be written as $`A_1=(\nu /2)\beta (𝒖,𝒘)`$. If only excluded–volume interactions are present, the second virial coefficient $`\beta `$ corresponding to the Onsager expression (4) is of the form $`\beta (𝒖,𝒘)=\beta ((𝒖𝒘)^2)`$, with $`\beta (x)=2bL^2\sqrt{1x}`$. Specifically, expanding $`\beta (x)`$ in a Taylor series and interchanging summation and averaging in this expansion, we get $`\beta (x)=_{n=0}^{\mathrm{}}a_nx^n`$, where $`a_n`$ are numerical coefficients. Each average $`x^n`$ can be represented in terms of cumulants $`x^k_c`$ of order $`kn`$. Resummation of the series leads to

$$\beta (x)=\beta (x)+\underset{m=1}{\overset{\mathrm{}}{}}\frac{1}{m!}\left(\frac{x^2_c}{2}\right)^m\beta ^{(2m)}(x)+\mathrm{},$$
(7)
where $`\beta ^{(2m)}`$ is the $`2m`$–th derivative of $`\beta `$ and ellipses denote terms including third or higher–order cumulants as factors. Therefore, the functional $`A_1`$ (4) can be split as $`A_1=A_1^{(1)}+B`$, where $`A_1^{(1)}`$ is the free energy, corresponding to the total neglect of second and higher–order cumulants in each term of the expansion,

$$A_1^{(1)}=\frac{U}{2}\sqrt{1𝒖𝒖:𝒖𝒖}.$$
(8)
In terms of the order parameter $`\mathrm{SS}`$, $`A_1^{(1)}`$ may be rewritten as $`A_1^{(1)}=(U/\sqrt{6})\sqrt{1(3/2)\mathrm{SS}:\mathrm{SS}}`$. By the mean value theorem, it is easy to see that $`A_1^{(1)}`$ gives an upper bound to $`A_1`$, $`A_1A_1^{(1)}`$, for the present case of excluded–volume interactions. The functional $`B`$ contains the higher–order cumulants. While all powers of the second cumulants are displayed in Eq. (7), in general it is not a priori clear whether it is more important to keep powers of the second cumulant or higher cumulants. However, we generally expect the linear term in the second order cumulant to be most important. The corresponding term $`m=1`$ in Eq. (7), $`A_1^{(2)}=(1/2)x^2_c\beta ^{\prime \prime }(x)`$, gives the first non-vanishing contribution to $`B`$,

$$A_1^{(2)}=\frac{U}{16}[(𝒖𝒘)^2𝒖𝒖:𝒖𝒖]^2(1𝒖𝒖:𝒖𝒖)^{3/2}.$$
(9)
Keeping only the first $`N`$ cumulants in the expansion (7), the functional (4) is approximated by non–linear functions of the first $`2N`$ moments of $`\psi `$. Inserting this approximation in the time evolution equation (1) amounts to an approximation of the self–consistent Onsager potential $`V`$ in terms of a polynomial of order $`2N`$ in $`𝒖`$ but with non–linear coefficients. In this approximation, the time evolution of the $`2n`$–th moment contains only the first $`(2n+2N)`$–th moments. With this, moment equations can be generated, that approach the original equations in a systematic way, thereby containing only a finite number of moments at each stage.

## IV Testing the approximation

Clearly, the above procedure is most valuable if the first terms, $`A_1^{(1)}`$, etc. , already provide a good approximation to the full expression $`A_1`$. While the general validity of $`A_1^{(1)}`$ as a good approximation to $`A_1`$ is a rather delicate problem, it should be mentioned that it is so at least in two limiting cases. Namely, for the isotropic state, the value of $`A_1^{(1)}`$ differs from $`A_1`$ for less than 5%, while in the fully ordered state the approximation (8) becomes exact. Moreover, on the submanifold of distribution functions of the form

$$\psi _\alpha (𝒖)=\frac{\alpha }{4\pi \mathrm{sinh}\alpha }\mathrm{cosh}(\alpha 𝒖𝒏),$$
(10)
where $`𝒏`$ is an arbitrary unit vector, and $`0\alpha \mathrm{}`$, the functional $`A_1^{(1)}`$ turns out to approximate $`A_1`$ very well for all values of the parameter $`\alpha `$ between the isotropic state, $`\alpha =0`$, and the fully ordered state, $`\alpha =\mathrm{}`$. To show this, we plot in Fig. 1 the functions $`A_1(\alpha )`$, $`A_1^{(1)}(\alpha )`$ and $`A_1^{(1)}(\alpha )+A_1^{(2)}(\alpha )`$, that result upon inserting the ansatz (10) into (4), (8), (9), respectively. Note, that $`A_1(\alpha )`$, $`A_1^{(1)}(\alpha )`$ and $`A_1^{(2)}(\alpha )`$ can be calculated analytically. For convenience, we plot the functions against the scalar order parameter, defined as $`S=\sqrt{(3/2)\mathrm{SS}:\mathrm{SS}}`$. Including $`A_1^{(2)}`$ does not only reduce the error of the approximate value of $`A_1`$ in the isotropic state to 1.5%, but improves the accuracy of the approximation over the whole range of $`S`$. For comparison, we included in Fig. 1 also the free energy $`A_1^{\mathrm{MS}}`$, corresponding to the Maier–Saupe expression (6), thereby choosing the undetermined constant so that the limit of the fully ordered state is matched correctly. Note, however, that in any case $`A_1^{\mathrm{MS}}`$ decays asymptotically like $`1/\alpha `$, for $`\alpha 1`$, whereas $`A_1`$ and $`A_1^{(1)}`$ behave like $`1/\sqrt{\alpha }`$ in this regime. We included in Fig. 1 also the derivative of the above functions, since they are related to the self–consistent potential $`V`$. Fig. 1 shows that also the derivative of $`A_1^{(1)}`$ provides a good approximation to the derivative of $`A_1`$, with correct limiting behavior near the isotropic and fully ordered state. Note that including the first correction $`A_1^{(2)}`$ yields excellent agreement to the true Onsager prediction. The Maier–Saupe potential captures the main features but, besides an undetermined constant, shows the wrong behavior near the fully ordered state. The ansatz (10), originally proposed by Onsager , is known to approximate the equilibrium distribution very well. Therefore we conclude that $`A_1^{(1)}`$ represents a good approximation to $`A_1`$, at least on a representative subset of distribution functions.

## V Thermodynamic consistency

It is worth mentioning that the presentation given so far can easily be cast into the recently developed GENERIC formalism of nonequilibrium thermodynamics . In the absence of potential forces, the example of rigid dumbbells, which are equivalent to the model of rigid rods, is formulated within the GENERIC formalism in Ref. . The mean field potentials considered above can be included in a straightforward manner, if we recognize that $`A_1=S_1`$, where $`S_1`$ is the entropic contribution per molecule to the free energy of interaction divided by $`k_\mathrm{B}`$. Formulating the original model as well as the approximations within the GENERIC formalism guarantees that our treatment is in accordance with the principles of nonequilibrium thermodynamics. This becomes especially important if the present model is considered in nonisothermal situations. For example, the structure of the GENERIC formalism requires the polymeric contribution to the elastic stress to be

$$\sigma _{\alpha \beta }^\mathrm{e}=3\nu k_\mathrm{B}TS_{\alpha \beta }\nu k_\mathrm{B}T(𝒖\times \frac{\delta S_1}{\delta \psi })_\alpha u_\beta .$$
(11)
Eq. (11) agrees with the result of Doi , obtained upon varying the free energy functional.

## VI The lowest order approximation

In the sequel, we will adopt the lowest order approximation $`A=A^{(1)}=A_0+A_1^{(1)}`$, where $`A_1^{(1)}`$ is given by Eq. (8), and $`A_0`$ is given by Eq. (3). This amounts to neglect of all higher order correlations in Eq. (7), or, equivalently, setting $`B=0`$. Substituting $`A^{(1)}`$ instead of $`A`$ into Eq. (1), we derive

$$_t\psi =\widehat{D}_\mathrm{r}\left[\psi \psi \left(\frac{U𝒖𝒖:𝒖𝒖}{2\sqrt{1𝒖𝒖:𝒖𝒖}}\right)\right].$$
(12)
It is now possible to identify the self–consistent potential as

$$V^{(1)}(𝒖)=\left(\frac{Uk_\mathrm{B}T}{2}\right)\frac{1𝒖𝒖:𝒖𝒖}{\sqrt{1𝒖𝒖:𝒖𝒖}},$$
(13)
which can be compared to the expression obtained from inserting the Maier–Saupe free energy (6) into Eq. (1)

$$V_{\mathrm{MS}}(𝒖)=a_2a_1Uk_\mathrm{B}T𝒖𝒖:𝒖𝒖,$$
(14)
where $`a_2`$ is an arbitrary constant. The normalized equilibrium solutions to the Eq. (12) are $`\psi _{\mathrm{eq}}^{(1)}=Z^1\mathrm{exp}[V^{(1)}/k_\mathrm{B}T]`$.

The rotational diffusivity (5) is related to the free energy of interaction, since $`\overline{D}_\mathrm{r}=D_\mathrm{r}\left[\frac{4}{\pi }A_1/(U/2)\right]^2`$. Substituting $`A_1=A_1^{(1)}`$ gives

$$\overline{D}_\mathrm{r}^{(1)}=(3\pi ^2/32)D_\mathrm{r}[1(3/2)\mathrm{SS}:\mathrm{SS}]^1.$$
(15)
The diffusion coefficient $`\overline{D}_\mathrm{r}^{(1)}`$ (15) is positive in the entire physically meaningful range of the order parameter $`\mathrm{SS}`$. Expression (15) should be compared with the Doi phenomenological result:

$$\overline{D}_{\mathrm{rD}}=D_\mathrm{r}[1(3/2)\mathrm{SS}:\mathrm{SS}]^2.$$
(16)
While we have not found an argument which of the two powers, $`1`$ or $`2`$, is more consistent, it should be stressed that our derivation of the diffusion coefficient does not need any further assumptions or adjustable parameters, while the derivation of Eq. (16)  requires the matching of $`\overline{D}_\mathrm{r}`$, resp. $`A_1^{\mathrm{MS}}`$ in both, the isotropic and the fully ordered state. Due to its relation to $`A_1^{(1)}`$, the diffusion coefficient $`\overline{D}_\mathrm{r}^{(1)}`$ (15) has a correct limit in the fully ordered state ($`D_\mathrm{r}/\overline{D}_\mathrm{r}^{(1)}=0`$ as soon as $`\mathrm{SS}:\mathrm{SS}=2/3`$ in the ordered state), while the opposite limit of the isotropic state ($`\overline{D}_\mathrm{r}=D_\mathrm{r}`$) is matched within $`8\%`$. Again, the first correction (9) reduces the error in this limit to less than $`3\%`$.

If we adopt (5) and approximate $`\overline{D}_\mathrm{r}`$ by $`\overline{D}_\mathrm{r}^{(1)}`$ (15), the time evolution of the order parameter $`\mathrm{SS}`$ can be derived from Eq. (12) by the so–called Prager procedure

$$_t\mathrm{SS}=6\overline{D}_\mathrm{r}^{(1)}\mathrm{SS}+6\overline{D}_\mathrm{r}^{(1)}\frac{U^{}}{\sqrt{1(3/2)\mathrm{SS}:\mathrm{SS}}}(\mathrm{SS}𝒖𝒖\mathrm{SS}:𝒖𝒖𝒖𝒖),$$
(17)
with $`U^{}=U/\sqrt{6}`$. This expression differs from the result of Doi  not only in the diffusion coefficient and in the reduced excluded–volume $`U`$ due to the undetermined constant in the Maier–Saupe potential (14), but contains a non–polynomial dependence on the order parameter $`\mathrm{SS}`$, which becomes important in the nematic state.

## VII Conclusion

We have presented a systematic procedure that allows to derive approximate moment equations for the Doi model of LCP, which contain only a finite number of higher order moments. The first approximation for the Onsager excluded–volume interaction results in an effective potential (13) proportional to $`𝒖`$$`𝒖`$, but different from the Maier–Saupe form (14) and without free parameters. We find indications, that (13) approximates the true Onsager potential better than the Maier–Saupe potential. For higher accuracy, the first correction seems to be the most important contribution. All these approximations are in accordance with nonequilibrium thermodynamics.

Note, that we have not addressed the problem of solving the resulting kinetic equations or “closing” the moment equations. This work is currently under preparation. Nevertheless, for comparing Eq. (17) to the corresponding equation with the Maier–Saupe potential, we follow Refs.  and consider the decoupling approximation $`\mathrm{SS}:𝒖𝒖𝒖𝒖=\mathrm{SS}:𝒖𝒖𝒖𝒖`$. If the order parameter is assumed to be of the form $`S_{\alpha \beta }=S(t)[n_\alpha n_\beta \delta _{\alpha \beta }]`$, the relaxation equation for the scalar parameter $`S`$ is found to be as follows:

$`_tS`$ $`=`$ $`6\overline{D}_\mathrm{r}^{(1)}{\displaystyle \frac{A^{(1)}(S,U^{})}{S}},`$ (19)
$`A^{(1)}(S,U^{})`$ $`=`$ $`{\displaystyle \frac{S^2}{2}}{\displaystyle \frac{U^{}}{9}}\sqrt{1S^2}\left(1{\displaystyle \frac{3S}{2}}+2S^2\right){\displaystyle \frac{U^{}}{6}}\mathrm{arcsin}(S),`$ (20)
where $`U^{}=U/\sqrt{6}`$. Due to a non–polynomial character of the function $`A^{(1)}`$ (20), the relaxation equation (19) differs formally from the Landau–de Gennes counterpart derived by Doi for the Maier–Saupe potential. Expansion of the function (20) around $`S=0`$ reproduces the result of Doi for the Maier–Saupe potential, subject to a renormalization of the strength of the excluded–volume potential, and a difference in the coefficient in front of the $`S^4`$ term. Moreover, the relaxation implied by $`A^{(1)}`$ (20) is qualitatively similar to the one given by Doi result and distinguishes the same three regimes. For $`U<U_1`$, $`A^{(1)}`$ has only one minimum at $`S=0`$, so that the system finally becomes isotropic. For $`U_1<U<U_2`$, a second local minimum occurs. The system either becomes isotropic or nematic depending on the initial value of $`S`$. Finally, for $`U>U_2`$, the isotropic state becomes unstable and the system always approaches a nematic state.

Due to the undetermined constant $`a_1`$ in the Maier–Saupe potential (14), the Doi theory predicts the values $`U_1`$ and $`U_2`$ also in terms of $`a_1`$. On the contrary, the self–consistent potential (13) contains no free parameters, so that $`U_1=3^{1/4}\sqrt{8}/(\sqrt{3}1)6.22`$ and $`U_2=3\sqrt{6}7.34`$ may directly be compared to the values in the Onsager theory $`U_1=8.38`$ and $`U_2=10.67`$. However, the values of the order parameter $`S_1=1/4`$ and $`S_2=1/2`$ at $`U_1`$ resp. $`U_2`$ in the Doi theory do not depend on $`a_1`$ and can therefore be compared to the values predicted by (13): $`S_1=(\sqrt{3}1)/20.37`$ and $`S_2=1/\sqrt{2}0.71`$.

In Fig. 2, the equilibrium order parameter $`S_{\mathrm{eq}}`$ is shown as a function of $`\nu /\nu _2`$. The lower solid line shows the prediction of the Doi theory, whereas the upper solid line corresponds to the approximation (13). For the Maier–Saupe potential, $`S_{\mathrm{eq}}0`$ is given by

$$S_{\mathrm{eq}}=\frac{1}{4}+\frac{3}{4}\sqrt{1\frac{8\nu _2}{3\nu }}.$$
(21)
For the free energy (20), $`S_{\mathrm{eq}}`$ is given implicitly as the solution to the algebraic equation

$$\frac{1+S_{\mathrm{eq}}2S_{\mathrm{eq}}^2}{\sqrt{1S_{\mathrm{eq}}^2}}=\frac{\nu _2}{\nu }.$$
(22)
For large values of $`\nu /\nu _2`$, the solution (22) approaches the value $`S_{\mathrm{eq}}=1`$ and asymptotically behaves like the solution of the Doi theory $`S_{\mathrm{eq}}1\nu _2/\nu `$, for large $`\nu `$. Note, that the decoupling approximation corrects the asymptotic behavior of the Maier–Saupe potential near the fully ordered state.

As is well–known, the detailed form of the interaction potential can have significant effect on the behavior of the order parameter in the nematic phase . Specifically, the amount of order at the transition is known to be much smaller in the Maier–Saupe theory than in the Onsager model. For comparison, we included in Fig. 2 the values of the order parameter obtained from minimizing the true Onsager free energy numerically , where $`\nu _2`$ now corresponds to the true nematic transition. Although the analysis of the phase transitions via the dynamical approach is affected by the use of the decoupling approximation, the prediction of the self–consistent approach is much closer to the true Onsager values than is the Maier–Saupe potential.

Finally, it should be mentioned that approximations to the Onsager potential like Eqs. (8) and (9) can also be used in the case of potential flows, following the approach of Thirumalai  without additional assumptions.

To summarize, we have developed a direct approach to the Doi model with the Onsager potential. We have demonstrated that the resulting kinetic equation has much in common with the Doi model with the phenomenological Maier–Saupe potential. Corrections to the approximation developed here can be found in a systematic way from Eqs. (7) to (9) by taking into account higher order correlations. The approach to derive self–consistent moment equations is applicable to other kinetic equations which can be cast into the form (1).

## Figure captions

Free energy of excluded–volume interaction for distribution functions (10) plotted against the scalar order parameter $`S(\alpha )`$. From top to bottom: approximation $`A_1^{(1)}`$ (8), with first correction $`A_1^{(1)}+A_1^{(2)}`$ (9), true Onsager expression $`A_1`$ (4) and the Maier–Saupe free energy $`A_1^{\mathrm{MS}}`$ (6), when the limits $`S=0`$ and $`S=1`$ are matched. In the inset, the derivative of the above functions is shown as a function of $`S(\alpha )`$. The order of the curves from top to bottom is the same.

The equilibrium order parameter $`S_{\mathrm{eq}}`$ as a function of $`\nu /\nu _2`$. The figure shows the behavior due to the Maier–Saupe potential (21), lower curve, and the solution of (22) corresponding to $`A^{(1)}`$, upper curve, in the decoupling approximation. Full circles indicate the order parameter for the true Onsager potential in the static case (from ).


text:
# Quantum phase transitions in a linear ion trap

## I Introduction

More than two decades ago, when quantum optics was young, the quantum dynamics of collective spin systems interacting with a single bosonic degree of freedom was a major research problem. The model arose as an attempt to describe the interaction between a collection of two level atoms and a single mode of the radiation field. Walls and co workers were among the first to realise that such models provided ideal examples of the role of quantum fluctuations in the nonlinear interaction between matter and light. Quantum fluctuations were shown to drastically change the predictions of semiclassical theory in such systems. This phenomenon has appeared more recently in the discovery of quantum phase transitions in quantum spin glasses and other many body quantum systems. While the collective spin models did not directly apply to achievable experiments at the time, they did provide insight that subsequently proved important for many other quantum optical experiments including anti-bunching, squeezing, and cavity QED. In this paper we show that the models of a collective spin interacting with one or more bosonic modes can now be experimentally realised in modern ion trap systems of the kind proposed for quantum computation. An enormous effort has gone into making such systems work at the quantum level, with little interference form classical sources of noise, and a number of such experiments exist today. It would thus appear worthwhile to reconsider the collective spin models, and the associated quantum many-body effects exhibited by such systems, with a view to direct experimental realisation.

In particular we consider the Tavis-Cummings (TC) model, which can be realised in a linear ion trap of $`N`$ ions with the bosonic degree of freedom appearing as the quantised collective centre-of-mass motion. If each ion is coupled to the vibrational motion using an identical external (classical) laser detuned to the first red-sideband transition, the symmetry is such that the electronic degree of freedom for the ions can be described as a collective spin ($`N`$) and the reversible dynamics is well described by the TC model. The TC model is known to exhibit important nonlinear quantum effects including a quantum phase transition in which the (zero temperature) ground state undergoes a morphological change as a parameter is varied and averages of intensive quantities undergo a bifurcation.

## II The Tavis-Cummings model

The interaction Hamiltonian for N ions interacting with the centre of mass vibrational mode can be controlled by using different kinds of Raman laser pulses. A considerable variety of interactions has already been achieved or proposed . In this paper we consider the first red-sideband transition. The ion is assumed to be in a three dimensional anisotropic harmonic potential. Two dimensions are very tightly bound and are neglected. In the remaining dimension, an external laser couples the electronic state to the vibrational motion. If the vibrational frequency is large enough and the Lamb-Dicke limit applies the motional sidebands of the absorption of the electronic transition can be resolved and a laser detuned below the electronic resonance by one unit of the trap frequency can excite the electronic transition by absorbing one vibrational phonon, the additional energy required being made up by the laser. We will assume that the laser ( or lasers if a Raman process is used) is sufficeintly strong that it can be treated classically. Under these assumptions the Hamiltonian, in the interaction picture, is

$$H_I=\mathrm{}\mathrm{\Omega }\underset{i=1}{\overset{N}{}}(a\sigma _+^{(i)}+a^{}\sigma _{}^{(i)})$$
(1)
where the coupling constant is $`\mathrm{\Omega }=\eta \mathrm{\Omega }_0`$ where $`\eta ^2=E_r/(\mathrm{}M\omega _0)`$ is the Lambe-Dicke parameter with $`E_r`$ the recoil kinetic energy of the atom, $`\omega _0`$ is the trap vibrational frequency, and $`M`$ is the effective mass for the centre-of-mass mode. The Lamb-Dicke limit assumes $`\eta <<1`$, which is easily achieved in practice. The frequency, $`\mathrm{\Omega }_0`$ is the effective Rabi frequency for the electronic transition involved. The raising and lowering operators for each ion are defined by $`\sigma _{}=|ge|`$ and $`\sigma _+=|eg|`$. This sideband transition can be used to efficiently cool the ions to the collective centre-of-mass ground state, thus preparing the system in the vibrational ground state.

If the external laser field on each ion is identical (in amplitude and phase) the interaction Hamiltonian is

$$H_I=\mathrm{}\mathrm{\Omega }(a\widehat{J}_++a^{}\widehat{J}_{})$$
(2)
where we have introduced the bosonic annihilation operator $`a`$ for the centre-of-mass vibrational mode and where we have used the definition of the collective spin operators,

$$\widehat{J}_\alpha =\underset{i=1}{\overset{N}{}}\sigma _\alpha ^{(i)}$$
(3)
where $`\alpha =x,y,z`$ . Identical laser fields could easily be obtained by splitting a single, stabilised laser into multiple beams. The interaction Hamiltonian in Eq (2) specifies the Tavis-Cummings model. This model first appeared in quantum optics where the bosonic mode is the quantised field in a cavity. However this realisation is difficult to achieve experimentally. In contrast the vibrational mode realisation should be readily achieved. The dynamics resulting from this Hamiltonian is quite rich. Collective spin models of this kind were considered many decades ago in quantum optics. In much of that work however the collective spin underwent an irreversible decay. In the case of an ion trap model however we can neglect such decays due to the long lifetimes of the excited states. On the other hand heating of the vibrational centre-of-mass mode can induce irreversible dynamics in the system in a manner that has not been previously considered, and that is reminiscent of thermal effects in condensed matter physics.

We are interested in the driven Tavis-Cummings model in which the vibrational mode is subject to a linear forcing term which can easily be achieved by a suitable combination of Raman laser pulses, or by appropriate AC voltages applied to the trap electrodes. In this case the Hamiltonian, in the interaction picture, is given by

$$H_I=\mathrm{}\mathrm{\Omega }(a\widehat{J}_++a^{}\widehat{J}_{})+\mathrm{}E(a+a^{})$$
(4)
This may be written in terms of the hermitian canonical oscillator variables $`\widehat{X}=(a+a^{})/\sqrt{2}`$, $`\widehat{Y}=i(aa^{})/\sqrt{2}`$, and the canonical angular momentum variables $`\widehat{J}_x=(\widehat{J}_++\widehat{J}_{})/2`$, $`\widehat{J}_y=i(\widehat{J}_+\widehat{J}_{})/2`$, $`\widehat{J}_z=[\widehat{J}_+,\widehat{J}_{}]/2`$. It takes the form

$$H=\widehat{X}\widehat{J}_x\widehat{Y}\widehat{J}_y+\chi \widehat{X}$$
(5)
with $`\chi =E/\mathrm{\Omega }`$ and we have scaled the Hamiltonian by $`HH/\sqrt{2}\mathrm{\Omega }`$. This indicates that time is measured in units of $`\frac{1}{\sqrt{2}\mathrm{\Omega }}`$.

Alsing has shown that the ground state of this system, for weak driving, is a product state in which the bosonic mode is squeezed and the electronic states are rotated in the angular momentum space. We provide a direct proof of this statement below. However it is first useful to consider the dynamics of the equivalent semiclassical model as many of the results in the quantum case can be interpreted in terms of the features of the semiclassical model.

### A Semiclassical Tavis-Cummings model

The Tavis-Cummings model represents an interaction between a simple harmonic oscillator and a linear top for which there is a classical model which we now define. We choose the classical model so that the equations of motion are of the same form as the Heisenberg equations of motion for the quantum model. The classical Hamiltonian is defined as

$$=X𝒥_xY𝒥_y+\chi EX$$
(6)
where $`X,Y`$ are respectively the canonical oscillator position and momentum variables with the canonical Poisson bracket $`\{X,Y\}=1`$, while $`𝒥_k`$ are the three components of angular momentum for a classical top with the canonical Poisson brackets $`\{𝒥_i,𝒥_k\}=_kϵ_{ijk}𝒥_k`$. The equations of motion for a canonical coordinate $`w`$ is given as usual by Poisson bracket with the Hamiltonian $`\dot{w}=\{w,H\}`$. The equations of motion are,

$`\dot{X}`$ $`=`$ $`𝒥_Y`$ (7)
$`\dot{Y}`$ $`=`$ $`𝒥_X\chi `$ (8)
$`\dot{𝒥_x}`$ $`=`$ $`Y𝒥_z`$ (9)
$`\dot{𝒥_y}`$ $`=`$ $`X𝒥_z`$ (10)
$`\dot{𝒥_z}`$ $`=`$ $`X𝒥_y+Y𝒥_x`$ (11)
Note that these equations have a conservation law $`𝒥_x^2+𝒥_y^2+𝒥_z^2=\text{constant}`$.

We now justify this choice of classical Hamiltonian by noting that the Heisenberg equations of motion for the Hamiltonian Eq(4) have the same form as the semiclassical equations of motion with all variables replaced by the corresponding operators. We thus see that the semiclassical equations result form taking moments of the Heisenberg equations and factorising all product moments. The factorisation assumptions ignores correlations which scale as $`1/N`$ for the scaled operators $`\widehat{J}_\sigma /N`$. The conservation law $`𝒥_x^2+𝒥_y^2+𝒥_z=constant`$ is a reflection of the operator relation

$$\widehat{J}^2=\frac{N}{2}\left(\frac{N}{2}+1\right)$$
(12)
which in the semiclassical limit indicates that $`𝒥_x^2+𝒥_y^2+𝒥_z=\frac{N^2}{4}`$.

The classical equations have one nontrivial fixed point at $`X^{}=Y^{}=𝒥_y^{}=0`$ and $`𝒥_x^{}=\chi `$, $`𝒥_z^{}=\sqrt{N^2/4\chi ^2}`$. However as the conservation law requires that $`|𝒥_x|N/2`$ we see that we must have

$$\frac{2E}{N\mathrm{\Omega }}1\text{(below threshold)}$$
(13)
which corresponds to an energy of $`=0`$. We will refer to this as the below threshold case. As $`E`$ is increased from zero, the fixed point for the angular momentum system rotates about the $`𝒥_y`$ direction eventually reaching the equatorial plane at $`𝒥_x=N/2`$ at the threshold condition. The oscillator system always has zero amplitude below threshold. If we linearise around this fixed point we discover that it is an unstable hyperbolic point with time constant proportional to $`\frac{1}{\sqrt{𝒥_z^{}}}`$. Note that this time constant goes to infinity as the fixed point is approached as is typical for a hyperbolic fixed point.

We now consider the above threshold case

$$\frac{2E}{N\mathrm{\Omega }}1\text{(above threshold)}$$
(14)
Clearly the value of $`|𝒥_x|`$ cannot increase above $`N/2`$. Indeed there is no fixed point above threshold. However there is a special solution curve that continuously joins to the below threshold case for phase curves with $`=0`$.

To see this we consider making a canonical transformation by a rotation in both the $`XY`$ plane and in the $`𝒥_x,𝒥_y`$ plane (see figure 1). The canonical transformations are

$`X`$ $`=`$ $`\overline{X}\mathrm{cos}\theta +\overline{Y}\mathrm{sin}\theta `$ (15)
$`Y`$ $`=`$ $`\overline{Y}\mathrm{cos}\theta \overline{X}\mathrm{sin}\theta `$ (16)
$`𝒥_x`$ $`=`$ $`\overline{𝒥}_x\mathrm{cos}\theta \overline{𝒥}_y\mathrm{sin}\theta `$ (17)
$`𝒥_y`$ $`=`$ $`\overline{𝒥}_x\mathrm{cos}\theta +\overline{𝒥}_y\mathrm{sin}\theta `$ (18)
The Hamiltonian then takes the form

$$=\overline{X}(\overline{𝒥}_x+\chi \mathrm{cos}\theta )\overline{Y}(\overline{𝒥}_y\chi \mathrm{sin}\theta )$$
(19)
The phase curves with $`=0`$ now correspond to either

$`\overline{X}`$ $`=`$ $`0;\overline{𝒥}_y=\chi \mathrm{sin}\theta `$ (20)
or

$`\overline{Y}`$ $`=`$ $`0;\overline{𝒥}_x=\chi \mathrm{cos}\theta `$ (21)
These phase curves smoothly join the fixed point at threshold if $`𝒥_x=N/2`$ which implies

$$\mathrm{cos}\theta =\frac{N\mathrm{\Omega }}{2E}$$
(22)
These solutions are illustrated in figure 1. Note that as $`E\mathrm{}`$ we have that $`\overline{𝒥}_x`$ eventually points in the direction of $`𝒥_y`$ while phase curve in the oscillator phase space points along the $`Y`$ axis, indicating that for large driving the system is essentially a particle in a linear potential which accelerates at constant rate. These results were first obtained by Alsing and Carmichael.

## III Quantum states

First note that the ground state when there is no driving is $`|j,j_z|0_v`$ with a zero eigenvalue. This ground state corresponds to the fixed point of the semiclassical model with zero oscillator amplitude and angular momentum pointing in the $`𝒥_z`$ direction. We postulate that as the driving is increased form zero the ground state of the Hamiltonian Eq(4) is given by

$$|_0=S(r)R(\theta )|j,j_z|0_v$$
(23)
where $`|j,j_z|0_v`$ corresponds to all ions in the ground state and the vibrational mode in the ground state. The operator $`S(r)`$ is a squeezing operator defined by

$$S^{}(r)aS(r)=\mu a+\nu a^{}$$
(24)
with $`\mu =\mathrm{cosh}r,\nu =\mathrm{sinh}r`$.

The rotation operator $`R(\theta )`$ is defined by

$$R(\theta )=e^{\theta (\widehat{J}_+\widehat{J}_{})}$$
(25)
and corresponds to a rotation of $`2\theta `$ around the $`\widehat{J}_y`$ axis. Consider now

$$H_I|_0=SR\left(R^{}S^{}HSR\right)|j,j_z|0_v$$
(26)
If we now transform the Hamiltonian and require that

$$R^{}S^{}HSR|j,j_z|0_v=0$$
(27)
we find the following conditions,

$`\nu (1+\mathrm{cos}2\theta )`$ $`=`$ $`\mu (1\mathrm{cos}2\theta )`$ (28)
$`\mathrm{\Omega }j\mathrm{sin}2\theta `$ $`=`$ $`E`$ (29)
which requires that

$$\mathrm{cos}2\theta =e^{2r}$$
(30)
and the ground state energy is taken to be $`_0=0`$. The ground state is thus a product of a squeezed state for the vibrational mode and a rotated angular momentum state, rotated about the $`\widehat{J}_y`$ axis.

The above results are consistent with the semiclassical approximation. The mean amplitude of a squeezed vacuum state is zero, corresponding to the semiclassical fixed point at $`\overline{X}=\overline{Y}=0`$ while the rotation around the $`\widehat{J}_y`$ axis corresponds to the semiclassical fixed point at $`\overline{𝒥}_x=\chi `$.

If we continue to increase $`E`$ above the threshold value the system adiabatically follows a zero energy state, although this is no longer a ground state. In fact the canonical transformation used in the semiclassical analysis can be applied to the quantum operator valued Hamiltonian. The result is the same as the semiclassical case, Eq (19) with all variables replaced with the corresponding operators. The zero energy state then corresponds to the zero energy eigenstate of $`\widehat{Y}\mathrm{cos}\theta \widehat{X}\mathrm{sin}\theta `$ with $`\mathrm{cos}\theta =N\mathrm{\Omega }/2E`$. This is of course just a rotated, infinitely squeezed state. The electronic state is likewise a angular momentum eigenstate rotated from $`|j,j`$ in the equatorial plane (orthogonal to $`\widehat{J}_z`$). Thus above threshold the zero energy eigenstate deforms continuously from the state at threshold.

Let us summarise these results. For no driving the ground state corresponds to the oscillator in the ground state and all ions in the ground state. As the driving is increased, but kept below threshold, this state deforms to a squeezed oscillator state while the collective spin system begins to rotate about the $`\widehat{J}_y`$ axis. Note that the mean oscillator amplitude $`a`$ remains zero as does the mean of the $`y`$-component of the collective spin. As the driving increases through the threshold value, this state changes its character so that a non zero value of $`\widehat{J}_y`$ is acquired and the oscillator is infinitely squeezed in a direction at an angle $`\mathrm{cos}\theta =N\mathrm{\Omega }/2E`$ to the below threshold squeezing. This morphological change of the state as the driving passes the semiclassical critical point is a quantum phase transition. The quantum phase transition can be seen in the mean value for $`\widehat{J}_y`$ and $`\widehat{J}_z`$ as shown in figure 2. Below threshold the scaled mean values are given by

$`{\displaystyle \frac{\widehat{J}_y}{N/2}}`$ $`=`$ $`0`$ (31)
$`{\displaystyle \frac{\widehat{J}_z}{N/2}}`$ $`=`$ $`\sqrt{1x^2}`$ (32)
and above threshold we have

$`{\displaystyle \frac{\widehat{J}_y}{N/2}}`$ $`=`$ $`\sqrt{1{\displaystyle \frac{1}{x^2}}}`$ (33)
$`{\displaystyle \frac{\widehat{J}_z}{N/2}}`$ $`=`$ $`0`$ (34)
where $`x=2E/N\mathrm{\Omega }`$.

What are the experimental manifestations of this transition ? Needless to say no one is ever going to observe an infinitely squeezed state in an experiment. So what does happens at $`2\theta =\pi /2`$ when the electronic state is the $`\widehat{J}_x`$ eigenstate $`|jj_x`$ and the vibrational mode appears to be infinitely squeezed ? Is such a state physically possible ? Suppose for example we begin in the ground state of the Hamiltonian with no driving ($`E=0`$) which is simply $`|j,j_z|0_v`$, and adiabatically increase the driving strength. It would appear that the system would then adiabatically evolve into the squeezed vibrational state described above. If we were ever able to reach the case $`2\theta =\pi /2`$ we would have reached an infinite energy state for the vibrational mode at a finite driving strength. Clearly this is not possible and to understand why it is useful to reconsider the semiclassical dynamics for this model. The adiabatic approximation requires that we vary the driving strength on a time scale slower than all other time scales in the system. The key time scale for the ground state variation is just the time scale associated with the hyperbolic unstable fixed point, $`(N^2/4\chi ^2)^{1/2}`$, which goes to infinity as we approach $`\frac{2E}{\mathrm{\Omega }N}=1`$. Thus the adiabatic increase of the driving must proceed infinitely slowly, that is it must be switched to the finite value $`E=\frac{\mathrm{\Omega }N}{2}`$ in an infinite amount of time. This pumps an infinite amount of energy into the system and results in infinite squeezing in the centre of mass vibrational mode. Obviously in practice this cannot be achieved so the totally squeezed ground state is not possible. However it will still be possible to achieve some squeezing of the vibrational mode at smaller values of the driving. This would make an interesting observation for current ion trap experiments even with only a few ions. The squeezing of the vibrational mode can be observed using the dynamical method of reference

In current ion trap experiments, laser cooling techniques allow the centre of mass mode to be prepared in the ground state. Unfortunately it does not stay there. Heating due to a variety of sources, including fluctuating linear potentials, lead to an irreversible evolution away from the ground state. If such heating is present during the coupling of the electronic and vibrational motions, irreversible dynamics will be spread to the collective spin degrees of freedom as well.

As an example we consider what happens if we use the Tavis-Cummings interaction (excitation on first red sideband) in the presence of strong heating. Heating of the centre-of-mass mode due to fluctuating liner potentials may be described in the interaction picture by the master equation,

$$\frac{dW}{dt}=i\mathrm{\Omega }[a\widehat{J}_++a^{}\widehat{J}_{},W]+\frac{\gamma }{2}\left(𝒟[a]+𝒟[a^{}]\right)W$$
(35)
where $`W`$ is the density operator for the spin and vibrational degrees of freedom and the superoperator $`𝒟`$ is defined by $`𝒟[A]\rho =2A\rho A^{}A^{}A\rho \rho A^{}A.`$ The irreversible term corresponds to two point processes in which phonons are removed or added from centre of mass mode at the rates $`\gamma a^{}a`$ and $`\gamma aa^{}`$ respectively. This does not change any first order moments, however it does lead to a diffusion in energy as $`\frac{da^{}a}{dt}=\gamma .`$ The effect of heating can be included in the semiclassical analysis by adding an appropriate stochastic term. In the Ito calculus the effect is to add to the equations for $`X,Y`$ terms of the form

$`dX`$ $`=`$ $`(\mathrm{})+\sqrt{\gamma }dW_x(t)`$ (36)
$`dY`$ $`=`$ $`(\mathrm{})+\sqrt{\gamma }dW_y(t)`$ (37)
where $`dW_i(t)`$ are independent Wiener processes. If the heating rate is small enough these terms can be neglected. However if they are large new steady states can occur in the semiclassical and quantum descriptions which will be described in a future publication.

We would like to thank Howard Carmichael for useful discussions.