text:
# A New Look to Massive Neutron Cores

## Introduction

We review in Section 1 the basic equations of the models being considered as well as the concept of principal transform of a 3-dimensional Riemannian metric which is at the core of our new point of view to understand these models.

Section 2 is devoted to lay down the fundamental system of equations to be integrated. We use spherical space coordinates of the quo-harmonic class which allow to implement $`𝒞^1`$ class smoothness across a sharp border when there is one.

In Section 3 we define the concept of proper mass $`M_p`$ to be used to define the binding energy $`E_b`$ of the models as the difference between the active gravitational mass $`M_a`$ and $`M_p`$. We define also the concept of proper mass density which is a fundamental hybrid concept related to $`M_p`$ and the principal transform of the quotient space metric. Its relation with the pressure $`p`$ characterizes the fluid source independently of the solution being considered.

Section 4 is devoted to establish that the binding energy $`E_b`$ can be obtained as the integral over all space of a localized energy density $`\sigma `$ which depends only on some of the gravitational potentials and its first radial derivatives.

In Section 5 we linearize the fundamental system of equations to obtain the linearized expression of the binding energy $`E_b`$, which coincides with the familiar Newtonian one, thus providing a partial justification to the definition of $`E_b`$ in the non linear regime. We obtain also the linearized expression of the energy density $`\sigma `$.

The last section contains our proposed application of the new point of view to the study of massive neutron cores. The equation of state is the usual one for a degenerate neutron gas, but both the density and the pressure are variants of those used by Oppenheimer and Volkoff . Our main result is that the maximum limit mass that a neutron core can have is approximately 3.8 $`M_{}`$ instead of 0.7 $`M_{}`$ in .

## 1 Static Spherically Symmetric Models

We shall be interested in this paper on global spherically symmetric models, which we shall write using a time adapted coordinate and polar-like space coordinates:

$$ds^2=A^2(r)c^2dt^2+d\widehat{s}^2$$
(1)
where:

$$d\widehat{s}^2=B^2(r)dr^2+B(r)C(r)r^2d\mathrm{\Omega }^2,d\mathrm{\Omega }^2=d\theta ^2+\mathrm{sin}^2\theta d\phi ^2,$$
(2)
solution of Einstein’s field equations:

$$S_{\alpha \beta }=\chi T_{\alpha \beta },\chi =8\pi G/c^4$$
(3)
where the r-h-s describes a compact fluid source, or with fast decreasing density, with two “flavours”: isotropic or a special kind of anisotropic pressure to be presented in a moment.

The quotient 3-dimensional metric (2) can be written using a variety of supplementary coordinate conditions belonging to two different types: algebraic or differential. The most often used is the curvature condition which uses a radial $`\stackrel{~}{r}`$ coordinate with the algebraic condition:

$$\stackrel{~}{C}=\stackrel{~}{B}^1$$
(4)
We shall refer also to any other quantity which assumes the use of this radial coordinate with a tilde overhead.

We shall use here almost exclusively the quo-harmonic condition which restricts the $`r`$ coordinate with the differential condition:

$$C^{}=\frac{2}{r}(BC)$$
(5)
where the prime means derivative with respect to r. The use of a differential condition makes possible the construction of global $`𝒞^1`$ models with a sharp boundary with vacuum, something which is not possible when using the coordinate $`\stackrel{~}{r}`$.

Notice that if $`B(r)`$ and $`C(r)`$ are known the curvature coordinate $`\stackrel{~}{r}`$ is simply the following function of $`r`$:

$$\stackrel{~}{r}=\sqrt{B(r)C(r)}r$$
(6)
while, on the contrary, to obtain explicitly the inverse function is in most cases of interest impossible or very cumbersome.

Despite the emphasis we put on the use of the quo-harmonic coordinate $`r`$ let us be clear from the beginning that the main conclusion of this paper will not owe anything to a particular choice of coordinates. It will owe instead all to a new concept: that of a principal transformation of a 3-dimensional Riemannian metric such as (2).

By definition, in the particular case we are considering, the principal transform of (2) is a new 3-dimensional Riemannian metric 111See for instance references  or :

$$d\overline{s}^2=\mathrm{\Phi }^2(r)B^2(r)dr^2+\mathrm{\Psi }^2(r)B(r)C(r)r^2d\mathrm{\Omega }^2$$
(7)
such that: i)

$$\overline{R}_{jkl}^i=0,i,j,k,l=1,2,3$$
(8)
and: ii)

$$(\widehat{\mathrm{\Gamma }}_{jk}^i\overline{\mathrm{\Gamma }}_{jk}^i)\widehat{g}^{jk}=0,$$
(9)
where, with otherwise obvious notations, the quantities with a hat overhead refer to the metric (2) and the quantities with a bar overhead refer to the metric (7). Notice that both conditions above being tensor conditions under any transformation of space coordinates the concept of Principal transformation is intrinsic to the Killing time congruence we are considering.

One of the practical conveniences of using polar quo-harmonic space coordinates from the outset in the very process of model-building is that, when appropriate boundary conditions are taken into account, the principal transform of (2) is just:

$$d\overline{s}^2=dr^2+r^2d\mathrm{\Omega }^2$$
(10)
and:

$$\mathrm{\Phi }(r)=1/B(r),\mathrm{\Psi }(r)=1/\sqrt{B(r)C(r)}$$
(11)
We shall consider two types of energy-momentum tensors: i) those describing standard perfect fluids with isotropic pressure:

$$T_{00}=\rho A^2,T_{ij}=p\widehat{g}_{ij},(x^0=ct)$$
(12)
and a new type of fluid where the isotropy of the pressure is meant in the sense of the principal transform (7):

$$T_{00}=\rho A^2,T_{ij}=p\overline{g}_{ij}$$
(13)
i.e. in the sense of (10) if polar quo-harmonic coordinates are used.

Notice that this second case, that we shall present in Section 3 as being the truly isotropic fluid, can be considered from the standard point of view as an anisotropic fluid of a particular type. Namely one for which the radial and tangential pressures are related to a single function $`p(r)`$ as follows:

$$p_r=pB^2,p_t=p(BC)^1$$
(14)
We can deal with both cases at once introducing in the field equations a two-valued flag $`f=0`$ or $`f=1`$ and the following expressions for $`p_r`$ and $`p_t`$:

$$p_r=(fB^2+1f)p,p_t=(f(BC)^1+1f)p$$
(15)
thus $`f=0`$ corresponds to the standard case and $`f=1`$ corresponds to the new case.

As usual we shall make use of a compressibility equation:

$$\rho =\rho (p)\text{or}p=p(\rho )$$
(16)
depending on convenience, to describe the physics of the source.

The boundary conditions at infinity will be:

$$\underset{r\mathrm{}}{lim}\mathrm{\Xi }=1,\underset{r\mathrm{}}{lim}r^3\mathrm{\Xi }^{}=0$$
(17)
where $`\mathrm{\Xi }=A,B`$ or $`C`$; and the regularity conditions at the centre of symmetry of the configuration will be:

$$\mathrm{\Xi }_0^{}=0,\mathrm{\Xi }_0^{\prime \prime }<\mathrm{}$$
(18)
Finally, in those cases where the model has a sharp boundary between an interior with $`\rho >0`$ and vacuum we shall require the continuity of both $`Xi`$ and $`X^{}`$, so that the space-time metric will be of global $`𝒞^1`$ class.

The system of units we shall use throughout will be such that:

$$\overline{c}=1,\overline{G}=1$$
(19)
where $`\overline{c}`$ is the speed of light in vacuum and $`\overline{G}`$ is Newton’s constant. A more specific system of units of this class will be chosen in Sect. 6.

## 2 Explicit field equations

Taking into account the coordinate condition (5) and also its derivative:

$$C^{\prime \prime }=2\frac{B^{}}{r}\frac{6}{r^2}(BC)$$
(20)
the field equation $`S_{00}=\rho A^2`$ can be written as follows:

$$B^{\prime \prime }=3\frac{BB^{}}{rC}\frac{B^2}{r^2C}+\frac{5}{4}\frac{B_{}^{}{}_{}{}^{2}}{B}+\frac{B^3}{r^2C^2}8\pi \rho B^3$$
(21)
The field equation $`S_{22}=p_tBCr^2`$ (or $`S_{33}=p_tBCr^2\mathrm{sin}^2\theta `$) can be written using the preceding Eq. (21) and the second Eq. (15) as:

$`A^{\prime \prime }=`$ $`{\displaystyle \frac{1}{rC}}\left({\displaystyle \frac{1}{2}}AB^{}BA^{}\right)+{\displaystyle \frac{1}{2}}{\displaystyle \frac{BA}{r^2C}}\left({\displaystyle \frac{B}{C}}1\right)+{\displaystyle \frac{1}{8}}{\displaystyle \frac{AB_{}^{}{}_{}{}^{2}}{B^2}}+{\displaystyle \frac{1}{2}}{\displaystyle \frac{A^{}B^{}}{B}}`$ (22)
$`+4\pi \rho B^2A+8\pi pAB\left(f/C+B(1f)\right)`$
The remaining equation we have to take care of is $`S_{11}=p_rB^2`$, or taking into account (5):

$$\frac{B}{r^2C}\frac{1}{4}\frac{B_{}^{}{}_{}{}^{2}}{B^2}\frac{B^2}{r^2C^2}\frac{A^{}B^{}}{BA}\frac{B^{}}{rC}2\frac{BA^{}}{rCA}=8\pi p(f+(1f)B^2)$$
(23)
On the other hand from the conservation equation:

$$_\alpha S_1^\alpha =0(x^1=r)$$
(24)
or:

$$S_1^1+\mathrm{\Gamma }_{10}^0(S_1^1S_0^0)+\mathrm{\Gamma }_{12}^2(S_1^1S_2^2)+\mathrm{\Gamma }_{13}^3(S_1^1S_3^3)=0$$
(25)
and Eqs. (15) we derive the equation:

$$p^{}=\frac{A^{}}{A}\left[\rho (fB^2+1f)+p\right]+f\left[2\frac{B^{}}{B}\left(\frac{B^{}}{B}+\frac{2}{r}\frac{B}{C}\right)\left(1\frac{B}{C}\right)\right]p$$
(26)
Taking into account the regularity conditions (18) at the centre of symmetry of the configuration in the preceding equations it is easy to see that they imply:

$$B_0=C_0,p_{t0}=p_{r0},S_{110}=B_0^2p_{t0}$$
(27)
Therefore the remaining field equation $`S_{11}=p_rB^2`$ is satisfied at the origin $`r=0`$ and from from (21), (22) and (25) it follows that it is satisfied everywhere.

What follows is the summary of this section and the preceding one: The models we are considering will be fully described by the field variables $`A,B,C`$ and the source variables $`\rho ,p`$; the latter being related by a compressibility equation of either type (16). This complete set of variables is constrained to satisfy the system of differential equations (5), (21), (22) and (26). Appropriate initial conditions will be $`A_0`$, $`B_0=C_0`$ and $`\rho _0>0`$ (or $`p_0>0`$).The initial values of $`A_0`$ and $`B_0`$ have to be chosen such that the asymptotic conditions (17) are satisfied. The boundary of the source will be defined by the first zero $`r=R`$ of the pressure $`p`$, and beyond the vacuum field equations will be required. The continuity of $`A`$, $`B`$ and $`C`$ and its first derivatives is automatically implemented across the boundary of the source.

## 3 Physical and geometrical interpretations

This will be the more difficult section of this paper, although it does not contain any calculation, because it deals about the meaning of words of common use.

When we look at Eqs. (3) as equations to be solved we all refer to the r-h-s as the source term. But this is not quite correct because the energy-momentum tensor depends on the coefficients of the unknown metric. In the case we are considering in this paper the real source variables are the so called density $`\rho `$ and pressure $`p`$ related by a compressibility equation (16). The meaning of these three ingredients deserve to be examined with some detail. A density by definition is a mass per unit volume, and a pressure is a force per unit surface. Therefore to be clear about them we must tell of what mass are we talking about and to what geometry of space are we referring when using the words volume and surface.

As we all know the concept of mass is tricky because it comes in three flavours: inertial mass, passive gravitational mass and active gravitational mass. Newtonian theory assumes the proportionality of the three masses and General relativity assumes the proportionality of inertial and passive gravitational masses of test bodies. Beyond that we have a few decisions to be taken.

To decide what geometry of space to use is also a tricky problem in relativity theory because we have to decide whether this geometry has to be known before we solve the field equations or will be known only after they have been solved, in which case the meaning of $`\rho `$, $`p`$ and the compressibility Eq. (16) will also be known only after the problem has been solved.

Of the three types of mass, only the meaning of active gravitational mass was settled very early by Tolman  identifying it with the Newtonian mass at infinity and proving that it can be calculated as the following integral over the source:

$$M_a=4\pi _0^{\mathrm{}}\rho \stackrel{~}{r}^2𝑑\stackrel{~}{r}$$
(28)
where we remind that $`\stackrel{~}{r}`$ is the curvature radial coordinate. Using (6) this formula can be written equivalently as:

$$M_a=4\pi _0^{\mathrm{}}\rho F(r)r^2𝑑r$$
(29)
where:

$$F(r)=\frac{1}{2}\sqrt{BC}(B^{}Cr+2B^2)$$
(30)
Bonnor  proposed to eliminate the ambiguities that remain defining the passive gravitational mass as:

$$M(B)=4\pi _0^{\mathrm{}}(\rho +p)\stackrel{~}{B}\stackrel{~}{r}^2𝑑\stackrel{~}{r}$$
(31)
This means that $`\rho +p`$ is interpreted as a density of inertial, or passive, gravitational mass and that the metric that gives a meaning to the word volume is (2). But this metric is known only once the problem has been solved and then it depends on the point of the body which is considered. This deprives the meaning of the variables $`\rho `$ and $`p`$ and the compressibility equation (16) of any a priori significance. Other difficulties that arise from this definition are discussed in Bonnor’s paper.

We do not believe that the consideration of a single body at rest, as we have been doing here, can say anything about its passive or inertial gravitational mass because this would require to know how it reacts to the presence of another comparable body. On the other hand we believe that we should be able to define its proper mass $`M_p`$ if we want to know what is the binding energy $`E_b`$ of any given configuration as defined by:

$$E_b=M_aM_p$$
(32)
More precisely, our point of view, along the lines of a long enduring effort to understand the concept of rigidity and establishing a theory of frames of reference in special and general relativity, consists in accepting the usual generalization of Schwarzschild’s “substantial mass” 222An english translation of Schwarzschild’s paper has been provided by S. Antoci in arXiv:physics/9912033 as proper mass:

$$M_p=4\pi _0^{\mathrm{}}\rho B^2Cr^2𝑑r,\sqrt{\widehat{g}}=B^2Cr^2\mathrm{sin}\theta $$
(33)
and defining at the same time a proper mass density $`\rho _p`$:

$$\rho _p=\rho B^2C$$
(34)
such that $`M_p`$ could be written:

$$M_p=4\pi _0^{\mathrm{}}\rho _pr^2𝑑r,\sqrt{\overline{g}}=r^2\mathrm{sin}\theta $$
(35)
This means then defining $`\rho _p`$ as a density of proper gravitational mass and interpreting the words volume and surface in the sense of the universal euclidian geometry (10) related to the quotient metric (2) by a principal transformation (7). This guarantees the independence of the meanings of $`\rho _p`$, $`p`$ and the compressibility equation (16) independently of the location of the element of the fluid in the object and independently of the solution of the field equations that one is considering. This guarantees also that $`M_p`$ can be identified with an appropriate number of identical samples of a fluid as weighted with a balance at the “shop store” before being assembled into the body.

The interpretation we have just given of the metric (10) implies also as a corollary that, as above-mentioned, a fluid with isotropic pressure should be described by an energy-momentum tensor as written in (13). And that for the compressibility equation to have a well defined a priori meaning it should be given as a relationship between $`\rho _p`$ and $`p`$:

$$\rho _p=\rho _p(p)\text{or}p=p(\rho _p)$$
(36)
## 4 Localized energy density

Following suit to the ideas of the preceding section we exhibit the quantity $`E_b`$ as an integral extended over all space of an energy density function depending only on $`r,B,C`$ and $`B^{},C^{}`$. From (29) and (35) it follows that:

$$E_b=4\pi _0^{\mathrm{}}\rho B^3Yr^2𝑑r$$
(37)
where:

$$Y=\frac{F}{B^3}\frac{C}{B}$$
(38)
while $`\rho B^3`$ can be obtained from (21) as:

$$\rho B^3=\frac{1}{8\pi }(XB^{\prime \prime })$$
(39)
with:

$$X=\frac{B(3rB^{}B+B^2/C)}{r^2C}+\frac{5}{4}\frac{B_{}^{}{}_{}{}^{2}}{B}$$
(40)
Therefore we have:

$$E_b=\frac{1}{2}_0^{\mathrm{}}XYr^2𝑑r\frac{1}{2}_0^{\mathrm{}}B^{\prime \prime }Yr^2𝑑r$$
(41)
Integrating by parts the second integral we obtain:

$$E_b=4\pi _0^{\mathrm{}}\sigma (r)r^2𝑑r4\pi \underset{r\mathrm{}}{lim}\sqrt{\frac{C}{B}}B^{}r^2\left(\frac{1}{4}\frac{B^{}Cr}{B^2}+1\sqrt{\frac{C}{B}}\right)$$
(42)
where:

$$\sigma (r)=\frac{1}{8\pi }(XY+Z)$$
(43)
with:

$$Z=\frac{B^{}}{8rB^3\sqrt{BC}}(5r^2C^2B_{}^{}{}_{}{}^{2}+2rBC(4\sqrt{BC}+C)B^{}8B^3(2\sqrt{BC}B+C))$$
(44)
depends only on $`r,B,C`$ and $`B^{}`$. From the asymptotic conditions (17) it follows that :

$$\underset{r\mathrm{}}{lim}B/C=1\text{and}\underset{r\mathrm{}}{lim}r^3B_{}^{}{}_{}{}^{2}=0$$
(45)
and therefore, the limit in (42) being zero, the final result is:

$$E_b=4\pi _0^{\mathrm{}}\sigma (r)r^2𝑑r$$
(46)
## 5 Linear approximation

We consider here the linear approximation of the models that we have been considering, to take a closer look to two of the de concepts that we have implemented in Sect. 3. Namely: the mass defect, or binding energy $`E_b`$ and the proper mass density $`\rho _p`$.

We assume that $`A`$, $`B`$ and $`C`$ can be written as:

$$A=1+A_1,B=1+B_1,C=1+C_1$$
(47)
where $`A_1`$, $`B_1`$ and $`C_1`$ are small quantities, of order $`ϵ`$ say. We assume also that $`\rho `$ is also of order $`ϵ`$ and that $`p`$ is of order $`ϵ^2`$ and can be ignored, as well as any other quantity of the same order or smaller, in the field equations.

The coordinate condition (5) and the field equations (21), (22) become then:

$`C_1^{}`$ $`=`$ $`{\displaystyle \frac{2}{r}}(B_1C_1)`$ (48)
$`B_1^{\prime \prime }`$ $`=`$ $`3{\displaystyle \frac{B_1^{}}{r}}+{\displaystyle \frac{1}{r^2}}(B_1C_1)8\pi \rho `$ (49)
$`A_1^{\prime \prime }`$ $`=`$ $`{\displaystyle \frac{1}{2r}}(B_1^{}2A_1^{})+{\displaystyle \frac{1}{2r^2}}(B_1C_1)+4\pi \rho `$ (50)
As our purpose is purely illustrative here we consider below the simplest case where the source is a spherical body of finite radius $`R`$ and constant $`\rho `$. The interior solution satisfying the regularity conditions (18) at the centre is:

$$A_1=\frac{2}{3}\pi \rho r^2+a_0,B_1=\frac{16}{15}\pi \rho r^2+b_0C_1=\frac{8}{15}\pi \rho r^2+b_0$$
(51)
where $`a_0`$ and $`b_0`$ are two allowed constants of integration; the exterior solution satisfying the asymptotic conditions (17) at infinity is:

$$A_1=\frac{b_1}{r},B_1=\frac{b_1}{r}+\frac{b_3}{r^3},C_1=\frac{2b_1}{r}\frac{2b_3}{r^3}$$
(52)
where $`a_1`$ and $`b_3`$ are two new constants of integration. Demanding the continuity of $`A_1`$, $`B_1`$ and $`C_1`$ and their first derivatives across the border $`r=R`$ fixes $`a_0`$ and all the $`b`$’s as follows:

$$a_0=2\pi \rho R^2,b_0=\frac{8}{3}\pi \rho R^2,b_1=\frac{4}{3}\pi \rho R^3,b_3=\frac{4}{15}\rho R^5$$
(53)
From the preceding results we can calculate the leading approximation, which is of order $`ϵ^2`$, of the localized energy density (43). For $`r<R`$ the result is using an arbitrary system of units:

$$\sigma =\frac{3}{20}\frac{G^2M^2r^2}{\pi R^6}$$
(54)
where at this approximation $`M`$ is either $`M_a`$ or $`M_p`$. For $`r>R`$ the result is:

$$\sigma =\frac{3}{160}\frac{G^2M^2(5r^4+6R^2r^23R^4)}{\pi r^8}$$
(55)
The binding energy can be calculated using (46), (54) and (55), or (32), (29) and (35) at the appropriate approximation. The result is, using arbitrary units, the familiar Newtonian amount:

$$E_b=\frac{3}{5}\frac{GM^2}{c^2R},$$
(56)
a result that can be obtained using a variety of other approaches 333See for instance reference

## 6 Massive Neutron Cores

Any particular model will be characterized by an equation of state and the value of its central density, or central pressure, or both in the important case in which one assumes that the density is constant. Taking into account the regularity conditions the initial conditions of the gravitational potentials $`A,B`$ and $`C`$ have to be chosen such that:

$$B_0=C_0,A_0^{}=B_0^{}=C_0^{}.$$
(57)
And taking into account the asymptotic conditions the values of $`A_0`$ and $`B_0`$ have to be chosen such that:

$$\underset{r\mathrm{}}{lim}A=\underset{r\mathrm{}}{lim}B=1$$
(58)
the condition:

$$\underset{r\mathrm{}}{lim}C=1$$
(59)
as, well as the remaining asymptotic conditions, being then automatically satisfied because the solution behaves as the exterior Schwarzschild one at infinity.

The numerical integration of the system of equations (5), (21), (22) and (26) where:

$$\rho =\rho _p(f(B^2C)^1+1f)$$
(60)
and $`\rho _p`$ is a known function of $`p`$, is a trial and error procedure. Arbitrary values of $`A_0,B_0`$ and $`p_0>0`$ have to be chosen; the integration has to proceed until $`p=0`$; then the equation of state has to be abandoned and $`\rho =p=0`$ has to be required; the integration has then to proceed to sufficiently large values of $`r`$ to check the asymptotic conditions (58). If the check is not satisfactory the whole process has to be started again with new values of $`A_0,B_0`$ and $`p_0>0`$.

As an important example we consider the equation of state of a degenerate neutron gas as it suits to a model of massive neutron cores. This was considered in a famous paper by Oppenheimer and Volkoff from the standard point of view which consists in putting $`f=0`$. Our point of view consist in using instead the value $`f=1`$.

The equation of state can be written in parametric form, including both points of view, as:

$`\rho `$ $`=`$ $`K(\mathrm{sinh}uu)(f(B^2C)^1+1f)`$ (61)
$`p`$ $`=`$ $`{\displaystyle \frac{1}{3}}K(\mathrm{sinh}u8\mathrm{sinh}{\displaystyle \frac{1}{2}}u+3u)`$ (62)
where,using arbitrary units:

$$K=\frac{m^4c^5}{4h^3}$$
(63)
$`m`$ being the mass of a neutron.

We recall below the values of $`M_a`$ obtained in  for several values of the inial values of $`u_0`$, and include the values of $`M_p`$ as calculated from (33):

| $`u_0`$ | $`M_a`$ | $`M_p`$ | $`E_b`$ |
| --- | --- | --- | --- |
| 1 | 0.033 | 0.033 | -0.000 |
| 2 | 0.066 | 0.071 | -0.005 |
| 3 | 0.077 | 0.088 | -0.011 |
| 4 | 0.070 | 0.085 | -0.015 |
| 5 | 0.060 | 0.074 | -0.014 |
| 6 | 0.049 | 0.062 | -0.013 |
| 7 | 0.042 | 0.054 | -0.012 |

The system of units that has been used is that satisfying (19) completed with the supplementary condition:

$$\overline{K}=1/4\pi $$
(64)
which is the choice made in  to define a unit of mass. The most notorious result is existence of a maximum mass corresponding approximately to $`u_0=3`$ whose value is $`M_a=0.078`$ which corresponds to $`M_a=0.71M_{}`$. Oppenheimer and Volkoff concluded also from a very crude non relativistic argument that above $`u_0=3`$ the equilibrium configurations were not stable.

For $`f=1`$ the results that we have obtained are:

| $`u_0`$ | $`M_a`$ | $`M_p`$ | $`E_b`$ |
| --- | --- | --- | --- |
| 1 | 0.04 | 0.04 | -0.00 |
| 2 | 0.12 | 0.13 | -0.01 |
| 3 | 0.24 | 0.28 | -0.04 |
| 4 | 0.36 | 0.47 | -0.11 |
| 5 | 0.41 | 0.56 | -0.15 |
| 6 | 0.36 | 0.51 | -0.15 |
| 7 | 0.33 | 0.48 | -0.15 |

Here also we obtain that there is a maximum mass $`M_a=0.41`$ which corresponds to $`M_a=3.8M_{}`$ but it is substantially larger than the value in  as well as the mass of some models with anisotropic pressure considered by Corchero in  . It is even somewhat larger than the limit value, $`M_a=3.2M_{}`$ obtained by Rhoades and Ruffini in  from very general considerations complying with the conventional point of view.

## Acknowledgments

I gratefully acknowledge the help provided by J. M. Aguirregabiria checking parts of this manuscript, the patience of A. Chamorro for a careful reading of it, and their comments as well as those of E. S. Corchero and J. Martín. I also gratefully acknowledge the position of visiting professor to the UPV/EHU that I have been holding while this paper was being prepared.


format: text

text:
# Potential resonant screening effects on stellar 12C+12C reaction rates

## Abstract

The 12C+12C fusion cross sections show resonant behavior down to the lowest energies accessible so far in the laboratory. If this tendency continues into the astrophysical energy range, the stellar 12C+12C reaction rates have to be corrected for resonant screening effects, in addition to the conventional screening corrections. We estimate the resonant screening effects in the weak electron screening limit for hydrostatic burning and white dwarf environments.

Nuclear reactions are the energy source for many astrophysical objects . When applied in the stellar environment, however, the nuclear reaction rate, determined from measurements for bare nuclei (for electron screening corrections to low-energy laboratory cross sections, see ), has to be corrected for plasma effects. The medium induces a screening potential $`V_{sc}(r)`$ which effectively reduces the Coulomb barrier between the two colliding nuclei and enhances the reaction cross section. This effective increase of the in-medium cross section over the laboratory cross section is conventionally described by an enhancement factor $`f_s`$. If $`V_{sc}(r)`$ is approximately constant during the penetration process, i.e. for radii smaller than the outer barrier turning point, the potential can be replaced by the screening energy $`U_0=V_{sc}(0)`$. Such a situation occurs in astrophysical sites in which the average Coulomb energy between neighboring ions is much less than the thermal energy of the plasma. As has been shown by Salpeter , in this weak screening limit Debye-Hückel theory is applicable leading to an enhancement of the thermally averaged nuclear reaction rate in the plasma. In cold ($`T0`$), but dense plasma fusion reactions are induced by density fluctuations. Again, the medium-induced screening potential strongly enhances the reaction cross sections. The appropriate formalism is developed for example in .

Besides this global screening enhancement an additional plasma correction to the nuclear reaction rate can occur if the rate is dominated by narrow resonances. The general theory has been outlined in , but no application has been identified so far. We will point out in this Note that this situation may occur for the 12C+12C fusion reaction, where the cross sections show noticeable resonant structures down to the lowest energies measured so far in the laboratory ($`E2.4`$ MeV, ). Until now, the resonant structure of the 12C+12C fusion cross section has been ignored in astrophysical applications and the respective reaction rate has been determined as a smooth average over the resonant contributions. Consistent with this procedure only the global plasma screening enhancement has been considered. However, if the resonant structure continues to even lower energies and the astrophysical reaction rate is due to the contributions of narrow resonances, one then has to consider that the entrance channel width of these resonances will be modified in the plasma. This resonant screening effect will reduce the reaction rate. Our discussion is based on the fact that for 12C+12C resonances far below the height of the Coulomb barrier, the entrance channel width ($`\mathrm{\Gamma }_a`$) is much smaller than the total resonance width. The latter (which is of order $`100`$ keV for the observed resonances above 2.4 MeV) is also noticeably smaller than the resonance energy. Hence we can treat the contribution of these resonances to the thermally averaged reaction rate by the formalism developed for narrow, isolated resonances (e.g. ), which for the 12C+12C system is:

$$\sigma v_R=\left(\frac{2\pi }{\mu kT}\right)^{3/2}\mathrm{}^22(2J+1)\mathrm{\Gamma }_a(E_r)\mathrm{exp}\left(\frac{E_r}{kT}\right)$$
(1)
where $`\mu `$ is the reduced mass and $`E_r`$ and $`J`$ are the energy and total angular momentum of the resonance.

In the medium the Coulomb potential $`V_c(r)`$ has to be replaced by $`V_c(r)+V_{sc}(r)`$ ($`V_{sc}(r)<0`$). This modification will change the resonance energy in the plasma to $`E_r^{}<E_r`$. Relatedly the resonance width in the plasma $`\mathrm{\Gamma }^{}(E_r^{})`$ is proportional to the penetration factor

$$\mathrm{\Gamma }^{}(E_r^{})exp\left(2\sqrt{\frac{2\mu }{\mathrm{}^2}}_{r_{in}}^{r_{out}}\sqrt{(V_c(r)+V_{sc}(r)E_r^{})}𝑑r\right)$$
(2)
where $`r_{in}`$ and $`r_{out}`$ are the inner and outer barrier turning points, respectively, and $`\mu `$ the reduced mass of the two colliding nuclei. The situation simplifies, if one has $`V_{sc}(r)=V_{sc}(0)=U_0`$ for radii $`r<r_{in}`$. Then, in a good approximation one has $`E_r^{}=E_rU_0`$. In the weak screening limit one even has $`V_{sc}(r)V_{sc}(0)`$ for $`r<r_{out}`$, resulting in $`\mathrm{\Gamma }^{}(E_r^{})=\mathrm{\Gamma }(E_r)`$ Such a situation holds in a massive star ($`20M_{}`$) where carbon core burning occurs at temperatures around $`10^9`$ K and for densities of order $`\rho =210^5`$ g/cm3 . Thus no resonant screening correction for the 12C+12C reaction cross section is expected in hydrostatic carbon burning.

In passing we mention that many nuclear reactions occuring in explosive hydrogen burning  are dominated by narrow resonances, e.g. the $`{}_{}{}^{15}O(\alpha ,\gamma )^{19}Ne`$ reaction which is usually believed to be a key for the break-out of the matter-flow from the hot CNO-cycle in novae or x-ray bursters. However, in both sites the weak screening conditions approximately hold and no resonant screening corrections have to be applied.

In general the magnitude of $`V_{sc}(r)`$ decreases with increasing radius; at large distances the medium-induced screening potential cancels the Coulomb repulsion between the colliding nuclei, i.e. $`V_c(r)+V_{sc}(r)=0`$ for large $`r`$. If $`V_{sc}(r)`$ already decreases for radii $`r`$ under the barrier ($`r_{in}<r<r_{out}`$) the width of the barrier to be penetrated in the medium is longer than the barrier width for bare nuclei. Correspondingly one has for the resonance width $`\mathrm{\Gamma }^{}(E_r^{})<\mathrm{\Gamma }(E_r)`$ which leads to an additional screening correction of the resonant cross section which we like to denote by $`f_r=\frac{\mathrm{\Gamma }^{}(E_r^{})}{\mathrm{\Gamma }(E_r)}`$. The total screening correction is then given by $`f_{tot}=f_sf_r`$. Obviously the resonant screening correction reduces the cross section in the medium, counteracting the overall enhancement by $`f_s`$.

For reactions in cold, dense plasma Salpeter and van Horn have derived an appropriate analytical expression for the resonant screening correction. This description is appropriate for a carbon white dwarf environment with $`T=510^7`$ K and $`\rho =210^9`$ K . Following  one then finds a resonant screening correction factor for the 12C+12C reaction if this reaction indeed proceeds via a narrow resonance at energy $`E_r`$:

$$f_r=exp\left(\frac{\pi }{\sqrt{\lambda ϵ}}\left[\frac{1.22}{ϵ^3}\frac{3.1}{ϵ^6}+\frac{75}{ϵ^9}\right]\right)$$
(3)
with the abbreviations

$$\lambda =\frac{1}{432}\left(\frac{\rho _{12}}{1.629}\right)^{1/3},ϵ=\frac{1}{36}\left(\frac{1.629}{\rho _{12}}\right)^{1/3}\frac{E_r}{49.6}$$
(4)
where $`\rho _{12}`$ measures the density in $`10^{12}`$ g/cm3 and the resonance energy $`E_r`$ is defined in keV. Obviously the correction becomes more important with increasing density and at lower resonance energies.

Although the 12C+12C resonance energies below $`E2.4`$ MeV are currently not known, we will now assume that the resonance structures observed in the fusion data at $`E2.4`$ MeV continue to lower energies. While this will probably not strongly affect the smooth energy dependence of the averaged astrophysical S-factor derived from the data at higher energies and used for the reaction rate in astrophysical applications so far, it will, however, change the screening corrections to the rate in the astrophysical environment. We will estimate the reduction of the 12C+12C reaction rate due to the resonant screening correction factor $`f_r`$ for carbon white dwarf environment hypothetically assuming a resonance in the energy interval $`0.42`$ MeV.

As is shown in the Fig. 1. the resonant screening correction can lead to a significant reduction of the reaction rate, which amounts to more than 11 orders of magnitude if the resonance energy is as low as $`E_r=400`$ keV. Such a change might influence the carbon ignition density in white dwarfs. We note, however, that despite such a potentially large correction global screening still dominates, which due to Table 1 of  amounts to an enhancement of the rate by 12 orders of magnitude. But this estimate might also change if the short-ranged potential which gives rise to the resonances is taken into account .

In summary, the 12C+12C reaction rate in a white dwarf environment might change significantly, if the reaction is dominated by narrow resonances at the relevant energies. This can only be established experimentally. If such a resonance is indeed being observed it is desirable to study the plasma enhancement factor due to many-nuclear correlations including a nuclear potential which supports resonances at low 12C+12C scattering energies.

Rob Cussons was partly supported by a Leeds-Aarhus student exchange program. Theodore Liolios acknowledges a Marie-Curie fellowship from the European Commission for a stay at the University of Aarhus. Our work is supported in parts by the Danish Research Council.


text:
# Dark energy and dissipation

## 1 General setting

The Friedmann equation plus the conservation equations for matter and dark energy in a Friedmann–Robertson–Walker universe dominated by these two components (non–interacting with one another), in terms of the density parameters, read

$`1`$ $`=`$ $`\mathrm{\Omega }_m+\mathrm{\Omega }_\varphi +\mathrm{\Omega }_k,`$ (1)
$`\dot{\mathrm{\Omega }}`$ $`=`$ $`(3\gamma 2)H(\mathrm{\Omega }1)\mathrm{\Omega },`$ (2)
$`\dot{\mathrm{\Omega }}_\varphi `$ $`=`$ $`[2+(3\gamma 2)\mathrm{\Omega }3\gamma _\varphi ]\mathrm{\Omega }_\varphi H,`$ (3)
where $`\mathrm{\Omega }\mathrm{\Omega }_m+\mathrm{\Omega }_\varphi `$, and $`\gamma `$ stands for the overall adiabatic index $`\gamma =(\gamma _m\mathrm{\Omega }_m+\gamma _\varphi \mathrm{\Omega }_\varphi )/\mathrm{\Omega }`$, with $`\gamma _{m,\varphi }1+(p_{m,\varphi }/\rho _{m,\varphi })`$, and such that $`1\gamma _m2`$ and $`0\gamma _\varphi <1`$ (it should be noted that in general $`\gamma _m`$ and $`\gamma _\varphi `$ may vary with time).

It is not our aim to propose a new dark energy model but rather to show that if one wishes to have late cosmic acceleration and simultaneosly solve the coincidence problem using general relativity, then the matter fluid component must be dissipative -in the sense explained below- irrespective of the potential of the scalar field .

From the above equations it is immediately seen that for $`\mathrm{\Omega }=1`$ Eq. (3) implies that $`\dot{\mathrm{\Omega }}_\varphi >0`$. Consequently, at large times $`\mathrm{\Omega }_\varphi 1`$ and $`\mathrm{\Omega }_m0`$, i.e., the accelerated expansion ($`q<0\gamma _\varphi <2/3`$) and the coincidence problem cannot be solved simultaneously within this approach. Moreover, for the solution $`\mathrm{\Omega }=1`$ to be stable the overall adiabatic index must comply with the upper bound $`\gamma <2/3`$ which is uncomfortably low.

A way out is to assume the matter fluid dissipative, i.e., endowed with a nonequilibrium scalar pressure $`\pi `$. This is always negative for expanding universes as required by the second law of thermodynamics , , . Further, the occurrence of $`\pi `$ is rather natural since (exception made of superfluids) all matter fluids found in Nature are dissipative -see e.g. , . Thereby Eqs. (2) and (3) generalize to

$`\dot{\mathrm{\Omega }}`$ $`=`$ $`\left[3\left(\gamma +{\displaystyle \frac{\pi }{\rho }}\right)2\right]H(\mathrm{\Omega }1)\mathrm{\Omega },`$ (4)
$`\dot{\mathrm{\Omega }}_\varphi `$ $`=`$ $`\left\{2+\left[3\left(\gamma +{\displaystyle \frac{\pi }{\rho }}\right)2\right]\mathrm{\Omega }3\gamma _\varphi \right\}H\mathrm{\Omega }_\varphi .`$ (5)
We now may have $`\dot{\mathrm{\Omega }}_\varphi <0`$ as well as $`\mathrm{\Omega }_m\mathrm{\Omega }_{m0}0`$ and $`\mathrm{\Omega }_\varphi \mathrm{\Omega }_{\varphi 0}0`$ for late time so long as the stationary condition

$$\gamma _m+\frac{\pi }{\rho _m}=\gamma _\varphi =2\frac{\dot{H}}{3H^2}$$
(6)
is satisfied. Besides, the constraint $`\gamma <2/3`$ is replaced by $`\gamma +(\pi /\rho )<2/3`$, which is somewhat easier to fulfill.

## 2 Stability analysis

For spatially flat FRW universes the asymptotic stability of the stationary solution $`\mathrm{\Omega }_{m0}`$ and $`\mathrm{\Omega }_{\varphi 0}`$ can be studied from equation (5). By slightly perturbing $`\mathrm{\Omega }_\varphi `$ it follows that the solution is stable (and therefore an attractor) provided the quantity $`\gamma _m+\frac{\pi }{\rho _m}\gamma _\varphi <0`$ and tends to zero as $`t\mathrm{}`$. This coincides with the stationary condition (6).

For $`\mathrm{\Omega }1`$ (i.e., when $`k0`$), it is expedient to introduce the ansatz $`ϵ=ϵ_0+\delta `$ in Eqs (4) & (5), where $`ϵ_0(\mathrm{\Omega }_m/\mathrm{\Omega }_\varphi )_0𝒪(1)`$ and $`\delta ϵ_0`$. One finds that

$$\dot{\delta }=\frac{3}{\mathrm{\Omega }_\varphi }\left(\frac{2}{3}\gamma _\varphi \right)\mathrm{\Omega }_kH(ϵ_0+\delta ).$$
(7)
As a consequence, the stationary solution will be stable for open FRW universes ($`\mathrm{\Omega }_k>0`$). For closed FRW universes one has to go beyond the linear analysis.

A realization of these ideas can be found in . There it is seen that the space parameter is ample enough that no fine tuning is required to have late acceleration together with the fact that both density parameters tend to constant values compatible with observation.

Moreover we have shown that, for a wide class of dissipative dark energy models, attractor solutions are themselves attracted towards a common asymptotic behavior. This “superattractor” regime provides a model of the recent universe that also exhibits an excellent fit to supernovae luminosity observations and no age conflict .

## 3 Concluding remarks

Although dissipative pressures have been repeteadly invoked when building cosmological models -see e.g. , , , , \- one may ask about its origin in the present context. On the one hand if the CDM particles are really self–interacting (as it has been recently argued ), then they are bound to arise. (CDM self–interaction with mean free path between about $`1`$ kpc and $`1`$ Mpc seems to be key to explain the structure of the halo of galaxies). On the other hand, it is well known that particle production somehow plays the role of a dissipative pressure . The latter is linked to the rate $`\mathrm{\Gamma }`$ of particle production by $`\pi =\gamma _m\rho _m\mathrm{\Gamma }/(3H)`$. Therefore, even a very small rate can produce a significant dissipative pressure .

According to the standard picture of cosmic structure formation a sufficiently long matter–dominated era must have taken place during which the observed structure grew from the density fluctuations measured by CMB anisotropy experiments -see e.g. . As a consequence the transition to the currently observed accelerated regime should have been quite recent in the evolution of the universe. This seems to imply that dissipative effects in dark matter were small until density inhomogeneities became large. If so, it is suggestive to think that the size of dissipative processes, as measured by the ratio $`\pi /\rho _m`$, increased with dark matter density and became large precisely because of the development of inhomogeneities . The observed smoothness of halo structure and the correlation between high energy cosmic ray production and dark matter clustering might provide support for this relationship .

Before closing, we would like to emphasize that our analysis does not apply to scenarios in which the dark energy component interacts with the matter component -see e.g. , . This is easy to understand as the interaction somehow mimics a dissipative pressure.

## Acknowledgments

This work has been partially supported by the Spanish Ministry of Science and Technology under grant BFM 2000-C-03-01 and 2000-1322, and the University of Buenos Aires under Project X223.


text:
# The Y2 Stellar Evolutionary Tracks

## 1 Introduction

The $`Y^2`$ group has recently released their updated isochrones (Yi et al., 2001; Kim et al., 2002). Theses isochrones have been tested against many high-quality observational data and proved reasonably accurate. Besides such successes, their wide coverage in metallicity and age, which is useful to population synthesis studies, triggered many requests for the corresponding stellar evolutionary tracks. Stellar tracks are the basic building blocks of isochrones, and thus it is important to ascertain their accuracy by all means before we feel confident about isochrones. A compelling test of such kinds is made on binary data which can yield mass and radius information on individual stars independently of the infamous uncertainty in distance.

## 2 Stellar Evolution Models

The critical pieces of input physics and model parameters are listed in Table 1, while greater details are available from Yi et al. (2001) and from Kim et al. (2002). The models were evolved from the pre-main-sequence (MS) stellar birthline to the onset of helium burning in the core at the red giant branch tip (RGBT). The mass range is approximately 0.4 $`M_{}`$ through 5 $`M_{}`$, which has been chosen for constructing full isochrones from MS to RGBT for ages 0.1 – 20 Gyr.

We release the entire database of the $`Y^2`$ stellar models in the form shown in Tables 2 and 3. Each stellar track is stored in two files. Track1 files hold information simply for 24 select young ages (1 – 80 Myr). In case of relatively low-mass models, this file represents a pre-MS track. Track2 files are for older ages extended to RGBT. A total of 150 grids, which some might want to call equal evolutionary points, are chosen according to a highly complicated scheme. The first 50 grids are based on the central helium abundance ($`Y_c`$) which grows monotonically and gradually as soon as central hydrogen burning begins. The next 100 grids are based on the core mass ($`M_c`$) that starts growing about the time $`Y_c`$ stops increasing after reaching its maximum. This transition may be regarded as the terminal-age main sequence. The grids in $`Y_c`$ and $`M_c`$ are not fixed but variable in different tracks because their maximum values and the ways they reach their maximum values are complex functions of mass, metallicity and so on. Therefore, in addition to the usual track information (age, temperature, luminosity), we present $`Y_c`$ and $`M_c`$ as well.

Figure 1 shows a sample set of stellar models. For clarity, pre-MS parts appear separately in dotted lines. For this separation, we arbitrarily chose the point where $`Y_c`$ is 0.2% enhanced from the initial value, which means that the star has just begun burning atomic hydrogen in its core. Note that the original tracks, as shown in this figure, have orders of magnitude finer evolutionary grids and thus are slightly more accurate but far less portable than the tracks ditributed through this paper. If high accuracy is required for good scientific reasons, the original-grid tracks can be made available for a set of original parameters ($`Z`$, $`M`$, \[$`\alpha `$/Fe\], overshoot) upon request to the authors.

## 3 Database and Interpolation Code

The full set of stellar models and a FORTRAN package that work for mass, metallicity, and $`\alpha `$-enhancement interpolation are available from the authors upon request or directly from our Web sites listed below, as well as from the electronic edition of the Astrophysical Journal.

* UK: www-astro.physics.ox.ac.uk/$``$yi/yystar.html
* US: www.astro.yale.edu/demarque/yystar.html
* Korea: csaweb.yonsei.ac.kr/$``$kim/yystar.html

The interpolation code adopts a polynomial interpolation scheme. We encourage the readers to use it when they find our model grids are not sufficiently small for their purposes. However, it should be emphasized that the interpolated tracks do not reproduce exactly the genuine tracks constructed using the stellar evolution code, as some properties evolve in much more complicated manners within the given grid space than a simple interpolation can mimic. The interpolation near the critical mass ($`M_{crit}^{conv}`$), above which a stellar core becomes convective on the MS (Yi et al., 2001), is a good example. Theory suggests a transition from a radiative core to a convective core as stellar mass reaches $`M_{crit}^{conv}`$, but determining $`M_{crit}^{conv}`$ is not trivial because the development of a convective core is gradual with respect to mass. In this regard, the accuracy of our critical masses listed in Table 2 of Yi et al. (2001) is no better than our typical mass grid, $`0.1M_{}`$. As a result, interpolated tracks based on the given models are unavoidably subject to such internal uncertainties. Besides, the abrupt change in the overshoot parameter near $`M_{crit}^{conv}`$ (overshoor is numerically turned on at $`M>M_{crit}^{conv}`$) introduces, perhaps inaccurately, discontinuity in the shape of the main-sequence turn-off region as a function of stellar mass, which makes the interpolations near $`M_{crit}^{conv}`$ difficult. Keeping these inherent limitations in mind, we purposely kept the interpolation scheme as simple as possible. We suggest that the readers should use the models and code with proper care.

This research has been supported by Korean Research Foundation Grant KRF-2002-070-C00045 (YCK), and received partial support from NASA grant NAG5-8406 (PD).


text:
# Temporal variation of coupling constants and nucleosynthesis

## 1 Introduction

In previous works, different constraints have been obtained for possible temporal variations of fundamental constants. Recent investigations involved primordial nucleosynthesis , cosmic microwave background (CMB) , quasar absorption lines , stellar nucleosynthesis , meteorites , the Oklo phenomenon  and atomic clocks . In most of these works temporal variations of the fine structure constant were investigated. Presently the only analysis showing a time variation definitely different from zero stems from the analysis of atomic multiplet spectra in quasar absorption lines through intervening interstellar clouds at $`0.5<z<3.5`$ .

In this work we investigate the triple-alpha process and the Oklo phenomenon to obtain cosmological constraints on the time variation of the electromagnetic fine structure constant, and the strength of the strong interaction, i.e., the nucleon and meson masses or the QCD scale parameter $`\mathrm{\Lambda }_{\mathrm{QCD}}`$.

## 2 Triple–Alpha process

In stellar nucleosynthesis one of the most sensitive reactions with respect to possible variations of the coupling constants is the triple–alpha process leading to the production of 12C . Stellar nucleosynthesis of carbon and oxygen in helium burning of ancient stars is therefore particularly of interest for two reasons: The process is very sensitive to variations in fundamental parameters and also took place at the same time ($`t_\mathrm{B}513`$ Gyr or $`z0.53`$) that absorption lines were being created in the light from quasars due to intervening gas clouds . The C/O abundance ratio and the absorption lines in the quasar spectra are therefore both sensitive to variations in coupling constants at the same epoch.

We compare the change of abundance ratios of C/O by variation of fundamental parameters like the fine structure constant $`\alpha _{\mathrm{em}}`$, and/or the nucleon $`m_\mathrm{N}`$ and meson masses $`m_\mathrm{B}`$. In the first step we investigated the changes of the 0$`{}_{}{}^{+}{}_{2}{}^{}`$–resonance in 12C using nuclear microscopic model calculations  (Table 1).

The observation of the C/O–abundance ratios in ancient stars and their computation in stellar models do not seem to vary by more than a factor 3 (see e.g., ). Such a change would be produced in helium burning of Red Giants by a variation of the fine structure constant of not more than $`\pm 0.6\%`$. This leads to the constraint $`\mathrm{\Delta }\alpha _{\mathrm{em}}/\alpha _{\mathrm{em}}\pm 6\times 10^3`$ shown in Table 2. The time variation of the fine structure constant through the analysis of atomic multiplet spectra in quasar absorption lines through intervening interstellar clouds of the order of $`\mathrm{\Delta }\alpha _{\mathrm{em}}`$/$`\alpha _{\mathrm{em}}10^5`$ is still much weaker than the constraint obtained from the triple–alpha process of the order of $`10^310^2`$ (see Table 2).

Almost the same value for the constraint is obtained when assuming identical linear scaling of the nucleon $`m_\mathrm{N}`$ and meson masses $`m_\mathrm{B}`$, i.e., $`\mathrm{\Delta }m_\mathrm{N}/m_\mathrm{N}=\mathrm{\Delta }m_\mathrm{B}/m_\mathrm{B}=\mathrm{\Delta }\mathrm{\Lambda }_{\mathrm{QCD}}/\mathrm{\Lambda }_{\mathrm{QCD}}`$. That the same value is obtained by varying either the fine-structure constant or simultaneously the nucleon and meson masses (the first and fourth line in Table 1) can also be verified in first–order perturbation theory. In this case, the corresponding constraint of the fine structure constant is obtained by assuming that changes in $`\alpha _{\mathrm{em}}`$ are of the order $`3060`$ larger than for $`\mathrm{\Lambda }_{\mathrm{QCD}}`$ . Therefore, the constraint obtained for $`\alpha _{\mathrm{em}}`$ is one to two orders of magnitude more stringent than the one by varying solely the fine structure constant.

## 3 The Oklo phenomenon

The Oklo phenomenon that occurred about 2 Gyr ago gives the most stringent constraints for variations of fundamental parameters like the fine structure constant as it is seen from Table 2. The Oklo phenomenon refers to a natural fission reactor which was operating at the Oklo uranium mine in Gabon. By examining the isotopic ratios of Sm in the Oklo reactor limits on cosmological time variations of fundamental parameters can be obtained . These bounds were derived by calculating the energy shift of a resonant state in 150Sm lying just 0.0937 eV above the threshold of the reaction $`{}_{}{}^{149}\mathrm{Sm}+\mathrm{n}^{150}\mathrm{Sm}+\gamma `$. The quantity of interest is the energy of that state $`E_\mathrm{R}=QE^{}=0.0937`$ eV, where $`Q`$ is the Q–value of the reaction and $`E^{}`$ is the excitation energy of the compound nucleus 150Sm.

We employ the Relativistic Mean–Field theory (RMFT)  to calculate the ground–state properties of 149Sm and 150Sm nuclei. We obtain the variation in the Coulomb energy difference between the two nuclei by varying the fine structure constant $`\alpha _{\mathrm{em}}`$ fractionally by 10-5, 10-4 and 10-3. The change in the Coulomb energy difference shows a linear relationship with the variation in $`\alpha _{\mathrm{em}}`$. We obtain for the Coulomb energy difference between 150Sm and 149Sm a value of 1.18 MeV, which is only slightly larger than the values obtained from the phenomenological Bethe–Weizsäcker formula. Therefore, the results claimed in recent work  remain unchanged in the light of RMFT–calculations.

We also varied the masses of the nucleon and the $`\sigma `$–, $`\omega `$– and $`\rho `$–mesons. Fractional changes were made of 10-5, 10-4 and 10-3 in the masses. The total energy difference between the isotopes 150Sm and 149Sm shows again a linear behaviour with the variation in the masses. The change in the total energy difference turns out to be about a factor of 10 larger than the above corresponding change in the Coulomb energy due to the variation in the fine structure constant. This would correspond to a constraint for the fine structure constant that is again as before at least two factors of magnitude more stringent than the one obtained by solely varying the fine structure constant.

In any case the Oklo phenomenon gives a more stringent constraint than the analysis of atomic multiplet spectra in quasar absorption lines as can be seen from Table 2. A possible explanation could be a non–linear temporal or even a non–monotonic variation of the fine structure constant as has been suggested by some recent papers, since the Oklo phenomenon occurred at a much later time than the absorption lines in the quasar spectra were created.

## Acknowledgments

Supported by OTKA-D32513/T037548, FKFP-0242-2000/0147-2001, and by the John Templeton Foundation (983-COS153).


text:
# Searching for dark matter halos in the Suprime-Cam 2 sq deg field

## 1 Introduction

Dark matter halos, which are the virialized systems formed via the gravitational amplification of initial density fluctuations, are one of the most important cosmological probes. Since their evolution is governed purely by gravity, their number density and the evolution of clustering properties can be precisely predicted analytically as well as numerically. Weak lensing is a unique observational technique to detect dark matter halos because it detects halos not by their light but by more essential physical parameter, the mass. It also allows us to directly estimate the halo mass bypassing a use of ambiguous empirical relations like the mass-temperature (or velocity dispersion) or mass-luminosity relations. Therefore, comparison of theories and weak lensing observations is a direct and reliable way to test the paradigm of the structure formation and to constrain the cosmological parameters.

Theoretical prescriptions for predicting statistics of weak lensing halos have been developed so far. Kruse & Schneider (1999) first calculated the expected number of dark matter halos detected by means of the aperture mass method (Schneider, 1996). Predicted number density of halos exceeds 10 per deg2 for cluster-normalized cosmologies if one obtains moderately deep images whose galaxy density is $`n_g`$30 arcmin-2. Bartelmann et al. (2001) pointed out that the number density is sensitive to the density profile of halos. On the other hand, Jain & van Waerbeke (2000) discussed the statistics of the peaks on convergence map. They claimed that the statistics can be an alternative approach to probe the mass function of dark matter halos and to constrain $`\mathrm{\Omega }_m`$.

Observationally, weak lensing has become reliable measures to map the dark matter distribution in clusters of galaxies (see a review by Mellier 1999). Recently, serendipitous detections of mass concentrations outside clusters on a reconstructed mass map have been reported by several authors (Erben et al., 2000; Umetsu & Futamase, 2000; Wittman et al., 2001; Miralles et al., 2002). All cases have no clear optical counterparts except the one by Wittman et al. (2001) which is identified as a cluster of galaxies at $`z=0.276`$.

These findings demonstrate that dark matter halo search via weak lensing is feasible and can detect even lower luminosity dark matter halos that are overlooked in conventional optical cluster surveys. The number of samples obtained so far is, however, limited. In order to make statistically significant discussions, weak lensing survey over a representative cosmic volume is crucial. The survey requires not only wide field but also high image quality to make reliable lensing analysis. Suprime-Cam on the 8 m Subaru telescope is the most ideal camera for the purpose. The median seeing in $`I_c`$ band monitored over a year and a half is reported to be $``$ 0.6 arcsec (Miyazaki et al., 2002). In this Letter, we report the first results from the Suprime-Cam weak lensing survey that covers contiguous 2.1 deg2 field. We examine statistical properties of local peaks on the reconstructed convergence map by comparing with model predictions.

## 2 Observation and data reduction

We set field size at 2.1 deg2, which is the largest size possible to be accomplished during the commissioning phase of the Suprime-Cam. The field we chose is centered at R.A. = 16h04m43s, decl. = +4312’19” (J2000.0). In order to compare light and mass we observed one of the seven fields that Gunn et al. (1986) made cluster surveys (GHO survey). The 2.1 deg2 field has twelve GHO clusters and one Abell cluster (A 2158). Note, however, that the number density of clusters in this field is comparable to the one averaged over the seven GHO fields.

We obtained $`R_c`$-band images on the nights of 2001 April 23-25. Suprime-Cam provides a field of view of 34’$`\times `$ 27’ with a scale 0”.202 pixel-1. Nine contiguous fields were observed in 3$`\times `$3 mosaic configuration. Each exposure on a given field was offset by 1’ from other exposures to remove cosmic rays and defects on the CCDs. Total exposure time on each field was 1800 sec (360 sec $`\times `$ 5).

The individual images were de-biased and then flattened using a median of all object frames taken during the observing run. Registration on each field was first performed assuming a simple geometrical model. The parameters of the model include (1) optical distortion, (2) displacement and rotation of each CCD from the fiducial position and (3) the offset of pointing between the exposures (Miyazaki et al., 2002). These parameters can be solved by minimizing position errors of control stars common to all exposures and used for the first transformation of the images to a common coordinate system. After the transformation, the residual of the position of the control stars is $``$ 0.5 pixel rms. The residual is due to several effects not considered in the simple model including atmospheric dispersion and asymmetric aberration of optics. We note, however, that the residual vector is smoothly changing over position and can be well modeled as the third order bi-linear polynomial function of position. This model then gives the fine correction of the transformation function. As a result the final residual decreased typically down to 0.07 pixel rms (14 miliarcsec). The seeing in the resulting image is 0”.68 FWHM and the scatter among the fields is quite small of 0”.04 rms.

We used the software suite imcat, an actual implementation of Kaiser et al. (1995), to carry out object detection, photometry and shape measurements of objects. Catalogs created for the nine fields are registered using stars on overlapping regions to result in a final catalog whose total field of view is 1.64 $`\times `$ 1.28. Differences in photometric zero point among the fields due to variation of sky condition were compensated at this stage but the adjustment was not significant ($``$ 0.05 mag). For weak lensing analysis, we adopt galaxies of $`22<R<25.5`$ whose signal-to-noise (S/N) ratio, $`nu`$, calculated in the imcat exceed 15 (Bacon et al., 2001). The number of objects in the catalog is 297547 (39.3 arcmin-2) and the count peaks at R = 25.2.

## 3 Galaxy shape analysis

The estimate of the shear, $`\stackrel{}{\gamma }`$, from the observed ellipticities of galaxies, $`\stackrel{}{e}=\{I_{11}I_{22},2I_{12}\}/(I_{11}+I_{22})`$ involves two steps where $`I_{ij}`$ are quadrupole moments of the surface brightness of objects. First, PSF anisotropy is corrected using the image of stars as references,

$$\stackrel{}{e^{}}=\stackrel{}{e}\frac{P_{sm}}{P_{sm}^{}}\stackrel{}{e^{}},$$
(1)
where $`P_{sm}`$ is the smear polarisability tensors and is mostly diagonal (Kaiser et al., 1995). $`(P_{sm}^{})^1\stackrel{}{e^{}}`$ was calculated for stars scattered over the field of view and we made fifth order bi-polynomial fit to the values as a function of the positions. This function is used in Eqn.( 1) to correct the ellipticities of faint galaxies. The nine fields are treated independently in this procedure. The rms value of ellipticities of the reference stars, $`|\stackrel{}{e^{}}|^2^{\frac{1}{2}}`$, is reduced from 2.7% to 1.3% as a result of the correction. Note that the rms before the correction is already small thanks to the superb image quality of Subaru telescope.

Luppino & Kaiser (1997) worked out how to correct the ellipticities for the effect of the seeing. The pre-seeing shear $`\stackrel{}{\gamma }`$ is described as

$$\stackrel{}{\gamma }=P_\gamma ^1\stackrel{}{e^{}},P_\gamma =P_{sh}\frac{P_{sh}^{}}{P_{sm}^{}}P_{sm},$$
(2)
where $`P_{sh}`$ is the shear polarisability tensor. The $`P_\gamma `$ of individual galaxies are, however, known to be noisy estimate and we thus adopted smoothing and weighting method developed by van Waerbeke et al.(2000; see also Erben et al. 2001 for a detail study of the smoothing scheme). For each object, twenty neighbors are first identified in the $`r_g`$-$`magnitude`$ plane where $`r_g`$ is a measure of size of objects which the imcat yields. A median value of $`P_\gamma `$ among the neighbors is adopted for a smoothed $`P_\gamma `$ of the object. The variance of raw $`\gamma `$ before the smoothing among the neighbors, $`\sigma _\gamma ^2`$, is used to estimate the weight of the object, $`w_n`$, as $`w_n=1/(\sigma _\gamma ^2+\alpha ^2)`$ where $`\alpha ^2`$ is the variance of all the objects in the catalog and $`\alpha 0.4`$ here. Under the weighting scheme, an averaged value of a certain observable, $`Q`$, is calculated as $`_{i=1}^Nw_{n,i}Q_i/_{i=1}^Nw_{n,i}`$ instead of $`_{i=1}^NQ_i/N`$. The $`\mathrm{tr}(P_\gamma )`$ for all objects is calculated as $`0.31`$ but $`P_\gamma `$ depends on $`r_g`$. For objects with $`r_g>3.5`$, $`\mathrm{tr}(P_\gamma )`$ is almost constant, $`0.65`$, and becomes small below the $`r_g`$ down to $`\mathrm{tr}(P_\gamma )0.2`$ for $`r_g=1.5`$.

## 4 Results and discussions

The convergence field, $`\kappa `$, is estimated using the original Kaiser & Squires (1993) inversion algorithm with a Gaussian window function $`W(\theta )=\mathrm{exp}(|𝜽|^2/\theta _G^2)/\pi \theta _G^2`$. In order to estimate the noise of $`\kappa `$ field, we randomized the orientations of the galaxies in the catalog and created a $`\kappa _{\mathrm{noise}}`$ map. The variance of the $`\kappa _{\mathrm{noise}}`$ depends on the smoothing scale $`\theta _G`$ and is represented as $`\sigma _{\mathrm{noise}}0.023\theta _G^1`$. This turned out to be consistent with theoretical estimate,

$$\sigma _{\mathrm{noise}}=\frac{\sigma _ϵ}{2\theta _G\sqrt{\pi n_g}},$$
(3)
where $`\sigma _ϵ=0.42`$ and $`n_g=39.3`$arcmin-2 are the observed variance of galaxy ellipticities and number density, respectively, yielding $`\sigma _{\mathrm{noise}}=0.020\theta _G^1`$.

We repeated the randomization 100 times and compute the rms among them at each grid where $`\kappa `$ is computed. Since $`\kappa `$ is Gaussian in its errors, this rms represents 1-$`\sigma `$ noise level on each grid, and thus the measured signal divided by the rms gives the S/N ratio, $`\nu `$, at that grid. Thick contour in Fig. 4 shows the S/N-map with a smoothing scale of $`\theta _G=1^{}`$. Threshold of the contour is set at 3.

One notes that numbers of dark halo candidates are detected. For reference, number density of moderately bright galaxies are super-imposed as thin contour. It is surprising that correlation of mass and light is clearly visible from this relatively simple analysis on the single band image. A halo that has the highest significance is designated (X) in Fig. 4. In this halo, the galaxy density also shows an excess as large as the other known clusters. In fact, a spectroscopic observation of two member galaxies implies that it is a cluster located at the redshift of 0.41. This object is not listed on the NASA/IPAC Extragalactic Database (NED). Other follow-up observations are underway and detailed comparison of mass and light will be discussed in subsequent papers.

Let us turn to statistical properties of the S/N-map, and we will compare the results with the model predictions. First, we examine the probability distribution function (PDF) of peaks. Peaks are identified by pixels having a higher or lower value of $`\nu `$ than the all eight surrounding pixels as suggested by Jain & van Waerbeke (2000). The PDF of the peak heights is shown as thick histograms of left two panels in Fig. 4. The thin histogram shows a noise PDF measured from $`\kappa _{\mathrm{noise}}`$ maps, the average of 100 randomizations is plotted and the error bars show the rms among them. The bimodal feature of the noise PDF is due to superposition of two contributions each of which is due to positive peaks and negative peaks (= troughs). Excess counts over the noise PDF is apparent on both positive ($`\nu 3`$) and negative ($`\nu 3`$) side of the histogram. The excess becomes more apparent as $`\theta _G`$ increases. The positive peaks imply the existence of dark matter halos whereas the natural explanation of the negative peaks seen is voids (Jain & van Waerbeke, 2000). This is the first observational evidence of the existence of voids in the matter distribution revealed by weak lensing.

Right panels of Fig. 4 show the peak PDF measured from numerical experiments of the gravitational lensing by large-scale structures. The experiments were performed using the ray-tracing technique combined with large N-body simulations and its details are described in Ménard et al. (2002). Briefly, $`N`$-body data from Very Large simulation carried out by the Virgo Consortium (Jenkins et al 2001; Yoshida, Sheth & Diaferio 2001) were used to generate the dark matter distribution. A $`\mathrm{\Lambda }`$CDM model ($`\mathrm{\Omega }_m=0.3`$, $`\mathrm{\Omega }_\lambda =0.7`$, $`h=0.7`$ and $`\sigma _8=0.9`$) is assumed. The multiple-lens plane ray-tracing algorithm was used to follow the light rays (Hamana & Mellier 2001, see also Jain, Seljak & White 2000 for technical detail). The lensing convergence and shear were computed for 10242 pixels with the pixel size of 0.25 arcmin for a single source plane of $`z_s=1`$. We generated the simulated noisy $`\kappa `$ map following the procedure described in Jain & van Waerbeke (2000). The noise is modeled as a Gaussian random field with the variance measured in the real data. The peak is identified as the same manner as for the observed data. Ten realizations of the numerical experiment were used to compute the peak PDF (the total field size is $`10\times 4.27^2`$ deg2, much larger than the observed data, and thus the peak distribution is very smooth).

The simulated PDF is very similar to the observed PDF. Jain & van Waerbeke (2000) pointed out that the peak PDF can be used to discriminate the cosmological model. They demonstrated using the ray-tracing experiments of a high density ($`\mathrm{\Omega }_m=1.0`$, $`\mathrm{\Omega }_\lambda =0`$, $`\sigma _8=0.9`$) and a low density($`\mathrm{\Omega }_m=0.3`$, $`\mathrm{\Omega }_\lambda =0`$, $`\sigma _8=0.9`$) models that the high density model has broader distribution than the low density model, especially on the negative side. It should be noticed that they adopted $`\sigma _ϵ=0.2`$ which is about half of our measured value, therefore one should not directly compare our real PDFs with theirs (Figure 3 of Jain & van Waerbeke (2000)).

In addition to the limited statistics of this work, theoretical prescription of the peak PDF is, however, still premature for placing constraint on the cosmological model in a quantitative manner. We therefore concentrate only on excess on the positive side here; the halo number count. We calculated number of the positive peaks, $`N_{obs}`$, in the observed S/N map exceeding $`\nu _{\mathrm{th}}`$. Halo number count is defined as $`N_{halo}=N_{obs}N_r\pm \sqrt{N_{obs}+\sigma _{N_r}^2}`$ where $`N_r`$ and $`\sigma _{N_r}`$ is average peak counts obtained on the randomized catalogs and their variances, respectively. The halo count is shown in Fig. 4 as well as theoretical predictions which are constructed following Kruse & Schneider (1999). We assumed a high positive peak in $`\kappa `$ map comes from a single halo whose number density is described by Press-Schechter mass function (Press & Schechter 1974; we used its modified model by Sheth & Tormen 1999). We further assumed that all halos have a universal density profile that Navarro et al. (1996) proposed (NFW); $`\rho (r)=\rho _sr^1(1+r)^2`$. We also compared the results with another halo profile, singular isothermal sphere model (SIS); $`\rho (r)=\rho _sr^2`$. For the relation between the scale length of the halo and halo mass we adopted the fitting form by Bullock et al (2001). The $`\sigma _{\mathrm{noise}}`$ measured from observed data is used to convert $`\kappa `$ to $`\nu `$. Note that Kruse & Schneider (1999) and Bartelmann et al. (2001) adopted $`\sigma _ϵ=0.2`$, thus their predicted count is much larger than our prediction.

First, if one considers only highly significant halos, say $`\nu 5`$, the predicted halo number count does not depends strongly on the background cosmology as is presented in Fig.3, where dashed line and solid line shows the SCDM model ($`\mathrm{\Omega }_m=1`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0`$, $`\sigma _8=0.6`$, $`h=0.5`$) and $`\mathrm{\Lambda }`$CDM ($`\mathrm{\Omega }_m=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, $`\sigma _8=0.9`$, $`h=0.7`$), respectively. The count, however, depends on the halo profile; the SIS predicts less counts as is shown in the dotted line in the figure. The observed count prefers NFW profile although the poisson error in $`\nu 6`$ is large in this work. It is demonstrated here that number count of highly significant halos can be an useful probe for the dark matter halo profile because of the independence of background cosmology. Meanwhile, the dependence becomes clearly visible in the counts of moderately significant haloes ($`\nu <5`$); the SCDM predicts more counts than LCDM does due mainly to larger $`\mathrm{\Omega }_m`$. Our observed result shows high number counts which favors high density universe. This result, however, cannot be taken at face value since the large excess is apparent in Fig.4 as a spike around $`\nu 3.5`$. Therefore, the excess count in the moderately significant halos is likely due to the cosmic variance. Comparison of data taken on multiple (say 5 to 10) uncorrelated fields of this size is necessary to evaluate the contribution from the cosmic variance and eventually to discuss cosmological models.

In summary, the present work clearly demonstrates that weak lensing technique is indeed a promising way to probe matter distribution purely through the mass. We detected not only the high positive peaks due to massive halos but also the excess in the peak PDF on the negative side which is the first observational evidence of the voids in dark matter distribution. Dark matter halo counts is proved to be the useful tool to test the paradigm of the structure formation and cosmological models. The cosmic variance that we encountered in this work should be evaluated by the data obtained by future high resolution wider field imaging survey.

We are grateful to L. van Waerbeke for helpful discussions about galaxy shape analysis, and M. Bartelmann for useful discussions. We thank Richard Ellis and Alexandre Refregier for their careful reading of the manuscript and the comments. T.H. F.N. and M.O. acknowledge support from Research Fellowships of the Japan Society for the Promotion of Science.


text:
# Abnormal phenomena in a one-dimensional periodic structure containing left-handed materials

## Abstract

The explicit dispersion equation for a one-dimensional periodic structure with alternative layers of left-handed material (LHM) and right-handed material (RHM) is given and analyzed. Some abnormal phenomena such as spurious modes with complex frequencies, discrete modes and photon tunnelling modes are observed in the band structure. The existence of spurious modes with complex frequencies is a common problem in the calculation of the band structure for such a photonic crystal. Physical explanation and significance are given for the discrete modes (with real values of wave number) and photon tunnelling propagation modes (with imaginary wave numbers in a limited region).

Left-handed materials (LHM) with negative permittivity and negative permeability, which were first suggested in veselago , have attracted much attention recently smith ; shelby ; pendry1  thanks to the experimental realization of such materials science  and the debate on the use of a LHM slab as a perfect lens to focus both propagating waves and evanescent waves pendry ; prl ; surface .

Multi-layered structures containing LHM have been investigated through calculating the transmittance or reflectance of the structure matrix ; bragg . Effects of photon tunnelling and Bragg diffraction have been observed in these works. In the present paper, we consider a one-dimensional (1D) photonic crystal formed by alternative layers of LHM and RHM. The band structure is investigated by analyzing the explicit dispersion equation. Some abnormal phenomena such as spurious modes with complex frequencies, discrete modes and photon tunnelling modes are found and explained in the present paper.

Consider a 1D periodic structure with alternating layers $`(ϵ_1,\mu _1)`$ and $`(ϵ_2,\mu _2)`$ as shown in Fig. 1. $`d_1`$ and $`d_2`$ are the widths of the two inclusion layers respectively, and $`a=d_1+d_2`$ is the period. Note that the conservation of energy requires that $`ϵ_i\mu _i>0`$, $`i=1,2`$. The corresponding refractive index is given by $`n_i=\pm \sqrt{ϵ_i\mu _i}`$ (negative sign for LHM). We consider an off-line propagation of monochromatic electromagnetic field (with time-dependence $`e^{i\omega t}`$) in a periodic structure with off-line wave vector $`\beta `$ along the $`x`$ axis. For the E-polarization case, one has

$`\{\begin{array}{ccc}E_{1y}=e^{i\beta x}(e^{ik_1z}+Ae^{ik_1z})\hfill & & \\ H_{1x}=\frac{k_1}{\omega \mu _1}e^{i\beta x}(e^{ik_1z}Ae^{ik_1z})\hfill & & \\ H_{1z}=\frac{\beta }{\omega \mu _1}e^{i\beta x}(e^{ik_1z}+Ae^{ik_1z})\hfill & & \end{array}\text{in region 1},`$ (4)
and

$`\{\begin{array}{ccc}E_{2y}=e^{i\beta x}(Be^{ik_2z}+Ce^{ik_2z})\hfill & & \\ H_{2x}=\frac{k_2}{\omega \mu _2}e^{i\beta x}(Be^{ik_2z}Ce^{ik_2z})\hfill & & \\ H_{2z}=\frac{\beta }{\omega \mu _2}e^{i\beta x}(Be^{ik_2z}+Ce^{ik_1z})\hfill & & \end{array}\text{in region 2},`$ (8)
where $`k_i`$ is the component of the wave vector along the $`z`$axis in region $`i`$ ($`i=1,2`$), i.e., $`k_i^2=\omega ^2ϵ_i\mu _i/c^2\beta ^2`$. Here $`c`$ is the speed of light in vacuum.

The tangential electric and magnetic fields should be continuous at $`z=0`$, i.e.,

$`\{\begin{array}{cc}E_{1y}(z=0^{})=E_{2y}(z=0^+)\hfill & \\ H_{1x}(z=0^{})=H_{2x}(z=0^+),\hfill & \end{array}`$ (11)
To obtain the dispersion relation for this 1D photonic crystal, we need to use the following periodic conditions according to the Bloch theorem,

$`\{\begin{array}{cc}E_{2y}(z=d_2)=E_{1y}(z=d_1)e^{iqa}\hfill & \\ H_{2x}(z=d_2)=H_{1x}(z=d_1)e^{iqa},\hfill & \end{array}`$ (14)
where $`q`$ is in the first Brillouin zone $`\frac{\pi }{a}q\frac{\pi }{a}`$ . Substituting systems (1) and (2) into systems (3) and (4), we obtain the following dispersion relation for the E-polarization case,

$`\mathrm{cos}(k_1d_1)\mathrm{cos}(\stackrel{~}{n}k_1d_2){\displaystyle \frac{\stackrel{~}{n}^2+\mu ^2}{2\stackrel{~}{n}\mu }}\mathrm{sin}(k_1d_1)\mathrm{sin}(\stackrel{~}{n}k_1d_2)=\mathrm{cos}(qa),`$ (15)
where $`\stackrel{~}{n}=k_2/k_1=\sqrt{n^2(k_1^2+\beta ^2)\beta ^2}/k_1`$, $`\mu =\mu _2/\mu _1`$, and $`n=n_2/n_1`$ is the index ratio of the two media.

Similarly, we can derive the following dispersion relation for the H-polarization case,

$`\mathrm{cos}(k_1d_1)\mathrm{cos}(\stackrel{~}{n}k_1d_2){\displaystyle \frac{\stackrel{~}{n}^2+ϵ^2}{2\stackrel{~}{n}ϵ}}\mathrm{sin}(k_1d_1)\mathrm{sin}(\stackrel{~}{n}k_1d_2)=\mathrm{cos}(qa),`$ (16)
where $`ϵ=ϵ_2/ϵ_1`$.

From analytical equations (5) and (6), we can find out how $`k_1`$ depends on $`q`$. The corresponding dispersion relation between $`\omega `$ and $`q`$ can then be obtained according to $`\omega ^2=c^2(k_1^2+\beta ^2)/(ϵ_1\mu _1)`$ for each fixed $`\beta `$ ($`\beta =0`$ is associated to the in-line propagation, i.e, the case when the propagation direction is along the normal of the material interfaces).

When the index ratio $`n<0`$, the solution to the dispersion equation and the associated band structure have quite different behaviors as compared to those for the usual case when $`n>0`$. For simplicity, we first illustrate this for the case when $`\beta =0`$. In this case, $`\frac{\stackrel{~}{n}^2+\mu ^2}{2\stackrel{~}{n}\mu }=\frac{\stackrel{~}{n}^2+ϵ^2}{2\stackrel{~}{n}ϵ}`$and the dispersion equations for the E-polarization and the H-polarization have the same form (as expected).

Spurious modes with complex $`\omega `$. When one of the constitutive media is of LHM (i.e., $`n<0`$), the analytical dispersion equation may have some complex solutions for $`k_1`$. These complex values of $`k_1`$ lead to complex $`\omega `$ since $`\omega ^2=c^2(k_1^2+\beta ^2)/(ϵ_1\mu _1)`$. Fig. 2 gives the in-line (i.e., $`\beta =0`$) band structure when $`n=2.5,\mu =2,d_1=d_2=0.5a`$. The real part of $`\omega `$ varies little as $`q`$ varies from $`\pi /a`$ to $`\pi /a`$. The problem of complex $`\omega `$ will also appear when $`\beta 0`$ (i.e., the case of off-line propagation). Obviously, these complex $`\omega `$ have no physical significance since physically the frequency $`\omega `$ must have a real value.

One would meet a similar problem of complex $`\omega `$ in the calculation of the band structure using other methods such as the plane wave expansion method (one needs to remove these spurious complex $`\omega `$ in the band structure). In the plane wave expansion method, the $`\omega `$ values are determined from the following eigenvalue equation, $`\frac{1}{\mu }\times (\frac{1}{ϵ}\times 𝐇)=\frac{\omega ^2}{c^2}𝐇`$. The existence of complex $`\omega `$ indicates that the operator $`\frac{1}{\mu }\times (\frac{1}{ϵ})`$ is not a Hermit operator or has no equivalent Hermit operator. One can also expect that these complex $`\omega `$ (associated with an exponential increasing factor of time) will make the finite-difference time domain (FDTD) method (with a periodic boundary condition) divergent in the calculation of the band structure of a photonic crystal consisting of alternate RHM and LHM inclusions (with constant permittivity and permeability) PRE ; chan ; wu . We have verified numerically with the plane wave expansion method and the FDTD method the existence of the spurious modes with complex $`\omega `$.

The fields associated with the spurious complex $`\omega `$ satisfy Maxwell’s equations, the boundary condition and the enforced periodic condition. Therefore, an additional requirement that the frequency $`\omega `$ must be real should be applied to make the solution significant in physics.

In the rest of the paper, we consider only situations when $`\omega `$ is real.

When $`\beta =0`$, a complex solution of $`k_1`$ always leads to a spurious complex $`\omega `$. However, when $`\beta 0`$ some complex solutions of $`k_1`$ may correspond to real $`\omega `$ and thus have physical significance. Since $`\omega ^2=c^2(k_i^2+\beta ^2)/(ϵ_i\mu _i)`$, a complex solution of $`k_ik_{iR}+jk_{iI}`$ (with $`k_{iI}0`$) should satisfies the following conditions in order to make $`\omega `$ real,

$`k_{iR}=0,k_{iI}^2<\beta ^2.`$ (17)
In other words, an imaginary solution of $`k_1`$ can have physical significance when $`\beta 0`$ (i.e., the off-line propagation case). Imaginary solutions of $`k_1`$ correspond to evanescent waves. In a periodic structure consisting of alternate RHM and LHM layers, the special property that the evanescent waves (in both directions) decrease in the RHM layer and increase in the LHM layer pendry  or vice versa permits the evanescent waves to propagate in the periodic structure (i.e., a propagation mode, for which the energy density does not decay over a period; we will call it a photon tunnelling mode and will discuss it in more details below). In the periodic structure composed solely by RHM or LHM, complex solutions of $`k_1`$ will not occur since evanescent waves can not form a propagation mode in these structures.

As discussed before, imaginary solutions of $`k_1`$ may have physical significance when $`\beta 0`$. This will be illustrated by the abnormal photon tunnelling modes in the band structure for a photonic crystal consisting of alternate RHM and LHM layers. Another abnormal (and interesting) phenomenon is that some discrete modes (with real values of wave number) are also found in the band structure.

To simplify the analysis, we consider a special case when $`n=1`$ for the E-polarization. Then $`\stackrel{~}{n}=1`$ (independent of $`\beta `$) and Eq. (15) is simplified to

$`\mathrm{cos}(k_1d_1)\mathrm{cos}(k_1d_2){\displaystyle \frac{1+\mu ^2}{2\mu }}\mathrm{sin}(k_1d_1)\mathrm{sin}(k_1d_2)=\mathrm{cos}(qa).`$ (18)
One notices that the above equation for $`k_1`$ is independent of $`\beta `$.We discuss situations when $`\mu 1`$ (so that the RHM and LHM layers have different values of impedance). It is particularly interesting when the two inclusion layers have the same width (i.e., $`d_1=d_2=0.5a`$) since in this case only discrete modes and photon tunnelling modes exist (conventional propagating band disappears in the band structure).

Discrete modes. Since $`\mathrm{cos}(k_1d_1)\mathrm{cos}(k_1d_2)+\frac{1+\mu ^2}{2|\mu |}\mathrm{sin}(k_1d_1)\mathrm{sin}(k_1d_2)1`$ (when $`k_1`$ is real), real solutions of $`k_1`$ to Eq. (8) exist only when $`q=0`$ (notice that $`\pi /aq\pi /a`$), which corresponds to discrete solutions $`k_1=2N\pi /a`$, $`N=0,\pm 1,\pm 2,\mathrm{}`$. Since these discrete solutions satisfy the transparent condition of $`k_1d_1=N\pi ,k_2d_2=N\pi `$ (one can easily show that each layer is transparent to neighboring layers at these discrete frequencies for a fixed $`\beta `$), the field with this wave vector will surely pass through the periodic structure and becomes a propagating mode. On the other hand, fields with any other wave vector (i.e., any other frequency for the fixed $`\beta `$) satisfies the Bragg reflective condition $`k_1d_1+k_2d_2=0`$ and can not propagate in this periodic structure. Therefore, only those discrete propagation modes exist. As an example, the discrete real solutions of $`k_1`$ to Eq. (8) are shown by the circles in Fig. 3 (with $`\mu =2`$). This abnormal phenomenon (existing for any $`\beta `$) can be utilized to make a very narrow filter (with no sidelope, which is different from any conventional type of filters).

Photon tunnelling modes. From Eq. (8) one sees that only complex $`k_1`$ can be the solutions if $`q0`$. Then according to condition (7) the solutions $`k_1`$ with physical significance must be imaginary. The expression on the left-hand side of Eq. (18) always decreases when $`k_{1I}`$ increases or decreases away from $`0`$ (a symmetric function of $`k_{1I}`$ with a maximum at $`k_{1I}=0`$). Therefore, only two conjugate imaginary solutions exist. The dashed lines in Fig. 3 show the imaginary solutions of $`k_1`$ to Eq. (8) for an example. When $`q(0)`$ varies in the first Brillouin zone, the two conjugate imaginary solutions vary continuously in a region from $`0.561j\times 2\pi /a`$ to $`0.561j\times 2\pi /a`$. Note that this figure is independent of $`\beta `$ for this special case. For off-line waves with $`\beta ^2>k_{iI}^2`$, the imaginary solutions correspond to special propagation modes$``$photon tunneling modes (as discussed before). For this example, evanescent modes with $`k_{1I}`$ outside the region $`[0.561j\times 2\pi /a,0.561j\times 2\pi /a]`$ can not propagate inside the periodic structure. From Fig. 3 one can also see that $`k_{1I}`$ approaches 0 when $`q`$ approaches 0. Thus, any $`\beta 0`$ can permit the existence of some photon tunnelling modes since $`k_{1I}`$ can be arbitrarily small to satisfy $`k_{1I}^2<\beta ^2`$. From the above discussion, one sees that for any $`\beta 0`$ there always exist some evanescent modes which can propagate inside the periodic structure while some other evanescent modes can not (one notices that in the case of a perfect lens pendry  all the evanescent waves can go through the lens). This limitation of the photon tunnelling modes is due to the enforced periodic condition. Thus, the problem of square non-integrable field mentioned in prl  will not occur in the present situation (this is also true for the next example when $`d_1d_2`$).

Both discrete modes and photon tunnelling modes exist when $`\beta 0`$. However, when $`\beta =0`$ (as in the above example), which corresponds to normally incident waves (i.e., the in-line propagation case), only discrete modes can exist (this means that only fields at discrete frequencies can exist in the structure).

In the above discussion the case when $`d_1=d_2`$ is considered. Fig. 4(a) shows the solutions of $`k_1`$ for an example when $`d_1d_2`$. From this figure one sees that non-discrete real $`k_1`$ solutions (corresponding to some continuous bands shown in Fig. 4(b) exist besides some discrete modes and photon tunnelling modes. The discrete solutions are at $`k_1=(2.5+5N)\times 2\pi /a`$ for $`q=\pm \pi /a`$ and at $`k_1=5N\times 2\pi /a`$ for $`q=0`$, $`N=0,\pm 1,\pm 2,`$. They are located in the forbidden band-gaps (see Fig. 4(b) for the corresponding band structure for a given $`\beta `$).

In the above we have discussed only the ideal situations when the periodic structure consists of ideal materials with lossless and non-dispersive negative $`ϵ`$ and negative $`\mu `$ ($`ϵ`$ and $`\mu `$ are the relative permittivity and permeability of the two inclusion media). In a non-ideal situation with material loss or dispersion the abnormal phenomena such as spurious modes with complex $`\omega `$, discrete modes and photon tunnelling modes can also be expected to occur for the following reasons.

If the inclusion media are dispersive, $`ϵ`$ and $`\mu `$ will be frequency-dependent. Since the solutions can be obtained from Eq. (5) in a point-by-point manner, they are correct locally (i.e., in a small frequency range where the frequency-dependence of $`ϵ`$ and $`\mu `$ can be neglected). Thus, the abnormal phenomena can still exist in certain small frequency region. Furthermore, the frequency in the band structure is scalable and thus we can scale it to a frequency at which the material parameters have desired values in order to observe these abnormal phenomena.

If the LHM inclusion is lossy the relative $`ϵ`$ or $`\mu `$ has a very small positive imaginary part (e.g., $`10^5j`$; as claimed in some papers that LHM must have some loss prl ). The abnormal phenomena can still occur since the dispersion equation (Eq. (5)) will not give any rapid change when a very small imaginary part is added to $`ϵ`$ or $`\mu `$. Note that for the photon tunnelling modes $`q`$ should be complex (corresponding to a pseudo-periodic field with a damping factor along the propagation direction) when $`ϵ`$ or $`\mu `$ is complex. Since the imaginary parts of these parameters are very small, Eq. (5) is still valid when all the parameters (including $`ϵ`$, $`\mu `$ and $`q`$) are replaced with their real parts.

In conclusion, we have found some abnormal phenomena such as spurious modes with complex $`\omega `$, discrete modes and photon tunnelling modes by analyzing the explicit dispersion equation for a 1D periodic structure with alternating LHM and RHM layers. One will meet the problem of complex $`\omega `$ in the calculation of the band structure using other methods such as the plane wave expansion method and the FDTD method. The physical significance of the discrete modes and photon tunnelling modes has been explained. The discrete modes can be utilized to make a very narrow filter with no sidelope. For an off-line wave there exist some photon tunnelling propagation modes while some other evanescent modes can not propagate inside the periodic structure.

Acknowledgement. The partial support of National Natural Science Foundation of China (under a key project grant; grant number 90101024) is gratefully acknowledged.