format: text

text:
# Quasiequilibrium sequences of synchronized and irrotational binary neutron stars in general relativity. II. Newtonian limits

## I Introduction

The research for equilibrium figures of binary systems has a long history. The pioneer work was done by Roche who treated a synchronously rotating, incompressible ellipsoid around a gravitating point source. Following this work, Darwin constructed a synchronized binary system composed of double incompressible ellipsoids (see ). For the non-synchronized case, Aizenman calculated an equilibrium figure of incompressible ellipsoid with internal motion orbiting a point source companion. In 1980’s, the improvement of computer performance made possible to construct equilibrium sequences of synchronized binary systems with compressible equation of state in Newtonian gravity.

Since the beginning of 1990’s, the studies on binary systems gain some importance in relation to astrophysical sources of gravitational waves. For example, coalescing binary neutron stars are expected to be one of the most promising sources of gravitational radiation that could be detected by the interferometric detectors currently in operation (TAMA300) or under construction (GEO600, LIGO, and VIRGO). Until now, numerous theoretical studies have been done in this research field. Among them there are (i) post-Newtonian analytical treatments (e.g. ), (ii) black hole perturbation, (iii) Newtonian and post-Newtonian (semi-)analytical treatments including hydrodynamical effects of stars, (iv) post-Newtonian hydrodynamical computations, (v) fully relativistic hydrodynamical treatments, pioneered by the works of Oohara & Nakamura (see e.g. ), and recently developed by Shibata, Oohara & Nakamura, and the Neutron Star Grand Challenge group. (vi) In parallel of the dynamical studies in (v), a quasiequilibrium formulation of the problem has been developed  and successfully implemented . The basic assumption underlying the quasiequilibirum calculations is that the timescale of the orbit shrinking is larger than that of the orbital revolution in the pre-coalescing state. Consequently the evolution of the binary system can be approximated by a succession of exactly circular orbits, hence the name quasiequilibrium.

The first quasiequilibrium configurations of binary neutron stars in general relativity have been obtained four years ago by Baumgarte et al., followed by Marronetti et al.. However these computations considered synchronized binaries. As far as coalescing binary neutron stars are concerned, this rotation state does not correspond to physical situations, since it has been shown that the gravitational-radiation driven evolution is too rapid for the viscous forces to synchronize the spin of each neutron star with the orbit as they do for ordinary stellar binaries. Rather, the viscosity is negligible and the fluid velocity circulation (with respect to some inertial frame) is conserved in these systems. Provided that the initial spins are not in the millisecond regime, this means that close binary configurations are well approximated by zero vorticity (i.e. irrotational) states. Irrotational configurations are more complicated to obtain because the fluid velocity does not vanish in the co-orbiting frame (as it does for synchronized binaries).

We have successfully developed a numerical method to tackle this problem and presented the first quasiequilibrium configurations of irrotational binary neutron stars elsewhere. The numerical technique relies on a multi-domain spectral method  within spherical coordinates. Since then, two other groups have obtained relativistic irrotational configurations: (i) Marronetti, Mathews & Wilson by means of single-domain finite difference method within Cartesian coordinates and (ii) Uryu & Eriguchi by means of multi-domain finite difference method within spherical coordinates.

Recently, we have presented our method in detail with numerous tests  (hereafter Paper I). Using this method, we have calculated (quasi)equilibrium sequences of binary systems with synchronized and irrotational rotation states. In the present article, we will show the results in Newtonian gravity. The results of relativistic calculations will be given in a forthcoming article .

The plan of the article is as follows. We start by presenting the equations governing binary systems in Newtonian gravity in Sec. II. Section III is devoted to the brief explanation about the improvement on the cases of stiff equation of state such as $`\gamma >2`$. In Sec. IV some tests for numerical code are performed. Then, we will show the results of evolutionary sequences of both synchronized and irrotational binary systems constructed by double identical stars with polytropic equation of state in Sec. V. Section VI contains the summary.

Throughout the present article, we use units of $`G=c=1`$ where $`G`$ and $`c`$ denote the gravitational constant and speed of light, respectively.

## II Formulation

Since the method which we use in the present article for getting equilibrium configurations has already been explained in Paper I, we will only briefly mention the basic equations and the solving procedure in this section.

### A Basic equations

The basic equations governing the problem are as follows (see also Sec. II.F of Paper I):

1. Equation of state:
We adopt a simple equation, i.e. a polytropic equation of state;

$$p(H)=\kappa n(H)^\gamma ,$$
(1)
$`p`$ being the fluid pressure, $`n`$ the fluid baryon number density, $`\kappa `$ a constant, $`\gamma `$ the adiabatic index, and $`H`$ the specific enthalpy. The relation between $`H`$ and $`n`$ is given by

$$n(H)=\left[\frac{\gamma 1}{\gamma }\frac{m_\mathrm{B}}{\kappa }H\right]^{1/(\gamma 1)},$$
(2)
where $`m_\mathrm{B}`$ denotes the mean baryon mass ($`m_\mathrm{B}=1.66\times 10^{27}`$ kg).
2. Euler equation:

$$\frac{𝐯}{t}+𝐯\stackrel{}{}𝐯=\frac{1}{m_\mathrm{B}n}\stackrel{}{}p\stackrel{}{}\nu ,$$
(3)
$`\nu `$ being the gravitational potential, and $`𝐯`$ the velocity field in the inertial frame which is expressed as

$$𝐯=𝐮+𝛀\times 𝐫$$
(4)
in terms of the velocity field in the corotating frame $`𝐮`$ and the orbital motion $`𝛀\times 𝐫`$. In the rigidly rotating case, since there is no motion in the corotating frame ($`𝐮=0`$), we obtain

$$𝐯=𝛀\times 𝐫.$$
(5)
On the other hand, in the irrotational case, the velocity field in the inertial frame is potential, i.e.

$$𝐯=\stackrel{}{}\mathrm{\Psi },$$
(6)
where $`\mathrm{\Psi }`$ is a scalar function.

Under the stationary condition the Euler equation (3) can be integrated as

$$H+\nu \frac{1}{2}(𝛀\times 𝐫)^2=\mathrm{const}$$
(7)
in the synchronized case and

$$H+\nu +\frac{1}{2}(\stackrel{}{}\mathrm{\Psi })^2(𝛀\times 𝐫)\stackrel{}{}\mathrm{\Psi }=\mathrm{const}$$
(8)
in the irrotational case.
3. Equation of continuity:

$$\frac{n}{t}+\stackrel{}{}(n𝐯)=0.$$
(9)
This equation is trivially satisfied by stationary rigid rotation. In the irrotational case, we can rewrite Eq. (9) as

$$n\mathrm{\Delta }\mathrm{\Psi }+\stackrel{}{}n\stackrel{}{}\mathrm{\Psi }=(𝛀\times 𝐫)\stackrel{}{}n.$$
(10)
4. Poisson equation for the gravitational field:

$$\mathrm{\Delta }\nu =4\pi m_\mathrm{B}n.$$
(11)

### B Solving procedure

* First of all, we prepare the equilibrium figures of two spherical stars. In the present computation, we treat binary systems composed of equal mass stars with the same equation of state. However, we symmetrize only with respect to the orbital plane. Although this treatment elongates the compute time, we do not symmetrize the stars because we are planning to study binary systems composed of different mass stars in the series of this research.
* The separation between the centers of the two stars is held fixed. Here we define the center as the point of the maximum enthalpy (or equivalently maximum density - see Eq. (2)).
* By setting the central value of the gradient of enthalpy to be zero, we calculate the orbital angular velocity $`\mathrm{\Omega }`$ (see Sec. IV.D.2 of Paper I for details).
* For irrotational configurations, the velocity potential $`\mathrm{\Psi }`$ is obtained by solving Eq. (10).
* Taking into account the gravitational potential from the companion star and the centrifugal force, we calculate the new enthalpy field from Eq. (7) or (8).
* Using the new enthalpy, we search for the location of $`H=0`$; This defines the stellar surface since $`H=0`$ is equivalent to $`p=0`$ (see Eqs. (1) and (2)). Having determined the new stellar surface, we change the position of the inner domain boundary to make it fit with the stellar surface.
* By substituting the new enthalpy in the new domain into Eq. (2), the new baryon density is obtained.
* Inserting the new baryon density into the source term of the Poisson equation (11), we can get the new gravitational potential.
* Finally, we compare the new enthalpy field with the old one. If the difference is smaller than some fixed threshold, the calculation is considered as having converged. If not, return to step (c), and continue.

During the calculation, we make the baryon mass to converge to a fixed value by varying the central enthalpy. Moreover, some relaxation is performed on the gravitational potential, the centrifugal force, and the orbital angular velocity, when such a treatment helps the iteration to converge. The explicit procedures are detailed in Paper I.

## III Improvement on the cases of stiff equation of state

### A Regularization technique

The numerical method which we use in this series of research  is a multi-domain spectral method. In general, spectral methods lose much of their accuracy when non-smooth functions are treated because of the so-called Gibbs phenomenon. This phenomenon is well known from the most familiar spectral method, namely, the theory of Fourier series: the Fourier coefficients ($`c_n`$) of a function $`f`$ which is of class $`𝒞^p`$ but not $`𝒞^{p+1}`$ decrease as only $`1/n^p`$ with increasing $`n`$. In particular, if the function has some discontinuity, its approximation by a Fourier series does not converge towards $`f`$ at the discontinuity point. The multi-domain spectral method circumvents the Gibbs phenomenon. The basic idea is to divide the space into domains chosen so that the physical discontinuities are located onto the boundaries between the domains.

However, even if the circumvention of the Gibbs phenomenon has been done by introducing multi-domains, there remains the Gibbs phenomenon at the boundaries between the domains when some physical field has infinite derivative at the boundary. For example, when we consider a star, its surface coincides with the boundary of a domain. If the star has a stiff equation of state, i.e. $`\gamma >2`$ in terms of the adiabatic index, the density decreases so rapidly that $`n/r`$ diverges at the surface. In order to recover high accuracy, we have developed a method to regularize the density profile by extracting its diverging part (see Sec. IV of Ref. ).

Here we briefly summarize this method. For a polytrope with adiabatic index $`\gamma `$, the matter density $`n`$ behaves as $`H^{1/(\gamma 1)}`$ (see Eq. (2)). In the steady state configuration, $`H`$ is Taylor expandable at the neighborhood of the stellar surface because the gravitational potential $`\nu `$ is Taylor expandable there (cf. Eq. (7) or (8)). Therefore $`H`$ vanishes as $`rR(\theta ,\phi )`$, where $`R(\theta ,\phi )`$ is the equation of the stellar surface. Consequently $`n`$ behaves as $`n[rR(\theta ,\phi )]^{1/(\gamma 1)}`$.

We introduce a known potential $`\mathrm{\Phi }_{\mathrm{div}}`$ such that $`n_{\mathrm{div}}:=\mathrm{\Delta }\mathrm{\Phi }_{\mathrm{div}}`$ has the same pathological behavior as $`n`$ and such that $`nn_{\mathrm{div}}`$ is a regular function (at least more regular than $`n`$). We then numerically solve

$$\mathrm{\Delta }\mathrm{\Phi }_{\mathrm{regu}}=nn_{\mathrm{div}},$$
(12)
where $`\mathrm{\Phi }_{\mathrm{regu}}:=\mathrm{\Phi }\mathrm{\Phi }_{\mathrm{div}}`$. Consider, for instance,

$$\mathrm{\Phi }_{\mathrm{div}}=F(\xi ,\theta ,\phi )(1\xi ^2)^{(\alpha +2)},$$
(13)
where $`\alpha =1/(\gamma 1)`$. $`F`$ is an arbitrary regular function and $`\xi `$ is a new radial variable such that $`\xi =1`$ at the surface of the star. It is easy to see that $`\mathrm{\Delta }\mathrm{\Phi }_{\mathrm{div}}`$ has a term vanishing at the surface with the same pathological behavior as $`n`$, i.e., $`(1\xi ^2)^\alpha `$. We have indeed

$`\stackrel{~}{\mathrm{\Delta }}\mathrm{\Phi }_{\mathrm{div}}`$ $`=`$ $`\stackrel{~}{\mathrm{\Delta }}F(1\xi ^2)^{(\alpha +2)}`$ (16)
$`4(\alpha +2)\xi (1\xi ^2)^{(\alpha +1)}_\xi F`$
$`+(\alpha +2)[6(1\xi ^2)^{(\alpha +1)}+4(\alpha +1)\xi ^2(1\xi ^2)^\alpha ]F,`$
where $`\stackrel{~}{\mathrm{\Delta }}`$ is the Laplacian computed with respect to $`(\xi ,\theta ,\phi )`$.

For calculational simplicity, we choose an harmonic function for $`F(\xi ,\theta ,\phi )`$. Then, we can write $`\mathrm{\Phi }_{\mathrm{div}}=_{l,m}a_{lm}\mathrm{\Phi }_{lm}`$, where

$$\mathrm{\Phi }_{lm}:=\xi ^l(1\xi ^2)^{(\alpha +2)}Y_l^m(\theta ,\phi ),$$
(17)
$`a_{lm}`$ being some numerical coefficients to be determined and $`Y_l^m`$ the standard spherical harmonics. We then obtain

$$n_{\mathrm{div}}=\underset{l,m}{}a_{lm}C_l(\xi )Y_l^m(\theta ,\phi ),$$
(18)
with

$`C_l(\xi )`$ $`=`$ $`(\alpha +2)[(4l+6)(1\xi ^2)^{(\alpha +1)}\xi ^l`$ (20)
$`+4(\alpha +1)\xi ^{(l+2)}(1\xi ^2)^\alpha ].`$
We now have to determine the values of the coefficients $`a_{lm}`$ which give the most regular function $`n_{\mathrm{regu}}:=nn_{\mathrm{div}}`$. The criterion which seems to give the best results is the following one. We expand both $`n`$ and $`n_{\mathrm{div}}`$ as truncated series of spherical harmonics $`Y_l^m(\theta ,\phi )`$ and Chebyshev polynomial $`T_i(\xi )`$,

$`n(\xi ,\theta ,\phi )={\displaystyle \underset{i,l,m=0}{\overset{I,L,M}{}}}n_{ilm}T_i(\xi )Y_l^m(\theta ,\phi ),`$ (21)
$`n_{\mathrm{div}}(\xi ,\theta ,\phi )={\displaystyle \underset{i,l,m=0}{\overset{I,L,M}{}}}a_{lm}C_{li}T_i(\xi )Y_l^m(\theta ,\phi ),`$ (22)
where we expand each of the function $`C_l(\xi )`$ in a Chebyshev series as

$$C_l(\xi )=\underset{i=0}{\overset{I}{}}C_{li}T_i(\xi ).$$
(23)
The values of $`a_{lm}`$ is computed in such a way that the $`I`$th coefficient of the truncated series of $`n_{\mathrm{regu}}`$ vanishes:

$$a_{lm}=\frac{n_{Ilm}}{C_{lI}}.$$
(24)
By means of the above procedure, we eliminate in $`n`$ the pathological term which vanishes as $`(1\xi ^2)^\alpha `$ but we introduce another pathological term proportional to $`(1\xi ^2)^{\alpha +1}`$. However, the divergence occurs in a higher order derivative of this term so that it has a much weaker effect on the accuracy of the results. The method can be improved by taking

$`\mathrm{\Phi }_{\mathrm{div}}`$ $`=`$ $`F(\xi ,\theta ,\phi )(1\xi ^2)^{(\alpha +2)}[a_1+a_2(1\xi ^2)`$ (26)
$`+a_3(1\xi ^2)^2+\mathrm{}+a_K(1\xi ^2)^{K1}]`$
instead of Eq. (13). The coefficients $`a_K`$ are chosen in such a way that the first, second, …, $`K`$th derivatives of $`n_{\mathrm{regu}}`$ vanish at $`\xi =1`$. Let us call $`K`$ the regularization degree of the procedure.

The explicit procedure to determine the coefficients $`a_K`$ is as follows: Let us consider the case of $`K=2`$ and choose an harmonic function for $`F(\xi ,\theta ,\phi )`$ for example. In this case, we need two diverging terms, $`n_{\mathrm{div}:k}(1kK)`$, where we define

$$n_{\mathrm{div}:k}:=\underset{l,m}{}a_{lm}^kC_l^k(\xi )Y_l^m(\theta ,\phi ),$$
(27)
with

$`C_l^k(\xi )`$ $`=`$ $`(\alpha +k+1)[(4l+6)(1\xi ^2)^{(\alpha +k)}\xi ^l`$ (29)
$`+4(\alpha +k)\xi ^{(l+2)}(1\xi ^2)^{\alpha +k1}].`$
First of all, we expand $`n`$, $`n_{\mathrm{div}:1}`$ and $`n_{\mathrm{div}:2}`$ by $`Y_l^m(\theta ,\phi )`$ and $`T_i(\xi )`$. Then, we can express $`C_l^k(\xi )`$ as an expansion in a Chebyshev series,

$$C_l^k(\xi )=\underset{i=0}{\overset{I}{}}C_{li}^kT_i(\xi ).$$
(30)
Finally, we solve the simultaneous equations

$`n_{Ilm}`$ $`=`$ $`a_{lm}^1C_{lI}^1+a_{lm}^2C_{lI}^2,`$ (31)
$`n_{I1lm}`$ $`=`$ $`a_{lm}^1C_{lI1}^1+a_{lm}^2C_{lI1}^2,`$ (32)
so that the $`I`$th and $`(I1)`$th coefficients of the truncated series of $`n_{\mathrm{regu}}`$ vanish.

### B Illustration and test

To validate the above regularization technique, we investigate the improvement in the relative error in the virial theorem for a spherical static star, which is defined as

$$\mathrm{Virial}\mathrm{error}=\frac{|W+3P|}{|W|},$$
(33)
$`W`$ being the gravitational potential energy and $`P`$ the volume integral of the fluid pressure. We give the values of the virial error in Fig. 1. In this figure we show the cases of (1) $`\gamma =3`$, (2) $`\gamma =2.5`$, (3) $`\gamma =2.25`$, and (4) $`\gamma =2`$, and in each small figure except for the case of $`\gamma =2`$ for which no Gibbs phenomenon occurs, we have drawn four lines: “no regularization”, regularization degree $`K=1`$, 2, and 3. It is found from these plots that the regularization dramatically improves the accuracy. However, the improvement by the regularization saturates at $`K=2`$. Note here that since there is no Gibbs phenomenon in the calculation of a spherical star in the case of $`\gamma =2`$, the virial error decreases rapidly (exponential decay - evanescent error) with the increase of the number of radial collocation point $`N_r`$, and reaches the minimum value due to the finite number of digits (15=double precision) used in the numerical computation.

Considering the saturation of the virial error, we will choose the regularization degree $`K=2`$ in the following numerical computations for $`\gamma >2`$. We note that since the matter density $`n`$ behaves as $`H^{1/(\gamma 1)}`$, the derivatives of $`n`$ of order higher than the second one diverge at the stellar surface even if $`\gamma <2`$. Therefore we will also use the regularization ($`K=2`$) for the cases of $`\gamma <2`$.

## IV Tests of the numerical code

Numerous tests have been already performed in the previous paper in the irrotational case (see Sec. V.B of Paper I). In particular a direct comparison with the numerical results of Uryu & Eriguchi  has been performed and the result presented in Table I of Paper I. Here we will focus on some tests in the synchronized case, not presented in Paper I. We will also show the relative error on the virial theorem along an evolutionary sequence in both cases, extending the test presented in Paper I (Fig. 7) to adiabatic indices different from 2.

### A Comparison with analytical solutions

In order to investigate the discrepancy between the results from the numerical code and those from Taniguchi & Nakamura’s analytic solution in the synchronized case with $`\gamma =2`$, we present the relative differences on global quantities as functions of the separation in a log-log plot in Fig. 2 (in the same way as Fig. 12 of Paper I for irrotational configurations). The relative differences are defined as follows:

$`{\displaystyle \frac{E_{\mathrm{num}}E_{\mathrm{ana}}}{GM^2/R_0}},`$ (34)
$`{\displaystyle \frac{J_{\mathrm{num}}J_{\mathrm{ana}}}{Md^2\mathrm{\Omega }_{\mathrm{Kep}}/2}},`$ (35)
$`{\displaystyle \frac{\mathrm{\Omega }_{\mathrm{num}}\mathrm{\Omega }_{\mathrm{ana}}}{\mathrm{\Omega }_{\mathrm{Kep}}}},`$ (36)
$`|\delta \rho _{c:\mathrm{num}}\delta \rho _{c:\mathrm{ana}}|,`$ (37)
where $`\mathrm{\Omega }_{\mathrm{Kep}}`$ is the Keplerian velocity for point mass particles

$$\mathrm{\Omega }_{\mathrm{Kep}}:=\left(\frac{2GM}{d^3}\right)^{1/2}.$$
(38)
In the Taniguchi & Nakamura’s analytic solution , the global quantities are expanded up to $`O[(R_0/d)^6]`$. After subtracting these analytic solutions from the numerical one, only the terms of order higher than $`O[(R_0/d)^6]`$ should remain. Indeed, we can see from Fig. 2 that the discrepancies between numerical and analytical solutions for the energy $`E`$ and the relative change in central density $`\delta \rho _c`$ are both higher than $`O[(R_0/d)^8]`$, and those for the angular momentum $`J`$ and the orbital angular velocity $`\mathrm{\Omega }`$ are both higher than $`O[(R_0/d)^7]`$. This means that the numerical solution and the analytical one agree up to the accuracy of this latter ($`O[(R_0/d)^6]`$).

### B Virial theorem

One of the best indicators of the accuracy of numerical solutions for binary systems in Newtonian gravity is the relative error in the virial theorem. This error is defined as

$$\mathrm{Virial}\mathrm{error}=\frac{|2T+W+3P|}{|W|},$$
(39)
where $`T`$, $`W`$, and $`P`$ denote respectively the kinetic energy of the binary system, its gravitational potential energy and the volume integral of the fluid pressure. By virtue of the virial theorem, the quantity defined by Eq. (39) should be zero for an exact solution. We show it in Fig. 3. This Figure can be considered as an extension to that presented in Paper I (Fig. 7) to adiabatic indices different from 2. It also contains the synchronized case. We can see from Fig. 3 that the code is quite accurate as $`\mathrm{Virial}\mathrm{error}<10^5`$ for $`d_G/R_03`$ even if $`\gamma =3`$, where $`d_G`$ denotes the separation between two centers of masses of each star.

## V Results

Using the numerical method explained in the above sections and in Paper I, we have constructed equilibrium sequences of both synchronized and irrotational binary systems in Newtonian gravity. Here, we consider only the case of binary systems composed of stars with equal mass and equation of state. We use 3 domains (one for the fluid interior) for each star and the following number of spectral coefficients: $`N_r\times N_\theta \times N_\phi =33\times 25\times 24`$ in each domain. A view of the configuration at the energy turning point (Sec. V C below) along a sequence with a $`\gamma =3`$ EOS is shown in Fig. 4 for synchronized binaries, and in Fig. 5 for irrotational one. We are interested in the turning point of total energy (and/or total angular momentum) because it corresponds to the onset of secular instability in the synchronized case and that of dynamical instability in the irrotational one at least for the ellipsoidal approximation.

### A Equilibrium sequences

Our results for equilibrium constant-mass sequences (evolutionary sequences) with adiabatic indices $`\gamma =3,2.5,2.25,2`$ and 1.8 are presented in Tables I and II. In these tables, $`d`$ denotes the separation between the centers of the two stars. Let us recall that the center of a star is defined as the point of maximum enthalpy (or equivalently maximum density). On the other hand, $`d_G`$ denotes the orbital separation between centers of masses of two stars. $`R_0`$, $`a_1`$, $`a_2`$, $`a_3`$, and $`a_{1,\mathrm{opp}}`$ are the radius of a spherical star of same mass, the radius parallel to x-axis toward the companion star, the radius parallel to y-axis, the radius parallel to z-axis, and the radius parallel to x-axis opposite to the companion star. The $`(x,y,z)`$ axes are the same as in Fig. 1 of Paper I. $`\rho _c`$ and $`\rho _{c0}`$ indicate the central density of a star and that of a spherical star of same mass. The normalized quantities $`\overline{\mathrm{\Omega }}`$, $`\overline{J}`$, and $`\overline{E}`$ are defined by

$`\overline{\mathrm{\Omega }}`$ $`:=`$ $`{\displaystyle \frac{\mathrm{\Omega }}{(\pi G\rho _0)^{1/2}}},`$ (40)
$`\overline{J}`$ $`:=`$ $`{\displaystyle \frac{J}{(GM^3R_0)^{1/2}}},`$ (41)
$`\overline{E}`$ $`:=`$ $`{\displaystyle \frac{E}{GM^2/R_0}},`$ (42)
where $`\mathrm{\Omega }`$, $`J`$, and $`E`$ denote respectively the orbital angular velocity, the total angular momentum, and the total energy, and $`\rho _0`$ is the averaged density of a spherical star of the same mass:

$$\rho _0:=\frac{3M}{4\pi R_0^3}.$$
(44)
Also listed in Tables I and II is the ratio

$$\chi :=\frac{(H/r)_{\mathrm{eq},\mathrm{comp}}}{(H/r)_{\mathrm{pole}}},$$
(45)
where $`(H/r)_{\mathrm{eq},\mathrm{comp}}`$ \[resp. $`(H/r)_{\mathrm{pole}}`$\] stands for the radial derivative of the enthalpy at the point on the stellar surface located in the orbital plane and looking toward the companion star \[resp. at the intersection between the surface and the axis perpendicular to the orbital plane and going through the stellar center ($`z`$ axis)\]. This quantity is useful because the mass shedding limit (“Roche limit”) corresponds to $`\chi =0`$ (cf. Sec. IV.E of Paper I). When $`\chi =0`$, an angular point (cusp) appears at the equator of the star in the direction to the companion.

An equilibrium sequence terminates either by a contact configuration, which corresponds to $`d/a_1=2`$ or by the mass shedding limit given by $`\chi =0`$. Note that these two conditions are not exclusive, as one can have $`\chi =0`$ at contact. By means of our numerical method, it is difficult to get exactly such configurations, so that the final points in Tables I and II are close to, but not exactly equal to the real end points $`d/a_1=2`$ or $`\chi =0`$.

For synchronized binaries, the real end point is the contact one . In that configuration, the surface of each star has a cusp at the contact point. However, our numerical multi-domain method assumes that the boundary between each domain is a differentiable surface. Therefore, for computing very close configurations, we stop the adaptation of the domain to the surface of the star when $`\chi <0.30.35`$ (see Sec. IV.E of Paper I). For synchronized binaries, the use of adapted domain is not essential, since no boundary condition is set at the stellar surface. The last lines for each $`\gamma `$ in Table I are obtained in this way<sup>*</sup><sup>*</sup>*Due to the overlapping of the external domains for contact configurations, it was not possible to compute these latter with high accuracy. This is why we stopped Table I slightly before $`d/a_1=2`$.

On the other hand, in the irrotational case, leaving the adaptation of the domain to the stellar surface, results in some numerical error since the solving method for the equation governing the velocity potential \[Eq. (10)\] assumes that the domain boundary coincide with the stellar surface (see Appendix B of Paper I). Therefore we stop the calculations at some points which are very close to the cusp points but slightly separated.

The symbol $``$ in Tables I and II indicates the points of minimum total energy (and total angular momentum) along the sequence, also called turning points. In the synchronized case, the minimum points exist for $`\gamma 2`$, and in the irrotational case, they do for $`\gamma 2.5`$. In both synchronized and irrotational cases, the minimum points of total energy and total angular momentum coincide with each other. It is worth to note here that we cannot exclude the possibility of existence of minimum points for $`\gamma <2`$ in the synchronized case and for $`\gamma <2.5`$ in the irrotational one, because we do not calculate up to contact points. However, we suspect that the critical values of the existence of minimum points are $`\gamma 1.8`$ in the synchronized case and $`\gamma 2.3`$ in the irrotational one. The detailed discussions are given later (see Sec. V C).

We show the total energy, the total angular momentum, the orbital angular velocity, and the relative change in central density along a constant mass sequence in Figs. 69. We can see from figures of total energy and total angular momentum that the values for synchronized systems are larger than those for irrotational ones. This is because the effect of spin of each star in the synchronized case is larger than that in the irrotational one. Such a spin produces not only the direct effects like spin energy or spin angular momentum but also larger deformations which results in larger quadrupole moments. We can see clearly a turning point in energy and angular momentum curves for large values of $`\gamma `$ on Figs. 68. This feature will be discussed in Sec. V C.

It is found from Fig. 9 that the central density decreases in both cases of synchronized and irrotational. The decrease of the central density in the synchronized case is about one order larger than that in the irrotational one. These behaviors are analytically known. For the synchronized case, Chandrasekhar obtained the lowest order change in the central density about 70 years ago and Taniguchi & Nakamura have calculated the higher order change, and Taniguchi & Nakamura have also shown it for the irrotational case.

In Figs. 1012, we show isocontours of baryon density in the synchronized case with $`\gamma =3`$, $`2`$ and $`1.8`$, respectively. These figures correspond to the configurations in the last lines (for each adiabatic index $`\gamma `$) in Table I. Note here that the small rough on the stellar surface of $`\gamma =3`$ in Fig. 10 is an artifact of the graphical software. The solved figure is a completely smooth one, thanks to the technique discussed in Sec. III A.

Isocontours of baryon density with $`\gamma =3`$, $`2`$ and $`1.8`$ for irrotational binaries are shown in Figs. 1315. These configurations correspond to the semi-final lines (for each adiabatic index $`\gamma `$) in Table II. Again note that the small rough on the stellar surface of $`\gamma =3`$ in Fig. 13 is an artifact of the graphical software and that the solved figure is a completely smooth one, thanks to the technique discussed in Sec. III A. We do not depict the configurations of the last lines of Table II for the following reason. As mentioned above and in Paper I, we make the boundary the inner domain fit with the stellar surface. This procedure is essential to accurately solve equilibrium figures in the irrotational case (Appendix B of Paper I). However, due to the apparition of the cusp on the stellar surface, we can no longer adapt the domain to that surface because all quantities are expressed by summation of a finite number of differentiable functions. For very close configurations, just prior to the apparition of the cusp, the surface is highly distorted so that there appear unphysical oscillations when using finite series of differentiable functions (Gibbs phenomenon). In the last lines for each $`\gamma `$ in Table II, such oscillations start to be seen, although they do not alter appreciably the global quantities. Therefore to avoid any misunderstanding, we choose not to plot the final lines in Table II, although we listed them.

Figures 1618 show some isocontours of the velocity potential, as well as the velocity field in the co-orbiting frame, for the polytropic indices $`\gamma =3`$, $`2`$ and $`1.8`$, respectively. In these figures, we show only one of the two star. The companion star is located at the position symmetric with respect to the $`yz`$ plane. The velocity potential $`\mathrm{\Psi }_0`$ is defined as

$$\mathrm{\Psi }_0:=\mathrm{\Psi }𝐖_0𝐫,$$
(46)
where $`𝐖_0`$ is the constant translational velocity field defined as the central value of $`𝐖:=𝛀\times 𝐫`$. Note that the vector field is tangent to the stellar surface, as it should be.

### B End points of sequences: contact vs. cusp

An equilibrium sequence terminates by the contact between the two stars $`(d/a_1=2)`$ or by a cusp at the onset of mass shedding $`(\chi =0)`$. In order to investigate which final fate occurs, it is helpful to display a sequence in the $`d/a_1\chi `$ plane. This is done in Fig. 19 where we compare the synchronized sequence with the irrotational one for $`\gamma =2`$. It is found from this figure that the value of $`\chi `$ in the irrotational case is larger than that in the synchronized case for large separations. However, when the separation decreases below $`d/a_13`$, $`\chi `$ in the irrotational case decreases rapidly and becomes smaller than that in the synchronized case for $`d/a_1<2.5`$. We magnify the region near the end of the sequences in Fig. 20. If we extrapolate the results up to the zero value of $`\chi `$ in the figure, we can speculate about the final fates of sequences. In the synchronized case, it seems that all the lines will reach $`\chi =0`$ at $`d/a_1=2`$. This means that the cusp does not appear before the two stars contact with each other. On the other hand, in the irrotational case, it seems that the lines may reach $`\chi =0`$ before $`d/a_1=2`$. The values of $`d/a_1`$ may be $`d/a_12.1`$ for $`\gamma =3`$ and $`d/a_12.25`$ for $`\gamma =1.8`$. This means that the cusp may be created before contact of stars in the irrotational case. This behavior in the irrotational case agrees with the results of Uryu.

It is worth to note here that the point where the cusp may be created is expressed by the orbital separation divided by the radius to the companion star ($`a_1`$). Then, it looks that the cusp appears later along equilibrium sequences for higher adiabatic index, i.e., stiffer equation of state. However, if we express the cusp point by the orbital separation divided by the radius of a spherical star with the same mass ($`R_0`$), the order of cusp appearance will be reversed. Although these behaviors are completely different from each other, the origin is the same. The fluid with lower adiabatic index is less affected by the tidal force and slightly deformed. It results in the larger $`d/a_1`$ with fixing the orbital separation $`d`$.

A simple explanation of the difference in $`\chi `$ between in the synchronized case and in the irrotational one is as follows. Using Eq. (46) and the decomposition of the orbital motion

$$𝛀\times 𝐫=𝐖_0+𝐖_s,$$
(47)
we can rewrite Eq. (8) as

$$H+\nu \frac{1}{2}(𝛀\times 𝐫)^2+\frac{1}{2}(\stackrel{}{}\mathrm{\Psi }_0𝐖_s)^2=\mathrm{const},$$
(48)
where $`𝐖_0`$ has been defined in Sec.V A \[see Eq. (46)\] and $`𝐖_s`$ is the spin part of the orbital motion defined sor star No.1 as

$$𝐖_s:=\mathrm{\Omega }(y_1,x_1,0),$$
(49)
where $`(x_1,y_1,z_1)`$ are the Cartesian coordinate centered on star 1 (see Fig. 1 in Paper I). In the following explanation, we pay particular attention to star 1. Note here that the last term on the left-hand side of Eq. (48) is the kinetic energy of the velocity field in the corotating frame (see Eq. (4)). Comparing Eq. (7) with Eq. (48), we see that the enthalpy fields for synchronized and irrotational binaries differ precisely by that kinetic energy.

Then $`\chi `$ for each case becomes

$`\chi ^{\mathrm{synch}}={\displaystyle \frac{\left({\displaystyle \frac{\nu }{x_1}}\right)_{\mathrm{eq},\mathrm{comp}}+\mathrm{\Omega }^2\left({\displaystyle \frac{d}{2}}+a_1\right)}{\left({\displaystyle \frac{\nu }{z_1}}\right)_{\mathrm{pole}}}},`$ (50)
$`\chi ^{\mathrm{irrot}}={\displaystyle \frac{\left({\displaystyle \frac{\nu }{x_1}}\right)_{\mathrm{eq},\mathrm{comp}}+\mathrm{\Omega }^2({\displaystyle \frac{d}{2}}+a_1){\displaystyle \frac{1}{2}}{\displaystyle \frac{}{x_1}}(\stackrel{}{}\mathrm{\Psi }_0𝐖_s)^2|_{\mathrm{eq},\mathrm{comp}}}{\left({\displaystyle \frac{\nu }{z_1}}\right)_{\mathrm{pole}}}}.`$ (51)
Of course, due to the difference of deformation of the stars, the gravitational potential $`\nu `$ for irrotational binaries is different from that for synchronized ones, even if we set two stars at the same orbital separation in the binary systems. However, the largest difference between $`\chi ^{\mathrm{synch}}`$ and $`\chi ^{\mathrm{irrot}}`$ is the last term in the numerator of Eq. (51). We show this behavior in the following. First, we divide $`\chi `$ into two parts;

$`\chi ^{\mathrm{synch}}=\chi _{\mathrm{pot}}^{\mathrm{synch}}+\chi _{\mathrm{flow}}^{\mathrm{synch}},`$ (52)
$`\chi ^{\mathrm{irrot}}=\chi _{\mathrm{pot}}^{\mathrm{irrot}}+\chi _{\mathrm{flow}}^{\mathrm{irrot}},`$ (53)
where

$`\chi _{\mathrm{pot}}^{\mathrm{synch}}={\displaystyle \frac{\left({\displaystyle \frac{\nu }{x_1}}\right)_{\mathrm{eq},\mathrm{comp}}+\mathrm{\Omega }^2\left({\displaystyle \frac{d}{2}}+a_1\right)}{\left({\displaystyle \frac{\nu }{z_1}}\right)_{\mathrm{pole}}}}(\mathrm{for}\mathrm{synchronized}\mathrm{binaries}),`$ (54)
$`\chi _{\mathrm{flow}}^{\mathrm{synch}}=0,`$ (55)
$`\chi _{\mathrm{pot}}^{\mathrm{irrot}}={\displaystyle \frac{\left({\displaystyle \frac{\nu }{x_1}}\right)_{\mathrm{eq},\mathrm{comp}}+\mathrm{\Omega }^2\left({\displaystyle \frac{d}{2}}+a_1\right)}{\left({\displaystyle \frac{\nu }{z_1}}\right)_{\mathrm{pole}}}}(\mathrm{for}\mathrm{irrotational}\mathrm{binaries}),`$ (56)
$`\chi _{\mathrm{flow}}^{\mathrm{irrot}}={\displaystyle \frac{{\displaystyle \frac{1}{2}}{\displaystyle \frac{}{x_1}}(\stackrel{}{}\mathrm{\Psi }_0𝐖_s)^2|_{\mathrm{eq},\mathrm{comp}}}{\left({\displaystyle \frac{\nu }{z_1}}\right)_{\mathrm{pole}}}}(\mathrm{for}\mathrm{irrotational}\mathrm{binaries}).`$ (57)
Then, we plot the differences $`\delta \chi _{\mathrm{pot}}=\chi _{\mathrm{pot}}^{\mathrm{irrot}}\chi _{\mathrm{pot}}^{\mathrm{synch}}`$ and $`\delta \chi _{\mathrm{flow}}=\chi _{\mathrm{flow}}^{\mathrm{irrot}}\chi _{\mathrm{flow}}^{\mathrm{synch}}`$ in Fig. 21 in the $`\gamma =2`$ case. Here, $`\delta \chi _{\mathrm{pot}}`$ denotes the difference in $`\chi `$ between synchronized and irrotational cases which includes the gravitational potential force plus centrifugal force, and $`\delta \chi _{\mathrm{flow}}`$ denotes one which has relation to the internal flow in the corotating frame.

We can confirm from this figure that, at small separation, the difference in the gravitational force plus centrifugal force is smaller than that in the force related internal flow. Therefore, we will consider only the last term in the numerator of Eq. (51) below.

Now, we assume the following form for $`\stackrel{}{}\mathrm{\Psi }_0`$:

$$\stackrel{}{}\mathrm{\Psi }_0=\mathrm{\Omega }(f(d,x_1,y_1,z_1),g(d,x_1,y_1,z_1),h(d,x_1,y_1,z_1)),$$
(58)
where $`f`$, $`g`$ and $`h`$ are some scalar functions of the orbital separation $`d`$ and the coordinates $`(x_1,y_1,z_1)`$. Then the last term in the numerator of Eq. (51) can be calculated as

$`F:=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \frac{}{x_1}}(\stackrel{}{}\mathrm{\Psi }_0𝐖_s)^2,`$ (59)
$`=`$ $`\mathrm{\Omega }^2\left[(f+y_1){\displaystyle \frac{f}{x_1}}+(gx_1)\left({\displaystyle \frac{g}{x_1}}1\right)+h{\displaystyle \frac{h}{x_1}}\right].`$ (60)
Accordingly, the surface value at $`(r_1,\theta _1,\phi _1)=(a_1,\pi /2,0)(x_1,y_1,z_1)=(a_1,0,0)`$ becomes

$$F_{\mathrm{eq},\mathrm{comp}}=\mathrm{\Omega }^2a_1\left[\frac{f(a_1)}{a_1}\frac{f}{x_1}(a_1)+\left(\frac{g(a_1)}{a_1}1\right)\left(\frac{g}{x_1}(a_1)1\right)\right],$$
(61)
where we have used $`h(d,x_1,y_1,0)=0`$ because the velocity field is antisymmetric with respect to the $`xy`$ plane. Note that we simplify the argument list of $`f`$ and $`g`$ from $`(a_1,0,0)`$ to $`(a_1)`$. Let us consider the dominant terms in $`f`$ and $`g`$ which are

$`f`$ $``$ $`\mathrm{\Lambda }(d)y_1,`$ (62)
$`g`$ $``$ $`\mathrm{\Lambda }(d)x_1.`$ (63)
Here $`\mathrm{\Lambda }`$ is a function of the separation $`d`$. The forms (62) and (63) can be justified from the studies by Lai, Rasio & Shapiro or Taniguchi & Nakamura. Indeed, for ellipsoidal models, $`\mathrm{\Lambda }`$ becomes

$$\mathrm{\Lambda }=\frac{a_1^2a_2^2}{a_1^2+a_2^2},$$
(64)
and its dependence of $`d`$ is $`O[(d/R_0)^3]`$ because the dominant effect in the deviation from a spherical star is produced by the tidal force. From Eqs. (62)-(63) it appears that the first term in the brackets on the right-hand side of Eq. (61) is negligible as compared with the second term, so that one is left with

$$F_{\mathrm{eq},\mathrm{comp}}\mathrm{\Omega }^2a_1\left(\frac{g(a_1)}{a_1}1\right)\left(\frac{g}{x_1}(a_1)1\right).$$
(65)
Taking into account Eqs. (63) and (64), we can have only two types of behaviors for the function $`g`$:

* $`{\displaystyle \frac{g(a_1)}{a_1}}<1`$ and $`{\displaystyle \frac{g}{x_1}}(a_1)1`$,
* $`{\displaystyle \frac{g(a_1)}{a_1}}<1`$ and $`{\displaystyle \frac{g}{x_1}}(a_1)>1`$.

Since $`g`$ tends to zero for very large orbital separations, the case (i) will occur in the earlier stage of the sequence. In this case, $`F_{\mathrm{eq},\mathrm{comp}}`$ becomes negative value so that it makes $`\chi ^{\mathrm{irrot}}`$ larger value than $`\chi ^{\mathrm{synch}}`$ It is worth to note that since the terms $`(\nu /x_1)`$, $`(\nu /z_1)`$, and $`\mathrm{\Omega }^2(d/2+x_1)`$ have negative values, $`\chi ^{\mathrm{irrot}}`$ becomes larger for the negative value of $`F`$.. However, if the case (ii) actually occur, $`F_{\mathrm{eq},\mathrm{comp}}`$ can take positive value, and $`\chi ^{\mathrm{irrot}}`$ can be smaller than that in the synchronized case.

In Fig. 22, the y-axis component of $`\stackrel{}{}\mathrm{\Psi }_0`$, i.e., the function $`g`$, and its derivative are shown as a function of the coordinate $`x_1`$ (normalized by the surface value $`a_1`$). This figure corresponds to the last line of the $`\gamma =2`$ case in Table II. We can see from this figure that the y-axis component of $`\stackrel{}{}\mathrm{\Psi }_0`$ remains smaller than that of $`𝐖_s`$ throughout the interior of the star, but the derivative near the stellar surface becomes larger than unity, i.e., the derivative of the y-axis component of $`𝐖_s`$. This means that the case (ii) really occurs for a very close configuration. It is worth to note here that the decrease of $`\chi ^{\mathrm{irrot}}`$ by the term $`F`$ occurs only near the stellar surface. On the contrary, the term $`F`$ increases $`\chi ^{\mathrm{irrot}}`$ around the center of the star.

Next, we show the surface value of y-axis component of $`\stackrel{}{}\mathrm{\Psi }_0`$ and its derivative along an evolutionary sequence in Fig. 23. It is found that the y-axis component of $`\stackrel{}{}\mathrm{\Psi }_0`$ is smaller than that of $`𝐖_s`$ throughout the sequence even if we extrapolate the lines to $`d/a_12.1`$ for $`\gamma =3`$ and to $`d/a_12.25`$ for $`\gamma =1.8`$. Moreover, the y-axis component of the derivative of $`\stackrel{}{}\mathrm{\Psi }_0`$ becomes larger than that of $`𝐖_s`$ for every $`\gamma `$ when the orbital separation decreases than $`d/a_1<2.5`$. This explains why $`\chi ^{\mathrm{irrot}}`$ becomes smaller than $`\chi ^{\mathrm{synch}}`$ for $`d/a_1<2.5`$ as we have found from Fig. 19. We can also see from Fig. 23 that the orbital separation where the derivative of $`\stackrel{}{}\mathrm{\Psi }_0`$ overcomes that of $`𝐖_s`$ is larger for smaller $`\gamma `$. This fact leads to the earlier appearance of a cusp for smaller $`\gamma `$.

Finally, let us discuss the dependence on the orbital separation in the y-axis component of $`\stackrel{}{}\mathrm{\Psi }_0`$ (and also its derivative). We present the y-axis component of $`\stackrel{}{}\mathrm{\Psi }_0`$ and its derivative along an evolutionary sequence in log-log plot in Fig. 24. It is found that both y-axis component of $`\stackrel{}{}\mathrm{\Psi }_0`$ and its derivative behave as $`O[(d/R_0)^3]`$ for the case of larger separation. For very close case as $`d/R_0<3`$, the y-axis component of $`\stackrel{}{}\mathrm{\Psi }_0`$ becomes as $`O[(d/R_0)^{12}]`$ and its derivative reaches $`O[(d/R_0)^{18}]`$.

### C Turning points of total energy

We show values of $`\chi `$ and separations $`d_G/R_0`$ at which the total energy (and/or total angular momentum) takes its minimum along a sequence as a function of the adiabatic index $`\gamma `$ in Fig. 25. We are interested in the turning point of total energy (and/or total angular momentum) because it corresponds to the onset of secular instability in the synchronized case and that of dynamical instability in the irrotational one at least for the ellipsoidal approximation. Note that the turning points in the total energy and total angular momentum coincide. As discussed in Sec. V B, $`\chi =0`$ denotes the end point of equilibrium sequences for both synchronized and irrotational cases. If we extrapolate the results up to the zero value of $`\chi `$, we can obtain the critical value of $`\gamma `$ below which the turning points of the total energy does no longer exist, and the value of the corresponding separation.

We can see from Fig. 25 that $`\chi `$ seems to be zero at $`\gamma =1.71.8`$ with $`d_G/R_0=2.72.75`$ in the synchronized case and at $`\gamma =2.22.3`$ with $`d_G/R_0=2.82.85`$ in the irrotational one.

Taking into account the appearance the cusp and the existence of the turning point, we can expect that the subsequent merger process stars from

* the turning point for $`\gamma 1.8`$ (plunge ?),
* the contact point for $`\gamma 1.8`$

in the synchronized case, and from

* the turning point for $`\gamma 2.3`$ (plunge ?),
* the cusp point with mass shedding $`\gamma 2.3`$

in the irrotational one.

## VI Summary

We have studied equilibrium sequences of both synchronized and irrotational binary systems in Newtonian gravity with adiabatic indices $`\gamma =3,2.5,2.25,2`$ and 1.8. Through the present article, we have understood the qualitative differences of physical quantities between synchronized binary systems and irrotational ones. The summary of the results is as follows:

* The two stars contact with each other in the synchronized case as an end point of equilibrium sequence.
* Irrotational sequences may terminate instead by a cusp point (mass-shedding) corresponding to a detached configuration.
* The turning points of total energy appear in the cases of $`\gamma 1.8`$ in the synchronized case.
* The turning points of total energy appear in the cases of $`\gamma 2.3`$ in the irrotational case.
* The turning points of total energy and total angular momentum coincide with each other.

###### Acknowledgements.

We thank Koji Uryu for providing us with unpublished results from his numerical code. The code development and the numerical computations have been performed on SGI workstations purchased thanks to a special grant from the C.N.R.S. We warmly thank the referee for his detailed comments which helped us to improve the quality of the article.


text:
# Fermion vacuum energies in brane world models

## Abstract

The fermion representations and boundary conditions in five dimensional anti de Sitter space are described in detail. In each case the one loop effective action is calculated for massless fermions. The possibility of topological or Wilson loop symmetry breaking is discussed.

Brane world models contain two types of fields, some restricted to four dimensional sheets and some living in the higher dimensional bulk . The quantum properties of those fields that exist in the bulk is a rich subject. In this paper we focus on a simple example, namely massless fermions in five dimensions, and consider the effective action and the possibility of topological symmetry breaking .

One of the issues which arises is the contribution that vacuum fluctuations make towards the stability of two parallel branes. This is important for the Randall-Sundrum scenario, which relates the mass hierarchy problem to brane separations in anti-de Sitter space . The effects of vacuum fluctuations have already been considered for scalar and fermion fields . However, we feel the need to clarify some statements that have been made about the fermion boundary conditions in these models and to calculate the effective action for a variety of boundary conditions.

The theory we consider is five dimensional. The fifth dimension is taken to be an orbifold $`S_1/Z_2`$, where the circle runs from $`y=a`$ to $`y=a`$ and the $`Z_2`$ acts by $`yy`$. The space is equivalent to a five dimensional spacetime with two four dimensional branes making up the boundary.

The choice of fermion representations in five dimensions is as wide as it is in four. To start with, we might require that the Lagrangian be invariant under the full five dimensional Lorentz group. The Lorentz symmetry in this case would then only be broken by the presence of the brane worlds. The fermions transform by a matrix $`S`$, related to a set of gamma matrices $`\mathrm{\Gamma }^a`$. For the Lorentz transformation $`yy`$,

$`S^1\mathrm{\Gamma }^5S`$ $`=`$ $`\mathrm{\Gamma }^5`$ (1)
$`S^1\mathrm{\Gamma }^\mu S`$ $`=`$ $`\mathrm{\Gamma }^\mu `$ (2)
The smallest representation of the gamma matrices which can satisfy these relations has eight component spinors. The existence of a matrix $`\mathrm{\Gamma }^6`$ which anticommutes with the other gamma matrices allows $`S=i\mathrm{\Gamma }^5\mathrm{\Gamma }^6`$. The benefits of using eight component spinors have also been emphasised in discussions of fermion boundary value problems .

The eight dimensional representation can be reduced to a real eight dimensional Majorana representation, where $`\psi =C\psi ^{}`$ and $`C`$ satisfies $`C^1\mathrm{\Gamma }^aC=\mathrm{\Gamma }^a`$. These are the fermion representations which arise from the reduction of supersymmetric theories . In the supersymmetry literature, the eight dimensional fermions are usually regarded as a pair of four dimensional fermions related by a symplectic transformation .

If, instead of the full Lorentz group, we require only symmetry under proper Lorentz transformations, then the representation reduces to four dimensional Weyl representations. These are the representations that have been considered hitherto in the context of the Randall-Sundrum scenario . A rule for transforming spinors under the transformation $`yy`$ is still required. For this purpose, we can allow the massless Dirac equation to transform as a pseudoscalar,

$`S^1\mathrm{\Gamma }^5S`$ $`=`$ $`\mathrm{\Gamma }^5`$ (3)
$`S^1\mathrm{\Gamma }^\mu S`$ $`=`$ $`\mathrm{\Gamma }^\mu `$ (4)
The solution is $`S=i\mathrm{\Gamma }^5=\gamma ^5`$, where $`\mathrm{\Gamma }^\mu =\gamma ^\mu `$ are the usual gamma matrices in four dimensions. We will refer to these fermions as ‘five dimensional Weyl fermions’. They have the useful property that they induce a chiral particle theory on the branes .

The boundary conditions used in this paper are chosen for consistency with orbifold reductions of the fifth dimension. The fermions carry a representation of the $`Z_2`$ symmetry, hence $`\psi (y)=\pm S\psi (y)`$. Therefore, if

$$P_\pm =\frac{1}{/}2(1\pm S),$$
(5)
we must impose one of the two equivalent conditions $`P_\pm \psi =0`$ at $`y=0`$.

The fact that $`y`$ lies on a circle would normally imply that $`\psi (a)=\psi (a)`$. However, if there is a symmetry it is possible to make the identification up to a symmetry transformation. The possibility $`\psi (a)=\psi (a)`$ has already been used as a mechanism for breaking supersymmetry, . More general possibilities which would allow gauge symmetry breaking are discussed later, but for the moment we have two inequivalent boundary conditions

$`y=0:P_{}\psi =0,`$ $`y=a:P_{}\psi =0`$ (6)
$`y=0:P_+\psi =0,`$ $`y=a:P_{}\psi =0`$ (7)
The first set might be regarded as untwisted and the second set twisted. The twisted case has been considered in flat five dimensional models by Antoniadis et al. . For Weyl fermions, the untwisted boudary condition agrees with the boundary conditions used by Grossman et al. and gives results similar to Garriga et al.. We will give results for both Dirac and Weyl fermions and for both twisted and untwisted cases.

We will take the metric

$$ds^2=e^{2\sigma }\eta _{\mu \nu }dx^\mu dx^\nu +dy^2.$$
(8)
Anti de Sitter space corresponds to $`\sigma =\kappa y`$, $`\kappa `$ being a constant. This metric is conformally flat,

$$ds^2=e^{2\sigma }(\eta _{\mu \nu }dx^\mu dx^\nu +d\tau ^2)$$
(9)
where $`0\tau \beta `$, and

$$\beta [\sigma ]=_0^ae^\sigma 𝑑y.$$
(10)
For anti de Sitter space, $`\beta =\kappa ^1(e^{\kappa a}1)`$.

When boundaries are present it is convenient to regard the Dirac operator $`D`$ mapping one set of fermions into an image set. The adjoint mapping is denoted by $`D^{}`$. The one loop contribution to the effective action is then

$$W=\frac{1}{/}2\mathrm{log}det(D^{}D).$$
(11)
The boundary conditions on the image fermions can be determined by the existence of $`D^{}`$. If $`P_{}\psi =0`$, this requires that the normal derivative of $`P_+\psi `$ should vanish.

In the massless case, $`D=D^{}=i\mathrm{\Gamma }^j_j`$ where $`j`$ runs from $`1`$ to $`5`$. The conformal transformation properties of the massless operator imply $`D=e^{3\sigma }D_0e^{2\sigma }`$, where $`D_0`$ is the Dirac operator in the strip of flat space $`0\tau \beta `$. The boundary conditions are also conformally invariant. We can relate the effective action to the result in flat space by

$$W=W_0+C[\sigma ]$$
(12)
where $`W_0=\frac{1}{/}2\mathrm{log}det(D_0^{}D_0)`$ and $`C[\sigma ]`$ is a correction term . We shall discuss the significance of this term later.

For untwisted Dirac fields we have $`P_{}\psi =0`$ and $`(P_+\psi )/\tau =0`$ on either boundary. The eigenvalues of $`D_0^{}D_0=^2`$ are then $`k^2+m_n^2`$, with $`m_n=\pi n/\beta `$, where $`n=0,1,2\mathrm{}`$. The degeneracy $`g=8`$ for each value of $`k`$. Twisted Dirac fields have similar eigenvalues with $`m_n=(n+\frac{1}{/}2)\pi /\beta `$. For Weyl fermions, the eigenfunctions and eigenvalues are unchanged, but now the degeneracy factors are $`g=4`$ rather than $`g=8`$.

The logarithms can be evaluated using $`\zeta `$-function regularisation . For untwisted fields,

$$W_0=d^4x\frac{3g}{128\pi ^2}\zeta _R(5)\beta ^4$$
(13)
where $`\zeta _R`$ is the Riemann $`\zeta `$-function. For twisted fields,

$$W_0=d^4x\frac{3g}{128\pi ^2}\frac{15}{16}\zeta _R(5)\beta ^4.$$
(14)
The results depend on the separation of the branes only through $`\beta `$ given in equation (10). The untwisted Weyl case gives the same result as that obtained by Garriga et. al. .

In this particular problem there is no dependence on the renormalisation scale. The same result can be obtained by dimensional regularisation where the absence of pole terms also indicates no dependence on the renormalisation scale. The situation changes when the branes are curved, and renormalisation scale dependent curvature terms arise .

The quantity $`W_0`$ is also the one loop correction to the effective action of an infinite set of particles in four dimensions with mass $`m_n`$. The difference between $`W`$ and $`W_0`$, namely the cocycle function $`C[\sigma ]`$, can be regarded as an anomaly in this reduction. Such anomalies have been recognised by Frolov et al. . In general, $`C[\sigma ]`$ will depend on the geometry of the branes, but in the present context the anomaly only contributes a constant term to the matter Lagrangians $`_v`$ and $`_h`$ on the ‘visible’ and ‘invisible’ branes.

In the Randall-Sundrum metric (8), the classical action reduces to

$$S=d^4xe^{4\kappa a}\left(_v\frac{3\kappa }{4\pi G_5}\right)d^4x\left(_h+\frac{3\kappa }{4\pi G_5}\right),$$
(15)
where $`G_5`$ is the gravitational constant in five dimensions. Junction conditions on the metric imply that both terms vanish in the vacuum. Adding the correction $`W`$ for untwisted fermions gives an effective action which now has a minimum for a particular separation $`a`$. However, the values obtained cannot give the correct mass hierarchy ($`\kappa a>30`$) without a considerable degree of fine tuning.

We can also derive the vacuum energy density from the effective action quite simply,

$$T_{\mu \nu }=\frac{2}{\sqrt{g}}\frac{\delta W}{\delta g^{\mu \nu }}=e^{5\sigma }\frac{dW}{d\beta }g_{\mu \nu }.$$
(16)
For the untwisted anti-deSitter case, $`\sigma =\kappa y`$, the energy density is

$$\frac{3g}{32\pi ^2}\zeta _R(5)\kappa ^4(e^{\kappa (ay)}e^{\kappa y})^5.$$
(17)
Since this is strongly concentrated near the visible brane, the back reaction of this energy modifies the junction conditions, resulting in a value for $`a`$ in agreement with the minimum of the effective action.

If there is a gauge symmetry, the possibility of topological or Wilson loop symmetry breaking arises. The boundary conditions can be generalised by inserting gauge transformations $`U_h`$ and $`U_v`$ on the two branes, so that now

$$P_{}=\frac{1}{/}2(1SU),$$
(18)
where $`U=U_h`$ or $`U=U_v`$. The condition $`P_{}^2=P_{}`$ requires $`U_h^2=U_v^2=I`$. The symmetry is broken, leaving the centraliser of $`U_h`$ and $`U_v`$, which preserves the boundary conditions, as the residual symmetry group.

For a simple, non-trivial, example, consider the group $`U(2)`$ with $`U_h=\sigma _3`$, the Pauli matrix. If $`U_v=I`$, the residual symmetry group is $`U(1)\times U(1)`$ and the fermions decompose into one twisted and one untwisted fermion. The combined one loop correction is therefore $`\frac{1}{/}16W_0`$, where $`W_0`$ is the untwisted result (13). If $`U_v=\pm \sigma _3`$, the residual symmetry group is the same but the fermions are both twisted or both untwisted, giving a correction $`2W_0`$ or $`\frac{15}{/}8W_0`$. The final case is represented by $`U_v=\sigma _1`$, and the residual symmetry group is $`U(1)`$. The eigenvalues are now of the form $`k^2+m_n^2`$, with $`m_n=(n\pm \frac{1}{/}4)\pi /\beta `$. The effective action can be calculated as before, and takes the value $`\frac{15}{/}256W_0`$.

For massless fermions, the separation between the two branes is only stable when the correction to the action is positive. Clearly, there is an interplay in this scenario between supersymmetry breaking, gauge symmetry breaking and the mass hierarchy problem. We are presently extending these results to massive fermions so that we can investigate the consequences in low energy superstring models.

###### Acknowledgements.

Antonino Flachi is supported by a University of Newcastle Upon Tyne Ridley Studentship.


text:
# 1 Introduction

## 1 Introduction

The subject of quantum field theory in curved spacetime studies the interrelation between quantum theory and spacetime geometry in the approximation that the geometry is remained as a classical background. This subject, however, contains some basic problems concerning inherent ambiguities in the definition of physical states associated with a quantum field. In the absence of a gravitational background, it is always possible to use the Poincare symmetries to obtain a physical vacuum (the state of lowest energy). One may then assume that the physically admissible states most naturally arise as local excitations of this state. In the presence of a curved spacetime, however, this procedure does not apply, because on a general curved sapcetime one may not find a global symmetry. In this case the problem concerning the determination of the local states and the role of the global features of spacetime is of obvious importance.
The Hadamard state-condition provides a framework in which we may improve our understanding in this context. In this framework one postulates that the short-distance singularity structure of the two-point function of a linear quantum field is represented by the Hadamard ansatz . In this prescription, however, there exist problems in the specification of the state-dependent part of the two-point function. For characterizing the physical states, it therefore seems to be essential to find out a suitable scheme for the treatment of these problems. In the present paper we shall consider this issue for a quantum scalar field conformaly coupled to a two-dimensional gravitational background.
The organization of this paper is as follows: In section 2, we use the Hadamard ansatz for the derivation of the local constraints imposed on the state-dependent part of the two-point function. In section 3, we develop a dynamical model in order to analyze these constraints. This dynamical model, which is a two-dimensional analog of a model considered in the previous publication , uses a conformal invariant C-number field to characterize the stress-tensor associated with the local states. The introduction of this scalar field into the analysis leads to a basic connection between the specification of the local states and the determination of a conformal frame. There is a general consistency requirement on the behavior of the C-number field which necessitates the use of a conformal frame in which the trace of the stress-tensor is measured by a cosmological constant. This leads to the establishment of a two-dimensional gravitational equation in this frame. In section 4, we offer some concluding remarks. Throughout the following we shall work in units in which $`\mathrm{}=c=1`$ and follow the sign conventions of Hawking and Ellis .

## 2 Hadamard states in two spacetime dimensions

We consider a free massless quantum scalar field $`\varphi (x)`$ propagating in a two-dimensional gravitational background with the action functional (In the following the semicolon and $``$ denote covariant differentiation)

$$S[\varphi ]=\frac{1}{2}d^2x\sqrt{g}g_{\alpha \beta }^\alpha \varphi ^\beta \varphi .$$
(1)
This action gives rise to the field equation

$$\mathrm{}\varphi (x)=0.$$
(2)
A state of $`\varphi (x)`$ is characterized by a hierarchy of wightman functions

$$\varphi (x_1),\mathrm{},\varphi (x_n).$$
(3)
We are primarily interested in quasi-free states, i.e., states for which the truncated n-point functions vanish for $`n>2`$. Such states may be characterized by their two-point functions. In a linear theory the antisymmetric part of the two-point function is common to all states in the same representation. It is just the universal commutator function. Thus, in our case all the relevant information about the state-dependent part of the two-point function are encoded in its symmetric part, denoted in the following by $`G^+(x,x^{})`$, which satisfies Eq.(2) in each argument. Equivalence principle suggests that the leading singularity of $`G^+(x,x^{})`$ should have a close correspondence to the singularity structure of the two-point function of a free massless field in a two-dimensional Minkowski spacetime. One may therefore assume that for $`x`$ sufficiently close to $`x^{}`$ the function $`G^+(x,x^{})`$ can be written as

$$G^+(x,x^{})=\frac{1}{4\pi }\mathrm{ln}\sigma (x,x^{})+F(x,x^{}),$$
(4)
where $`\sigma (x,x^{})`$ is one-half the square of the geodesic distance between $`x`$ and $`x^{}`$, and $`F(x,x^{})`$ is a regular function. This may be viewed as a two-dimensional analog of the Hadamard ansatz. It should be remarked that in general there is a missing mass scale in the expression (4) emerging from the fact that the argument of the logarithm must be dimensionless. We shall return to this point in the next section.
The function $`F(x,x^{})`$ satisfies a general constraint obtaining from the symmetry condition of $`G^+(x,x^{})`$ and requiring that the expression (4) satisfies the wave equation (2). The study of this constraint has obviously a particular significance for analyzing the state-dependent part of the two-point function. We therefore present its derivation as follows: Applying Eq.(2) on (4) leads to

$$\sigma (x,x^{})\mathrm{}F(x,x^{})=\frac{1}{4\pi }(\sigma _{;\alpha }^{;\alpha }2),$$
(5)
in which we have used the definition of $`\sigma (x,x^{})`$, namely

$$\sigma (x,x^{})=\frac{1}{2}g_{\alpha \beta }\sigma ^{;\alpha }\sigma ^{;\beta }.$$
(6)
The explicit form of the right hand side of (5) can be specified by writing the covariant expansion of $`\sigma ^{;\alpha \beta }`$

$$\sigma ^{;\alpha \beta }=g^{\alpha \beta }\frac{1}{3}R_{\gamma \delta }^{\alpha \beta }\sigma ^{;\gamma }\sigma ^{;\delta }+\frac{1}{12}R_{\gamma \delta ;\rho }^{\alpha \beta }\sigma ^{;\gamma }\sigma ^{;\delta }\sigma ^{;\rho }+O(\sigma ^2).$$
(7)
Taking the trace of this relation and using the fact that in a two-dimensional spacetime the Einstein tensor identically vanishes

$$R_{\alpha \beta }=\frac{1}{2}g_{\alpha \beta }R,$$
(8)
leads to

$$\sigma _{;\alpha }^{;\alpha }=2\frac{1}{6}g_{\gamma \delta }R\sigma ^{;\gamma }\sigma ^{;\delta }+\frac{1}{24}g_{\gamma \delta }R_{;\rho }\sigma ^{;\gamma }\sigma ^{;\delta }\sigma ^{;\rho }+O(\sigma ^2).$$
(9)
Combining (5) and (9), yields

$$\sigma (x,x^{})\mathrm{}F(x,x^{})=\frac{1}{24\pi }g_{\gamma \delta }R\sigma ^{;\gamma }\sigma ^{;\delta }+\frac{1}{96\pi }g_{\gamma \delta }R_{;\rho }\sigma ^{;\gamma }\sigma ^{;\delta }\sigma ^{;\rho }+O(\sigma ^2).$$
(10)
Substituting (6) and the covariant expansion of $`F(x,x^{})`$ , namely

$$F(x,x^{})=F(x)\frac{1}{2}F_{;\alpha }(x)\sigma ^{;\alpha }+\frac{1}{2}F_{\alpha \beta }(x)\sigma ^{;\alpha }\sigma ^{;\beta }+\frac{1}{4}\{\frac{1}{6}F_{;\alpha \beta \gamma }(x)F_{\alpha \beta ;\gamma }(x)\}\sigma ^{;\alpha }\sigma ^{;\beta }\sigma ^{;\gamma }+O(\sigma ^2),$$
(11)
into (10), one finds

$$\frac{1}{2}F_\eta ^\eta (x)g_{\gamma \delta }\sigma ^{;\gamma }\sigma ^{;\delta }+\frac{1}{2}g_{\gamma \delta }[F_{\alpha \beta }^{;\alpha }(x)\frac{1}{2}F_{\eta ;\beta }^\eta (x)+\frac{1}{12}(\mathrm{}F(x))_{;\beta }\frac{1}{3}\mathrm{}(F_{;\beta }(x))\frac{1}{12}RF_{;\beta }(x)]\sigma ^{;\gamma }\sigma ^{;\delta }\sigma ^{;\beta }+O(\sigma ^2)$$

$$=\frac{1}{24\pi }g_{\gamma \delta }R\sigma ^{;\gamma }\sigma ^{;\delta }+\frac{1}{96\pi }g_{\gamma \delta }R_{;\beta }\sigma ^{;\gamma }\sigma ^{;\delta }\sigma ^{;\beta }+O(\sigma ^2).$$
(12)
We may now compare term by term up to the third order of $`\sigma ^{;\alpha }`$ to obtain the desired relations

$$F_\eta ^\eta (x)=\frac{1}{12\pi }R,$$
(13)
$$F_{\alpha \beta }^{;\alpha }(x)\frac{1}{2}F_{\eta ;\beta }^\eta (x)+\frac{1}{12}(\mathrm{}F(x))_{;\beta }\frac{1}{3}\mathrm{}(F_{;\beta }(x))\frac{1}{12}RF_{;\beta }(x)=\frac{1}{48\pi }R_{;\beta }.$$
(14)
The equations (13) and (14) are to be regarded as two independent constraints imposed on the state-dependent part of the two-point function. We may combine them to obtain

$$F_{\alpha \beta }^{;\alpha }(x)+\frac{1}{12}(\mathrm{}F(x))_{;\beta }\frac{1}{3}\mathrm{}(F_{;\beta }(x))\frac{1}{12}RF_{;\beta }(x)=\frac{1}{48\pi }R_{;\beta }.$$
(15)
It should be noted that in the derivation of this constraint we have used the covariant expansion of $`F(x,x^{})`$ and $`\sigma ^{;\alpha \beta }`$ only up to the second order terms. In general there exist additional constraints on the higher order expansion terms. However in our analysis we shall neglect these higher order constraints and concentrate our attention on (15).

## 3 The large scale behavior of local states

In the first step of analyzing the constraint (15) we note that it can be written as a total divergence. In fact we may use the relation (8) and the differential identity

$$\mathrm{}(F_{;\beta }(x))=(\mathrm{}F(x))_{;\beta }+R_{\alpha \beta }F^{;\alpha }(x),$$
(16)
to get the equation

$$^\alpha \mathrm{\Sigma }_{\alpha \beta }=0,$$
(17)
where

$$\mathrm{\Sigma }_{\alpha \beta }=\mathrm{\Sigma }_{\alpha \beta }^{(0)}+\mathrm{\Sigma }_{\alpha \beta }^{(1)},$$
(18)
and

$$\mathrm{\Sigma }_{\alpha \beta }^{(0)}=\frac{1}{2}(F_{;\alpha \beta }(x)\frac{1}{2}g_{\alpha \beta }\mathrm{}F(x))(F_{\alpha \beta }+\frac{1}{24\pi }g_{\alpha \beta }R),$$
(19)
$$\mathrm{\Sigma }_{\alpha \beta }^{(1)}=\frac{1}{48\pi }g_{\alpha \beta }R.$$
(20)
The tensor $`\mathrm{\Sigma }_{\alpha \beta }`$ is decomposed into a traceless part, which is denoted by $`\mathrm{\Sigma }_{\alpha \beta }^{(0)}`$, and the tensor $`\mathrm{\Sigma }_{\alpha \beta }^{(1)}`$ which leads to the trace anomaly in two-dimensions. The principal strategy now is to relate the tensor $`\mathrm{\Sigma }_{\alpha \beta }^{(0)}`$ to the traceless stress-tensor of a conformal invariant scalar field which is given by

$$T_{\alpha \beta }=_\alpha \psi _\beta \psi \frac{1}{2}g_{\alpha \beta }_\gamma \psi ^\gamma \psi ,$$
(21)
where $`\psi (x)`$ is a C-number scalar field satisfying the dynamical equation

$$\mathrm{}\psi (x)=0.$$
(22)
In this case the tensor $`\mathrm{\Sigma }_{\alpha \beta }`$ takes the form

$$\mathrm{\Sigma }_{\alpha \beta }=T_{\alpha \beta }+\frac{1}{48\pi }g_{\alpha \beta }R.$$
(23)
This relation indicates how the local states are related to the background geometry. Taking its trace we obtain

$$\mathrm{\Sigma }_\alpha ^\alpha =\frac{1}{24\pi }R,$$
(24)
which is the trace anomaly of a two-dimensional quantum stress-tensor of a conformal invariant scalar field. This raises the possibility of regarding $`\mathrm{\Sigma }_{\alpha \beta }`$ as the quantum stress-tensor induced by the two-point function. However by comparing (23) with (17) we observe that the existence of the trace anomaly on the background metric is in conflict with the interpretation of $`T_{\alpha \beta }`$ as a conserved stress-tensor. This observation reveals that the connection between the local properties of the state-dependent part of the two-point function with those of a C-number scalar field, $`\psi (x)`$, can not be consistently established on the background metric. This apparent discrepancy, however, may be removed by appealing to the conformal symmetry of the above model. The conformal invariance of Eq.(22) implies that at dynamical level it is not possible to single out a particular configuration for $`\psi (x)`$ among many different conformally related configurations. In order to determine a configuration one has to choose a particular conformal frame. This therefore raises the question of which of these conformal frames corresponds to the physical one. In general the choice of a conformal frame depends on the physical conditions one wishes to impose on such a frame and these conditions are suggested by the problem under consideration. Here we shall choose the conformal frame by demanding that $`T_{\alpha \beta }`$ shall be consistent with the conservation property of a stress-tensor. We first consider a conformal transformation

$$g_{\alpha \beta }\overline{g}_{\alpha \beta }=\mathrm{\Omega }^2g_{\alpha \beta },$$
(25)
$$\psi (x)\psi (x),$$
(26)
where $`\mathrm{\Omega }`$ is a nonvanishing, smooth spacetime function. Then we write Eq.(23) in a conformal frame describing by the metric $`\overline{g}_{\alpha \beta }`$ so that

$$\overline{\mathrm{\Sigma }}_{\alpha \beta }=\overline{T}_{\alpha \beta }+\frac{1}{48\pi }\overline{g}_{\alpha \beta }\mathrm{\Lambda },$$
(27)
where we have set

$$\overline{R}\mathrm{\Lambda }=0,$$
(28)
or, equivalently,

$$_\gamma \mathrm{ln}\mathrm{\Omega }^\gamma \mathrm{ln}\mathrm{\Omega }=\frac{1}{2}\mathrm{\Lambda }\mathrm{\Omega }^2\frac{1}{2}R+\frac{\mathrm{}\mathrm{\Omega }}{\mathrm{\Omega }}.$$
(29)
This is a differential equation which determines a particular conformal factor characterizing a conformal frame in which the relation (27) holds. Equation (28) corresponds to a two-dimensional analog of the vacuum Einstein field equations with a cosmological constant. This clearly emphasizes the significance of the gravitational coupling of local states in the choice of a conformal frame as a physical one. We should also remark that $`\mathrm{\Lambda }`$ provides now the missing mass-scale mentioned previously in connection with the Hadamard ansatz. This mass-scale presents itself through the root of an effective cosmological constant in a two-dimensional gravitational theory.
To make a closer look at the relation (29), let us investigate it for a particular case in which the background geometry is described by a two-dimensional Schwarzschild metric

$$ds^2=(1\frac{2GM}{r})dt^2+(1\frac{2GM}{r})^1dr^2.$$
(30)
In this case the trace anomaly is $`\mathrm{\Sigma }_\alpha ^\alpha =\frac{GM}{6\pi r^3}`$. To solve Eq.(29) in this particular case one needs some boundary conditions. One may intuitively expect that for a sufficiently small cosmological constant the conformal frame, characterized by (29), and the asymptotic regions of the background metric should have the same properties . We therefore choose the boundary conditions so that an asymptotic correspondence can be established between the conformal frame and the background frame in the Schwarzschild coordinate system. Since the trace anomaly vanishes when $`r\mathrm{}`$, Eq.(29) reduces to

$$_\gamma \mathrm{ln}\mathrm{\Omega }^\gamma \mathrm{ln}\mathrm{\Omega }=\frac{1}{2}\mathrm{\Lambda }\mathrm{\Omega }^2+\frac{\mathrm{}\mathrm{\Omega }}{\mathrm{\Omega }}.$$
(31)
For a sufficiently small $`\mathrm{\Lambda }`$ and in a static case, in which $`\mathrm{\Omega }`$ is only a function of $`r`$, this equation has a simple solution, namely $`\mathrm{\Omega }(r)=e^{(ar+b)}`$ with $`a`$ and $`b`$ being arbitrary constants. Applying this conformal factor to (30) gives in the asymptotic regions

$$ds^2=e^{2(ar+b)}(dt^2+dr^2).$$
(32)
The case $`a=b=0`$ corresponds to the Minkowski spacetime and the conformal frame meets the background metric for $`r\mathrm{}`$ in the Schwarzschild coordinate system. Thus $`a=b=0`$ characterizes the physically admissible conformal factor. When $`a0`$, a coordinate transformation

$$\eta =\frac{1}{a}e^{(ar+b)}\mathrm{sinh}(at),$$

$$\xi =\frac{1}{a}e^{(ar+b)}\mathrm{cosh}(at),$$
(33)
would change (32) to

$$ds^2=d\eta ^2+d\xi ^2.$$
(34)
demonstrating that the coordinates $`(t,r)`$ are the Rindler-type coordinates  associated with an accelerated observer in a flat spacetime. In this case the asymptotic correspondence between the two frames can not be established in the Schwarzschild coordinate system.

## 4 Concluding remarks

We have analyzed an analog of the Hadamard ansatz for the specification of local states associated with a linear quantum field conformally coupled to a two-dimensional gravitational background. For studying the constraints emerging from this ansatz we have employed a conformal invariant classical scalar field. The conformal invariance permits us to introduce a conformal frame in which the gravitational coupling of local states is established. Such a frame is constructed by relating the trace of the quantum stress-tensor to a cosmological constant. We have shown that the gravitational coupling of the local states leads to a two-dimensional gravitational theory.
Acknowledgment
The authors would like to acknowledge the financial support of the Office of Scientific Research of Shahid Beheshti University.


text:
# Spin injection through the depletion layer: a theory of spin-polarized p-n junctions and solar cells

\[

## Abstract

A drift-diffusion model for spin-charge transport in spin-polarized p-n junctions is developed and solved numerically for a realistic set of material parameters based on GaAs. It is demonstrated that spin polarization can be injected through the depletion layer by both minority and majority carriers, making all-semiconductor devices such as spin-polarized solar cells and bipolar transistors feasible. Spin-polarized p-n junctions allow for spin-polarized current generation, spin amplification, voltage control of spin polarization, and a significant extension of spin diffusion range.

\]

Spintronics has played an important role in defining novel applications that are either not feasible or ineffective with traditional semiconductor electronics. Spintronic devices have found their niche in industries for magnetic read heads and nonvolatile memory cells. Here we propose and demonstrate a new scheme for spintronics, a spin-polarized p-n junction, which amplifies spin density, significantly extends the range of spin diffusion, electronically tailors spin polarization, and, in combination with light as a spin-polarized solar sell, generates spin-polarized currents with tunable spin polarization. We prove these concepts by solving drift-diffusion equations for a realistic device model based on GaAs, which demonstrates that spin polarization can be injected through the depletion layer. Possibility of injecting spin polarization through a transistor is also discussed.

The electrical injection of spin-polarized carriers within all-semiconductor structures (from a magnetic into a nonmagnetic semiconductor) was recently realized experimentally  (the scheme proposed in ). Optical injection of spin-polarized carriers (both minority  and majority ) has been known for some time. In addition, the relatively long spin diffusion lengths, coherent spin transport across semiconductor interfaces, a successful fabrication of a magnetic/nonmagnetic p-n junction based on the novel (Ga,Mn)As material, and the recent demonstration of a gate-voltage control of magnetization in (In,Mn)As, make semiconductors promising materials for spintronic applications .

In this paper we investigate the spin-charge transport in semiconductors under the conditions of inhomogeneous bipolar doping (Flatte and Vignale have recently made an interesting proposal for spin diodes and transistors in unipolar semiconductor heterostructures–a very different case from ours); we are not concerned with spin injection per se. Our model device is a spin-polarized p-n junction with spin polarization induced (either optically–in which case we get a spin-polarized solar cell–or electronically) to minority or majority carriers. By studying spin-charge transport numerically across the depletion layer, we observe novel phenomena, all resulting from the fact that spin polarization is transferred (what we call injected) through the depletion layer.

We first introduce a drift-diffusion model for spin-charge transport in an inhomogeneously doped semiconductor illuminated with circularly polarized light (and later solve the model for GaAs). In addition to the approximations used in deriving the usual (unpolarized) equations , we assume that all the dopants are fully ionized, the carrier populations nondegenerate and varying only in one ($`x`$) direction. Further, we assume the spin polarization to be carried solely by conduction electrons (that is, consider holes unpolarized), as is the case of III-V semiconductors like GaAs, best candidates for photo-spintronics.

The following parameters describe our model: acceptor (donor) densities $`N_A`$ ($`N_D`$); electron (hole) densities $`n`$ ($`p`$), in equilibrium $`n_0`$ ($`p_0`$); intrinsic carrier density $`n_i`$ ($`n_i^2=n_0p_0`$); electron (hole) number current densities $`J_n`$ ($`J_p`$); electron (hole) mobilities and diffusivities $`\mu _n`$ ($`\mu _p`$) and $`D_n`$ ($`D_p`$); bipolar photoexcitation rate $`G`$; intrinsic generation-recombination rate constant $`w`$ and spin-relaxation time $`T_1`$. Electron parameters $`n`$, $`J_n`$, and $`G`$ will carry spin index $`\lambda `$ ($`\lambda =1,`$ for spin up and $`\lambda =1,`$ for spin down): $`n_{}+n_{}=n`$, $`J_n+J_n=J_n`$, and $`G_{}+G_{}=G`$. We also define spin-related quantities: spin density $`s=n_{}n_{}`$, spin polarization $`\alpha =s/n`$, and spin photoexcitation rate $`G^s=G_{}G_{}`$.

The time evolution and spatial distribution of carrier and spin densities is described by three sets of equations. (i) Poisson’s equation $`d^2\varphi /dx^2=\rho /ϵ`$, where $`\varphi `$ is the electrostatic potential (related to the electric field $`E=d\varphi /dx`$ in the $`x`$ direction), $`ϵ`$ is the sample dielectric permittivity, and $`\rho =e(N_DN_An+p)`$ is the local charge density with the elementary charge $`e`$. (ii) The linear response equations for the particle currents, $`J_{n\lambda }=\mu _nn_\lambda ED_n(dn_\lambda /dx)`$ and $`J_p=+\mu _ppED_p(dp/dx)`$. The mobilities and conductivities are connected via Einstein’s relation $`eD=k_BT\mu `$; $`k_B`$ is the Boltzmann constant, and $`T`$ the absolute temperature. (iii) The continuity equations

$`{\displaystyle \frac{dn_\lambda }{dt}}+{\displaystyle \frac{dJ_{n\lambda }}{dx}}`$ $`=`$ $`w(n_\lambda pn_0p_0/2){\displaystyle \frac{n_\lambda n_\lambda }{2T_1}}+G_\lambda ,`$ (1)
$`{\displaystyle \frac{dp}{dt}}+{\displaystyle \frac{dJ_p}{dx}}`$ $`=`$ $`w(npn_0p_0)+G,`$ (2)
expressing the particle generation and recombination as well as spin relaxation. Equation sets (i-iii), together with appropriate boundary conditions drawn from the actual physical situation, fully determine the steady state ($`dn_\lambda /dt=dp/dt=0`$) distribution of carrier and spin densities.

In a homogeneous case and steady state Eqs. 1 and 2 become $`w(npn_0p_0)=G`$, and $`wsp+s/T_1=G^s`$. Note that the rate at which electron spin relaxes is not $`1/T_1`$, but rather $`wp+1/T_1`$, reflecting the disappearance of spin by recombination (spin is effectively transferred to holes which then lose it). This is more pronounced in $`p`$-doped samples. Let us see if our equations recover what is already known about spin polarization in semiconductors. For polarization we get $`\alpha =\alpha _0(1+n_0p_0/np)/(1+1/T_1wp)`$, where $`\alpha _0=G^s/G`$ is the polarization at time of pair creation. Take a $`p`$-doped sample. Then $`pp_0`$, $`nn_0`$, and $`\alpha =\alpha _0/(1+\tau /T_1)`$, with electron lifetime $`\tau =1/wp_0`$, expresses spin orientation: if $`\tau T_1`$, which is usually the case, $`\alpha \alpha _0`$ and electron spins are effectively oriented throughout the sample. In an $`n`$-doped sample $`wpG/n_0`$, and the spin polarization is $`\alpha _0/(1+n_0/GT_1)`$. This is optical spin pumping: spin polarization depends on the illumination intensity $`G`$ and is noticeable at $`Gn_0/T_1`$.

Our prototype model device is a $`2`$ $`\mu `$m long GaAs sample, doped with $`N_A=3\times 10^{15}`$ cm-3 acceptors on the left and with $`N_D=5\times 10^{15}`$ cm-3 donors on the right along the $`x`$-axis (the doping profile is shown in Fig. 1). The left surface of the sample, $`x=0`$, is illuminated by circularly polarized light which creates electrons and holes with photodensities $`\delta n=\delta p=3\times 10^{14}`$ cm-3, and induces electronic spin polarization $`\alpha _0=0.5`$ (the value given by the band-structure symmetry) at the surface. Since holes lose their spin orientation faster than they gain it (also a band structure effect), their polarization is effectively zero and need not be considered. The task is to find the steady-state distribution of electron and hole densities $`n`$ and $`p`$, as well as electron spin density $`s=n_{}n_{}`$ and polarization $`\alpha =s/n`$, as a function of $`x`$.

We solve the equation sets (i-iii) numerically, in the steady-state regime, with the boundary conditions: $`p(0)=p_0(0)+\delta p`$ (where $`p_0(0)=0.5(N_A+(N_A^2+4n_i^2)^{0.5})`$), $`n(0)=p_0(0)/n_i^2+\delta n`$, and $`\alpha (0)=0.5`$, reflecting bipolar photoexcitation and spin orientation processes at the illuminated surface (there is no illumination in the bulk: $`G=0`$); $`n(2)=n_0(2)`$ (where $`n_0(2)=0.5(N_D+(N_D^2+4n_i^2)^{0.5})`$), $`p(2)=n_0(2)/n_i^2`$, and $`\alpha (2)=0`$, maintaining equilibrium at the right surface. For electrostatic potential $`\varphi `$ we use $`\varphi (0)=(k_BT/e)\mathrm{ln}(p_0(0)/n_i)`$ and $`\varphi (2)=(k_BT/e)\mathrm{ln}(n_0(2))/n_i)V`$, where $`V`$ is the applied bias, conventionally defined with respect to the dark built-in value. The room-temperature parameters of GaAs used in our model are : the intrinsic carrier density $`n_i=1.8\times 10^6`$ cm-3, electron and hole mobilities $`\mu _n=10\mu _p=4000`$ $`\mathrm{cm}^2\mathrm{V}^1\mathrm{s}^1`$, electron and hole diffusivities (from Einstein’s relation) $`D_n=10D_p=103.6`$ $`\mathrm{cm}^2\mathrm{s}^1`$, and dielectric permittivity $`ϵ=13.1ϵ_0`$ (where $`ϵ_0`$ is the permittivity of free space). As for the generation-recombination rate $`w`$ (assumed to be a constant independent of $`n`$ and $`p`$), we take the value of $`(1/3)\times 10^5`$ $`\mathrm{cm}^3\mathrm{s}^1`$. The electron lifetime in the $`p`$-region is then $`\tau 1/(wN_A)0.1`$ ns. Finally, spin-relaxation time $`T_1`$ is taken to be $`0.2`$ ns, a reasonable value for GaAs .

We first discuss the case of $`V=0`$ (the applied reverse bias from an external battery cancels the forward voltage due to photoexcitation), which is in Fig. 1. Electronic density starts off with the value of $`n\delta n`$ at $`x=0`$, decreases somewhat in the depletion layer, then rises by more than a decade to reach its equilibrium value of $`n_0N_D`$ in the $`n`$ region. The hole density sharply decreases from approximately $`N_A`$ through the depletion layer until it becomes $`n_i^2/N_D`$, the equilibrium value on the $`n`$ side. Spin density $`s`$ is $`0.5n`$ at the illuminated surface, and essentially follows the spatial dependence of $`n`$ through the depletion layer, but once in the $`n`$ region it decays towards zero. The surprising feature is the increase of $`s`$ through the depletion layer. In effect, the magnetization of the sample increases by more than an order of magnitude as a result of spin injection. The polarization $`\alpha `$ (Fig. 1, inset) starts at 0.5 at the illuminated surface, stays almost constant through the transition region, then decreases to zero at the right boundary. We checked that these results are robust against changes in $`T_1`$ and $`\tau `$ by up to two decades, as long as $`T_1\tau `$ (so that appreciable spin polarization can be induced in the $`p`$ region in the first place).

Spin polarization clearly survives the depletion layer. This is not an ordinary spin injection in which certain number of spin-polarized electrons tunnel through a contact and the spin density is equal on both sides. What we have, rather, is a spin pumping (leading to spin density amplification) mechanism. Indeed, in optical spin pumping  circularly polarized light creates spin polarization of majority carriers (electrons in $`n`$ region) by intensive illumination which generates spins at a faster rate than $`1/T_1`$. Here we illuminate the $`p`$ (not $`n`$) region, so we do not have optical spin pumping. The physics is the following: light induces spin polarization of minority carriers (electrons in the $`p`$ region) through optical spin orientation (see Ref. ), which diffuse towards the transition region where they are swiftly pushed by the built-in field into their native $`n`$ region. What we have is spin pumping through the minority channel: spin-polarized minority electrons bring spin into the $`n`$ region faster than the spin there relaxes or diffuses away.

In a sense, minority electrons play the role of circularly polarized light. The evidence for that is in Fig. 1. In the $`n`$ region, spin imbalance is present even though the carrier populations are well relaxed in equilibrium. As a result, spin diffusion is controlled by the majority diffusivity constant $`D_n`$ rather than the minority constant $`D_p`$, as would be the case of carrier diffusion. Spin decay beyond the depletion layer is therefore described by the single exponential $`s\mathrm{exp}(x/L_s)`$, where $`L_s^n=(D_nT_1)^{0.5}1.4`$ $`\mu `$m is the spin diffusion length of electrons in the $`n`$ region. In the $`p`$-region spin diffusion length is $`L_s^p=(D_p\tau _s)^{0.5}`$; electrons diffuse bipolarly with hole diffusivity $`D_p`$ for a time $`\tau _s=\tau T_1/(\tau +T_1)0.067`$ ns, so that $`L_s^p0.26`$ $`\mu `$m is of the order of the electron diffusion length in the $`p`$-region. Since $`L_s^n/L_s^p=(D_n/D_p)^{0.5}\times (1+T_1/\tau )^{0.5}`$ ($`5`$ in our model), the effective range of spin diffusion is extended far beyond the expected value. The depletion layer acts to extend the range of spin diffusion, with the effect most pronounced for $`T_1\tau `$, which is typically the case (especially at lower temperatures).

Our model scheme above constitutes a solar cell, as light illuminates the surface within the electron diffusion length from the depletion layer. As the electron diffusion length almost coincides with $`L_s^p`$, electrons arrive at the transition region spin polarized. The built-in electric field then sweeps the electrons in the $`n`$-region (and holes back into $`p`$). Because of the spin amplification in the $`n`$-region, the resulting electrical current is spin polarized and can be used for spintronic applications (in a combination with ferromagnetic semiconductors or metals). Our spin-polarized solar cell has the usual I-V characteristics, with 1.03 V open-circuit voltage.

How does spin polarization $`\alpha (x)`$ change with the applied bias $`V`$? The calculated profiles for our model device are in Fig. 2 (top). There is a strong monotonic decrease of $`\alpha `$ with increasing $`V`$. To quantify this dependence we consider the change of the total spin accumulated in the cell (integral of the spin density $`s`$ from $`x=0`$ to 2); the spin accumulates almost exclusively in the $`n`$ region (where electron density $`n`$ is large). The result is in Fig. 3. Total spin changes by almost 20 times when increasing $`V`$ from -0.2 V to 1.2 V. By (loose) analogy with junction capacitance, we call this effect spin capacitance. The spin accumulation in the $`p`$ region essentially follows the nonequilibrium density of electrons (charges) there, but in the $`n`$ region, where spin diffusion length is much greater than carrier diffusion length, nonequilibrium spin accumulates to a much greater distance; spin capacitance is not trivially connected with nonequilibrium charges. Clearly, after switching off the light the spin is lost, so to “store” spins in the solar cell one needs to supply energy (the fact that $`T_1\tau `$ will lead to special transient effects related to spin recovery). The reason why $`\alpha `$ depends so strongly on $`V`$ is that $`V`$ changes the extent of the depletion layer: as $`V`$ increases, the width of the depletion layer decreases. Since we are illuminating the same point $`x=0`$ in all cases, the amount of spin polarization that reaches the depletion layer decreases as the width of the layer decreases (so that the distance from the surface to the center of the layer increases). At large forward voltages (say, $`V=1.2`$ V, in Fig. 2), the injected spin essentially decays within $`L_s^p`$, which is shorter than the distance from the illumination surface to the depletion layer. Another effect affecting the dependence of $`\alpha `$ on $`V`$ is the value of the built-in field, which also decreases with increasing $`V`$. This electronic control of spin polarization could be measured by observing luminescence of electrons in the $`n`$ region (by, say, forcing them to recombine with holes in heterostructures attached to the right surface). The above effects are not limited to optical spin injection. The only requirement is that there is an electronic spin and carrier imbalance (which can be also created electronically) in the p region.

Finally, consider the case of spin polarization injected in the majority (here $`n`$) region (say, by ferromagnetic electrodes). There is no light illuminating the sample, and all the parameters and boundary conditions remain as above, except that at $`x=0`$ all the charge densities are the equilibrium ones, and $`\alpha (2)=0.5`$ (a quite favorable case). The calculated $`\alpha (x)`$ for different $`V`$ are in Fig. 2 (bottom). The polarization, as before, is injected through the depletion layer (into minority electrons), but now $`V`$ does not affect $`\alpha (x)`$ that much. Consider what would happen if an $`n`$ region was now attached to the sample from the left, that is, we have an n-p-n transistor. The spin-polarization in the base ($`p`$) would be injected and amplified in the $`n`$-collector (by the spin-pumping from the minority channel we introduced above). We conclude that spin polarization can be injected all the way through a transistor, from emitter into collector.

We thank Paul Crowell for useful discussions. This work was supported by DARPA and the U.S. O.N.R.