format: text

text:
# 1 Status of Observational Cosmology

## 1 Status of Observational Cosmology

The ten most significant parameters to be determined in observational cosmology are:

1. Age of the universe: $`t_0`$
2. The Hubble constant $`H_0`$ from the Hubble relation: $`v=H_0d`$
3. Density parameter : $`\mathrm{\Omega }_m=\frac{\rho _m}{3H_0^2/8\pi G}`$
4. Deceleration parameter: $`q_0=\frac{\ddot{R}(t_0)R(t_0)}{\dot{R}^2(t_0)}`$
5. The baryon density $`\mathrm{\Omega }_B`$ and the vacuum density $`\mathrm{\Omega }_\mathrm{\Lambda }`$
6. The parameters associated with microwave backgound fluctuations: $`n,\sigma _8,T/S,N_T`$

The strongest lower limit for $`t_0`$ is determined from studies of the stellar populations of globular clusters. The main error in the globular clusters age estimate comes from the uncertain distance to the globular clusters. A 0.25 magnitude error in the distance translates into a $`22\%`$ error in the cluster age . Independent age limits come from the cooling of white dwarfs. The best estimates give $`t_013`$ Gyr, with a lower limit of $`11`$ Gyr. For $`t_0>13`$ Gyr, we have $`h0.50`$ (where $`h`$ is defined by the Hubble parameter: $`H_0=100`$ h km $`s^1\mathrm{Mpc}^1`$) for matter density $`\mathrm{\Omega }_m=1`$, and for $`h0.73`$ we have $`\mathrm{\Omega }_m0.3`$ in spatially flat cosmologies with $`\mathrm{\Omega }_m+\mathrm{\Omega }_\mathrm{\Lambda }=1`$.

The Hubble parameter is now better determined (it used to be known to within a factor of two). Most measurements are now consistent with a value: $`h=0.65\pm 0.08`$ . It is remarkable that data obtained from several different methods for determing $`H_0`$ lead to similar results, which gives hope for an ultimate convergence of measurements. For $`\mathrm{\Omega }_m=0.4`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }=0.6`$ and $`h=0.65\pm 0.08`$, the age of the universe would be $`t_0=13\pm 2`$ Gyr, in agreement with globular cluster age estimates. This result is one of the strong arguments for a low matter density $`\mathrm{\Omega }_m0.3`$ and a non-zero cosmological constant $`\mathrm{\Omega }_\mathrm{\Lambda }0.7`$.

A most promising new way of measuring $`\mathrm{\Omega }_m`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }`$ on cosmological scales is to use small-angle anisotropies in the CMB radiation and high-redshift Type Ia supernovae (SNe Ia). The Supernovae Cosmology Project (CSP) and the High-Z Supernovae team have found a significant number of Type Ia supernovae . The more recent larger SCP data set of 42 high redshift data gives for the flat case $`\mathrm{\Omega }_m=0.28_{0.080.04}^{+0.09+0.05}`$ . The High-Z Supernovae group has also measured $`\mathrm{\Omega }_m`$ giving in the flat case $`\mathrm{\Omega }_m=0.4\pm 0.3`$. Two possible sources of problems are the dimming by dust and the assumption made that evolution for nearby and far supernovae is uniform.

In CMB anisotropy studies, the location of the first acoustic Doppler peak at angular wave number $`l250`$ is a strong indication of a flat universe $`\mathrm{\Omega }_m+\mathrm{\Omega }_\mathrm{\Lambda }=1`$. The MAXIMA and BOOMERANG balloon flights seem to confirm this result, and the existence of a second and possible third peak would appear to be consistent with the predictions of simple inflation models. New data from the NASDA Microwave Anisotropy Probe satellite will hopefully strengthen these results.

We can summarise the main observational results:

1. Age of universe$`t_0=916\mathrm{Gyr}`$ (from globular clusters) $`=917\mathrm{Gyr}^{}`$
2. Hubble parameter$`H_0=100hs^1Mpc^1,h=0.65\pm 0.08`$.
3. Baryon density$`\mathrm{\Omega }_bh^2=0.019\pm 0.001`$ (from $`D/H`$) $`>0.015`$ from Ly$`\alpha `$ forest opacity\*
4. Matter density$`\mathrm{\Omega }_m=0.4\pm 0.2`$ (from cluster baryons) $`=0.34\pm 0.1`$ from Ly$`\alpha `$ forest $`P(k)^{}(P(k)=Ak^n\mathrm{with}n=1`$ for the Harrison-Zel’dovich spectrum) $`=0.4\pm 0.2`$ from cluster evolution\* $`>\frac{3}{4}\mathrm{\Omega }_\mathrm{\Lambda }\frac{1}{4}\pm \frac{1}{8}`$ from SN Ia $`>0.3(2.4\sigma `$ from flows)
5. Total density$`\mathrm{\Omega }_m+\mathrm{\Omega }_\mathrm{\Lambda }=1\pm 0.3`$ (from CMB peak location)
6. Dark vacuum energy density$`\mathrm{\Omega }_\mathrm{\Lambda }=0.8\pm 0.3`$ (from last two lines)
7. Neutrino density$`\mathrm{\Omega }_\nu 0.001`$ (from Superkamiokande) $`0.1^{}`$

Here, the cosmological parameters with * are obtained by assuming $`\mathrm{\Lambda }CDM`$ i.e cold dark matter models with non-zero cosmological constant.

The Type Ia reshift measurements have indicated the remarkable result that the cosmic expansion is presently undergoing an acceleration.

## 2 Bimetric Gravity Theory

A new kind of vector-tensor and scalar-tensor theory of gravity, which exhibits a bimetric structure and contains two or more light cones \[6, Clayton+Moffat:1999a, 7\], has been introduced, recently. This type of model has attracted some attention , and similar effects have been noted elsewhere . The motivation for considering these models is derived form earlier work , which provided a scenario in which some of the outstanding issues in cosmology can be resolved. These models provide a fundamental dynamical mechanism for varying speed of light theories and generate a new mechanism for an inflationary epoch that could solve the initial value problems of early universe cosmology. In the following, we shall review some of the main features of this new kind of gravity theory and its application to cosmology.

In these models matter that satisfies the strong energy condition can nevertheless contribute to the cosmic acceleration. Our cosmological model can be mapped to a model with varying fundamental constants , albeit not uniquely and requiring some care in the interpretation of the varying constants that appear.

It is hoped that the models can shed some light on the new observational data that suggests the expansion of the universe at present is undergoing an acceleration . Although there has been some success in understanding the latter problem by the inclusion of a class of very particular scalar field potentials , it is fair to say that not all issues have been resolved. Using the scalar field version of the model, we expect that not only will we be able to generate sufficient inflation, but that a quintessence-like solution should be achieveable.

We shall be considering models wth an action of the form

$$S=\overline{S}_{\mathrm{gr}}[\overline{g}]+S[g,\psi ]+\widehat{S}[\widehat{g},\widehat{\varphi }^I].$$
(1)
The first term is the usual Einstein-Hilbert action for general relativity constructed from a metric $`\overline{g}_{\mu \nu }`$, and the final term is the contribution from the non-gravitational (matter) fields in spacetime $`\widehat{\varphi }^I`$, and is built from a different but related metric $`\widehat{g}_{\mu \nu }`$.

The contribution $`S[g,\psi ]`$ is constructed from a metric $`g_{\mu \nu }`$ and includes kinetic terms for a field or fields (unspecified as yet) $`\psi `$ that may be considered to be part of the gravitational sector, modifying the reaction of spacetime to the presence of the matter fields in $`\widehat{S}[\widehat{g},\widehat{\varphi }^I]`$. The manner in which $`\psi `$ accomplishes this is by modifying the metric that appears in each of the actions. For example, in  $`\psi `$ was a vector field, $`\overline{g}_{\mu \nu }=g_{\mu \nu }`$ and $`\widehat{g}_{\mu \nu }=g_{\mu \nu }+b\psi _\mu \psi _\nu `$, whereas in \[Clayton+Moffat:1999a\] $`\psi `$ was a scalar field, $`\overline{g}_{\mu \nu }=g_{\mu \nu }`$ and $`\widehat{g}_{\mu \nu }=g_{\mu \nu }+b_\mu \psi _\nu \psi `$. These relations imply that matter and gravitational fields propagate at different velocities if $`\psi `$ is non-vanishing.

Since the matter action $`\widehat{S}`$ is built using only $`\widehat{g}_{\mu \nu }`$, it is the null surfaces of $`\widehat{g}_{\mu \nu }`$ along which matter fields propagate. If we assume that other than the presence of a “composite” metric the matter action is otherwise a conventional form (perfect fluid, scalar field, Maxwell, etc.), then variation of the matter action yields the matter energy-momentum tensor $`\widehat{T}^{\mu \nu }`$, which will be conserved

$$\widehat{}_\nu \widehat{T}^{\mu \nu }=0,$$
(2)
by virtue of the matter field equations $`\widehat{F}_I=0`$. Throughout we will write, for example, $`\widehat{}_\nu `$ for the covariant derivative constructed from the Levi-Civita connection of $`\widehat{g}_{\mu \nu }`$. Since we also assume that the matter fields satisfy the dominant energy condition, we therefore know (assuming appropriate smoothness of $`\widehat{g}_{\mu \nu }`$) that matter fields cannot travel faster than the speed of light as determined by $`\widehat{g}_{\mu \nu }`$.

The gravitational action is written

$$\overline{S}_{\mathrm{gr}}[\overline{g}]=\frac{1}{\kappa }𝑑\overline{\mu }\overline{R},$$
(3)
where we use a metric with ($`+`$$``$$``$$``$) signature and have defined $`\kappa =16\pi G/c^4`$. We will denote the metric densities by, e.g., $`\overline{\mu }=\sqrt{det(\overline{g}_{\mu \nu })}`$ and in addition write $`d\overline{\mu }=\overline{\mu }dtd^3x`$. We will not consider a cosmological constant, since it can easily be included later. We can identify the metric $`\overline{g}_{\mu \nu }`$ as providing the light cone for the gravitational system.

We consider a Proca model with arbitrary potential

$$S[g,\psi ]=\frac{1}{\kappa }𝑑\mu \left(\frac{1}{4}B^2V(X)\right),$$
(4)
where we will use the definition

$$X=\frac{1}{2}\psi ^2,$$
(5)
and $`V^{}(X)=V(X)/X`$. We will also use $`B_{\mu \nu }=_\mu \psi _\nu _\nu \psi _\mu `$, $`\psi ^2=g^{\mu \nu }\psi _\mu \psi _\nu `$ and $`B^2=g^{\alpha \mu }g^{\beta \nu }B_{\alpha \beta }B_{\mu \nu }`$. We will assume that as $`\psi _\mu 0`$ we have $`V(X)m^2X`$ and therefore the linearized (in $`\psi _\mu `$) limit of our model is identical to Einstein-Proca field equations coupled to matter. The standard energy-momentum tensor for the vector field is

$$T^{\mu \nu }=B^{\mu \alpha }B_{}^{\nu }{}_{\alpha }{}^{}+\frac{1}{4}g^{\mu \nu }B^2+V^{}\psi ^\mu \psi ^\nu Vg^{\mu \nu }.$$
(6)
Although there exists a more general class of models, we will limit ourselves to the choice

$$\widehat{g}_{\mu \nu }=g_{\mu \nu }+b\psi _\mu \psi _\nu ,\overline{g}_{\mu \nu }=g_{\mu \nu }+g\psi _\mu \psi _\nu ,$$
(7)
where $`b`$ and $`g`$ are constants, so that the variations of $`\widehat{g}_{\mu \nu }`$ and $`\overline{g}_{\mu \nu }`$ are related to those of $`g_{\mu \nu }`$ and $`\psi _\mu `$.

The field equations are given by

$$_\mu B^{\mu \nu }+V^{}\psi ^\nu +gT^{\mu \nu }\psi _\mu +\kappa \frac{\widehat{\mu }}{\mu }(gb)\widehat{T}^{\nu \mu }\psi _\mu =0,$$
(8)
$$\overline{\mu }\overline{G}^{\mu \nu }=\frac{1}{2}\mu T^{\mu \nu }+\frac{1}{2}\kappa \widehat{\mu }\widehat{T}^{\mu \nu }.$$
(9)
It is clear that $`\widehat{g}_{\mu \nu }`$ and $`\overline{g}_{\mu \nu }`$ provide the characteristic surfaces for matter and gravitational fields, respectively.

We can prove that any matter model that conserves energy-momentum with respect to $`\widehat{g}_{\mu \nu }`$ is consistent with the gravitational structure that we have introduced .

The “most physical” metric is clearly $`\widehat{g}_{\mu \nu }`$, since it describes the geometry on which matter propagates and interacts. Because all matter fields are coupled to the same metric $`\widehat{g}_{\mu \nu }`$ in exactly the same way, the weak equivalence principle is satisfied. Furthermore, because one can work in a local Lorentz frame of $`\widehat{g}_{\mu \nu }`$, in which non-gravitational physics takes on its special relativistic form, the Einstein equivalence principle is also satisfied. However, because $`\widehat{g}_{\mu \nu }`$ does not couple to matter in the same way as in general relativity unless $`\psi _\mu =0`$, the strong equivalence principle will be violated.

The main motivation for considering these theories is that they should have something to say about the horizon problem in the early universe. If $`\psi _\mu 0`$, then if we choose $`b>g`$, matter fields will propagate outside the light cone of the gravitational field. As $`\psi _\mu 0`$ the matter light cone will ‘contract’ and matter and gravitational disturbances will eventually propagate at the same velocity. If one considers a frame in which gravitational waves propagate at a constant speed, then as the light cone of matter contracts, the universe will appear to material observers to expand acausally.

## 3 Cosmology

Implicit in the idea of a varying light speed is that the speed of light is changing with respect to some fixed frame of reference. If one introduces a fundamental frame for this, then it is perhaps sensible to introduce a function $`c(t,x)`$ to describe this variability . The models that we are considering are based on the idea that the speed of light can be changing with respect to the speed of gravitational disturbances, and therefore any indication of the speed of light as a function of spacetime is frame-dependent. In particular, we will see that a frame in which the speed of light is constant and the speed of gravitational disturbances is changing is connected via a diffeomorphism to a frame where the speed of gravitational disturbances is constant, and the speed of light is changing. Quantities of interest such as the local light cone, horizons, etc. are derived directly from the relevant metric, thereby avoiding any guesswork as to which ‘speed of light’ to use—the gravitational or electromagnetic . The constant $`c`$ is fixed in the present universe by making measurements of the electromagnetic field.

In a homogeneous and isotropic (FRW) universe, the vector field $`\psi _\mu `$ has components $`\psi _\mu =(c\psi _0(\tau ),0,0,0)`$. We will begin with the metric $`g_{\mu \nu }`$ in comoving form

$$g_{\mu \nu }dx^\mu dx^\nu =c^2d\tau d\tau R^2(\tau )\gamma _{ij}dx^idx^j,$$
(10)
and therefore

$$\widehat{g}_{\mu \nu }dx^\mu dx^\nu =\widehat{\mathrm{\Theta }}^2(\tau )c^2d\tau d\tau R^2(\tau )\gamma _{ij}dx^idx^j,$$

$$\overline{g}_{\mu \nu }dx^\mu dx^\nu =\overline{\mathrm{\Theta }}^2(\tau )c^2d\tau d\tau R^2(\tau )\gamma _{ij}dx^idx^j.$$
(11)
The spatial metric in spherical coordinates has the standard form

$$\gamma _{ij}=\mathrm{diag}(1/(1\mathrm{kr}^2),\mathrm{r}^2,\mathrm{r}^2\mathrm{sin}^2\theta ),$$
(12)
and we have defined

$$\widehat{\mathrm{\Theta }}=\sqrt{1+2bX},\mathrm{and}\overline{\mathrm{\Theta }}=\sqrt{1+2gX},$$
(13)
where from (5) we have $`X=\frac{1}{2}\psi _0^2`$.

Although we begin with the choice (10), once we have derived the field equations, we will make a coordinate transformation in order to put $`\widehat{g}_{\mu \nu }`$ in comoving form and thereby make a comparison with the standard cosmological results a simpler matter. Note that we are reversing the definitions of $`t`$ and $`\tau `$ as used in our previous article .

The matter energy-momentum tensor will have a perfect fluid form

$$\widehat{T}^{\mu \nu }=\left(\rho +\frac{p}{c^2}\right)\widehat{u}^\mu \widehat{u}^\nu p\widehat{g}^{\mu \nu },$$
(14)
where we have written the velocity field as $`\widehat{u}^\mu `$ to emphasize that it is normalized using the metric $`\widehat{g}_{\mu \nu }`$, so that

$$\widehat{g}_{\mu \nu }\widehat{u}^\mu \widehat{u}^\nu =c^2.$$
(15)
The matter conservation laws (2) lead to the usual relation

$$_\tau \rho +3\frac{_\tau R}{R}\left(\rho +\frac{p}{c^2}\right)=0.$$
(16)
The Friedmann equations take the form

$$\left(\frac{_\tau R}{R}\right)^2+\frac{kc^2\overline{\mathrm{\Theta }}^2}{R^2}=\frac{\kappa c^4}{6}\overline{\mathrm{\Theta }}^3\left[\frac{\rho }{\widehat{\mathrm{\Theta }}}+\frac{1}{\kappa c^2}(2XV^{}V)\right],$$
(17)
$$2\frac{_\tau ^2R}{R}+\left(\frac{_\tau R}{R}\right)^2+\frac{kc^2\overline{\mathrm{\Theta }}^2}{R^2}2\frac{_\tau R}{R}\frac{_\tau \overline{\mathrm{\Theta }}}{\overline{\mathrm{\Theta }}}=\frac{\kappa c^2}{2}\overline{\mathrm{\Theta }}\left[\widehat{\mathrm{\Theta }}p+\frac{1}{\kappa }V\right].$$
(18)
The single remaining Proca field equation from (8) is

$$\frac{1}{c\widehat{\mathrm{\Theta }}}\psi _0\left[\widehat{\mathrm{\Theta }}\left(\overline{\mathrm{\Theta }}^2V^{}gV\right)\kappa (bg)c^2\rho \right]=0.$$
(19)
We now perform the coordinate transformation

$$dt=\widehat{\mathrm{\Theta }}d\tau ,$$
(20)
and defining

$$\eta =\frac{\overline{\mathrm{\Theta }}}{\widehat{\mathrm{\Theta }}}=\sqrt{\frac{1+2gX}{1+2bX}},$$
(21)
we see that the metric $`\widehat{g}_{\mu \nu }`$ is put into comoving form

$$\widehat{g}_{\mu \nu }dx^\mu dx^\nu =c^2dtdtR^2(t)\gamma _{ij}dx^idx^j,$$

$$\overline{g}_{\mu \nu }dx^\mu dx^\nu =\eta ^2(t)c^2dtdtR^2(t)\gamma _{ij}dx^idx^j.$$
(22)
We have

$$\left(\frac{\dot{R}}{R}\right)^2+\frac{kc^2\eta ^2}{R^2}=\eta ^2\frac{\kappa c^4}{6}\rho _{\mathrm{eff}},$$
(23)
$$2\frac{\ddot{R}}{R}+\left(\frac{\dot{R}}{R}\right)^2+\frac{kc^2\eta ^2}{R^2}2\frac{\dot{R}}{R}\frac{\dot{\eta }}{\eta }=\eta ^2\frac{\kappa c^2}{2}p_{\mathrm{eff}},$$
(24)
where $`\dot{\rho }=_t\rho `$. In (23), we have defined the effective energy and pressure densities as

$$\rho _{\mathrm{eff}}=\eta \left(\rho +\frac{1}{\kappa c^2}\widehat{\mathrm{\Theta }}(2XV^{}V)\right),p_{\mathrm{eff}}=\frac{1}{\eta }\left(p+\frac{1}{\kappa \widehat{\mathrm{\Theta }}}V\right).$$
(25)
The reason for making these definitions is that (23) has exactly the form of the Friedmann equations for the metric $`\overline{g}_{\mu \nu }`$, and therefore these effective energy and momentum densities will also satisfy the conservation laws (16).

The function $`R(t)`$ is written in comoving coordinates and, therefore, the speed of light is constant. This emphasizes that having a ‘varying speed of light’ is a frame-dependent statement. In a frame where the speed of matter propagation (including electromagnetic fields) is constant, the speed of gravitational waves will be changing. In a frame where the speed of gravitational waves is constant, the speed of matter propagation will be changing. This, of course, is as it should be, since we have not introduced any nondynamical preferred frame into our model.

In the following we will specialize to a model where the vector field potential is a simple mass term:

$$V=m^2X,V^{}=m^2.$$
(26)
In this case (25) becomes

$$\rho _{\mathrm{eff}}=\eta \left(\rho +(bg)\rho _{pt}\widehat{\mathrm{\Theta }}X\right),$$

$$p_{\mathrm{eff}}=\frac{1}{\eta }\left(p+c^2(bg)\rho _{pt}\frac{X}{\widehat{\mathrm{\Theta }}}\right).$$
(27)
The nontrivial solution ($`\psi _00`$) of the field equation (19) leads to

$$\rho =\rho _{pt}\widehat{\mathrm{\Theta }}(1+gX),$$
(28)
where

$$\rho _{pt}=\frac{m^2}{\kappa c^2(bg)},H_{pt}=\sqrt{\frac{c^2m^2}{6(bg)}},$$
(29)
are the density at which $`\psi _0^2=0`$ is reached, and the inverse Hubble time at which this occurs (assuming that $`k=0`$).

We can now write the acceleration parameter as observed by material observers from (23) as

$$\widehat{q}=\frac{\ddot{R}}{H^2R}=\frac{\kappa c^4}{12}\frac{\eta ^2}{H^2}\left(\rho _{\mathrm{eff}}+\frac{3}{c^2}p_{\mathrm{eff}}\right)\frac{\dot{\eta }}{H\eta },$$
(30)
where we have defined the Hubble function $`H=\dot{R}/R`$. We have

$$\frac{\dot{\eta }}{H\eta }=3\frac{(bg)}{\rho _{pt}\overline{\mathrm{\Theta }}^2\widehat{\mathrm{\Theta }}(g+b+3bgX)}\left(\rho +\frac{p}{c^2}\right).$$
(31)
## 4 The Very Early Universe

For very short times following the initial singularity, we expect that $`\psi _0`$ is large, and if we assume that $`gX1`$ and $`bX1`$, then from (28) we find that

$$\rho =\rho _{pt}\sqrt{2b}gX^{3/2}.$$
(32)
This results in the Friedmann equation

$$\left(\frac{\dot{R}}{R}\right)^2+\frac{k\overline{c}^2}{R^2}=\frac{\overline{\kappa }\overline{c}^4}{6}\rho ,$$
(33)
where

$$\overline{c}=c\sqrt{\frac{g}{b}},\overline{G}=G\sqrt{\frac{g}{b}},\overline{\kappa }=\frac{16\pi \overline{G}}{\overline{c}^4}.$$
(34)
Although the behaviour of the solutions are well-known, it is worth pointing out that the ‘effective’ constants $`\overline{c}`$ and $`\overline{G}`$ are not interpretable as the effective speed of light and gravitational constant, rather they are effective constants that dictate how the gravitational field reacts to the presence of matter. Matter fields continue to propagate with speed $`c`$ consistent with (3). It is the gravitational field perturbations that propagate with speed $`\overline{c}`$, which is the justification for the notation.

During this phase there is clearly no inflation, but the horizon scales of the gravitational field and matter fields are related by

$$\overline{d}_H(t)=\frac{\overline{c}}{c}\widehat{d}_H(t)=\sqrt{\frac{g}{b}}\widehat{d}_H(t),\mathrm{where}\widehat{d}_H(t)=cR(t)_0^t\frac{ds}{R(s)},$$
(35)
with a similar definition for $`\overline{d}_H(t)`$ using the metric $`\overline{g}_{\mu \nu }`$. Because we have $`g<b`$ we expect that not only is the speed of gravitational disturbances slower than that of matter, but also that the coupling between matter and the gravitational sector is also lessened.

What we have here is very close to what was originally envisaged by one of us in . This is part of the motivation for including the $`g0`$ possibility, the other is that the approach to the initial singularity in this phase follows the same path as in ordinary GR+matter models, with a re-interpretation of the parameters. In this case we have a model that interpolates between this initial period where $`\overline{c}>c`$ and the later universe where $`\overline{c}=c`$.

## 5 Inflation and Light Cone Contraction

As $`\psi _0`$ decreases towards the point where $`gX1`$ the solution will no longer be a good approximation. If we now consider the solution when $`gX1`$, from (28) we have

$$\widehat{\mathrm{\Theta }}=\frac{\rho }{\rho _{pt}},\mathrm{or}X=\frac{1}{2b}\left[\left(\frac{\rho }{\rho _{pt}}\right)^21\right],$$
(36)
and the Friedmann equation (23) becomes

$$\left(\frac{\dot{R}}{R}\right)^2+\frac{kc^2\eta ^2}{R^2}\left(\frac{\rho _{pt}}{\rho }\right)^2=\frac{\kappa c^4}{12}\rho _{pt}\left[1+\left(\frac{\rho _{pt}}{\rho }\right)^2\right].$$
(37)
In this limit

$$\rho _{\mathrm{eff}}+\frac{3}{c^2}p_{\mathrm{eff}}=\frac{1}{\rho _{pt}}\left[\rho \left(\rho +\frac{3}{c^2}p\right)+\rho ^2\rho _{pt}^2\right],$$
(38)
which is greater than zero if the strong energy condition is satisfied, since $`\rho \rho _{pt}`$, and (30) reduces to

$$\widehat{q}=\frac{\kappa c^4\rho _{pt}}{12H^2}\left[\frac{1}{\rho }\left(\rho +\frac{3}{c^2}p\right)+1\left(\frac{\rho _{pt}}{\rho }\right)^2\right]\frac{3}{\rho }\left(\rho +\frac{p}{c^2}\right).$$
(39)
Since we expect that $`H^2`$ is large in the early universe (we can arrange that $`\rho _{pt}\rho _c`$ where $`\rho _c=12H^2/(\kappa c^4)`$, it is clear from (39) that even if matter satisfies the strong energy condition, the final term will dominate and $`\widehat{q}<0`$ (unless, perhaps, the weak energy condition is also violated). This is the expansion of the universe as seen by material observers. The acceleration of the gravitational geometry $`\overline{q}`$ would lack the final term and therefore $`\overline{q}>0`$.

That we get inflation was demonstrated previously , where an exact solution for $`k=0`$ and $`g=0`$ was found. Although we discovered that we could not get enough expansion to solve the horizon problem with pure radiation, a slowly rolling scalar field could provide the necessary negative pressure. The role that the extra structure of our model plays is that the fine-tuning that is required in a simple scalar field, potential-driven model is alleviated.

The flatness problem requires a bit more explanation. Dividing (23) by $`H^2`$ and defining

$$\overline{ϵ}=\frac{kc^2\eta ^2}{(\dot{R})^2},$$
(40)
we find a differential equation that $`\overline{ϵ}`$ satisfies by taking a derivative and using (23) to give

$$\dot{\overline{ϵ}}=\frac{\kappa c^4\eta ^2}{6H}\overline{ϵ}\left(\rho _{\mathrm{eff}}+\frac{3}{c^2}p_{\mathrm{eff}}\right).$$
(41)
Therefore, since $`\overline{ϵ}>0`$ and $`H>0`$ in the early universe, the only way for $`\overline{ϵ}=0`$ to be an attractor for (23) is for $`\rho _{\mathrm{eff}}+\frac{3}{c^2}p_{\mathrm{eff}}<0`$ at least for part of the history of the universe. What is not so obvious is whether the quantity $`\overline{ϵ}`$ as defined in (40) is of physical relevance.

The quantity of geometrical importance is the $`3`$-curvature of the spacelike slices, $`6k/R^2`$, which suggests that the physically meaningful quantity to examine would be

$$\widehat{ϵ}=\frac{kc^2}{(\dot{R})^2},$$
(42)
which has the equation of motion

$$\dot{\widehat{ϵ}}=2\widehat{ϵ}\widehat{q}.$$
(43)
Another way of stating this is that the curvature radius defines $`\mathrm{\Omega }`$ through

$$R_{\mathrm{curv}}=\frac{R}{|k|^{1/2}}=\frac{c}{H(|\mathrm{\Omega }1|)^{1/2}},$$
(44)
and so $`\widehat{ϵ}=|\mathrm{\Omega }1|`$. Since we found from (39) that $`\widehat{q}<0`$ in the early universe, clearly $`\widehat{ϵ}=0`$ is an attractor for (23), and since it is most-likely the quantity of physical importance for matter physics, we can also claim to have solved the flatness problem once the horizon problem is solved.

## 6 Conclusions

In our bimetric model, the universe generically accelerates ($`\widehat{q}<0`$) during some period in the early universe, and in the same period the physical importance of spatial curvature diminishes ($`|\mathrm{\Omega }1|`$ is decreasing). This can occur even when the matter fields satisfy the strong energy condition.

The model that we have considered generalizes that which appeared in  in a way that more closely follows the scenario discussed in . In the very early universe, matter and gravitational fields propagate with different and approximately constant velocities. During a period in which the matter light cone, originally much larger than the light cone of gravity, contracts, material observers will see an acausal expansion of the universe similar to inflation. Because the light cone of gravity does not undergo the same contraction, we expect there to be an observable difference in the scalar versus tensor contributions to the cosmic microwave background anisotropies.

## 7 Acknowledgements

This work was supported by the Natural Sciences and Engineering Research Council of Canada.


text:
# The dynamical evolution of the fragmented, bipolar dust shell around the carbon star IRC +10 216 based on observations performed with the 6 m telecope at the Special Astrophysical Observatory, Russia

## 1 Introduction

IRC +10 216 (CW Leo) is the nearest and best–studied carbon star and one of the brightest infrared sources. It experiences a strong mass loss at a rate of $`\dot{M}25\times 10^5`$M$`_{}`$yr-1 (see e.g. Loup et al. LoupForveilleEtAl93 (1993)). The central star of IRC +10 216 is a long–period variable star (LPV) with a period of approximately 649 days (Le Bertre LeBertre92 (1992)). Recent distance estimates of 110 pc to 135 pc (Groenewegen Groenewegen97 (1997)) and 150 pc (Crosas & Menten CrosasMenten97 (1997)) were reported. IRC +10 216’s initial mass can be expected to be close to 4 M (Guelin et al. GuelinForestiniEtAl95 (1995), Weigelt et al. WeigeltBalegaEtAl98 (1998)). The bipolar appearance of the nebula around this object was already reported by Christou et al. (ChristouRidgwayEtAl90 (1990)) and Kastner & Weintraub (KastnerWeintraub94 (1994)). The non-spherical structure is consistent with the conjecture that IRC +10 216 is in a phase immediately before entering the protoplanetary nebula stage. High–resolution observations of this object and its circumstellar dust shell were reported by McCarthy et al. (McCarthyMcLeodEtAl90 (1990)), Christou et al. (ChristouRidgwayEtAl90 (1990)), Dyck et al. (DyckBensonEtAl91 (1991)), Danchi et al. (DanchiBesterEtAl94 (1994)), Osterbart et al. (OsterbartBalegaEtAl97 (1997)), Weigelt et al. (WeigeltBalegaEtAl97 (1997), WeigeltBalegaEtAl98 (1998) WeigeltBloeckerEtAl99 (1999)), Skinner et al. (SkinnerMeixnerEtAl98 (1998)), and Haniff & Buscher (HaniffBuscher98 (1998)). The results of Dyck et al. (DyckBensonEtAl91 (1991)) and Haniff & Buscher (HaniffBuscher98 (1998)) showed that the structure of the dust shell of IRC +10 216 has been changing for some years. Detailed radiative transfer calculations for IRC +10 216 were recently performed by Ivezić & Elitzur (IvezicElitzur96 (1996)), Crosas & Menten (CrosasMenten97 (1997)), and Groenewegen (Groenewegen97 (1997)) using a large amount of spectroscopic and visibility data. The aim of this paper is to discuss the properties of the inner dust shell of IRC +10 216 on the basis of a series of high–resolution observations. In a second paper (Men’shchikov et al., in prep.) we will present a detailed two–dimensional radiative transfer model for this object.

## 2 Observational results

The IRC +10 216 speckle interferograms were obtained with the 6 m telescope at the Special Astrophysical Observatory in Russia and our Nicmos 3 camera at four epochs and with our Hawaii–array camera at one epoch (Nov. 1998). Table 1 lists the observational parameters.

### 2.1 $`J`$-, $`H`$-, and $`K`$-band reconstructions

Figs. 1 and 2 show the $`K`$ and $`H`$ images, respectively, of the central region of IRC +10 216 for all epochs. The high–resolution images were reconstructed from the speckle interferograms using the bispectrum speckle interferometry method (Weigelt Weigelt77 (1977), Lohmann et al. LohmannWeigeltEtAl83 (1983), Weigelt Weigelt91 (1991)). The achieved resolution of the images depends slightly on the data quality (seeing and number of recorded interferograms) and is indicated in the figure captions. The object power spectra were determined with the speckle interferometry method (Labeyrie Labeyrie70 (1970)). The speckle transfer functions were derived from speckle interferograms of the unresolved stars mentioned in Table 1.

We denote the resolved components in the center of the nebula as A, B, C, and D (see Figs. 1g and 2b) in the order of decreasing peak intensity (based on the $`K`$ band results from 1996). In addition, Fig. 2b shows three fainter components denoted with E, F, and G. In Fig. 2c, cuts through the centers of components A and B of the $`H`$ and $`K`$ images from 1997 are shown to illustrate the differences in the relative intensities of A and B for different wavelengths.

Fig. 3 shows images of the faint nebula around IRC +10 216, which seems to have a bipolar structure. The resolution of these images was reduced to 149 mas ($`J`$) and 95 mas ($`H`$ and $`K`$) to increase the signal–to–noise ratio in the outer parts of the nebula. The faint granular structure of the images ($`0.5\%`$ of the peak intensity) is partly caused by speckle noise.

Note that the faint extended feature at position angle PA $``$340° in the $`J`$ image corresponds quite well to a very faint component (denoted with E in Figs. 1g and 2b) visible in all images in Figs. 1 and 2 (assuming that the brightest component in the $`J`$ image is roughly coinciding with component A in the $`H`$ and $`K`$ images). Two other faint features at PA $``$20° and PA $``$50° (F and G in Fig. 2b) are only visible in some of the images.

### 2.2 Polarimetry

Fig. 4 shows the results of polarimetric observations with the HST Nicmos camera at a wavelength of 1.1 $`\mu `$m (raw data retrieved from the Hubble Data Archive, STScI). The data were obtained on April 30, 1997 at a photometric phase of $`\mathrm{\Phi }=0.76`$. Data for three polarization angles were available (0°, 120°, and 240°). It has been carefully checked that the three images were correctly registered. The images have been obtained without new pointing of the telescope so that no correction was necessary (see Voit Voit97 (1997)). From the data the total intensity (contours in Fig. 4a and b), the polarized intensity (Fig. 4c), the degree and the position angle of the polarization have been derived (see e.g. Voit Voit97 (1997), Fischer et al. FischerHenningEtAl94 (1994)). Fig. 4b shows the derived polarization map superimposed on the contours of the total intensity. In Fig. 4a the two images were superimposed by centering the brightest peaks onto each other. Looking at the southern tails of A and at the northern arms this centering seems to be appropriate. The structure of the 1.1 $`\mu `$m HST image at the position of B is affected by the diffraction pattern of A. Nevertheless, the contours indicate that a structure corresponding to B might be present in the 1.1 $`\mu `$m image.

### 2.3 The core structure and the faint nebula

The multi–component structure of the IRC +10 216 $`K`$–band image (of the innermost 300 mas $`\times `$ 300 mas) has already been reported by Weigelt et al. (WeigeltBalegaEtAl98 (1998)) and Haniff & Buscher (HaniffBuscher98 (1998)). This bright inner region is surrounded by a larger faint nebula (with $`1\%`$ of the peak brightness of A; see Fig. 3b and c). Three arms of the nebula can be seen at position angles of $``$30° (NE), 340° (NW), and 210° (SW) with respect to component A. At 160° (counterside of the NW-arm) the nebula is much fainter. The direction from component A to component B (PA$``$20°) can be taken as the direction of the main axis (see also Kastner & Weintraub KastnerWeintraub94 (1994)). This direction fits well to the main axis of the H13CN($`J=10`$) emission (Dayal & Bieging DayalBieging95 (1995)) which is weakly elongated on a scale of about 10′′. On larger scales asymmetries are not observed (Dayal & Bieging DayalBieging95 (1995), Groenewegen et al. GroenewegenVanDerVeenEtAl97 (1997), Mauron & Huggins MauronHuggins99 (1999)).

### 2.4 The bipolar structure in $`J`$ and at shorter wavelengths

The $`J`$–band image (Fig. 3a) and the 0.79 $`\mu `$m and 1.06 $`\mu `$m HST images (Haniff & Buscher HaniffBuscher98 (1998), Skinner et al. SkinnerMeixnerEtAl98 (1998), see also Fig. 4a for the total intensity of the HST 1.1 $`\mu `$m polarization data) show a bipolar shape of the nebula. The southern lobe has a cometary or fan–shaped structure whereas the northern area of the images shows two arms reminiscent to (but weaker than) the northern X–shaped structure of the Red Rectangle (see, e.g., Men’shchikov et al. MenshchikovBalegaEtAl98 (1998)). However, the fact that even in the polarized intensity the nebula is very faint on the southeastern side suggests the main axis to be at PA$``$20° to 30°.

### 2.5 $`HK`$ color image

For the January 1997 data, the integral intensities in the full fields of view of our camera (5$`\stackrel{}{.}`$1 in $`H`$ and 7$`\stackrel{}{.}`$8 in $`K`$) were compared to those of the photometric standard stars HIP 71284 (BS 5447) and HD 106965 (cf. Elias et al. EliasFrogelEtAl82 (1982)). The resulting IRC +10 216 magnitudes are $`K=2.3`$ and $`H=5.4`$. According to an extrapolation of the fitted light curve of Le Bertre (LeBertre92 (1992)), IRC +10 216 was very close to its light minimum in January 1997. Our magnitudes are in good agreement with this prediction. In a square aperture of 1$`\stackrel{}{.}`$6 the magnitudes can be determined to be $`K=2.5\pm 0.1`$ and $`H=5.7\pm 0.1`$. The integral color in this field is thus $`HK3.2`$. The resolved two–dimensional $`HK`$ color image is shown in Fig. 5 together with the contours of the $`H`$–band image. The $`H`$ and $`K`$ images used for the determination of the $`HK`$ color were reconstructed with a common resolution of 95 mas which is also the resolution of the color image. The dependence of the $`HK`$ color image on relative shifts between the two images was investigated because it is not a priori known whether the positions of the intensity maxima in $`H`$ and $`K`$ coincide. The relative position of the components A and B is very similar in the $`H`$ and $`K`$ images from 1997 so that a solution could be found where both components A and B are almost coinciding (within a few mas) for the two images. We found that the color image is not very sensitive to relative shifts within the realistic uncertainties of some milli–arcseconds.

### 2.6 Separation of components A and B

Fig. 6 shows the separation of the components A and B in the $`K`$–band images for different epochs. Phase 0 in Fig. 6 corresponds to JD=2449430. This date of the photometric maximum was derived from the results of Le Bertre (LeBertre92 (1992)). The separations are: 191 mas, 201 mas, 214 mas, 245 mas, and 265 mas, for the 5 epochs from 1995 to 1998 shown in Fig. 1. The linear regression fit gives a value of 23 mas/yr for the average increase in the apparent separation of the components. Interpreting this increase as a real motion would lead to 14 km/s within the plane of the sky (for a distance of $`D=130`$ pc).

The apparent relative motion of the nebula components is obviously not simply related to the stellar variability which has a period of approximately 649 days (Le Bertre LeBertre92 (1992)). It may thus be related to either an overall expansion of an inhomogeneous circumstellar dust medium or a variability of the dust shell with a period significantly larger than the stellar pulsation period (cf. Fleischer et al. FleischerGaugerEtAl92 (1992), Winters et al. WintersFleischerEtAl94 (1994), WintersFleischerEtAl95 (1995)).

Table 2 lists the apparent relative velocities within the plane of the sky of the components A to E with respect to either A or B or the center positions defined by A and B or by A, B, C, and D. The values were determined from linear regression fits.

### 2.7 Structural changes within the nebula.

The images by Haniff & Buscher (HaniffBuscher98 (1998)) taken in 1989 and 1997, as well as the one–dimensional data collected by Dyck et al. (DyckBensonEtAl91 (1991)) covering a larger range of epochs, already show that the structure of the envelope changes on time scales of some years. Our $`K`$–band observations now allow us to study the changes within approximately 2 stellar pulsation cycles in more detail. Besides the apparent motion of the components, Fig. 1 shows that these components change their shapes and relative fluxes. The brightest component A appears to be rather symmetric in 1995. At the later epochs it approximately keeps its size perpendicular to the axis A–B ($``$20°) but becomes narrower along this axis. The peak–to–peak intensity ratio of B and A is approximately constant from 1995 to 1997. Thereafter the component B started fading. At the same time the other components become brighter and detached from A. Note that the photometric phases of IRC +10 216 in January 1997 and November 1998 are almost identical. In fact, the integral $`K`$ magnitudes were the same ($`K=2.3`$). Again we find that the time scale for the changes seen in our images is significantly different from the period of the stellar pulsation.

## 3 Discussion

The observations presented in the previous section as well as other high–resolution observations (e.g. Haniff & Buscher HaniffBuscher98 (1998)) provide detailed information on the variable structure in IRC +10 216. The question of where, behind all the dust, the central star is located, is of specific interest to understand the physical properties of the nebula.

At short wavelengths the nebula shows a bipolar structure. A comparison of the observed structure with other bipolar nebulae like the Red Rectangle (Men’shchikov et al. MenshchikovBalegaEtAl98 (1998)) suggests that the X–like arms originate mainly from scattering of stellar light on the surfaces of cavities. The star is then at least partially obscured by an optically thick dust shell or torus. Haniff & Buscher (HaniffBuscher98 (1998)) argued convincingly that the main axis of the object is tilted with its southern side towards the observer.

### 3.1 Is component A the star?

Motion of the components. From molecular line observations the terminal velocity of the expansion of the circumstellar envelope was found to be $`14`$ to 15 km/s (Dayal & Bieging DayalBieging95 (1995), Gensheimer & Snyder GensheimerSnyder97 (1997)). The velocities of B and D (14 km/s and 13 km/s, resp.) with respect to A can, therefore, only be understood if the direction of motion is approximately within the plane of the sky. However, the structure of the bipolar nebula at short wavelengths shows that we are looking neither pole–on nor edge–on at the nebula but at an intermediate viewing angle (e.g. 50° to 60°). The assumption that A is the star is thus not satisfactorily consistent with our observation.

Polarimetry. The synthetic polarization maps of Fischer et al. (FischerHenningEtAl96 (1996)) show that a significant polarization at the position of the star is only present for nearly edge–on configurations. The high degree of polarization of A ($`14\%`$) is thus not in agreement with A being the star. At larger separations from A the polarization pattern is centrosymmetric. The center of such a pattern is thought to be at the position of the illuminating source (Fischer et al. FischerHenningEtAl94 (1994), FischerHenningEtAl96 (1996)). The fact that in Fig. 4b this center does not coincide with A but is located significantly north of it, supports our interpretation that A is not the star.

### 3.2 Is the star at or near the position of B?

In the following (a to d) we will argue that the assumption that the star is at or near the position of B is consistent with the observations.

(a) The cometary shapes of A in the $`H`$ and the $`J`$ images and the 0.79 $`\mu `$m and 1.06 $`\mu `$m HST images (Haniff & Buscher HaniffBuscher98 (1998)), as well as the polarization data strongly suggest that A is part of a scattering lobe within a bipolar structure. Consistently, A and its southern tails are relatively blue ($`HK`$ ranging from 2 to 3.2 in Fig. 5) compared to the integral color.

(b) The northern components B, C, and D are, on the other hand, rather red ($`HK4.2`$) in comparison with the integral color. This suggests that these structures are strongly obscured and reddened by circumstellar dust.

(c) The brightest northern component in the $`J`$ image (Fig. 3a) can hardly be seen in the $`H`$ image and is thus very blue. It has a separation of $``$500 mas from A at a position angle of $``$27°. In the $`K`$–band image from the same epoch (Fig. 1b) component B is at a separation of $``$200 mas and at a position angle of 21°. This means that B is almost in the center between the northern X–shaped $`J`$-band component and the southern cometary component A (see Sect. 2.4). This morphology suggests that the latter components are the opposite lobes of a bipolar structure around a central star approximately located at the position of B.

(d) The polarization map (Fig. 4b) fits well with the picture that the star is at B. The region between the two $`J`$–band lobes is only weakly polarized with polarization null points in the east and northwest. Such null points indicate a transition from an optically thick region with multiple scattering to a region with predominantly single scattering of photons (Piirola et al. PiirolaScaltritiEtAl92 (1994), Fischer et al. FischerHenningEtAl94 (1994), FischerHenningEtAl96 (1996)). The center of the centrosymmetric polarization pattern is located between the two lobes. Because of the scatter in the direction of the polarization vectors it is not possible to determine very precisely if the center coincides with B, but the polarization pattern is consistent with this assumption. We note that at the position of B the polarized intensity (Fig. 4c) may be slightly contaminated by the wings of the HST diffraction pattern associated with the dominant component A. The polarimetry of A itself, however, is not affected. Since the polarimetric features discussed above and in Sec. 3.1 are only slightly influenced by this contamination, the conclusions drawn so far do not change.

Changes in the mass loss rate. The change of the shape of component A and the fading of B can be attributed to an increasing mass loss which is accompanied by a gradual increase of the optical depth of the dust shell. This is most obvious for the later observation epochs, suggesting an enhanced mass loss since 1997. A strongly variable mass loss has, in fact, been predicted by theoretical models treating the dust formation mechanism in the envelopes of long–period variable carbon stars (Fleischer et al. FleischerGaugerEtAl92 (1992), Winters et al. WintersFleischerEtAl94 (1994), WintersFleischerEtAl95 (1995)). Periods of this mechanism may be significantly longer than the stellar pulsation period (Winters et al. WintersFleischerEtAl95 (1995)). An increasing optical depth of the inner dust shell also results in an increasing apparent separation of the dust shell structures. The apparent motion of the components is thus not solely determined by the velocity of the dust particles but also by the changes of the optical properties of the circumstellar dust.

Alternative model: the star between A and B. From the present observational data it is not possible to exclude the possibility that the star may be located between A and B, close to B, or in the center of A, B, C, and D. In particular, the precision of the polarization map is not sufficient to conclude whether the star is at the position of B or only close to it.

Radiative transfer modeling. An answer to the question of where the star is located requires radiative transfer calculations. In a second paper (Men’shchikov et al., in preparation), we will present the results of our two–dimensional radiative transfer modeling showing that the shapes of A and B cannot be reproduced when assuming a position of the star between A and B. On the other hand, it was possible to reproduce these shapes and the intensity ratio of A and B in the case where the star is assumed to be at the position of B. Clear preference is thus given to the latter model.

### 3.3 Stellar evolution and bipolar structure

IRC +10 216 is without doubt in a very advanced stage of its AGB evolution due to its long pulsational period, high mass-loss rate, and carbon-rich dust-shell chemistry indicating that a significant number of thermal pulses have already taken place. The star’s initial mass can be estimated to be $`4`$M$`{}_{}{}^{}\pm 1`$M due to the observed isotopic ratios of C, N and O in the dust shell (Guelin et al. GuelinForestiniEtAl95 (1995)) and the luminosity of the central star (Weigelt et al. WeigeltBalegaEtAl98 (1998)). Accordingly, the core mass should be $`0.7`$ to $`0.8`$M with corresponding thermal-pulse cycle times of $`1\mathrm{3\hspace{0.17em}10}^4`$yr (Blöcker Bloecker95 (1995)). Introducing the mean observed mass-loss rate to these thermal-pulse periods shows that the present stellar wind leads to a very effective erosion of the envelope per thermal pulse cycle, possibly as high as $`1`$M/cycle. Consequently, the whole envelope may be lost during the next few thermal pulses leading to the termination of the AGB evolution. Thus, it is not unlikely to assume that IRC +10 216 has entered a phase immediately before moving off the AGB. This is strongly supported by the non-spherical appearence of its dust shell showing even bipolar structures. Unlike AGB stars, post-AGB objects as protoplanetary nebulae often expose prominent features of asphericities, in particular in axisymmetric geometry (e.g. Olofsson Oloff93 (1993), Oloff96 (1996)). Accordingly, IRC +10 216 can be thought to be an object in transition. It is noteworthy that the establishment of bipolar structures, i.e. the metamorphosis into a protoplanetary nebulae, obviously already begins during the (very end of) AGB evolution. The clumpiness within the bipolar shape is probably due to small scale fluctuations of the dust condensation radius which, in turn, might be influenced by, e.g., giant surface convection cells (Schwarzschild Schwarzschild75 (1975)). The formation of giant convection cells can be assumed to be a common phenomenon in red giants.

The shaping of planetary nebulae can successfully be described by interacting stellar wind models (Kwok et al. Kwoketal78 (1978), Kahn & West KahnWest85 (1985)). Within this scenario a fast (spherical) wind from the central star interacts with the slow wind of the preceding AGB evolution. The slow AGB wind is asssumed to be non-spherical (axisymmetric) which leads to the observed morphology of planetary nebulae (Mellema Mel96 (1996)).

The cause of an aspherical AGB mass loss is still a matter of debate. Different mechanisms to provide the required equatorial density enhancements are discussed (cf. Livio Liv93 (1993)). Among these, binarity is one channel including common envelope evolution and spin up of the AGB star due to the interaction with its companion (Morris Mor81 (1981), Bond & Livio BondLiv90 (1990)). Not only stellar companions are found to be able to spin up the AGB star but also substellar ones as brown dwarfs and planets, most effectively by evaporation in the AGB star’s envelope (Harpaz & Soker SokHar92 (1992), Soker Soker97 (1997)). Currently there is no observational evidence for a possible binary nature of IRC +10 216. The fact that the polarization pattern in the southeastern part of the nebula at 1.1 $`\mu `$m has a different orientation than in the rest of the nebula might be an indication of a second illuminating source.

Mechanisms inherent to the star include rotation (Dorfi & Höfner DorHoef96 (1996), Garcia-Segura et al. GarBer99 (1999)), non-radial pulsations (Soker & Harpaz SokHar92 (1992)) and magnetic fields (Pacoli et al. Pacetal92 (1992), Garcia-Segura et al. GarBer99 (1999)). Both non-radial $`p`$-modes and magnetic fields appear to only be important for significant rotation rates. Often spin-up agents due to binarity are assumed. For instance, Groenewegen (Groenewegen96 (1996)) favours non-radial pulsation or an as yet unidentified companion which spun up the central star as the most likely explanation for the non-spherical shape of the dust shell of IRC +10 216.

AGB stars are known to be slow rotators. Stars with initial masses below $`1.3`$M can be expected to lose almost their entire angular momentum during the main sequence phase due to magnetic braking operating in their convective envelopes. Consequently, they are not believed to devolop non-spherical mass-loss due to rotation. Stars with larger initial masses are spun down due to expansion and mass loss in the course of evolution, but may achieve sufficiently high rotation rates at the end of their AGB stage (Garcia-Segura et al. GarBer99 (1999)). Even small rotation rates influence dust-driven winds considerably, yielding a mass loss preferentially driven in the equatorial plane (Dorfi & Höfner DorHoef96 (1996)). Furthermore, for supergiants leaving the Hayashi line, Heger & Langer (HegLan98 (1998)) found that significant spin up of the surface layers may take place. Thus, on second glance, inherent rotation might be able to support axisymmetric mass loss during the transition to the proto-planetary nebula phase for more massive AGB stars as IRC +10 216.

## 4 Conclusions

We have presented high-resolution $`J`$–, $`H`$–, and $`K`$–band observations of IRC +10 216 with the highest resolution so far at $`H`$ of 70 mas. A series of $`K`$–band images from five epochs between October 1995 and November 1998 shows that the inner nebula is non-stationary. The separations of the four dominant resolved components increased within the 3 years by up to $`35\%`$. For the two brightest components a relative velocity within the plane of the sky of about 23 mas/yr or 14 km/s was found. Within these 3 years the rather faint components C and D became brighter whereas component B faded. The general geometry of the nebula seems to be bipolar.

We find that the most promising model to explain the structures and changes in the inner nebula is to assume that the star is at or near the position of component B. The star is then strongly but not totally obscured at $`H`$ and $`K`$. Consistently component B is very red in the $`HK`$ color while A and the northern $`J`$–band components are relatively blue. The polarization pattern with strong polarization in the northern arms and also a significant polarization in the peak supports this model. The inner nebula and the apparent motions seem to be rather symmetric around this position and the observed changes are consistent with the assumption of an enhanced mass loss becoming apparent at least in 1997.

IRC +10 216 is without doubt in a very advanced stage of its AGB evolution. The observed bipolarity of its dust shell even reveals that it has possibly entered the phase of transformation into a protoplanetary nebula.

###### Acknowledgements.

We thank the referee for valuable comments. This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France. The HST images for the polarization analysis have been retrieved from the Hubble Data Archive operated at the STScI, Baltimore, USA.


text:
# Thermal diffusion of sine-Gordon solitons

## 1 Introduction

As a key subject within nonlinear science, the dynamics of emergent, coherent structures (solitons, vortices, etc) has been a research topic that has attracted very much attention in the past quarter century scott . One question, extensively investigated in the literature bass ; kiv ; oldanx ; kvbook ; yuripr ; anx  is the following: Is, and if so, how is the motion and the shape of those excitations modified by the presence of small perturbations? Indeed, when applied to physical situations of interest, nonlinear models must incorporate additional terms, such as damping, constant or periodic external forces, or noise, to name a few. Among those, stochastic perturbations are very much of interest in view of their highly non trivial effects on nonlinear systems gsbook , and a great deal of work has been devoted to them bass ; oldanx ; kvbook . In particular, of the very many nonlinear models applied to physical problems, the sine-Gordon (sG) equation has been considered in much detail in this context, as it applies to, e.g., propagation of ultra-short optical pulses in resonant laser media lamb , a unitary theory of elementary particles skyrme1 ; skyrme2 ; enz ; raja , propagation of magnetic flux in Josephson junctions Barone , transmission of ferromagnetic waves feld , epitaxial growth of thin films CW ; krug ; us , motion of dislocations in crystals frenkel ; Nabarro , flux-line unlocking in type II superconductors nuevo , or DNA dynamics eng ; DNA ; yaku , situations in which noise (of different origins) can play, and often does, a crucial role. As an example, let us mention the recent work on long Josephson junctions reported in exp , where the authors calculated the escape rate from the zero-voltage state induced by thermal fluctuations, obtaining very satisfactory results compared with the experimental ones.

Specifically, this work is devoted to the study of the diffusive dynamics of sG kink solitons subjected to a thermal bath, as given by the stochastically perturbed, damped sG equation

$$\varphi _{tt}\varphi _{xx}+\mathrm{sin}(\varphi )=\alpha \varphi _t+f(x,t,\varphi _,\mathrm{}),$$
(1)
where $`\alpha \varphi _t`$ is a damping term with a dissipation coefficient $`\alpha `$, and $`f(x,t,\varphi _,\mathrm{})`$ is a thermal (gaussian) noise term fulfilling

$$\begin{array}{c}f(x,t,\varphi _,\mathrm{})=\sqrt{D}\eta (x,t),\eta (x,t)=0,\hfill \\ \\ \eta (x,t)\eta (x^{},t^{})=\delta (xx^{})\delta (tt^{}),\hfill \end{array}$$
(2)
where $`\sqrt{D}`$ is related to temperature through the fluctuation-dissipation theorem $`D=2\alpha k_bT`$, $`k_b`$ being the Boltzmann constant and $`T`$ the temperature.

To our knowledge, the first results on problems directly related to the one we deal with here were obtained by Joergensen et al. joer , who performed experiments on Josephson junctions and presented a derivation of the diffusion constant for kinks. Subsequently, Kaup and Osman kaup  studied, in a more rigorous way, the motion of damped sG kinks, driven by a constant force, in the presence of thermal fluctuations by using a singular perturbation expansion. They analyzed the temperature effect on the mean velocity of the kink and also the changes in the shape of the kink. In addition, they calculated the diffusion coefficient of the kink up to first-order in temperature and the energy values corresponding to the translational ($`E_T=k_bT/2`$) and radiational ($`E_R=k_bT`$) modes. These values of the energy have been also obtained by Marchesoni march , who applied the McLaughlin and Scott approach McL  to investigate kink motion under thermal fluctuations (see bass ; kiv ; kvbook  for reviews).

For the sake of completeness, let us mention work done along a different line, namely that devoted to the diffusive motion of the kink in equilibrium with phonons in isolated sG systems (possibly perturbed) yuripr ; Wsch ; March2 ; March . In this case, the kink diffusive motion is characterized by two diffusion coefficients. The first one of them is proportional to $`T^2`$ and is related to the anomalous diffusion, that arises from the phase shifts of kinks colliding with phonons and takes place on a short time scale in which the collision among kink and wave packet is elastic; the kink retains the same velocity after the collision (non-dissipative diffusion) and suffers only a spatial shift. However, for large times and in slightly perturbed sG systems, this interaction is nonlinear and becomes inelastic and the velocity of the kink changes after the collisions Wada . This diffusive regime is called viscous and has a diffusion coefficient proportional to $`T^1`$. The diffusion of the kink when the low energy excitations are represented by breathers has also been studied, and in Niko  it has been demonstrated that both descriptions (breathers or phonons) are equivalent and give rise to the same diffusion coefficient in the anomalous regime.

In any event, we want to stress that, although in this type of diffusion problem there are many open questions Wada , we will concern ourselves with the other kind of diffusion problem, in which the phonons appear as a consequence of an external heat bath, represented by white noise correlated in space and in time and the damping is included explicitly à la Langevin. The main aim of this work is to extend a previous study of ours about the overdamped limit of sine-Gordon kink diffusion nos  to the more physical and general case of the underdamped dynamics (i.e., with finite dissipation coefficient). As we will see below, the general perturbative approach raja  we resorted to in the overdamped case can also be applied, albeit with more difficulties, to the underdamped problem. The corresponding theoretical analysis is presented in Sec. 2, where we obtain explicit expressions for the long-time diffusive dynamics of kinks up to second-order in temperature, thus going beyond the currently available knowledge. The accuracy and importance of the new terms is assessed by numerical simulations in Sec. 3: we will see there that the quadratic corrections are in good agreement with the simulations and, most importantly, that they must be taken into account even for not so large temperatures. Finally, in the conclusions we summarize the main results of this work, comparing the underdamped and overdamped dynamics of the sG equation and discussing other related questions.

## 2 Analytical results

We begin by briefly reviewing the basic results we need for our analytical approach. We will concern ourselves with the perturbation effect on the kink solutions of the unperturbed sG equation, whose static form is

$$\varphi _0(x,t)=4\text{arctan}[\mathrm{exp}(x)].$$
(3)
Small perturbations over this equation can be treated by calculating the spectrum of linear excitations around the kink solution bis : To this end, we write

$$\varphi (x,t)=\varphi _0(x)+\psi (x,t),\psi (x,t)<<\varphi _0(x),$$
(4)
substitute in (1) (with $`\alpha =D=0`$) and linearize around $`\varphi _0(x)`$, arriving at the following equation for $`\psi (x,t)`$:

$$\psi _{tt}=\psi _{xx}\left[1\frac{2}{\mathrm{cosh}^2(x)}\right]\psi .$$
(5)
Then, assuming that the solution of (5) has the form

$$\psi (x,t)=f(x)\mathrm{exp}\left(i\omega t\right)$$
(6)
we find the eigenvalue problem for $`f(x)`$,

$$\frac{^2f}{x^2}+\left[1\frac{2}{\mathrm{cosh}^2(x)}\right]f=\omega ^2f.$$
(7)
This equation admits the following eigenfunctions with their respective eigenvalues

$`f_T(x)={\displaystyle \frac{2}{\mathrm{cosh}(x)}},\omega _T^2`$ $`=`$ $`0,`$ (8)
$`f_k(x)={\displaystyle \frac{\mathrm{exp}(ikx)[k+i\text{tanh}(x)]}{\sqrt{2\pi }\omega _k}},\omega _k^2`$ $`=`$ $`1+k^2,`$ (9)
which represent, respectively, the translation (Goldstone) mode and the radiation modes. Importantly, the functions $`f_T(x)`$ and $`f_k(x)`$ form a complete set with the orthogonality relations

$`{\displaystyle _{\mathrm{}}^+\mathrm{}}f_T^2(x)𝑑x`$ $`=`$ $`8,`$ (10)
$`{\displaystyle _{\mathrm{}}^+\mathrm{}}f_T(x)f_k(x)𝑑x`$ $`=`$ $`0,`$ (11)
$`{\displaystyle _{\mathrm{}}^+\mathrm{}}f_k(x)f_k^{}^{}(x)𝑑x`$ $`=`$ $`\delta (kk^{}).`$ (12)
We can now proceed with our problem: In order to tackle Eq. (1), with noise as given in (2), we use the same Ansatz proposed for the overdamped case in nos  (or for the general Klein-Gordon system in raja ): We assume that the solution of Eq. (1) is

$$\varphi (x,t)=\varphi _0[xX(t)]+_{\mathrm{}}^+\mathrm{}𝑑kA_k(t)f_k[xX(t)],$$
(13)
where $`X(t)`$ is the kink position. We now insert (13) in (1) and use the orthogonality of $`f_k`$ and $`f_T`$ bis , obtaining the following system of differential equation for $`X(t)`$ and $`A_k(t)`$:

$`\ddot{X}(t)+\alpha \dot{X}(t)`$ $`=`$ $`{\displaystyle \frac{\alpha }{8}}\dot{X}(t){\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑kA_k(t)I_1(k){\displaystyle \frac{1}{16}}{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑k{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑k^{}A_k(t)A_k^{}(t)R_3(k,k^{})+`$ (14)
$`+`$ $`{\displaystyle \frac{\sqrt{D}}{8}}{\displaystyle _{\mathrm{}}^+\mathrm{}}f_T[xX(t)]\eta (x,t)𝑑x`$
$``$ $`{\displaystyle \frac{1}{48}}{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑k{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑k_1{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑k_2A_k(t)A_{k_1}(t)A_{k_2}(t)R_6(k,k_1,k_2)`$
$``$ $`{\displaystyle \frac{\dot{X}(t)}{4}}{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑k{\displaystyle \frac{A_k}{t}}I_1(k){\displaystyle \frac{\ddot{X}(t)}{8}}{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑kA_k(t)I_1(k)+`$
$`+`$ $`{\displaystyle \frac{\dot{X}^2(t)}{8}}{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑k{\displaystyle \frac{A_k}{t}}I_2(k),`$
$`{\displaystyle \frac{^2A_k}{t^2}}+\alpha {\displaystyle \frac{A_k}{t}}+\omega _k^2A_k(t)`$ $`=`$ $`\alpha \dot{X}(t){\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑k^{}A_k^{}(t)I_3(k^{},k)+`$ (15)
$`+`$ $`{\displaystyle \frac{1}{2}}{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑k{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑k^{}A_k(t)A_k^{}(t)R_4(k,k^{})`$
$``$ $`\sqrt{D}{\displaystyle _{\mathrm{}}^+\mathrm{}}f_k^{}[xX(t)]\eta (x,t)𝑑x+`$
$`+`$ $`{\displaystyle \frac{1}{6}}{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑k^{}{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑k_1{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑k_2A_k^{}(t)A_{k_1}(t)A_{k_2}(t)R_7(k^{},k,k_1,k_2)+`$
$`+`$ $`2\dot{X}(t){\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑k^{}{\displaystyle \frac{A_k^{}}{t}}I_3(k^{},k)+`$
$`+`$ $`\ddot{X}(t){\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑k^{}A_k^{}(t)I_3(k^{},k)+\dot{X}^2(t)I_1(k),`$
where

$`I_1(k)`$ $`=`$ $`{\displaystyle _{\mathrm{}}^+\mathrm{}}{\displaystyle \frac{f_k}{\theta }}f_T(\theta )𝑑\theta ={\displaystyle \frac{i\pi \omega _k}{\sqrt{2\pi }\mathrm{cosh}\left({\displaystyle \frac{\pi k}{2}}\right)}},`$
$`I_2(k)`$ $`=`$ $`{\displaystyle _{\mathrm{}}^+\mathrm{}}{\displaystyle \frac{^2f_k}{\theta ^2}}f_T(\theta )𝑑\theta ,`$
$`R_3(k,k^{})`$ $`=`$ $`{\displaystyle _{\mathrm{}}^+\mathrm{}}f_T(\theta ){\displaystyle \frac{f_T}{\theta }}f_k(\theta )f_k^{}^{}(\theta )𝑑\theta ={\displaystyle \frac{i(\omega _k^2\omega _k^{}^2)^2}{4\omega _k\omega _k^{}\mathrm{sinh}\left({\displaystyle \frac{\pi \mathrm{\Delta }k}{2}}\right)}},\mathrm{\Delta }k=k^{}k,`$
$`I_3(k,k^{})`$ $`=`$ $`{\displaystyle _{\mathrm{}}^+\mathrm{}}{\displaystyle \frac{f_k}{\theta }}f_k^{}^{}(\theta )𝑑\theta ,`$
$`R_4(k,k^{})`$ $`=`$ $`{\displaystyle _{\mathrm{}}^+\mathrm{}}[f_k^{}^{}(\theta )]^2{\displaystyle \frac{f_T}{\theta }}f_k(\theta )𝑑\theta ,R_4(k,k)={\displaystyle \frac{3i\omega _k}{8\sqrt{2\pi }\mathrm{cosh}\left({\displaystyle \frac{\pi k}{2}}\right)}},`$
$`R_6(k,k_1,k_2)`$ $`=`$ $`{\displaystyle _{\mathrm{}}^+\mathrm{}}{\displaystyle \frac{^2f_T}{\theta ^2}}f_k(\theta )f_{k_1}^{}(\theta )f_{k_2}(\theta )𝑑\theta ,`$
$`R_7(k,k^{},k_1,k_2)`$ $`=`$ $`{\displaystyle _{\mathrm{}}^+\mathrm{}}\mathrm{cos}(\varphi _0)f_k^{}^{}(\theta )f_k(\theta )f_{k_1}^{}(\theta )f_{k_2}(\theta )𝑑\theta .`$ (16)
It goes without saying that these equations can not be solved. Therefore, in order to extract information from them, we resort to a perturbative approach assuming the noise term is small, or equivalently, that the temperature and the dissipation are not too large (this is not a serious restriction since our single-kink approach does not apply to high temperatures, when kink-antikink pairs are thermally generated butt ). We then expand $`X(t)`$ and $`A_k(t)`$ in powers of $`\sqrt{D}`$, i.e., $`X(t)=_{n=1}^{\mathrm{}}(\sqrt{D})^nX_n(t)`$ and $`A_k(t)=_{n=1}^{\mathrm{}}(\sqrt{D})^nA_k^n(t)`$, since when $`\sqrt{D}=0`$ and $`\alpha =0`$ we recover the static kink solution (in this case initially centered at the origin) of the sG equation. By substituting these expansions in (14) and (15) we find a set of linear equations for the coefficients of these series. The first members of this hierarchy correspond to order $`O(\sqrt{D})`$:

$$\ddot{X}_1(t)+\alpha \dot{X}_1(t)=\frac{1}{8}_{\mathrm{}}^+\mathrm{}f_T[xX(t)]\eta (x,t)𝑑xϵ_1(t),$$
(17)
from where we obtain the statistical properties of $`ϵ_1(t)`$,

$`ϵ_1(t)=0,ϵ_1(t)ϵ_1(t^{})={\displaystyle \frac{1}{8}}\delta (tt^{}),`$ (18)
and

$$\frac{^2A_k^1}{t^2}(t)+\alpha \frac{A_k^1}{t}(t)+\omega _k^2A_k^1(t)=_{\mathrm{}}^+\mathrm{}f_k^{}[xX(t)]\eta (x,t)𝑑x\xi _k(t),$$
(19)
which in turn leads to

$`\xi _k(t)=0,\xi _k(t)\xi _k^{}(t^{})=\delta (tt^{})\delta (kk^{}).`$ (20)
Equations (17)–(20) have been obtained in sal  by using a similar, but more restrictive perturbative approach. By integrating these two equations we obtain the first-order terms, $`X_1(t)`$ and $`A_k^1`$:

$$\begin{array}{c}X_1(t)=_0^te^{\alpha t^{}}_0^t^{}e^{\alpha \tau }ϵ_1(\tau )𝑑\tau 𝑑t^{},A_k^1(t)=e^{{\scriptscriptstyle \frac{\alpha t}{2}}}\left\{C_1(t)\mathrm{sin}\omega t+C_2(t)\mathrm{cos}\omega t\right\},\hfill \\ \\ C_1(t)=\frac{1}{\omega }_0^t\xi _k(\tau )e^{{\scriptscriptstyle \frac{\alpha \tau }{2}}}\mathrm{cos}\omega \tau d\tau ,C_2(t)=\frac{1}{\omega }_0^t\xi _k(\tau )e^{{\scriptscriptstyle \frac{\alpha \tau }{2}}}\mathrm{sin}\omega \tau d\tau ,\hfill \end{array}$$
(21)
where $`\omega ^2=\omega _k^2(\alpha ^2/4)`$. From these relations we can calculate the mean values and correlation functions up to first order in $`\sqrt{D}`$:

$`X_1(t)=0,X(t)X(t^{})`$ $`=`$ $`DX_1(t)X_1(t^{})={\displaystyle \frac{D}{16\alpha ^3}}[e^{\alpha M}e^{\alpha |\mathrm{\Delta }t|}+e^{\alpha M\alpha |\mathrm{\Delta }t|}`$ (22)
$``$ $`e^{\alpha (t+t^{})}+e^{\alpha t}+e^{\alpha t^{}}+2(\alpha M1)],`$
$`\dot{X}_1(t)`$ $`=`$ $`0,\dot{X}(t)\dot{X}(t^{})=D\dot{X}_1(t)\dot{X}_1(t^{})={\displaystyle \frac{D}{16\alpha }}\left[e^{\alpha |\mathrm{\Delta }t|}e^{\alpha (t+t^{})}\right],`$ (23)
$`A_k^1(t)=0,A_k(t)A_k(t^{})`$ $`=`$ $`DA_k^1(t)A_k^1(t^{})={\displaystyle \frac{D}{\omega ^2}}e^{\alpha (t+t^{})/2}[{\displaystyle \frac{e^{\alpha M}1}{2\alpha }}\mathrm{cos}\omega \mathrm{\Delta }t`$ (24)
$``$ $`{\displaystyle \frac{\alpha e^{\alpha M}}{8\omega _k^2}}\mathrm{cos}\omega \mathrm{\Delta }t{\displaystyle \frac{\omega e^{\alpha M}}{4\omega _k^2}}\mathrm{sin}\omega |\mathrm{\Delta }t|+{\displaystyle \frac{\alpha }{8\omega _k^2}}\mathrm{cos}\omega (t+t^{})`$
$``$ $`{\displaystyle \frac{\omega }{4\omega _k^2}}\mathrm{sin}\omega (t+t^{})],`$
where $`\mathrm{\Delta }t=tt^{}`$, and $`M=\mathrm{min}(t,t^{})`$. Of course, for $`t^{}=t`$ in Eq. (22) we recover the result in kaup  for $`[X(t)]^2`$.

We now turn to the main point of our work, namely obtaining the next-order corrections for the position and the velocity of the center of the kink. This requires the calculation of the next two contributions to $`X(t)`$ as well as the second order in the radiation terms, which are:

$`O(D)`$

$`\ddot{X}_2(t)+\alpha \dot{X}_2(t)`$ $`=`$ $`ϵ_2(t),`$ (25)
$`ϵ_2(t)`$ $``$ $`{\displaystyle \frac{ϵ_1(t)}{8}}{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑kA_k^1(t)I_1(k){\displaystyle \frac{\dot{X}_1(t)}{4}}{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑k{\displaystyle \frac{A_k^1}{t}}(t)I_1(k)`$ (26)
$``$ $`{\displaystyle \frac{1}{16}}{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑k^{}{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑kA_k^1(t)A_k^{}^1(t)R_3(k,k^{}),`$
$`{\displaystyle \frac{^2A_k^2}{t^2}}(t)+\alpha {\displaystyle \frac{A_k^2}{t}}(t)+\omega _k^2A_k^2(t)`$ $`=`$ $`ϵ_1(t){\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑k^{}A_k^{}^1(t)I_3(k^{},k)+`$ (27)
$`+`$ $`{\displaystyle \frac{1}{2}}{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑k{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑k^{}A_k^1(t)A_k^{}^1(t)R_4(k,k^{})`$
$``$ $`\dot{X}_1^2(t)I_1(k)+2\dot{X}_1(t){\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑k^{}{\displaystyle \frac{A_k^{}^1}{t}}(t)I_3(k^{},k);`$
$`O([\sqrt{D}]^3)`$

$`\ddot{X}_3(t)+\alpha \dot{X}_3(t)`$ $`=`$ $`ϵ_3(t),`$ (28)
$`ϵ_3(t)`$ $``$ $`{\displaystyle \frac{ϵ_1(t)}{8}}{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑kA_k^2(t)I_1(k){\displaystyle \frac{ϵ_2(t)}{8}}{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑kA_k^1(t)I_1(k)`$ (29)
$``$ $`{\displaystyle \frac{1}{16}}{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑k{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑k^{}A_k^2(t)A_k^{}^1(t)R_3(k,k^{})`$
$``$ $`{\displaystyle \frac{1}{16}}{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑k{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑k^{}A_k^1(t)A_k^{}^2(t)R_3(k,k^{})`$
$``$ $`{\displaystyle \frac{1}{48}}{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑k{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑k_1{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑k_2A_k^1(t)A_{k_1}^1(t)A_{k_2}^1(t)R_6(k,k_1,k_2)+`$
$`+`$ $`{\displaystyle \frac{\dot{X}_1^2(t)}{8}}{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑k{\displaystyle \frac{A_k^1}{t}}(t)I_2(k){\displaystyle \frac{\dot{X}_1(t)}{4}}{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑k{\displaystyle \frac{A_k^2}{t}}(t)I_1(k)`$
$``$ $`{\displaystyle \frac{\dot{X}_2(t)}{4}}{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑k{\displaystyle \frac{A_k^1}{t}}(t)I_1(k).`$
Analogously to what we have done for the first-order corrections, from the solutions of Eqs. (25) and (28) we find that

$`X_2(t)`$ $`=`$ $`0,\dot{X}_2(t)=0,`$ (30)
$`X_3(t)`$ $`=`$ $`0,\dot{X}_3(t)=0.`$ (31)
As for higher moments, taking into account that the cross-correlation function of $`X_1(t)`$ and $`X_3(t^{})`$ is of the same order as $`X_2(t)X_2(t^{})`$, and also that $`X_1(t)X_2(t^{})=0`$ we obtain that

$`[X(t)]^2`$ $`=`$ $`D[X_1(t)]^2+`$ (32)
$`+`$ $`D^2([X_2(t)]^2+2X_1(t)X_3(t))+\mathrm{},`$
$`[\dot{X}(t)]^2`$ $`=`$ $`D[\dot{X}_1(t)]^2+`$ (33)
$`+`$ $`D^2([\dot{X}_2(t)]^2+2\dot{X}_1(t)\dot{X}_3(t))+\mathrm{}`$
The expressions for the functions $`[X_2(t)]^2`$, $`[\dot{X}_2(t)]^2`$, $`X_1(t)X_3(t)`$, and $`\dot{X}_1(t)\dot{X}_3(t)`$ can be obtained after a lengthy calculation, and are very cumbersome indeed. We therefore do not include them here. However, for large time ($`t>>1/\alpha `$) these relations can be simplified, yielding, as $`t\mathrm{}`$

$`[X(t)]^2`$ $`=`$ $`{\displaystyle \frac{k_bTt}{4\alpha }}\left\{1+{\displaystyle \frac{k_bT}{32}}\left(1+{\displaystyle \frac{9\sigma ^2}{4}}\right)\right\},`$ (34)
$`[\dot{X}(t)]^2`$ $`=`$ $`{\displaystyle \frac{k_bT}{8}}\left\{1+{\displaystyle \frac{3k_bT}{128}}\left(12+\sigma ^2\right)\right\},`$ (35)
with

$$\sigma =_{\mathrm{}}^+\mathrm{}\frac{dk}{\omega _k\mathrm{cosh}\left({\displaystyle \frac{\pi k}{2}}\right)}=1.62386.$$
(36)
To complete the characterization of the kink diffusion, we can now compute in a straightforward way the average value of the wave function $`\varphi (x,t)`$, defined as

$`\varphi (x,t)`$ $`=`$ $`\varphi _0[x\sqrt{D}X_1(t)]+O(D)=`$ (37)
$`=`$ $`{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑X_1p(X_1)\varphi _0[x\sqrt{D}X_1(t)],`$
where $`p(X_1)`$ is the probability distribution function for $`X_1`$. To find explicitly this function we note that, if we rewrite Eq. (17) as a system of two differential equations,

$`\dot{X}_1`$ $`=`$ $`V,`$
$`\dot{V}`$ $`=`$ $`\alpha V+ϵ_1(t),`$ (38)
the last equation represents an Ornstein-Uhlenbeck process for the velocity, and its distribution function is given by

$$p(V)=\sqrt{\frac{1}{2\pi V^2}}\mathrm{exp}\left(\frac{V^2}{2V^2}\right),$$
(39)
(see noise ). Subsequently, by integrating the first equation of (38), we obtain that $`X_1=_0^tV(\tau )𝑑\tau `$. Since $`V`$ has a Gaussian distribution function, $`X_1`$ has also a Gaussian distribution function, given by (recall that $`X_1(t)=0`$)

$$p(X_1)=\sqrt{\frac{1}{2\pi [X(t)]^2}}\mathrm{exp}\left(\frac{1}{2}\frac{X_1^2}{[X(t)]^2}\right),$$
(40)
where the first and second moments of $`X_1`$ were obtained before, see Eq. (22). With this result, the integral (37), can be evaluated numerically taking into account Eqs. (40) and (22). In the next section we will compare this result with the mean value of the wave function as obtained from simulations of the full partial differential equation (1).

## 3 Numerical simulations

In order to test the approximate theory developed in the previous section, we have simulated numerically Eq. (1) by using the Heun method maxra . In our simulations we begin with a kink, initially at rest, with free boundary conditions. For the damping coefficient we choose $`\alpha =0.1`$, which is not too small because from (34) we can see that $`[X(t)]^2`$ is proportional to $`1/\alpha `$. This means that if $`\alpha `$ is too small the kink can move in a much larger region, forcing us to increase the length of our simulated system in the simulations, already quite time consuming. Furthermore, the relation (34) is only valid for large times ($`t>>1/\alpha `$). Again, for too small $`\alpha `$ we would need to simulate our equation for very long times and, as $`[X(t)]^2`$ increases linearly with time \[see Eq. (34)\], the system length would once more have to be large. The other parameters are $`\mathrm{\Delta }x=0.2`$, $`\mathrm{\Delta }t=0.001`$ and the length of the system $`L=200`$. We have calculated all average values over 1000 realizations up to a final time $`400`$. It is important to point out that, this system being inertial, the accuracy of the averages is considerably less than for overdamped problems, this being the reason why we have to use such large ensembles of trajectories to obtain reasonably good results.

An important, nontrivial issue is the question as to how can we find the center of the kink. We solve this problem by finding all the discrete lattice points $`x_i`$ and $`x_{i+1}`$ such that $`\varphi _i\pi `$ and $`\varphi _{i+1}\pi `$ or vice versa, and then estimating the corresponding points $`\stackrel{~}{x}_i`$ where $`\varphi =\pi `$ by linear interpolation. Afterwards, among the $`n`$ such points $`\stackrel{~}{x}_n`$, we choose to be the center of the kink the value $`\stackrel{~}{x}=\stackrel{~}{x}_n`$, which minimizes $`_{i=1}^{L/\mathrm{\Delta }x}[\varphi _i(t)\varphi _0(x\stackrel{~}{x}_n,t)]^2`$, i.e., the discrete version of the integral of the square of the difference of $`\varphi `$ and $`\varphi _0`$. It has to be realized that this involves an assumption, namely that individual realizations of the kink have a shape similar to that of the unperturbed kink. As can be seen from Fig. 1, where the individual realizations are compared with the initial condition (an exact kink), this is indeed the case and our procedure is truly sensible. Therefore, we are sure that this method to compute the kink center avoids the false centers, which can appear for higher temperatures due to fluctuations introducing a systematic difference between the numerical and the theoretical results (see nos ). With the procedure we have just summarized, that works even for relatively large temperatures, we believe we find a very accurate approximation to the actual center of the kink. We will come back to Fig. 1 below.

As an example of the comparison of the numerical simulations of Eq. (1) with the theoretical results obtained in the previous section and valid for large times, Fig. 2 shows the numerically computed variance of the center of the kink, $`[X(t)]^2X(t)^2`$, as well as the first- and second-order analytical expressions. The plot clearly evidences that the numerical variance asymptotically coincides with the second-order expression: Note that to compare the different curves one has to look at the slopes at times $`t>>1/\alpha `$ (in this case, $`t100`$, for instance, as $`\alpha =0.1`$); the theoretical result is not valid at early times and therefore there is a bias between analytics and numerics coming from that. The small, irregular oscillations in the numerical curve arise from the difficulty in accurately computing averages in an underdamped system like this mentioned above; however, we believe that the present accuracy is enough to confirm the validity of our approach. We have observed the same agreement for other values of temperatures ($`k_bT=0.2,0.6,0.8`$, not shown). In all cases, we have computed the diffusion coefficient for large times as the slope of the variance of $`X(t)`$ again for $`t100`$, the regime in which we expect our analytical approximation to be valid. Summarizing our results, these numerical values of the diffusion coefficient are plotted in Fig. 3 together with the theoretical results. It is clear that for large temperatures the quadratic behavior in $`k_bT`$ of the diffusion coefficient becomes important. For higher values of the temperature, such as $`k_bT=0.8`$, the numerical value of the diffusion coefficient is not so close to the predicted one. This effect arises because of the large diffusivity of the kink in that range: Indeed, for this and higher temperatures the kink performs very long excursions away from the center, reaching the boundaries of the numerical integration interval; it is clear that when this occurs, the diffusion of the kink is not in free space anymore and hence those realizations spoil the quality of the averages. The way to solve this problem would be to resort to much larger numerical systems, but within our present computing capabilities this would necessitate a simultaneous decrease in the number of realizations in the average, leading again to poorer results. However, it is important to realize that this boundary effect leads to an underestimation of the diffusion coefficient (as the boundary prevents the kink from travelling as far as it should) and therefore the point in Fig. 3 for $`k_bT=0.8`$ is a lower bound for the diffusion coefficient, with the actual one lying even closer to our second order prediction.

Finally, there is one last question that deserves discussion, namely that of the physical significance of the mean value of the wave function $`\varphi `$. In Fig. 1 we can clearly see that, whereas individual realizations of kinks look very similar to the unperturbed ones, the mean value of $`\varphi `$ is a much wider excitation, not even close to the original kink. Figure 1 clearly shows that this does not mean that the width of individual kinks increases; indeed, much as we discussed regarding the overdamped problem nos , we have verified numerically that the mean wave function $`\varphi (x,t_{fix})`$ increases due to the variance of the kink center of individual realizations, and hence it should not be interpreted as the typical deformation of the shape of kinks: Indeed, the widening of the mean value of $`\varphi `$ arises from the contributions of the stochastically moving, but mostly undistorted kinks whose center positions have the distribution of a rigid, diffusing particle. To further check this interpretation, we can look at Fig. 4, where we have represented the mean value of the wave function for two fixed times $`t_{fix}=100,300`$, obtained from the numerical simulation of the full partial differential equation (1), for $`k_bT=0.6`$ and $`\alpha =0.1`$. The overimposed points, computed by using the Gaussian distribution function $`p(X_1)`$ \[Eq. (40)\] of the kink center $`X(t)=\sqrt{D}X_1`$ found in the last section, show the excellent agreement between our theory and the simulation. Of course, there is a small discrepancy that is likely to disappear if one would go to a next order calculation, but for the present purpose of understanding the mean wave function $`\varphi `$ the first order calculation is enough. In addition, we have plotted the initial kink (at rest) in order to see that the mean value of the wave function increases with time.

## 4 Discussion and conclusions

To summarize, in this work we have studied the diffusive dynamics of sine-Gordon kinks subjected to thermal fluctuations. We have analytically calculated expressions, valid up to second order in temperature, for the average position and variance of the kink center, as well as for the mean shape of the kink. We have numerically checked the validity of these results up to temperatures of the order of $`k_bT=0.8`$ (in dimensionless units, equivalent to about a 10% of the kink rest mass), already close to the temperature at which kink-antikink nucleation becomes a likely event. Therefore, our first conclusion is that the second-order theory developed here is the proper one, meaning it is accurate and higher order terms are negligible, to describe the thermal diffusion of sine-Gordon kinks in the single kink propagation regime. Interestingly, our calculation pinpoints the fact that the second-order correction in $`k_bT`$ comes from the interaction between kink and phonons. This implies that the physics behind this contribution is basically the same as for the case of anomalous diffusion in an isolated chain mentioned in the introduction yuripr ; Wsch ; March2 ; March . Note that we do not expect $`T^1`$ contributions in our analytical calculations, as they are carried out in a continuum sG equation yuripr  and, in any case, they would show up in simulations only for very low temperatures. Apart from that, it is also interesting to note that, according to Eq. (35), the second order term implies an increase of the energy carried by the kink beyond the $`k_bT/2`$ predicted by statistical mechanics (recall that the kink mass is 8 in our units). This can be interpreted in the following way: The kink is dressed by phonons which increase its mass. Thus, the kink energy would be $`M(T)\dot{X}^2/2`$, with a temperature dependent mass $`M(T)`$ whose expression can be easily found from Eq. (35). In order to confirm this interpretation, one could compute the energy carried by the phonons which dress the kink, but we believe that it is not necessary because, on the one hand, it would be a rather involved calculation (far beyond the scope of this work) and, on the other hand, we do not think that there is any other possible interpretation of this result.

A second relevant point of this study relates to the numerical simulation and center location procedures. As this is an underdamped (inertial) system, the thermal mobility of kinks is quite large, the larger the higher the temperature. Because of this, we have not been able to obtain very precise numerical averages at the top of the temperature range studied, since the lengths of the systems and the number of realizations required are very large and consequently time consuming. However, we believe that the results presented here are enough to verify our theory. This is reinforced by the very good agreement between analytics and numerics regarding the mean shape of the field, even for temperatures as large as $`k_bT=0.6`$ (see Fig. 3), which shows that our approach indeed captures the physics of the diffusion process. In addition, we want to emphasize that, to our knowledge, we have designed a new algorithm to detect the kink center which gives very good results even for the highest temperatures studied, where previous researches, such as nos , had found problems arising from the many false centers detected.

Another important issue is the comparison of the present analysis to that in nos  for the overdamped problem. We have found that the diffusion coefficient given by (34) for the present case practically coincides with that obtained in nos  for the overdamped limit of the equation: the difference in the second order is approximately $`0.06k_bT`$, i.e., very small compared to the magnitude of the quadratic contribution itself. Furthermore, the width of the mean value of the wave function increases with time for the overdamped case nos  in the same manner as that reported here. Therefore, we can conclude that for large times the dynamic of the underdamped sG kinks is very similar to the overdamped case. This is an important point, because in principle one can expect similar results for other kink-bearing systems such as the $`\varphi ^4`$ equation, for instance, whose overdamped diffusive dynamics is known (see dzi  for the $`\varphi ^4`$ case), thus avoiding the much more involved calculation of the underdamped case.

Finally, we want to mention the relevance of this work to experimental systems, such as long Josephson junctions. As has been shown in exp , the thermal sG equation (1) is a good description of the physics of in-line Josephson junctions (although different boundary conditions are needed in that case). The work in exp  compared the predictions from the sG model to experimentally measured escape rates from the zero voltage state. Therefore, it should be possible to design similar experiments in order to test our results and, specifically, the increased (quadratic) diffusivity of kinks at higher temperatures vs the linear behavior at lower ones. We hope that our theoretical work stimulates further experimental research in that direction.

## Acknowledgement

We are grateful to A. R. Bishop for his comments and Esteban Moro for his suggestions. Work at GISC (Leganés) has been supported by CICyT (Spain) grant MAT95-0325 and DGES (Spain) grant PB96-0119. Travel between Bayreuth and Madrid is supported by “Acciones Integradas Hispano-Alemanas”, a joint program of DAAD (Az. 314-AI) and DGES. This research is part of a project supported by NATO grant CRG 971090.