text:
# Gauge covariant heat kernel expansion at finite temperature

## I Introduction

Finite temperature quantum field theory has gained widespread applications in different areas of physics like e.g. for the description of the early universe or the study of heavy-ion collisions. This work is devoted to a finite temperature generalization of a non-perturbative technique which is very well established for the vacuum case. The heat kernel expansion Schwinger (1951); DeWitt (1965); Seeley (1967) (for a review see e.g. Ball (1989)) is frequently used in the context of the path integral method to integrate out degrees of freedom in a non-perturbative way. The outcome is a derivative expansion of the fields which correspond to the remaining degrees of freedom which were not integrated out. In other words the heat kernel expansion yields an effective field theory (cf. e.g. Ball (1989); Espriu et al. (1990); Gasser and Leutwyler (1984)). A typical application which we also address here is to integrate out fermions from the action of a gauge theory which is bilinear in the fermion fields. The obtained fermionic functional determinant is processed further by the heat kernel expansion. In this way the ultraviolet divergent terms which require renormalization can be properly isolated e.g. using the $`\zeta `$-function method Ball (1989). It may appear that part of the symmetries of the original action are spoiled by the need for renormalization. Such an appearance of anomalies can also be well accounted for with the heat kernel method. In the following we shall not consider anomalies explicitly. We keep the number of dimensions as well as the gauge group completely general. We restrict ourselves, however, to a flat space-time structure.

Naively generalizing the heat kernel expansion to finite temperature we will encounter problems concerning the gauge symmetries as we will see below. This, however, is not a manifestation of an anomaly, as it appears for any gauge symmetry and only at finite temperature. We will solve that problem subsequently by presenting an improved heat kernel expansion which preserves the gauge symmetries (apart from potential anomalies of course). Throughout this work we will use the Matsubara formalism Matsubara (1955); Das (1997) to represent fields at finite temperature. It amounts to a Euclidean metric with the boson (fermion) fields obeying (anti-)periodic boundary conditions in the Euclidean time direction with periodicity $`\beta =1/T`$ where $`T`$ is the temperature of the heat bath the fields are coupled to.

The task to extend the heat kernel expansion to finite temperature has already been addressed in Boschi-Filho et al. (1992). It has been claimed there that the ratio of the finite temperature heat kernel and the vacuum one is just a temperature dependent constant independent of the considered fields. In contradiction to that statement the authors of Actor et al. (1999) presented a counter-example where the mentioned ratio does depend on the gauge fields. It is the purpose of the present work to clarify this contradiction and to present a valid generalization of the heat kernel expansion to the finite temperature case. In vacuum the heat kernel expansion is an expansion in powers of the field amplitudes and their derivatives. For gauge fields this amounts to an expansion in powers of the field strengths and their gauge covariant derivatives. All the other information present in the vector potentials is due to the gauge freedom and therefore of no physical significance. At finite temperature, however, the field strengths and their derivatives are no longer sufficient to characterize the gauge fields. The crucial point is that due to the periodicity requirement mentioned above the Euclidean time direction is compactified. This leads to the appearance of an Aharonov-Bohm effect Aharonov and Bohm (1959). The Polyakov loop Polyakov (1978) supplements the information contained in the field strengths. It is defined by

$$L(x):=P\mathrm{exp}\left[\underset{x_0}{\overset{x_0+\beta }{}}𝑑tA_0(t,\stackrel{}{x})\right]$$
(1)
where $`A_0`$ is the zeroth component of the gauge potential and $`P`$ denotes path ordering. The latter is required for non-abelian gauge theories to obtain a gauge covariant result. Obviously the Polyakov loop constitutes a non-local object in terms of the gauge fields. It is easy to check that a Taylor expansion of $`L`$ in powers of the gauge potential and its derivatives e.g. around the point $`x`$ would yield objects which are ill-defined with respect to gauge transformations. Therefore, the Polyakov loop should be treated on equal footing with the field strength without any further expansion in terms of the amplitudes of the gauge potential. On the other hand, a heat kernel expansion is just an expansion in terms of the field amplitudes. Therefore to prevent an erroneous expansion of the Polyakov loop, the latter has to be isolated right from the start to obtain a heat kernel expansion which is appropriate for the finite temperature case.

The paper is organized in the following way: In the next section we take the usual method to perform a heat kernel expansion and naively generalize it to the finite temperature case. To some extent we follow here the work of Boschi-Filho et al. (1992). At a certain stage, however, we have to point out differences between our findings and the ones of Boschi-Filho et al. (1992). Finally we have to face a problem concerning the gauge covariance of the results obtained with this naive generalization of the heat kernel expansion technique. In Section III we present an improved method to perform the heat kernel expansion at finite temperature preserving gauge covariance in every single step. Gauge covariant Seeley-DeWitt coefficients are calculated in Section IV. It is also shown in this section that the results reduce to the already known results for the vacuum case. In addition, the general results for the Seeley-DeWitt coefficients are applied to the case of two-dimensional QED with constant fields. For this special case the heat kernel and thus also all its Seeley-DeWitt coefficients can be calculated exactly Actor et al. (1999). This provides the possibility to cross-check our results with the ones presented in Actor et al. (1999). Finally Section V summarizes our findings. Three Appendices are added containing technical aspects of the presented formalism.

## II Naive representation for the heat kernel

Consider the differential operator in Euclidean space

$$\mathrm{\Delta }_x=D_\mu ^xD_\mu ^x+X(x)$$
(2)
with the gauge covariant derivative

$$D_\mu ^x=_\mu ^x+A_\mu (x).$$
(3)
Here the functions $`A_\mu `$ and $`X`$ are no operators (but might be matrix-valued due to the group structure of the gauge group under consideration). Under gauge transformations the introduced quantities transform as

$$X(x)V(x)X(x)V^1(x),D_\mu ^xV(x)D_\mu ^xV^1(x).$$
(4)
The heat kernel is given by

$$H(x,x^{};\tau ):=\left(x|e^{\tau \mathrm{\Delta }}|x^{}\right)=e^{\tau \mathrm{\Delta }_x}\delta (xx^{}).$$
(5)
We shall assume that the differential operator (2) originates from a fermionic determinant which we want to evaluate by the heat kernel technique Ball (1989). A generalization of the presented formalism to the bosonic case is in principle straightforward. The representation for the heat kernel on the right hand side of (5) is completely appropriate for the vacuum case. In this section we will generalize this heat kernel to the finite temperature case. We will point out the problem which appears when this heat kernel is used at finite temperature. In the next section we will present an improved representation for the heat kernel which is better suited for the treatment at finite temperature. Hence some of the formulae which we present in this section should not be regarded as valuable results but as illustrations of the problems which appear at finite temperature.

For finite $`T`$ the $`\delta `$-function which appears in (5) is given by

$$\delta (xx^{})=\underset{n=\mathrm{}}{\overset{+\mathrm{}}{}}T\frac{d^{D1}k}{(2\pi )^{D1}}e^{i\stackrel{}{k}(\stackrel{}{x}\stackrel{}{x}^{})}e^{i\omega _n(x_0x_0^{})}=:\frac{d^Dk}{(2\pi )^D}e^{ik(xx^{})}$$
(6)
where $`D`$ is the number of space-time dimensions. Being interested in the calculation of a fermionic determinant via the heat kernel method the discrete energies appearing in (6) are defined as Matsubara (1955); Das (1997)

$$\omega _n=2\pi T\left(n+\frac{1}{2}\right).$$
(7)
It is straightforward to show that the diagonal entries of $`H`$ can be written as Nepomechie (1985)

$$H(x,x;\tau )={\displaystyle \frac{d^Dk}{(2\pi )^D}e^{\tau (k^2+M^2)}\underset{m=0}{\overset{\mathrm{}}{}}\frac{1}{m!}\tau ^m(2ikD_x\mathrm{\Delta }_x+M^2)^mf(x)}|_{f=1}.$$
(8)
We have introduced an arbitrary mass term $`M`$ to keep the formalism as general as possible. All derivatives in (8) act on an arbitrary function $`f`$ which has to be put to 1 at the very end. We want to emphasize already here that the precise meaning of “covariant derivative acting on $`f`$” is important. As $`f`$ basically substitutes a fermionic field the meaning of $`D_\mu ^xf(x)`$ really is $`_\mu ^xf(x)+A_\mu (x)f(x)`$ and not $`[D_\mu ^x,f(x)]`$. Therefore eventually putting $`f`$ to 1, a term like $`A_\mu (x)`$ might remain! This is in contrast to the claim in Boschi-Filho et al. (1992) that $`D_\mu f`$-terms should be dropped in the end (cf. the statement between eqs. (3.17) and (3.18) in Boschi-Filho et al. (1992)). This will have important consequences as we will see below.

Next we introduce generalized Seeley-DeWitt coefficients Seeley (1967); DeWitt (1965); Nepomechie (1985); Ball (1989) as follows:

$$H(x,x;\tau )=\left(\frac{1}{4\pi \tau }\right)^{D/2}e^{\tau M^2}\underset{n=0}{\overset{\mathrm{}}{}}g_n(x;\sqrt{\tau }T)\tau ^n.$$
(9)
Note that at finite temperature there is a remaining $`\tau `$ dependence of the coefficients $`g_n`$ but only via the specific combination $`\sqrt{\tau }T`$ Boschi-Filho et al. (1992). Hence the $`g_n`$’s become independent of $`\tau `$ for the vacuum case $`T=0`$. To calculate e.g. the first two coefficients, $`g_0`$ and $`g_1`$, we have to evaluate the sum over $`m`$ in (8) up to $`m=2`$. We note that all terms with odd numbers of $`k_\mu `$ drop out by the $`k`$ summation/integration. We first calculate

$$c_m:=(2ikD_x\mathrm{\Delta }_x+M^2)^mf(x)|_{f=1}|_{\text{drop odd powers of }k}.$$
(10)
We get

$`c_0`$ $`=`$ $`1,`$ (11a)
$`c_1`$ $`=`$ $`D_\mu A_\mu \stackrel{~}{X},`$ (11b)
$`c_2`$ $`=`$ $`4k_\mu k_\nu D_\mu A_\nu +D_\nu D_\nu D_\mu A_\mu D_\mu D_\mu \stackrel{~}{X}\stackrel{~}{X}D_\mu A_\mu +\stackrel{~}{X}^2,`$ (11c)
where we have introduced

$$\stackrel{~}{X}:=XM^2.$$
(12)
Note that the momenta scale as $`1/\sqrt{\tau }`$ due to the $`\mathrm{exp}(\tau k^2)`$ term in (8). Therefore it is easy to see that $`c_0`$ contributes to $`g_0`$, $`c_1`$ and the first term of $`c_2`$ contribute to $`g_1`$, and the remaining terms of $`c_2`$ contribute to $`g_2`$. Hence

$`g_0`$ $`=`$ $`(4\pi \tau )^{D/2}{\displaystyle \frac{d^Dk}{(2\pi )^D}e^{\tau k^2}c_0}=I_0(\sqrt{\tau }T),`$ (13a)
$`g_1`$ $`=`$ $`(4\pi \tau )^{D/2}{\displaystyle \frac{d^Dk}{(2\pi )^D}e^{\tau k^2}\left(c_1\frac{1}{2}\tau \mathrm{\hspace{0.17em}4}k_\mu k_\nu D_\mu A_\nu \right)}=(D_\mu A_\mu \stackrel{~}{X})I_0(\sqrt{\tau }T)D_0A_0I_{2t}(\sqrt{\tau }T)D_iA_iI_{2s}(\sqrt{\tau }T)`$ (13b)
$`=`$ $`\stackrel{~}{X}I_0(\sqrt{\tau }T)+D_0A_0\left(I_0(\sqrt{\tau }T)I_{2t}(\sqrt{\tau }T)\right),`$
where the relevant sum-integrals are collected in Appendix A.

The important point now is that the term $`D_0A_0`$ does not transform homogeneously under gauge transformations in contrast to $`\stackrel{~}{X}`$.111Note that $`\stackrel{~}{X}`$ transforms in the same way as $`X`$. Therefore this cannot be a proper result for the generalized Seeley-DeWitt coefficient. This problem does not show up at zero temperature since then both integrals $`I_0`$ and $`I_{2t}`$ become identical (cf. Appendix A). Hence the dubious term drops out in this case. We note that this $`D_0A_0`$-term was overlooked in Boschi-Filho et al. (1992). It was still present in eq. (3.15) of Boschi-Filho et al. (1992) but was disregarded afterwards by arguing that $`D_0f`$ vanishes in the end. As already pointed out, this would be true for $`[D_0,f]`$ but not for $`D_0f`$. Ignoring this subtle difference the authors of Boschi-Filho et al. (1992) did not encounter any problem concerning gauge invariance. However, as shown above there is a problem with gauge invariance and actually there should be one on physical grounds: By construction the heat kernel expansion (8) is an expansion in powers of field amplitudes and their derivatives. Any non-local object is expanded in terms of local ones. In vacuum this is a procedure which preserves gauge invariance — apart from well-defined anomalies which should be there on physical grounds. The reason why this procedure works for the vacuum case can be traced back to the fact that the field strengths (and their derivatives) are sufficient to characterize the physical content of the gauge fields. The field strengths, however, are local objects. The situation is completely different at finite temperature as already mentioned in the introduction. In this case the Polyakov loop carries additional information. This object is manifestly non-local. An expansion of the Polyakov loop in powers of the field amplitudes would yield misleading gauge dependent results. Hence the Polyakov loop has to be isolated as an additional degree of freedom before processing the heat kernel expansion.

Having established the physical reason for the erroneous result (13b) we next have to analyze what went wrong on the technical level when proceeding from (5) to (8). Concerning the behavior with respect to gauge transformations we observe that formally the heat kernel (5) transforms as a non-local object,

$$H(x,x^{};\tau )V(x)e^{\tau \mathrm{\Delta }_x}V^1(x)\delta (xx^{})=V(x)e^{\tau \mathrm{\Delta }_x}\delta (xx^{})V^1(x^{})=V(x)H(x,x^{};\tau )V^1(x^{}),$$
(14)
on account of the transformation properties (4). It is illustrative to decompose this transformation property in the following way:

$$e^{\tau \mathrm{\Delta }_x}V(x)e^{\tau \mathrm{\Delta }_x}V^1(x)\text{and}\delta (xx^{})V(x)\delta (xx^{})V^1(x^{}).$$
(15)
Correspondingly the left hand side of (8) transforms gauge covariantly

$$H(x,x;\tau )V(x)H(x,x;\tau )V^1(x).$$
(16)
In contrast, the single Fourier modes on the right hand side of (8) do not show that behavior with respect to gauge transformations. In general, only after summing up all Fourier modes and performing the sum over $`m`$ one obtains the heat kernel with its well-defined transformation properties. For the vacuum case it turns out that already the summation over the Fourier modes is sufficient; each single Seeley-DeWitt coefficient obtained in that way transforms gauge covariantly. As we have demonstrated above this is no longer true at finite temperature. Therefore it is more appropriate to use an approach which preserves in each step as much of the good transformation properties as possible. To be specific, we are looking for a Fourier type representation of the $`\delta `$-function which preserves the transformation property (15) for each single Fourier mode. In this way we will never encounter the problem that only a sum of objects transforms properly while the single terms do not. We stress once more that such an approach is not mandatory for the vacuum case. After the summation of the Fourier modes is performed the results transform properly with respect to gauge transformations.

## III Gauge covariant heat kernel

To obtain gauge covariant results for the Seeley-DeWitt coefficients it might be not mandatory but it is obviously more save to start with a Fourier type representation of the $`\delta `$-function which transforms gauge covariantly in every single Fourier mode. An important ingredient for such an approach is the gauge link (see e.g. Yang (1974); Bralic (1980); Elze et al. (1986))

$$U[A](x;x^{}):=P\mathrm{exp}\left[\underset{x^{}}{\overset{x}{}}𝑑z_\mu A_\mu (z)\right]$$
(17)
where $`P`$ denotes path ordering. For our purposes it is useful to choose the following path parametrized by $`z`$: the path connects with straight lines the sequence of points $`x^{}=(x_0^{},\stackrel{}{x}^{})(\tau _0,\stackrel{}{x}^{})(\tau _0,\stackrel{}{x})x=(x_0,\stackrel{}{x})`$. Here $`\tau _0`$ is an arbitrary time. It will turn out that none of the final results will depend on $`\tau _0`$. From the technical point of view this aspect can be seen as a consistency check of the whole approach. For the chosen path we find

$$U[A](x,x)=1.$$
(18)
Another important property of gauge links is their additivity Yang (1974); Bralic (1980); Elze et al. (1986):

$$U[A](x;x^{})U[A](x^{};x^{\prime \prime })=U[A](x;x^{\prime \prime })$$
(19)
valid for arbitrary $`x^{}`$ on the contour connecting $`x`$ and $`x^{\prime \prime }`$. Finally under gauge transformations $`U`$ transforms as Yang (1974); Bralic (1980); Elze et al. (1986)

$$U[A](x;x^{})V(x)U[A](x;x^{})V^1(x^{}).$$
(20)
Therefore it is tempting to write

$$\delta (xx^{})=\delta (xx^{})U[A](x,x^{})={\displaystyle \frac{d^Dk}{(2\pi )^D}e^{ik(xx^{})}U[A](x,x^{})}.$$
(21)
Indeed, now the modified single Fourier modes transform gauge covariantly, i.e. in the same way as the $`\delta `$-function in (15):

$$e^{ik(xx^{})}U[A](x,x^{})V(x)e^{ik(xx^{})}U[A](x;x^{})V^1(x^{}).$$
(22)
However this cannot be the full story: At finite temperature the $`\delta `$-function which corresponds to fermions must be antiperiodic in Euclidean time with periodicity $`\beta =1/T`$:

$$\delta (x_0+\beta x_0^{})=\delta (x_0x_0^{}).$$
(23)
This is also true for the single Fourier modes

$$e^{ik_0(x_0+\beta x_0^{})}=e^{ik_0(x_0x_0^{})}$$
(24)
on account of (7). The modified Fourier modes displayed on the left hand side of (22) do not share that important feature as the gauge link $`U`$ is not periodic in time:

$$U[A](x_0+\beta ,\stackrel{}{x};x_0^{},\stackrel{}{x}^{})=\underset{=L(x)}{\underset{}{U[A](x_0+\beta ,\stackrel{}{x};x_0,\stackrel{}{x})}}U[A](x_0,\stackrel{}{x};x_0^{},\stackrel{}{x}^{})U[A](x_0,\stackrel{}{x};x_0^{},\stackrel{}{x}^{})$$
(25)
where we have used (19). Note that the term which spoils the periodicity is nothing but the Polyakov loop defined in (1).

We now construct a modified gauge link which both has the transformation property (20) and is periodic in time. The key point is to isolate the (generalized) zero mode of the time component of the gauge potential which is closely related to the Polyakov loop. We first define the generalized zero mode

$$a_0(x):=T\mathrm{log}L(x)=T\mathrm{log}U[A](x_0+\beta ,\stackrel{}{x};x_0,\stackrel{}{x}).$$
(26)
The properties of $`a_0`$ are collected in Appendix B. In general already the log-function of complex numbers is multiple valued. This is of course also true for the log-function of matrices. In the following we use the log-function defined via its Taylor expansion (cf. Appendix B)

$$\mathrm{log}(1+𝒪):=\underset{n=1}{\overset{\mathrm{}}{}}(1)^{n+1}\frac{1}{n}𝒪^n.$$
(27)
Hence we restrict ourselves to Polyakov loops with values close to the unity matrix. We postpone the further discussion of that issue to Section V. We introduce the modified gauge link

$$\stackrel{~}{U}[A](x;x^{}):=e^{x_0a_0(x)}U[A](x;x^{})e^{x_0^{}a_0(x^{})}.$$
(28)
Using the properties of $`a_0`$ collected in Appendix B it is straightforward to prove the following properties of $`\stackrel{~}{U}`$:

$`\stackrel{~}{U}[A](x;x)`$ $`=`$ $`1,`$ (29a)
$`\stackrel{~}{U}[A](x_0+\beta ,\stackrel{}{x};x_0^{},\stackrel{}{x}^{})`$ $`=`$ $`\stackrel{~}{U}[A](x_0,\stackrel{}{x};x_0^{},\stackrel{}{x}^{}),`$ (29b)
$`\stackrel{~}{U}[A](x;x^{})`$ $``$ $`V(x)\stackrel{~}{U}[A](x;x^{})V^1(x^{}).`$ (29c)
With these developments at hand the heat kernel can be represented as

$$H(x,x^{};\tau )=e^{\tau \mathrm{\Delta }_x}\left(\delta (xx^{})\stackrel{~}{U}[A](x;x^{})\right).$$
(30)
For the diagonal part one gets

$$H(x,x;\tau )={\displaystyle \frac{d^Dk}{(2\pi )^D}e^{\tau (k^2+M^2)}\underset{m=0}{\overset{\mathrm{}}{}}\frac{1}{m!}\tau ^m(2ikD_x\mathrm{\Delta }_x+M^2)^m\stackrel{~}{U}[A](x;x^{})}|_{x^{}=x}.$$
(31)
## IV Gauge covariant Seeley-DeWitt coefficients

From the representation (31) it is straightforward to calculate the generalized Seeley-DeWitt coefficients defined in (9). We present in the following the calculation of the first three generalized Seeley-DeWitt coefficients. For this task we introduce

$$\stackrel{~}{c}_m:=(2ikD_x\mathrm{\Delta }_x+M^2)^m\stackrel{~}{U}[A](x;x^{})|_{x^{}=x}|_{\text{drop odd powers of }k}.$$
(32)
Note that by construction every $`\stackrel{~}{c}_m`$ transforms gauge covariantly in contrast to $`c_m`$ as defined in (10). For the lowest order terms we get

$`\stackrel{~}{c}_0`$ $`=`$ $`1,`$ (33a)
$`\stackrel{~}{c}_1`$ $`=`$ $`(D_\mu D_\mu \stackrel{~}{X})\stackrel{~}{U}[A](x;x^{})|_{x^{}=x},`$ (33b)
$`\stackrel{~}{c}_2`$ $`=`$ $`(4k_\mu k_\nu D_\mu D_\nu +D_\mu D_\mu D_\nu D_\nu \stackrel{~}{X}D_\mu D_\mu D_\mu D_\mu \stackrel{~}{X}+\stackrel{~}{X}^2)\stackrel{~}{U}[A](x;x^{})|_{x^{}=x},`$ (33c)
$`\stackrel{~}{c}_3`$ $`=`$ $`4k_\mu k_\nu (D_\mu D_\nu D_\alpha D_\alpha D_\mu D_\alpha D_\alpha D_\nu D_\alpha D_\alpha D_\mu D_\nu +D_\mu D_\nu \stackrel{~}{X}+D_\mu \stackrel{~}{X}D_\nu +\stackrel{~}{X}D_\mu D_\nu )\stackrel{~}{U}[A](x;x^{})|_{x^{}=x}`$ (33d)
$`+\text{terms which contribute to }g_3,`$
$`\stackrel{~}{c}_4`$ $`=`$ $`16k_\mu k_\nu k_\alpha k_\beta D_\mu D_\nu D_\alpha D_\beta \stackrel{~}{U}[A](x;x^{})|_{x^{}=x}+\text{terms which contribute to }g_3\text{}g_4.`$ (33e)
It is lengthy but straightforward to determine the corresponding generalized Seeley-DeWitt coefficients. In the course of the calculation there appears the following type of $`n`$-fold derivatives:

$$D_{x,0}^n\stackrel{~}{U}[A](x;x^{})|_{x^{}=x}=a_0^n(x).$$
(34)
The proof for this relation is given in Appendix C. The formulae collected in Appendix A serve to perform the $`k_i`$ integrations and to rewrite the finite temperature sums. The lowest order generalized Seeley-DeWitt coefficients turn out to be

$`g_0`$ $`=`$ $`I_0(\sqrt{\tau }T),`$ (35a)
$`g_1`$ $`=`$ $`\stackrel{~}{X}I_0(\sqrt{\tau }T)+a_0^2\left(I_0(\sqrt{\tau }T)I_{2t}(\sqrt{\tau }T)\right),`$ (35b)
$`g_2`$ $`=`$ $`{\displaystyle \frac{1}{2}}\stackrel{~}{X}^2I_0(\sqrt{\tau }T){\displaystyle \frac{1}{6}}[D_i,[D_i,\stackrel{~}{X}]]I_0(\sqrt{\tau }T){\displaystyle \frac{1}{6}}[D_0,[D_0,\stackrel{~}{X}]]\left(3I_0(\sqrt{\tau }T)2I_{2t}(\sqrt{\tau }T)\right)`$ (35c)
$`\left([D_0,\stackrel{~}{X}]a_0+\stackrel{~}{X}a_0^2\right)\left(I_0(\sqrt{\tau }T)I_{2t}(\sqrt{\tau }T)\right)+{\displaystyle \frac{1}{12}}[D_i,D_j]^2I_0(\sqrt{\tau }T)+{\displaystyle \frac{1}{6}}[D_i,D_0]^2\left(2I_0(\sqrt{\tau }T)I_{2t}(\sqrt{\tau }T)\right)`$
$`+{\displaystyle \frac{1}{6}}[D_i,[D_0,[D_i,D_0]]]\left(I_0(\sqrt{\tau }T)I_{2t}(\sqrt{\tau }T)\right)+{\displaystyle \frac{1}{3}}[D_i,[D_i,D_0]]a_0\left(I_0(\sqrt{\tau }T)I_{2t}(\sqrt{\tau }T)\right)`$
$`+{\displaystyle \frac{1}{2}}a_0^4\left(I_0(\sqrt{\tau }T)2I_{2t}(\sqrt{\tau }T)+I_{4t}(\sqrt{\tau }T)\right).`$
This constitutes the main result of the present work. Of course, by construction every coefficient transforms gauge covariantly. Obviously the method can be extended to calculate higher order coefficients. The expressions for $`I_{\mathrm{}}`$ are collected in Appendix A. Note that the coefficients $`g_n`$ do not depend on the number of space-time dimensions $`D`$. The dependence of $`g_n`$ on the temperature is non-analytic as can be seen by inspecting (55) and (56). This leads to a very smooth transition of the vacuum to finite temperature results. In the following we shall apply these general results to two special cases. The purpose is to check the consistency of our results.

In vacuum all $`I_{\mathrm{}}`$-functions become unity (see Appendix A). One obtains the standard results for this limiting case Ball (1989):

$`g_0^{\mathrm{vac}}`$ $`=`$ $`1,`$ (36a)
$`g_1^{\mathrm{vac}}`$ $`=`$ $`\stackrel{~}{X},`$ (36b)
$`g_2^{\mathrm{vac}}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\stackrel{~}{X}^2{\displaystyle \frac{1}{6}}[D_i,[D_i,\stackrel{~}{X}]]{\displaystyle \frac{1}{6}}[D_0,[D_0,\stackrel{~}{X}]]+{\displaystyle \frac{1}{12}}[D_i,D_j]^2+{\displaystyle \frac{1}{6}}[D_i,D_0]^2`$ (36c)
$`=`$ $`{\displaystyle \frac{1}{2}}\stackrel{~}{X}^2{\displaystyle \frac{1}{6}}[D_\mu ,[D_\mu ,\stackrel{~}{X}]]+{\displaystyle \frac{1}{12}}[D_\mu ,D_\nu ]^2.`$
As a further check of our results (35) we compare them to the special case discussed in Actor et al. (1999): two-dimensional massless QED with a constant electric field strength $`E`$ and in addition a constant zero mode. The (trace of the) heat kernel for such a system is calculated exactly in Actor et al. (1999):

$$\mathrm{tr}H(x,x;\tau )=\frac{E}{2\pi }\frac{1}{\mathrm{tanh}E\tau }\left\{1+2\underset{n=1}{\overset{\mathrm{}}{}}(1)^n\mathrm{cos}\left[n\frac{Ex_1+2\pi Ta}{T}\right]\mathrm{exp}\left(\frac{n^2E}{4T^2\mathrm{tanh}E\tau }\right)\right\}.$$
(37)
Below we will expand the exact result to obtain the corresponding generalized Seeley-DeWitt coefficients. The latter can be compared to our results which we will apply now to the case at hand.

The vector potentials are given by

$$A_0=i(Ex_1+2\pi Ta),A_1=0.$$
(38)
For the zero mode defined in (26) we find222The already mentioned ambiguity concerning the logarithm will be discussed below.

$$a_0=A_0=i(Ex_1+2\pi Ta).$$
(39)
The task is to calculate (approximately) the heat kernel for the operator

$$\overline{)}D^2=D_\mu D_\mu X$$
(40)
with

$$X=\frac{1}{4}[\gamma _\mu ,\gamma _\nu ][D_\mu ,D_\nu ]=E\sigma _3$$
(41)
where we have used the convention Actor et al. (1999)

$$\gamma _0=\sigma _1,\gamma _1=\sigma _2$$
(42)
with the Pauli matrices $`\sigma _i`$.

For this case the generalized Seeley-DeWitt coefficients are (of course we use $`M=0`$):

$`g_0`$ $`=`$ $`I_0(\sqrt{\tau }T),`$ (43a)
$`g_1`$ $`=`$ $`E\sigma _3I_0(\sqrt{\tau }T)(Ex_1+2\pi Ta)^2\left(I_0(\sqrt{\tau }T)I_{2t}(\sqrt{\tau }T)\right),`$ (43b)
$`g_2`$ $`=`$ $`E\sigma _3(Ex_1+2\pi Ta)^2\left(I_0(\sqrt{\tau }T)I_{2t}(\sqrt{\tau }T)\right)+{\displaystyle \frac{1}{6}}E^2\left(I_0(\sqrt{\tau }T)+I_{2t}(\sqrt{\tau }T)\right)`$ (43c)
$`+{\displaystyle \frac{1}{2}}(Ex_1+2\pi Ta)^4\left(I_0(\sqrt{\tau }T)2I_{2t}(\sqrt{\tau }T)+I_{4t}(\sqrt{\tau }T)\right).`$
The expansion of (37) has to be performed in the following way: Expand in powers of $`\tau `$, but with $`T^2\tau `$ fixed. There is, however, one exception to this rule, namely the $`2\pi Ta`$ term. The temperature appearing there does not originate from integrating out the thermal fermions but was put in by hand in Actor et al. (1999) to make $`a`$ dimensionless. To account for that fact we use (39) and rewrite (37) as

$$\mathrm{tr}H(x,x;\tau )=\frac{1}{2\pi \tau }\frac{E\tau }{\mathrm{tanh}E\tau }\left\{1+2\underset{n=1}{\overset{\mathrm{}}{}}(1)^n\mathrm{cos}\left[\sqrt{\tau }n\frac{ia_0}{\sqrt{\tau }T}\right]\mathrm{exp}\left(\frac{n^2}{4\tau T^2}\frac{E\tau }{\mathrm{tanh}E\tau }\right)\right\}.$$
(44)
Now we apply our expansion rule to (44):

$`\mathrm{tr}H(x,x;\tau )`$ $``$ $`{\displaystyle \frac{1}{2\pi \tau }}\left(1+{\displaystyle \frac{1}{3}}\tau ^2E^2\right)`$ (45)
$`\times \left\{1+2{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}(1)^n\left(1+\tau {\displaystyle \frac{n^2a_0^2}{2\tau T^2}}+\tau ^2{\displaystyle \frac{n^4a_0^4}{24\tau ^2T^4}}\right)\left(1{\displaystyle \frac{1}{3}}\tau ^2E^2{\displaystyle \frac{n^2}{4\tau T^2}}\right)\mathrm{exp}\left({\displaystyle \frac{n^2}{4\tau T^2}}\right)\right\}`$
$``$ $`{\displaystyle \frac{1}{2\pi \tau }}\left(1+{\displaystyle \frac{1}{3}}\tau ^2E^2\right)\left\{I_0+\left(\tau a_0^2{\displaystyle \frac{1}{6}}\tau ^2E^2\right)(I_0I_{2t})+{\displaystyle \frac{1}{2}}\tau ^2a_0^4(I_02I_{2t}+I_{4t})\right\}`$
$``$ $`{\displaystyle \frac{1}{2\pi \tau }}\left[\left(1+{\displaystyle \frac{1}{3}}\tau ^2E^2\right)I_0+\left(\tau a_0^2{\displaystyle \frac{1}{6}}\tau ^2E^2\right)(I_0I_{2t})+{\displaystyle \frac{1}{2}}\tau ^2a_0^4(I_02I_{2t}+I_{4t})\right]`$
where we have used the relations (55) and (56) to rewrite the finite temperature sums. Finally recalling (9) we can read off the (trace of the) generalized Seeley-DeWitt coefficients:

$`\mathrm{tr}g_0`$ $`=`$ $`2I_0,`$ (46a)
$`\mathrm{tr}g_1`$ $`=`$ $`2a_0^2(I_0I_{2t}),`$ (46b)
$`\mathrm{tr}g_2`$ $`=`$ $`{\displaystyle \frac{1}{3}}E^2(I_0+I_{2t})+a_0^4(I_02I_{2t}+I_{4t}).`$ (46c)
These results are in complete agreement with (43).

## V Summary

We have presented a generalization of the heat kernel expansion technique to finite temperature. This task has already been envisaged in Boschi-Filho et al. (1992). Technically we followed the approach of Boschi-Filho et al. (1992) to some extent. There are deviations, however, concerning the formal as well as the physics point of view - which are of course intertwined. Physically we expect to find (at least) two aspects which make the finite temperature case distinct from the vacuum case. First, we expect the appearance of terms which “break” Lorentz invariance like e.g. a different coefficient in front of the electric field strength (squared) as compared to the magnetic field strength. Such differences have not been found in Boschi-Filho et al. (1992) whereas our calculation shows that they are indeed present (cf. e.g. the fifth and sixth term in (35c)). Second, as already pointed out in the introduction at finite temperature the field strengths are not sufficient to characterize the gauge fields. Due to the compactness of the Euclidean time direction the Polyakov loop carries additional information — a manifestation of the Aharonov-Bohm effect. Hence we expect the appearance of a corresponding object in the heat kernel expansion. Indeed in our approach the generalized zero mode which is the logarithm of the Polyakov loop appears. Again this object does not show up in the approach of Boschi-Filho et al. (1992).

Concerning the technical aspects of our approach as compared to Boschi-Filho et al. (1992) we have taken care of the correct meaning of the gauge covariant derivative acting eventually on the arbitrary function $`f`$. The latter appears in Nepomechie’s constructive method to obtain the heat kernel expansion; cf. (8). Only the correct treatment of the covariant derivative revealed the problem of the naively calculated Seeley-DeWitt coefficients concerning gauge covariance. Subsequently we developed an expansion method which preserves gauge covariance in every single step. While this is not necessary for the vacuum case — albeit it does not hurt — it turned out to be very useful for the finite temperature case. The requirement to simultaneously preserve gauge covariance and (anti-)periodicity led to the introduction of the generalized zero mode. In this way the Polyakov loop entered our formalism.

Formally the results of Boschi-Filho et al. (1992) could be reproduced if we replaced $`I_{2t}`$ and $`I_{4t}`$ by $`I_0`$ in (35). Obviously a large number of terms would then drop out. These terms have not been considered in Boschi-Filho et al. (1992). By inspecting (55) it is clear that such a replacement is only correct at zero temperature. Consequently the terms which are neglected in Boschi-Filho et al. (1992) are the ones which are not Lorentz invariant. However, the appearance of such terms is expected as the heat bath “breaks” Lorentz invariance. On the other hand, rotational invariance is still preserved. Obviously this requirement is fulfilled by our result (35).

Further support for the soundness of our results comes from Actor et al. (1999) where the case of two-dimensional QED at finite temperature with constant fields is addressed. The possible appearance of a zero mode which cannot be gauged away at finite temperature constitutes the discrepancy between Boschi-Filho et al. (1992) and Actor et al. (1999) as noted in the latter reference. As shown above our results are in complete agreement with the special case of two-dimensional QED with constant fields discussed in Actor et al. (1999). One might regard our results as a generalization of the results of Actor et al. (1999) to arbitrary dimensions and non-constant fields. On the other hand, in contrast to Actor et al. (1999) we cannot offer an exact solution here but only an expansion in terms of field amplitudes and their derivatives.

Concerning gauge fields at finite temperature the crucial point is that this expansion is not only an expansion in powers of the field strength but in addition in powers of the generalized zero mode. As already emphasized in Actor et al. (1999) even a constant zero mode cannot be gauged away at finite temperature. It is important to stress one aspect which is different in our approach as compared to the exact result presented in Actor et al. (1999). Inspecting (37) reveals a periodicity of the heat kernel with respect to

$$aa+1.$$
(47)
Our approximate result (43) seems not to reproduce that behavior of the exact solution. Of course, an expansion in powers of the zero mode amplitude might not be appropriate any more for large amplitudes. On the other hand, due to the periodicity property (47) the zero mode can be made arbitrarily large. Actually this problem is just another manifestation of the fact that the logarithm introduced in (26) is in principle multiple valued. Indeed the Polyakov loop which corresponds to the problem studied in Actor et al. (1999) is given by

$$L(x)=e^{i(\beta Ex_1+2\pi a)}.$$
(48)
Obviously at this stage the periodicity property is still present. The problem appears with the calculation of the zero mode according to (26) which apparently yields (39). This, however, is not fully true: If the logarithm is defined via its Taylor expansion (27) the outcome of (26) is

$$a_0(x)=i[Ex_1+2\pi T(am)]$$
(49)
where $`m`$ is an integer chosen such that $`\varphi :=|\beta Ex_1+2\pi (am)|`$ is minimal. And even this is only true if $`\varphi `$ is small enough, i.e. $`|e^{i\varphi }1|<1`$. Otherwise the Taylor expansion (27) for the logarithm would not converge. With the correct result (49) for the zero mode we see that the periodicity property (47) is still present. On the other hand, by considering the convergence properties of the Taylor expansion for the logarithm it becomes obvious that for an expansion of the heat kernel small amplitudes are mandatory right from the start. Strictly speaking so far our formalism can only be applied to Polyakov loops with values close to unity. Of course, nothing is wrong with an expansion in powers of the field amplitudes which works when the amplitudes are small but breaks down when they are large. On the other hand, it clearly would be desirable to develop a heat kernel expansion in which the Polyakov loop enters directly and not only via the generalized zero mode. At the present stage it is not clear whether this is possible.

## Appendix A Useful sum-integrals

The following sum-integrals show up in the calculation of the lowest order terms of the heat kernel expansion:

$`(4\pi \tau )^{D/2}{\displaystyle \frac{d^Dk}{(2\pi )^D}e^{\tau k^2}}`$ $`=:`$ $`I_0(\sqrt{\tau }T),`$ (50a)
$`(4\pi \tau )^{D/2}{\displaystyle \frac{d^Dk}{(2\pi )^D}e^{\tau k^2}\mathrm{\hspace{0.17em}2}\tau k_0^2}`$ $`=:`$ $`I_{2t}(\sqrt{\tau }T),`$ (50b)
$`(4\pi \tau )^{D/2}{\displaystyle \frac{d^Dk}{(2\pi )^D}e^{\tau k^2}\mathrm{\hspace{0.17em}2}\tau k_ik_j}`$ $`=:`$ $`I_{2s}(\sqrt{\tau }T)\delta _{ij},`$ (50c)
$`(4\pi \tau )^{D/2}{\displaystyle \frac{d^Dk}{(2\pi )^D}e^{\tau k^2}\frac{4}{3}\tau ^2k_0^4}`$ $`=:`$ $`I_{4t}(\sqrt{\tau }T),`$ (50d)
$`(4\pi \tau )^{D/2}{\displaystyle \frac{d^Dk}{(2\pi )^D}e^{\tau k^2}\mathrm{\hspace{0.17em}4}\tau ^2k_0^2k_ik_j}`$ $`=:`$ $`I_{4m}(\sqrt{\tau }T)\delta _{ij},`$ (50e)
$`(4\pi \tau )^{D/2}{\displaystyle \frac{d^Dk}{(2\pi )^D}e^{\tau k^2}\mathrm{\hspace{0.17em}4}\tau ^2k_ik_jk_lk_m}`$ $`=:`$ $`I_{4s}(\sqrt{\tau }T)\left(\delta _{ij}\delta _{lm}+\delta _{il}\delta _{jm}+\delta _{im}\delta _{jl}\right).`$ (50f)
Note that sum-integrals with an odd number of $`k`$ factors vanish. The spatial $`k`$-integrations can be readily evaluated using

$`{\displaystyle \frac{d^Dk}{\pi ^D}e^{\stackrel{}{k}^2}}`$ $`=`$ $`1`$ (51a)
$`{\displaystyle \frac{d^Dk}{\pi ^D}e^{\stackrel{}{k}^2}k_ik_j}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\delta _{ij},`$ (51b)
$`{\displaystyle \frac{d^Dk}{\pi ^D}e^{\stackrel{}{k}^2}k_ik_jk_lk_m}`$ $`=`$ $`{\displaystyle \frac{1}{4}}\left(\delta _{ij}\delta _{lm}+\delta _{il}\delta _{jm}+\delta _{im}\delta _{jl}\right).`$ (51c)
One gets

$`I_0(\sqrt{\tau }T)=I_{2s}(\sqrt{\tau }T)=I_{4s}(\sqrt{\tau }T)`$ $`=`$ $`(4\pi \tau )^{1/2}T{\displaystyle \underset{n=\mathrm{}}{\overset{+\mathrm{}}{}}}\mathrm{exp}(\tau \omega _n^2),`$ (52a)
$`I_{2t}(\sqrt{\tau }T)=I_{4m}(\sqrt{\tau }T)`$ $`=`$ $`(4\pi \tau )^{1/2}T{\displaystyle \underset{n=\mathrm{}}{\overset{+\mathrm{}}{}}}\mathrm{exp}(\tau \omega _n^2)\mathrm{\hspace{0.17em}2}\tau \omega _n^2,`$ (52b)
$`I_{4t}(\sqrt{\tau }T)`$ $`=`$ $`(4\pi \tau )^{1/2}T{\displaystyle \underset{n=\mathrm{}}{\overset{+\mathrm{}}{}}}\mathrm{exp}(\tau \omega _n^2){\displaystyle \frac{4}{3}}\tau ^2\omega _n^4.`$ (52c)
These finite-temperature sums can be rewritten using a generalized theta function transformation. We adopt formula (3.11) from Boschi-Filho et al. (1992) (with a misprint corrected):

$$\underset{n=\mathrm{}}{\overset{+\mathrm{}}{}}\mathrm{exp}\left(\gamma (n+1/2)^2\right)\left[\gamma (n+1/2)^2\right]^p=\left(\frac{\pi }{\gamma }\right)^{1/2}\frac{(2p1)!!}{2^p}\left[1+2\underset{n=1}{\overset{\mathrm{}}{}}(1)^n\mathrm{exp}(\pi ^2n^2/\gamma )\underset{l=0}{\overset{p}{}}\left(\frac{\pi ^2n^2}{\gamma }\right)^lC_l^p\right]$$
(53)
with

$$C_l^p:=\frac{2^lp!}{l!(pl)!(2l1)!!},$$
(54)
$`p`$ a non-negative integer and $`(2p1)!!:=13\mathrm{}(2p1)`$. We apply (53) to (52) with $`\gamma =(2\pi T)^2\tau `$:

$`I_0`$ $`=`$ $`1+2{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}(1)^n\mathrm{exp}\left({\displaystyle \frac{n^2}{4\tau T^2}}\right),`$ (55a)
$`I_{2t}`$ $`=`$ $`1+2{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}(1)^n\mathrm{exp}\left({\displaystyle \frac{n^2}{4\tau T^2}}\right)\left(12{\displaystyle \frac{n^2}{4\tau T^2}}\right),`$ (55b)
$`I_{4t}`$ $`=`$ $`1+2{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}(1)^n\mathrm{exp}\left({\displaystyle \frac{n^2}{4\tau T^2}}\right)\left[14{\displaystyle \frac{n^2}{4\tau T^2}}+{\displaystyle \frac{4}{3}}\left({\displaystyle \frac{n^2}{4\tau T^2}}\right)^2\right].`$ (55c)
We also present some combinations of $`I_{\mathrm{}}`$-functions which appear in (35):

$`I_0I_{2t}`$ $`=`$ $`{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}(1)^n\mathrm{exp}\left({\displaystyle \frac{n^2}{4\tau T^2}}\right){\displaystyle \frac{n^2}{\tau T^2}},`$ (56a)
$`I_02I_{2t}+I_{4t}`$ $`=`$ $`{\displaystyle \frac{1}{6}}{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}(1)^n\mathrm{exp}\left({\displaystyle \frac{n^2}{4\tau T^2}}\right)\left({\displaystyle \frac{n^2}{\tau T^2}}\right)^2.`$ (56b)
## Appendix B Properties of the generalized zero mode

In this appendix we collect and prove some important properties for the generalized zero mode defined in (26). We first note that for an abelian gauge theory $`a_0`$ is just the zero mode of the Fourier decomposition of the time component of the vector potential:

$$A_0(x_0,\stackrel{}{x})=\underset{n=\mathrm{}}{\overset{+\mathrm{}}{}}a_n(\stackrel{}{x})e^{i2n\pi Tx_0}.$$
(57)
For a non-abelian gauge theory such an easy relation cannot be obtained. The abelian case, however, justifies our name choice “generalized zero mode”. The relations presented in the following hold for arbitrary gauge theories.

Periodicity:

$$a_0(x_0+\beta ,\stackrel{}{x})=T\mathrm{log}P\mathrm{exp}\left[\underset{x_0+\beta }{\overset{x_0+2\beta }{}}𝑑tA_0(t,\stackrel{}{x})\right]=T\mathrm{log}P\mathrm{exp}\left[\underset{x_0}{\overset{x_0+\beta }{}}𝑑t\underset{=A_0(t,\stackrel{}{x})}{\underset{}{A_0(t+\beta ,\stackrel{}{x})}}\right]=a_0(x_0,\stackrel{}{x}).$$
(58)
Behavior under gauge transformations:
We restrict ourselves to gauge transformations which are periodic in $`x_0`$. These are the only ones appropriate for the fermion fields which are integrated out via the heat kernel method. On account of (20) we find

$$U[A](x_0+\beta ,\stackrel{}{x};x_0,\stackrel{}{x})V(x_0+\beta ,\stackrel{}{x})U[A](x_0+\beta ,\stackrel{}{x};x_0,\stackrel{}{x})V^1(x_0,\stackrel{}{x})=V(x_0,\stackrel{}{x})U[A](x_0+\beta ,\stackrel{}{x};x_0,\stackrel{}{x})V^1(x_0,\stackrel{}{x}).$$
(59)
To obtain the transformation property of $`a_0`$ we decompose $`U`$ as

$$U[A](x_0+\beta ,\stackrel{}{x};x_0,\stackrel{}{x})=:1+𝒪$$
(60)
and use a Taylor expansion of the log-function:

$`\mathrm{log}U[A](x_0+\beta ,\stackrel{}{x};x_0,\stackrel{}{x})`$ $``$ $`\mathrm{log}(VUV^1)=\mathrm{log}(V(1+𝒪)V^1)=\mathrm{log}(1+V𝒪V^1)={\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}(1)^{n+1}{\displaystyle \frac{1}{n}}(V𝒪V^1)^n`$ (61)
$`=V\mathrm{log}(1+𝒪)V^1=V(\mathrm{log}U)V^1.`$
Thus

$$a_0(x)V(x)a_0(x)V^1(x).$$
(62)
Finally one can also show using the Taylor expansion of the exponential function:

$$\mathrm{exp}\left(\alpha a_0(x)\right)\mathrm{exp}\left(\alpha V(x)a_0(x)V^1(x)\right)=V(x)\mathrm{exp}\left(\alpha a_0(x)\right)V^1(x)$$
(63)
for arbitrary numbers $`\alpha `$.

## Appendix C Proof of equation (34)

In Elze et al. (1986) the derivative of a straight line gauge link $`U_S[A](b,a)`$ with respect to its endpoints $`a(x)`$, $`b(x)`$ has been considered. One gets

$`_\mu ^xU_S[A](b,a)`$ $`=`$ $`{\displaystyle \frac{b_\nu }{x_\mu }}\left[A_\nu (b)(ba)_\alpha {\displaystyle \underset{0}{\overset{1}{}}}𝑑ssU_S[A](b,z)F_{\alpha \nu }(z)U_S[A](z,b)\right]U_S[A](b,a)`$ (64)
$`+U_S[A](b,a)\left[A_\nu (a)(ba)_\alpha {\displaystyle \underset{0}{\overset{1}{}}}𝑑s(s1)U_S[A](a,z)F_{\alpha \nu }(z)U_S[A](z,a)\right]{\displaystyle \frac{a_\nu }{x_\mu }}`$
with the straight line path connecting $`a`$ and $`b`$,

$$z=z(s)=a+s(ba),$$
(65)
and the field strength

$$F_{\mu \nu }:=[D_\mu ,D_\nu ].$$
(66)
This implies for the gauge links of interest

$$[D_0^x,U[A](x_0+\beta ,\stackrel{}{x};x_0,\stackrel{}{x})]=0$$
(67)
and

$$D_0^xU[A](x;x^{})=U[A](x;x^{})_0^x.$$
(68)
We make use of the decomposition of $`U`$ given in (60). Using (67) we find

$$[D_0,𝒪]=0$$
(69)
and therefore

$$[D_0,𝒪^n]=0.$$
(70)
This yields

$$[D_0,a_0]=T[D_0,\mathrm{log}(1+𝒪)]=T\underset{n=1}{\overset{\mathrm{}}{}}(1)^{n+1}\frac{1}{n}[D_0,𝒪^n]=0.$$
(71)
Hence we obtain

$`[D_0,\mathrm{exp}\left(x_0a_0(x)\right)]`$ $`=`$ $`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{n!}}[D_0,x_0^na_0^n]={\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{n!}}[D_0,x_0^n]a_0^n+{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{n!}}x_0^n[D_0,a_0^n]={\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{n!}}nx_0^{n1}a_0^n`$ (72)
$`=`$ $`a_0\mathrm{exp}\left(x_0a_0(x)\right).`$
Finally we get

$`D_0^x\stackrel{~}{U}[A](x;x^{})`$ $`=`$ $`D_0^xe^{x_0a_0(x)}U[A](x;x^{})e^{x_0^{}a_0(x^{})}`$ (73)
$`=`$ $`[D_0^x,e^{x_0a_0(x)}]U[A](x;x^{})e^{x_0^{}a_0(x^{})}+e^{x_0a_0(x)}D_0^xU[A](x;x^{})e^{x_0^{}a_0(x^{})}`$
$`=`$ $`a_0(x)\stackrel{~}{U}[A](x;x^{})`$
where we have used (68) in the last step. Using (73), (71) and (29a) leads to (34).


format: text

text:
# The High-Energy Polarization-Limiting Radius of Neutron Star Magnetospheres I – Slowly Rotating Neutron Stars

## 1 Introduction

The thermal radiation of isolated NS stars has the potential of teaching us much about the properties of NSs. Its advantage over non-thermal emission (in radio, optical, X-rays, and gamma-rays) is that the theory behind the emission is significantly better understood and the radiation actually comes from the surface of the compact object. Because the thermal emission is expected to be intrinsically polarized, more information could potentially be learned by the detection and analysis of polarization measurements. Recent observations with the ROSAT, ASCA, Chandra and XMM-Newton missions have shown that some of these sources are bright enough to be potential candidates for X-ray polarimetry in future missions. Moreover, the thermal radiation of some of the isolated NSs can even be detected in optical wavelengths. Thus, it is worthwhile understanding how this polarization is generated and conserved, and what additional information can actually be extracted from its measurement.

In the presence of strong magnetic fields, the opacity of ionized matter to the transfer of photons becomes polarization dependent (Lodenqual et al. 1974 ; Ventura 1979 ). This is chiefly because it is easier to scatter electrons in the direction along the magnetic field than it is in a perpendicular direction. Thus, the opacity of light rays with their electric vector polarized perpendicular to the magnetic field would be significantly reduced. A typical photon with this polarization is emitted or scattered last deeper in the atmosphere than one in the other mode. The deeper regions of the atmosphere are hotter, so more flux emerges in this polarization state. Nearly complete polarization can result for the thermal emission (Kanno 1975 ; Pavlov & Shibanov 1978 ; Ventura, Nagel & Meszaros 1979 ).

An observer will see photons originating from the entire surface of the NS hemisphere facing her. Thus, the different polarizations should be added together appropriately. If nothing happens to the photons and their polarization as they propagate from the NS surface, then the polarizations observed at infinity can be added rather simply, as was done by Kanno 1975  in a simple model for the atmosphere. Pavlov & Zavlin 2000  used a more realistic atmosphere and calculated the net observed polarization while taking the effects that GR has on the magnetic field and on light ray bending. In both cases, net polarizations of order 5% to 30% are obtained because the polarizations of the radiation arriving from different regions of the surface tends to cancel each other.

The above analyses, however, did not consider the effects of QED-induced vacuum birefringence. When QED is coupled to strong magnetic fields, several interesting consequences are obtained. For example, it was shown by Pavlov & Shibanov 1978  and Meszaros & Ventura 1979  that QED has to be taken into account when calculating the appropriate opacity, especially near the cyclotron resonance. This is the case even though a priori it appears that the plasma effects should dominate. More relevant to us is the fact that QED turns the vacuum into a birefringent medium when strong magnetic fields are present (Tsai & Erber 1975 ; Heyl & Hernquist 1997 ).

In a series of recent papers, we have examined several consequences of vacuum birefringence in neutron-star magnetospheres. When the fields are significantly stronger than the critical QED field of $`B_{QED}=4.4\times 10^{13}`$ G then the index of refraction of one polarization state can be significantly different from unity and magnetic lensing can result (Shaviv, Heyl & Lithwick 1999 ). The main result of this lensing effect is that the effective surface area of the NS as measured by the two polarization states is different.

When weaker fields are present, the birefringence can still have interesting implications. At a particular frequency, the vacuum will only decouple the polarization modes out to a particular distance from the surface of the star. Up to this radius, radiation polarized perpendicular to the magnetic field will remain perpendicular to the local direction of the magnetic field even if the direction of the field changes along the path. If the modes only begin to mix at a significant fraction of the distance to the light cylinder, even if the intrinsic polarization at the surface is constant over energy, photons of different energies will exhibit different directions of polarization after passing through the magnetosphere (Heyl & Shaviv 2000 ) and a circular component of the polarization will develop. Similar effects arise at lower frequencies, when plasma birefringence is considered (Cheng & Ruderman 1979 ; Barnard 1986 ).

In Heyl & Shaviv 2002 , we showed that QED birefringence is also important for the polarization evolution close to the NS. When it is properly taken into account, a very large net polarization should be observed. This is counter to previous predictions (e.g. Pavlov & Zavlin 2000 ). As a result, larger polarization signals will be observable which will allow more information to be extracted from the observation of the thermal radiation. Measurement of the high polarization will also serve as the first direct evidence of the birefringence of the magnetized vacuum due to QED and a direct probe of the behavior of the vacuum at magnetic fields of order of and above the critical QED field of $`B_{QED}=4.4\times 10^{13}`$ G. This should be contrasted with the decades of Earth based experiments which have not succeeded thus far in detecting the vacuum birefringence induced by strong magnetic fields (Iacoponi & Zavattini 1979 ; Bakalov et al. 1998 ; Rizzo 1998 ; Nezrick 1999 ).

We begin in §2 with the description of the physics needed to calculate the polarization to be observed at infinity. In §3, we elaborate the results presented in Heyl & Shaviv 2002  and build upon them to understand the observational signatures of vacuum polarization in rotating neutron stars. §3.3.2 estimates the strength of the polarized signal averaged over the rotation of the star for the subset of radio pulsars for which we know the geometry of the dipole field. We end in §4 with a discussion of the ramifications of this effect.

## 2 Calculations

Several ingredients are needed for the calculation of the net polarization to be observed at infinity. First, the structure of the magnetic field must be specified. We assume that the magnetic field is a centered dipole. Second, we need a model for the intrinsic polarization emitted by a magnetized atmosphere. For simplicity, we assume here that the atmospheres emit completely (linearly) polarized radiation. In a large frequency range, it is more than an adequate approximation because the effective temperature for the two polarizations will be markedly different if high magnetic fields are present. For example, at photon energies $`E_\gamma `$ much below the electron rest energy and cyclotron energy $`E_{cyc,e}`$, but much above the ion cyclotron energy, the typical degree of linear polarization $`p_L`$ should be $`1p_L𝒪(E_\gamma /E_{\text{cyc,e}})^2`$, unless the angle between the magnetic field $`𝐁`$ and the photon wavevector $`𝐤`$ is very small (Pavlov & Zavlin 2000 ).

To calculate the observed polarization at infinity, we need to calculate the trajectories of the light rays, which due to GR are bent. Along these trajectories, we have to solve for the evolution of the polarization. This will be dominated by the vacuum birefringence.

### 2.1 Photon Trajectories

In calculating the trajectories of the photons we assume that the field is sufficiently weak such that the index of refraction for both modes is approximately unity throughout the magnetosphere (c.f. Shaviv, Heyl & Lithwick 1999 ). We also neglect the effect of the rotation of the star on the spacetime surrounding it.

Without rotation, all planes that pass through the center of the star are equivalent, so we can integrate the equations of motion for a photon in the equatorial plane of the Schwarzschild metric. The trajectory is determined uniquely by the impact parameter $`b`$. To integrate the polarization, we require the position of the photon as a function of the proper length along its path. Misner, Thorne & Wheeler 1973  give the differential equations for the trajectory

$`{\displaystyle \frac{\mathrm{d}t}{\mathrm{d}\lambda }}`$ $`=`$ $`\left(1{\displaystyle \frac{2M}{r}}\right)^1`$ (1)
$`{\displaystyle \frac{\mathrm{d}r}{\mathrm{d}\lambda }}`$ $`=`$ $`\left[1{\displaystyle \frac{b^2}{r^2}}\left(1{\displaystyle \frac{2M}{r}}\right)\right]^{1/2}`$ (2)
$`{\displaystyle \frac{\mathrm{d}\theta }{\mathrm{d}\lambda }}`$ $`=`$ $`0`$ (3)
$`{\displaystyle \frac{\mathrm{d}\varphi }{\mathrm{d}\lambda }}`$ $`=`$ $`{\displaystyle \frac{b}{r^2}},`$ (4)
and also

$$\frac{\mathrm{d}l}{\mathrm{d}\lambda }=\left(1\frac{2M}{r}\right)^{1/2}$$
(5)
where $`M`$ is the gravitational mass of the star. Combining Equation (2) and Equation (4) yields a separable equation for $`\varphi (r)`$ which is useful for quickly determining to which part of the star a region of the image corresponds (Page 1995 )

$$\varphi \varphi _0=x_0^{M/R}\left[\left(1\frac{2M}{R}\right)\left(\frac{M}{R}\right)^2\left(12u\right)u^2x^2\right]^{1/2}du,$$
(6)
where $`xb/R_{\mathrm{}}`$ and $`R_{\mathrm{}}R\left(12M/R\right)^{1/2}`$ and $`R`$ is the circumferential radius of the star.

### 2.2 Polarization Trajectories

The polarization lies in the plane perpendicular to the trajectory of the photon. Unlike in flat spacetime, because the photon travels along a curved path, the orientation of this plane with respect to a distant observer necessarily varies along the path. Under the assumption that space surrounding the neutron star is devoid of material and nongravitational fields, the polarization is constant, if it is defined in a basis consisting a vector in the plane of the trajectory and one perpendicular to that plane (Pineault 1977 ).

To naturally include the standard general relativistic result, we calculate the evolution of the polarization due to the birefringence of the vacuum in the aforementioned basis. Figure 1 shows an example trajectory along with the polarization basis at a particular point. $`\alpha `$ is the angle between the magnetic dipole ($`\stackrel{}{m}`$) and the line of sight ($`\stackrel{}{O}`$), and $`\beta `$ is the angle between the trajectory plane and the $`\stackrel{}{m}\stackrel{}{O}`$plane.

Kubo & Nagata 1983  find that the evolution of the polarization of a wave traveling through a birefringent and dichroic medium in the limit of geometric optics is given by

$$\frac{𝐬}{l}=\widehat{𝛀}\times 𝐬+\left(\widehat{𝐓}\times 𝐬\right)\times 𝐬,$$
(7)
where $`l`$ is the proper distance along the trajectory, $`𝐬`$ is the normalized Stokes vector (Jackson 1975 ), and $`\widehat{𝛀}`$ and $`\widehat{𝐓}`$ are the birefringent and dichroic vectors. The Stokes vector consists of the four Stokes parameters, $`S_0,S_1,S_2`$ and $`S_3`$. The vector $`𝐬`$ consists of $`S_1/S_0,S_2/S_0`$ and $`S_3/S_0`$. The result was found by Kubo & Nagata 1983  for any dielectric medium and it was extended for a medium that is both dielectric and permeable by Heyl & Shaviv 2000 .

As Heyl & Shaviv 2000  argue, QED decouples the polarization states in the vacuum for sufficiently strong fields. Here we will restrict ourselves to fields substantially less than $`B_{\text{QED}}4.4\times 10^{13}`$ G. A field of $`10^{12}`$ G is sufficient to decouple the polarization states at the surface of a neutron star for $`\nu 10^{12}`$ Hz. A plasma with the Goldreich-Julian density decouples the polarization states of photons with $`\nu 10^{14}`$ Hz. Here we will focus on ultraviolet through X-ray radiation, so the plasma contribution to the index of refraction may be neglected. We will consider photons with $`\nu 10^{14}`$ Hz (i.e. for which the plasma would be important) but we will also neglect the plasma contribution so we can connect our results with those of Pavlov & Zavlin 2000  who neglect the birefringence of the magnetosphere entirely.

If one neglects the plasma and takes the weak-field limit, the dichroic vector vanishes and the magnitude of the birefringent vector is

$$\left|\widehat{𝛀}\right|=\frac{2}{15}\frac{\alpha _{\text{QED}}}{4\pi }\frac{\omega }{c}\left(\frac{B_{}}{B_{\text{QED}}}\right)^2,$$
(8)
and it points in the direction of the projection of the magnetic field onto the Poincaré sphere. $`B_{}`$ is the strength of the magnetic field perpendicular to the direction of the photon’s propagation. We assume that the magnetic field is a centered dipole and we neglect the distortion of the magnetic field due to general relativity. For the masses and radii that we are considering, the perturbation to the field strength is at most a factor of two and the change in the direction of the field is less than 5 throughout (Ginzburg & Ozernoĭ 1965 ). Both the value of the polarization-limiting radius and the emergent flux depend weakly on the strength of the magnetic field at the surface – both increase as $`B^{0.4}`$; therefore, this simplification does not have an important effect on the results.

Heyl & Shaviv 2000  find that the polarization states are decoupled as long as the gradient of the index of refraction is not too large, specifically

$$\left|\widehat{𝛀}\left(\frac{1}{|\widehat{𝛀}|}\frac{|\widehat{𝛀}|}{l}\right)^1\right|0.5.$$
(9)
For radial trajectories this yields

$$rr_{\text{pl}}\left(\frac{\alpha _{\text{QED}}}{45}\frac{\nu }{c}\right)^{1/5}\left(\frac{\mu }{B_{\text{QED}}}\mathrm{sin}\alpha \right)^{2/5}1.2\times 10^7\mu _{30}^{2/5}\nu _{17}^{1/5}\left(\mathrm{sin}\alpha \right)^{2/5}\text{cm}.$$
(10)
where $`\mu `$ is the magnetic dipole moment of the star and $`\nu `$ is the frequency of the photon. Also, $`\mu _{30}=\mu /(10^{30}\text{ G cm}^3)`$ and $`\nu _{17}=\nu /10^{17}\text{Hz}`$.

If this condition is not met, the polarization remains constant, i.e. the polarization modes are coupled. Figure 2 depicts the radii within which either the plasma or the vacuum effectively decouples the polarization modes in the magnetosphere – we have assumed that the plasma density is given by the Goldreich & Julian 1969  result. At distances closer to the star than the polarization-limiting radius ($`r_{\text{pl}}`$), the polarization of the radiation remains in one of the two polarization modes of the strongly magnetized plasma or vacuum.

If the entire surface emits in one polarization mode, i.e. the surface emission is initially fully polarized, the radiation will remain in that mode until the polarization-limiting radius, so one can estimate the observed extent of the polarization geometrically by calculating the solid angle subtended by the image of the surface of the star at the polarization-limiting radius. As the angular size of the image at $`r_{\text{pl}}`$ vanishes, the polarized fraction approaches unity. Numerical integration of the photon paths bears this out.

## 3 Results

These results assume that $`BB_{QED}`$ at the decoupling radius. This is always the case for the photon energies of interest, because even if the field on the surface is greater than the QED value, as is the case in magnetars, the decoupling takes place far enough from the surface such that the $`r^3`$ term will make the field significantly sub-critical.

### 3.1 Effects along a trajectory

We begin by integrating the photon light trajectories for specific photons leaving the NS surface and following the evolution of their polarization. Because the QED vacuum is not dichroic by itself, eq. (7) dictates that the amplitude of $`𝐬`$ does not change along a ray. In the course of evolution, however, the linear component of $`𝐬`$, i.e. , the 1-2 components, may change direction, and the amount of circular component $`s_3`$ can change as well. The top part of figure 3 depicts the evolution of the angle of $`𝐬`$ in the 1-2 plane of the Poincaré Space, together with the angle of the birefringent vector $`\widehat{𝛀}`$ (determined by the magnetic field component perpendicular to the ray).

As the particular photon leaves the surface, the magnetic field orientation rotates by about 2.2 radians which corresponds to 4.4 radians in the 1-2 plane of the Poincaré Space. Because the coupling is weaker at lower frequencies, the lower frequency photons follow the direction of the magnetic field up to a smaller distance. Beyond the polarization-limiting distance, the polarization direction freezes. Its direction roughly corresponds to the direction of the magnetic field where the modes couple. Because modes couple gradually, the direction of the magnetic field can change during the coupling process. If the change in direction is rapid, a large circular polarization results. This is seen in the bottom part of figure 3. For the low frequency and high frequency photons, the coupling takes place before the magnetic field can significantly change or after it has stopped changing, so for these frequencies, the circular component obtained is small. For the intermediate frequency, which for $`10^{12}`$ G NSs corresponds to the optical or UV region, the coupling takes place while the direction of the magnetic field is changing and a large circular component is generated.

At high frequencies (hard X-ray to $`\gamma `$-rays), the coupling can take place at a significant fraction of the light cylinder radius. In this case, it was shown by Heyl & Shaviv 2000  that one should take into account the rotation of the NS. Because the modes of higher frequencies couple farther from the NS, the directions of the rotating magnetic field at the polarization-limiting radii are different for different photon energies, such that phase leads between different wave bands can result. Because coupling takes place while the NS is rotating, circular polarization can again result.

### 3.2 The Polarization “Image” of a NS

The next step is to construct a polarization “image” of a NS. The apparent surface is projected onto a surface perpendicular to the NS-observer direction. The image is then divided into elements of equal solid angle. Next, light rays are followed from each element of the apparent surface to the observer taking into account GR light bending (eq. 6) and polarization evolution (eq. 7), as is described in §3.1.

Typical results are portrayed in figure 4, which shows the polarization observed at infinity overlaying the GR lensed image of the NS. The typical polarization is an ellipse. The major axis describes the direction of the linear polarization (in real space) while the ratio of the minor to major axis gives the amount of circular polarization. Not given in the figure is the sense of rotation of the circular polarization. From symmetry, one obtains that the circular components in the top half of the images are opposite to those in the bottom half. The same is true for the $`s_2`$ polarization which describes the $`\pm 45^{}`$ polarization directions in real space. Unlike the $`s_3`$ antisymmetry however, the $`s_2`$ antisymmetry is apparent in the images.

The anti-symmetry of the $`s_2`$ and $`s_3`$ components implies that when the polarizations from all the star will be added together, only a net $`s_1`$ component will result. Namely, the net polarization from the NS will be either in the direction of the magnetic dipole axis or perpendicular to it. This statement will not be true if the cylindrical symmetry is broken, either by the magnetic field, by rotation or by the atmospheric emission.

### 3.3 The Net Polarization of a NS

Once a polarization “image” is calculated, the net polarization seen by an observer is found by integrating the intensity contributed by each of the normal modes of the atmosphere to each of the observed Stokes’s parameters. Here we will treat the simple case where the intensity in one mode vanishes ($`I_O=0`$), and the intensity in the other mode is isotropic ($`I_X=`$ constant). In this case, the value of $`S_1/S_0`$ is simply the mean value of $`s_1`$ evaluated over the observed polarization field (e.g. as depicted in Figure 4). For this simple model, we denote $`S_1/S_0`$ summed over the entire image by $`\overline{s}_1`$.

The results for $`\overline{s}_1`$ as a function of the magnetic field strength and frequency $`\nu \mu _{30}^2`$ is given in the left panel of figure 5, for two inclination angles and three different NS radii, 6 km, 10 km and 18 km. The right panel depicts the net polarization $`\overline{s}_1`$ as a function of the angle between the line of sight and the magnetic dipole moment for the three NS radii and two frequencies.

Figure 5 depicts several important trends:

1. Higher frequency radiation is more strongly polarized.
2. More strongly magnetized stars exhibit stronger polarization.
3. As the line of sight approaches the direction of the dipole, the net polarization vanishes.
4. At high frequencies, the emission from larger stars is less polarized. The trend is reversed at low frequencies.

Equation (10) predicts the first three of these trends directly. As the frequency of the photon or the strength of the dipole moment increase, the polarization-limiting radius increases. Also as $`\mathrm{sin}\alpha `$ increases, the polarization-limiting radius increases. A larger value of $`r_{\text{pl}}`$ results in a larger polarized fraction because the solid angle subtended by the bundle of rays that eventually reach the detector decreases with distance from the star. Over successively smaller solid angles, the magnetic field geometry appears successively more uniform, and the polarization from different regions of the star is added more coherently.

The final trend requires a two-part explanation. If $`\nu \mu _{30}^2\mathrm{sin}^2\alpha 10^{12}`$ Hz, then $`r_{\text{pl}}R`$, so the net polarization depends almost entirely on the angular size of the ray bundle. Far from the surface of the star, the linear radius of the bundle is $`b_{\text{max}}`$. For a given mass, $`b_{\text{max}}`$ decreases with the radius of the star until it reaches a constant value for $`R<3M`$; consequently, smaller stars have smaller bundles and larger net polarizations.

For $`\nu \mu _{30}^2\mathrm{sin}^2\alpha 10^{12}`$ Hz, the polarization-limiting radius is comparable or smaller than the radius of the star. In this regime, magnetospheric birefringence has little effect on the polarized image; therefore, in this regime, the results of Pavlov & Zavlin 2000  are obtained. We see a larger fraction of the surface of more compact stars so the net polarization will decrease as $`M/R`$ increases, because the polarization is then added mostly incoherently. For $`R<3.5M`$ we see the entire surface and for $`R3M`$ we see an infinite number of images of the surface (e.g. Page 1995 ).

The paragraphs that follow examine the ramifications of these trends in more detail. In particular, we calculate of the polarized light curve of a neutron star and its average. We will continue in a subsequent publication with predictions of the net polarization for a realistic model of the emission from the surface of a neutron star.

#### 3.3.1 Polarization light curve of a NS

When an observer measures the polarization of a rotating NS, the amount and angle of polarization will generally vary because the rotation axis and magnetic axis are usually not aligned together. We define the angular separation between the magnetic and rotational axes as $`\gamma `$. The magnetic inclination angle $`i`$ ($`\pi /2\alpha `$) can be related to the inclination above the rotational equator $`i_r`$ and the rotational phase $`\varphi `$ between the last time the two axes coincided in the observer’s meridional plane. The relation is

$$\mathrm{sin}i=\mathrm{sin}i_r\mathrm{cos}\gamma +\mathrm{cos}i_r\mathrm{sin}\gamma \mathrm{cos}\varphi .$$
(11)
If we work in a coordinate system aligned with the rotational $`z`$-axis, and a $`y`$-axis that is perpendicular to the plane containing the line of sight and the $`z`$-axis, then the observer’s direction is:

$$\widehat{𝐨}=\mathrm{cos}i_r\widehat{𝐱}+\mathrm{sin}i_r\widehat{𝐳}.$$
(12)
If we use the rotational phase $`\varphi `$ and the separation $`\gamma `$ between the two axes, the direction of the magnetic axis is:

$$\widehat{𝐦}=\mathrm{sin}\gamma \mathrm{cos}\varphi \widehat{𝐱}+\mathrm{sin}\gamma \mathrm{sin}\varphi \widehat{𝐲}+\mathrm{cos}\gamma \widehat{𝐳}.$$
(13)
Using these relations, we can calculate the cosine and sine of twice the apparent angle $`\psi `$ that the magnetic axis makes with the $`y`$-axis. These are needed if we wish to known the direction of linear polarization. To do so, we project the polarization state $`|S_1`$, in which the net polarization will be in (e.g., §3.2) onto the polarization states $`|S_{O,1}`$ and $`|S_{O,2}`$ of the observer. $`|S_{O,1}`$ describes linear polarization in the observers $`y`$-axis and $`|S_{O,2}`$ describes polarization in a direction rotated by $`45^{}`$. The projections are

$`S_{O,1}|S_1`$ $`=`$ $`\mathrm{cos}2\psi ={\displaystyle \frac{2\left((\widehat{𝐦}(\widehat{𝐦}\widehat{𝐨})\widehat{𝐨})\widehat{𝐲}\right)^2}{|\widehat{𝐦}(\widehat{𝐦}\widehat{𝐨})\widehat{𝐨})|^2}}1`$ (14)
$`=`$ $`{\displaystyle \frac{2\left(\mathrm{sin}\gamma \mathrm{sin}\varphi \right)^2}{1\left(\mathrm{cos}\gamma \mathrm{sin}i_r+\mathrm{cos}i_r\mathrm{sin}\gamma \mathrm{cos}\varphi \right)^2}}1p_1(\gamma ,i_r,\varphi ),`$ (15)
and

$`S_{O,2}|S_1`$ $`=`$ $`\mathrm{sin}2\psi ={\displaystyle \frac{2\left((\widehat{𝐦}(\widehat{𝐦}\widehat{𝐨})\widehat{𝐨})\widehat{𝐲}\right)\left((\widehat{𝐦}(\widehat{𝐦}\widehat{𝐨})\widehat{𝐨})\widehat{𝐱}\right)}{|\widehat{𝐦}(\widehat{𝐦}\widehat{𝐨})\widehat{𝐨})|^2}}`$ (16)
$`=`$ $`{\displaystyle \frac{2\mathrm{sin}\gamma \mathrm{sin}i_r\mathrm{sin}\varphi \left(\mathrm{cos}\varphi \mathrm{sin}\gamma \mathrm{sin}i_r\mathrm{cos}\gamma \mathrm{cos}i_r\right)}{1\left(\mathrm{cos}\gamma \mathrm{sin}i_r+\mathrm{cos}i_r\mathrm{sin}\gamma \mathrm{cos}\varphi \right)^2}}p_2(\gamma ,i_r,\varphi ).`$ (17)
If the total net polarization observed at infinity $`|\overline{s}|`$ is calculated, the functions $`p_{1,2}(\gamma ,i_r,\varphi )`$ can now be used to find the actual polarization that will be measured if our polarizers are aligned with the $`z`$ axis or $`45^{}`$ to it respectively.

Results for the “polarization light-curve” for several frequencies and several NS angles are depicted in figures 6-8. The solid curve is the total polarization. It could be calculated if the measurement yields both $`S_{O,1}`$ and $`S_{O,2}`$. If however a polarimeter just measures one axis, and it is aligned with $`y`$ then its measurements will follow the dashed line. If rotated by $`45^{}`$, the data will follow the dotted line. Clearly, the best possible measurement is that of the time behavior of the polarization in two axes which yields the large net polarization $`|\overline{S}|`$ and the various angles in the system. The latter include the angle separating the axes $`\gamma `$, the observer’s inclination above the rotational equator $`i_r`$, as well as the direction of the rotational axis in the sky.

#### 3.3.2 Time averaged polarization of a NS

Although the best polarization measurement possible should be time resolved, often is it hard to do so. It is easier to measure the polarization averaged over the spin period. By looking at figures 6-8, we see that if our polarimeter is aligned with the rotational axis, a net polarization signal is obtained though it is typically significantly smaller than the absolute polarization. If the polarimeter is rotated by $`45^{}`$, symmetry dictates a null average signal.

The average polarization depends on $`\gamma `$ – the separation between the axis, and $`i_r`$, the inclination above the rotational equator. Although for the general population of NS, the two should not be correlated, this is not the case if we wish to study NSs for which their geometry is already known. NS geometry is known for some pulsars from linear polarization swing measurements in the radio. Because the objects have to be pulsars with beams passing close to the line of sight, there is a selection effect which chooses objects with only $`i_r\pi /2\gamma `$.

Figure 9 describes the average expected polarization for the different pulsars for which $`\gamma `$ and $`i_r`$ are known. The data are taken from Mitra & Deshpande 1999 . We find that the expected time averaged polarization for the thermal radiation of $`10^{12}`$ G type pulsars is going to be on average 5-7 times larger if QED effects are properly taken into account and measurement is done in the optical or X-ray. One example is PSR 0656+14 which has a measured thermal spectrum (Pavlov, Welty & Cordova 1997 ). The prediction is that the time average of its polarization is going to be about 25% times the typical intrinsic polarization of an average surface element. This should be compared with a 5% prediction times the typical average intrinsic polarization, if polarization dragging and aligning does not take place. Note that even if the surface elements were to emit completely polarized radiation, then the maximum possible time average polarization that can be obtained for any geometrical configuration is 12.5% if QED is neglected with typical values being significantly smaller. Thus, the time average measurement of the polarization of PSR 0656+14 is sufficient to prove the effects of QED on aligning the polarization. Because the pulsar was detected in optical and UV, polarimetry can in principle already be done with very long observations.

We must be careful not to overstate the observability of this effect in the optical and near ultraviolet. Typically, if thermal emission from the surface of the star dominates in the optical, the sources are exceptionally faint and would require approximately one night of observing time on a ten-meter-class telescope to detect the intrinsic polarization of the source. Furthermore, contamination by non-thermal emission is typically important. For example, even in the optical $`30\%`$ of the emission from PSR 0656+14 is non-thermal (Pavlov, Welty & Cordova 1997 ). Disentangling these two emission mechanisms in the optical is difficult but possible in principle. In the X-rays the signal is much stronger and non-thermal emission plays a lesser role. However, we do not now have instruments measure the polarization of X-ray radiation from astrophysical sources.

## 4 Discussion & Summary

It is well known that the intrinsic polarization of the thermal radiation emanating from any NS surface element should be highly polarized. This is a direct result of the effects that the magnetic field has on photon propagation. However, it was thought until recently that because each surface element has a different magnetic field orientation, the combined emission for all the different surface elements would result with a low net polarization for the integrated light. This conclusion, however, rests on the assumption that nothing special happens to the polarization angles along the way. Here, we have shown that QED does have a major effect on the polarization angles in the magnetosphere.

In the presence of QED, the vacuum becomes birefringent. If the differences in the indices of refraction are large enough, ‘adiabatic evolution’ of the polarization will evolve each polarization state separately up to a distance $`r_{\text{pl}}`$, the polarization-limiting radius. If the states change slowly because the magnetic field orientation changes, the direction of polarization will change as well. This phenomenon is known to be important for radio waves due to plasma birefringence (Cheng & Ruderman 1979 ). The main differences between the two effects are first that plasma birefringence becomes progressively more important for long wavelengths, as opposed to the shorter wavelengths in which vacuum birefringence becomes progressively more important. Second, the amount of actual plasma birefringence is hard to predict accurately, because the amount of plasma present varies according to the pulsar model adopted. Vacuum birefringence depends only on the magnetic field of the NS.

If the polarization limiting radius $`r_{\text{pl}}`$ is far from the surface of the NS, the adiabatic evolution arising from QED birefringence has a very interesting effect—it aligns the polarization angles such that large net polarizations are obtained. The further from the surface that the coupling of the modes takes place (where adiabatic evolution fails), the better is the alignment of rays originating from different surface elements. For typical magnetic fields of $`10^{12}`$ G, the alignment is already important for optical and UV photons. And it should be almost complete in X-rays. In stronger fields, as are predicted to exist on magnetars, the alignment should be almost complete even in the optical and the polarization would be very high. This should be compared with the predictions neglecting polarization alignment which always result with significantly lower polarizations.

If the magnetic and rotational axes are misaligned, as is generally the case, the direction of linear polarization changes with the rotation phase. As a result, the time-averaged measurements generally yield smaller net polarizations than time-resolved measurements. The latter are therefore much more preferable. However, they require a more elaborate measuring technique which for the very faint thermal signal of PSRs is highly nontrivial.

Generally, a circular polarization component along one ray arises when the polarization limiting radius is not orders of magnitude larger than the radius of the NS. This is due to the fact that while coupling of the states occur at $`rr_{\text{pl}}`$ the magnetic field is changing its orientation relative to the ray. However, when summed over the image, the circular component vanishes by symmetry. Therefore, any measurement of a non vanishing circular component in the thermal radiation would imply that the system has broken its symmetry between the apparent sides of the magnetic axis. This can happen for example if the magnetic field has a non symmetric component (rotation, higher multipoles, offset dipole, etc.). It can happen if the temperature (and therefore emission) is not only a function of magnetic latitude (e.g., if there are ‘hot spots’). It can also arise because a rotating NS will Doppler boost the radiation from one apparent side of the rotation axis to the blue and the other side to the red. A circular component was also shown to arise when taking the effects that rotation have on the decoupling process itself (Heyl & Shaviv 2000 ). The main difference between the two types of circular components is that the latter type increases with frequency, while the circular component that arises from asymmetries is largest for optical or UV (for $`10^{12}`$ G), when the polarization limiting radius is comparable to the radius of the NS.

Polarization measurements of the thermal radiation will clearly be very beneficial. First, the measurement of polarization will verify the birefringence induced by a magnetic field predicted by QED. Magnetic vacuum birefringence has not yet been detected. Moreover, measurement of polarization often elucidates the geometry of the systems. In this case however, it could also give information on actual physical parameters. For example, if the magnetic dipole moment $`\mu `$ is known (e.g. from spin-down rate measurement) then $`R`$ can be extracted from the polarization measurement which indicate how much alignment has taken place. The more alignment observed, the smaller the radius has to be because the NS apparent solid angle at $`r_{\text{pl}}`$ is then smaller. In AXPs, it could be used to verify their extreme magnetic nature.

Support for this work was provided by the National Aeronautics and Space Administration through Chandra Postdoctoral Fellowship Award Number PF0-10015 issued by the Chandra X-ray Observatory Center, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of NASA under contract NAS8-39073. NS wishes to thank CITA for the fellowship which supported him.


text:
# Effect of non–Gaussian noise sources in a noise induced transition

## 1 INTRODUCTION

During the last three decades a wealth of research results on fluctuations or noise have lead us to the recognition that in many situations noise can actually play a constructive role that induces new ordering phenomena. Some examples are stochastic resonance in zero–dimensional  and extended systems , noise induced transitions , noise induced phase transitions , noise induced transport , noise sustained patterns , etc.

Most of the studies on the noise induced phenomena indicated above assume that the noise source has a Gaussian distribution (either white or colored). In addition to the intrinsic interest in the study of non–Gaussian noises, there are some experimental evidences, particularly in sensory and biological systems , indicating that at least in some of these phenomena the noise sources could be non–Gaussian. The use of non–Gaussian noises in studies on noise induced phenomena is scarce mainly due to the mathematical difficulties. This is in contrast with the existence of some analytical tools when working with Gaussian (particularly white) noises.

Here we present some results on one of those noise induced phenomena when driven by a noise source taken as colored and non–Gaussian. It corresponds to a noise induced transition like those discussed in . The problem we discuss here corresponds the so called genetic model which follows from a particular biological modeling or from a chemical reaction . In Ref.  the effect of a colored Gaussian (Ornstein–Uhlenbeck) noise on that model was analyzed in detail and a novel reentrance phenomenon in the phase diagram was found. In order to study the effect that the non–Gaussian nature of the noise can have in the transition, we have used here a particular distribution  whose departure from the Gaussian behavior is governed by a parameter $`q`$: while $`q=1`$ corresponds to Gaussian noise, the value $`q>1`$ produces a long tail distribution and $`q<1`$ gives a cut-off distribution. As discussed for other noise induced phenomena in , we have found that the departure from the Gaussian distribution for the noise produces a strong effect. In this case such an effect corresponds to an anticipation of the transition when the distribution has long tails, and to a retardation when it is cut-off. This effect could have interesting consequences in the case of noise-induced-phase transitions .

In the next section we briefly discuss the form and properties of the non–Gaussian noise source. After that, we present the model we analyze here and a simple analytical approximation to treat the non–Gaussian noise that allows us to determine some of the general features of the phase diagram. We continue discussing numerical simulations that support and extend the theoretical results. In the last section we draw some conclusions.

## 2 NON–GAUSSIAN NOISE AND ITS PROPERTIES

We start considering the following general form of a Langevin equation

$$\dot{x}=f(x)+g(x)\eta (t).$$
(1)
As described in the Introduction, we assume that the noise term $`\eta (t)`$ has a non–Gaussian distribution. More precisely, we consider that $`\eta (t)`$ is a Markovian process generated as the solution of the following Langevin equation

$$\dot{\eta }=\frac{1}{\tau }\frac{d}{d\eta }V_q(\eta )+\frac{\sqrt{2D}}{\tau }\xi (t),$$
(2)
being $`\xi (t)`$ a standard Gaussian white noise of zero mean and correlation $`\xi (t)\xi (t^{})=\delta (tt^{})`$, and

$$V_q(\eta )=\frac{D}{\tau (q1)}\mathrm{ln}\left[1+\frac{\tau }{D}(q1)\frac{\eta ^2}{2}\right].$$
(3)
Although we believe that our results are quite general, such a particular form for the noise $`\eta (t)`$ allows us to easily control the departure from the Gaussian behavior by changing a single parameter $`q`$. $`D`$ and $`\tau `$ are noise parameters related to the noise intensity and the correlation time, as we now detail. The stationary properties of the noise $`\eta `$, including the time correlation function, have been studied in  and here we summarize the main results. The stationary probability distribution is given by

$$P_q^{st}(\eta )=\frac{1}{Z_q}\left[1+\frac{\tau }{D}(q1)\frac{\eta ^2}{2}\right]^{\frac{1}{q1}},$$
(4)
where $`Z_q`$ is the normalization factor. This distribution can be normalized only for $`q<3`$. The first moment, $`\eta =0`$, is always equal to zero, and the second moment,

$$\eta ^2=\frac{2D}{\tau (53q)},$$
(5)
is finite only for $`q<5/3`$. Furthermore, for $`q<1`$, the distribution has a cut-off and it is only defined for $`|\eta |<\eta _c\sqrt{\frac{2D}{\tau (1q)}}`$. Finally, the correlation time $`\tau _q`$ of the stationary regime of the process $`\eta (t)`$ diverges near $`q=5/3`$ and it can be approximated over the whole range of values of $`q`$ as $`\tau _q2\tau /(53q)`$. Clearly, when $`q1`$ we recover the limit of $`\eta `$ being a Gaussian colored noise, namely the Ornstein–Uhlenbeck process, $`\xi _{ou}(t)`$, with correlations $`\xi _{ou}(t)\xi _{ou}(t^{})=\frac{D}{\tau }\mathrm{e}^{|tt^{}|/\tau }`$ and probability distribution $`P^{st}(\xi _{ou})=Z^1\mathrm{e}^{\frac{\tau }{2D}\xi _{ou}^2}`$.

In Ref. , an effective Markovian approximation to the process $`\eta `$ via a path integral procedure was obtained, that is analogous to the “unified colored noise approximation” obtained in  for the case of Ornstein–Uhlenbeck noise. Such an approximation allowed us to get quasi-analytical results for the mean-first-passage time or transition rate. These results and its dependence on the different parameters in the case of a double well potential were compared with extensive numerical simulations with excellent agreement. That approximation was also exploited to study stochastic resonance with theoretical results  that are in very good agreement with experimental ones .

## 3 $`q1`$ APPROXIMATION FOR THE NON–GAUSSIAN NOISE

The “effective Markovian approximation” introduced in Ref.  and indicated above has been shown to have some drawbacks. As discussed in , such an approximation renders a non normalizable stationary probability distribution function, $`P_q^{st}(\eta )`$. This seems to be a general characteristic of such an approximation that, at least in the case discussed in , was solved introducing an adequate and controlled cutoff. However, it is not clear that such a methodology could be always used, and particularly for the models we can analyze in relation to noise induced transitions. Hence, for the kind of models we use in order to study the effect of non–Gaussian noises on a noise induced transition we have resorted to a more simple, albeit physically reasonable, approximation. Our aim is at least to obtain some of the main trends in the region $`|q1|1`$ (both for $`q<1`$ and $`q>1`$). The idea is that in such a region the non–Gaussian noise will only slightly depart from the Gaussian behavior. Hence in Eq. (2) we adopt

$`{\displaystyle \frac{1}{\tau }}{\displaystyle \frac{V_q(\eta )}{\eta }}`$ $`=`$ $`{\displaystyle \frac{\eta }{\tau }}\left[1+{\displaystyle \frac{\tau }{D}}(q1){\displaystyle \frac{\eta ^2}{2}}\right]^1`$
$``$ $`{\displaystyle \frac{\eta }{\tau }}\left[1+{\displaystyle \frac{\tau }{D}}(q1){\displaystyle \frac{\eta ^2}{2}}\right]^1a(\tau )\eta ,`$
with $`a(\tau )={\displaystyle \frac{53q}{2\tau (2q)}}`$. This implies a kind of renormalized Ornstein–Uhlenbeck (Gaussian) noise. Using this form of the noise we proceed to study the problem in Eq.(1).

The genetic model is described by the following equation

$$\dot{x}=\frac{1}{2}x+\lambda x(1x)+x(1x)\eta (t),$$
(7)
and, according to the discussion above, we consider in this work that $`\eta (t)`$ is a non–Gaussian noise described by Eqs.(2-3).

The analytical treatment simplifies if we make the change of variables

$$y=\mathrm{ln}\left(\frac{x}{1x}\right),$$
(8)
so that the original multiplicative stochastic differential equation (7) transforms into a new one but now with an additive noise

$$\dot{y}=\mathrm{sinh}(y)+\lambda +\eta (t).$$
(9)
From now on, as in Ref., and in order to simplify the algebra, in what follows we adopt $`\lambda =0`$.

The use of the approximation Eq. (3) (that we reiterate could be valid only for $`|q1|1`$) together with the “unified colored noise approximation” , yields

$$\dot{y}f(y)\frac{1}{[1\frac{f^{}(y)}{a(\tau )}]}+\frac{1}{\tau a(\tau )}\frac{1}{[1\frac{f^{}(y)}{a(\tau )}]}\xi (t),$$
(10)
with $`f(y)=\mathrm{sinh}(y)`$. This result can be obtained in two ways: (a) Applying a direct adiabatic-like elimination procedure to the Langevin system given by Eq. (9) and Eqs. (2,3) . That is, taking the derivative of Eq. (9) with respect to $`t`$ and setting $`\ddot{y}=0`$. (b) A formal one, through the application of the path-integral formalism to the indicated non Markovian Langevin equations, supplemented with a consistent Markovian approximation scheme .

Since, by the use of these “drastic” approximations, the problem has been reduced to Eq. (10) which contains a white, Gaussian, noise, it is easy now, using standard techniques, to obtain the stationary probability density function. For the original variable $`x`$, it reads

$$P_{st}(x)=N\mathrm{exp}[F(x)],$$
(11)
where $`N`$ is a normalization constant and

$`F(x)=`$ $``$ $`{\displaystyle \frac{1}{8D}}\left({\displaystyle \frac{53q}{2q}}\right)^2\left({\displaystyle \frac{12x+2x^2}{x(1x)}}\right)`$
$`\times \left(1+\left({\displaystyle \frac{2q}{53q}}\right)\left({\displaystyle \frac{\tau }{2}}\right)\left({\displaystyle \frac{12x+2x^2}{x(1x)}}\right)\right)`$
$`+`$ $`\mathrm{ln}\left(1+\left({\displaystyle \frac{2q}{53q}}\right)\tau \left({\displaystyle \frac{12x+2x^2}{x(1x)}}\right)\right)\mathrm{ln}(x(1x)).`$
A noise induced transition appears when this density function changes from being unimodal to bimodal. More specifically, the point $`x=1/2`$ changes from being a maximum to being a minimum. The condition of extremum at $`x=1/2`$ and the condition for finding the boundary separating the regions where $`x=1/2`$ is a maximum or a minimum yield

$$D(\tau )=\frac{1}{2}\left(\frac{53q}{2q}\right)^2\left(1+\frac{2q}{53q}2\tau \right)^2\left(1+\frac{2q}{53q}6\tau \right)^1.$$
(13)
such that for a fixed $`\tau `$, the distribution is bimodal for $`D>D(\tau )`$ and it is unimodal for $`D<D(\tau )`$. In the limit $`q1`$, this expression reduces to the result found in  for $`q=1`$:

$$D_0(\tau )=2\frac{[1+\tau ]^2}{1+3\tau }.$$
(14)
Figure 1 plots the $`D(\tau )`$ boundary (the “transition line”) for different values of $`q`$ in the $`(D,\tau )`$ plane. Several features can be noticed in this figure. It is clear that a minimum value of the noise intensity $`D`$ is always needed to induce the transition. The dependence with the parameter $`q`$ appears as an advancement of the transition to smaller values of $`D`$ when $`q>1`$ and a retardation when $`q<1`$. We believe that this will be a general behavior for long tail (power–law) and cut–off noise distributions, respectively. A closer look shows the persistence of the reentrance effect previously found in : for some fixed values of the noise intensity $`D`$ it is possible to obtain a transition, first from unimodal to bimodal and then back to unimodal, just by varying the time correlation $`\tau `$ of the noise.

These main features, as obtained from our very simplifying assumptions are found to be in qualitative agreement with the results coming from numerical simulations of the (not simplified) set of equations (2,3,7) (details of the numerical integration scheme are given in the appendix). In figure 2 we show that the probability distribution can switch from unimodal to bimodal just by increasing the parameter $`q`$, for fixed values of $`D`$ and $`\tau `$. The reentrance effect is illustrated, for $`q=1`$ in figure 3.

The numerical simulations show that there is a crossing of the transition lines for $`q<1`$ and $`q1`$, an aspect which is not reproduced by the approximated theoretical results. Another aspect observed in figure 1 is that for large $`\tau `$, the theoretically predicted transition line, has a smaller linear slope for increasing values of $`q`$. In fact, according to Eq. (13) we have $`D\frac{53q}{3(2q)}\tau `$ in that limit. There seems to be not such a strong $`q`$ dependence of the slope of the curves for large $`\tau `$ in the numerical results. Otherwise, the numerical and theoretical results are in good qualitative agreement.

## 4 FINAL REMARKS

We have presented some preliminary results on a noise induced phenomenon driven by a colored and non–Gaussian noise source that was generated by a $`q`$-distribution. The phenomenon studied here is a noise induced transition like those discussed in , and corresponds the so called genetic model. We have resorted to a simple approximation –that avoids the difficulties of the effective Markovian approximation in – introduced with the idea that for $`|q1|1`$, the non–Gaussian noise will only slightly depart from the Gaussian one, and its effect being only a correction to the Gaussian behavior. The resulting Eq. (13) shows some of the right tendencies with $`q>1`$ or $`q<1`$, although some features observed in the simulations (crossing of the transition lines and slight dependence of the linear slope for large $`\tau `$) are not reproduced by the simplified theoretical treatment. In addition, the stationary probability density function defined in Eqs.(11) and (3) clearly shows, for constant values of $`D`$ and $`\tau `$ and as a function of $`q`$, the transition from unimodal to bimodal behavior. Again, these general results and tendencies, together with the reentrance effect, have been confirmed by extensive numerical simulations as shown in figures 1, 2 and 3. However, the agreement between theory and simulations is only qualitative, as the first are unable to reproduce some of the details revealed by the second. For instance, from the numerical simulations a more complex behavior of the transition line for large $`\tau `$ is apparent, while it is not present in the theoretical predictions (see Fig. 1).

These results are in the same line that a series of previous studies on noise induced phenomena where the noise source was non–Gaussian. In all the cases, we have found that the response of the system results to be enhanced or to have a strong dependence for values of the parameter $`q`$ corresponding to a departure from Gaussian behavior. Some of such studies, including stochastic resonance , gated traps , Brownian motors , have been reviewed in .

Even though our results are so far only preliminary and further studies are required, we expect that the phenomenon discussed here will have remarkable effects on the noise induced phase transitions  as well as in the problem of coupled ratchets . For instance, we can expect changes in the phase diagram of the first, or in the hysteresis cycle of the second, that could have interesting consequences in technological applications. Such studies will be the subject of further work.

## ACKNOWLEDGMENTS

This works is supported by MCyT (Spain) and FEDER, projects BFM2001-0341-C02-01 and BFM2000-1108. HSW acknowledges the partial support from CONICET and ANPCyT, both Argentine agencies, and thanks the MECyD, Spain, for an award within the Programa de Sabáticos para Profesores Visitantes, and the Universitat de les Illes Balears for the kind hospitality extended to him.

## APPENDIX: NUMERICAL SIMULATIONS

In order to test our theoretical predictions,we have resorted to numerical simulations of the system indicated by Eqs. (7) and (2,3). For the numerical integration of those equations we have used a second order stochastic Runge–Kutta type algorithm . For the general problem

$`\dot{x}`$ $`=`$ $`f(x)+g(x)\eta (t)`$ (15)
$`\dot{\eta }`$ $`=`$ $`q(\eta )+ϵ\xi (t),`$ (16)
where $`\xi (t)`$ is a white noise of mean equal to zero and correlations $`\xi (t)\xi (t^{})=\delta (tt^{})`$, the algorithm is as following: time is discretized in steps of size $`h`$ and the value of the process $`\eta (t)`$ at the discretized times is obtained by the recursion relation

$`k(t)`$ $`=`$ $`hq(\eta (t))`$ (17)
$`l(t)`$ $`=`$ $`ϵh^{1/2}G(t)`$ (18)
$`\eta (t+h)`$ $`=`$ $`\eta (t)+{\displaystyle \frac{h}{2}}\left[q(\eta (t))+q(\eta (t)+k(t)+l(t))\right]+l(t),`$ (19)
with the initial value $`\eta (0)=0`$. Here the numbers $`G(t)`$ for different times $`t`$ are independent Gaussian random numbers of mean equal to zero and variance equal to one. They have been obtained by using a particularly efficient generator. In the case $`q<1`$, it might occur that the generated value of $`G(t)`$ is such that $`\eta (t+h)`$ is outside the cut-off interval $`(\eta _c,\eta _c)`$. If this happens, that value is discarded and a new $`G(t)`$ is generated until the condition $`|\eta (t+h)|<\eta _c`$ is fulfilled. We have found that the percentage of discarded $`G(t)`$ values is less than $`5\times 10^4`$ in the worst case analyzed here corresponding to $`q=0.5,\tau =4`$. No such a problem ever occurs for $`q1`$.

Once the values of $`\eta (t)`$ at the discretized times have been obtained, the integration of (15) proceeds with a standard Runge–Kutta second order method:

$`k_x(t)`$ $`=`$ $`f(x(t))`$ (20)
$`x(t+h)`$ $`=`$ $`x(t)+{\displaystyle \frac{h}{2}}[f(x(t))+f(x(t)+k_x(t))+`$ (22)
$`g(x(t))\eta (t)+g(x(t)+k_x(t))\eta (t+h)].`$
We have used a time step $`h=0.01`$ and the histogram for $`P^{st}(x)`$ has been obtained from $`10^8`$ values of $`x`$ separated by a time $`t=1`$ (corresponding, hence, to $`100`$ integration steps).