text:
# BPS Branes in Cosmology

## I Introduction

During the last few years there were many attempts to construct a consistent brane cosmology, see e.g.  and references therein. One of the most interesting possibilities is to use supersymmetric BPS branes in cosmology. Studies of this idea developed from the no-go theorems for nonsingular supersymmetric domain walls  to the construction of supergravity in singular spaces  where the bulk and brane actions are supersymmetric. Investigation of the BPS brane cosmology, i.e. the theory of interacting and moving near-BPS branes, has brought an additional level of complexity, both on the technical and on the conceptual level.

One of the most challenging recent attempts to construct a consistent cosmology based on a picture of colliding BPS branes is ekpyrotic scenario . In this paper we will discuss some of the general issues of the BPS brane cosmology by critically analyzing the ekpyrotic scenario.

It was claimed in  that the ekpyrotic scenario is based on the Hor̆ava-Witten (HW) phenomenology, and it solves all major cosmological problems without using inflation . However, in  it was argued that these claims were overly optimistic. First of all, the standard HW phenomenology (including most of its versions with non-standard embedding) is based on the assumption that we live on the positive tension brane . Meanwhile, the model proposed in  was based on an unconventional approach to HW theory, assuming that we live on the negative tension brane. This required a substantial reformulation of HW phenomenology . In particular, one must revise the standard assumption  that the gauge coupling on the hidden brane is large. In  a different version of this scenario was proposed, assuming that we live on the positive tension brane, in accordance with . It was called pyrotechnic scenario.In the pyrotechnic scenario, unlike in the ekpyrotic scenario, we do not make any attempts to avoid inflation. In fact we have argued in  that avoiding inflation in this scenario may require additional fine-tuning. We will discuss this in Section II of our paper and explain that the relevant issue is not the standard versus non-standard embedding, but Hor̆ava-Witten phenomenology  versus Benakli-Lalak-Pokorski-Thomas  phenomenology.

Other concerns include extreme fine-tuning of the brane potential $`V(Y)`$ required in the ekpyrotic scenario. In particular, the potential must be extremely small (suppressed by a factor at least as small as $`10^{50}`$) near the hidden brane. This makes it very difficult to understand how this scenario could be made consistent with the brane stabilization . To solve the homogeneity problem in this scenario, one would need the branes from the very beginning to be parallel to each other with an accuracy better than $`10^{60}`$ on a scale $`10^{30}`$ times greater than the distance between the branes. In our opinion, these problems, as well as several other problems of the ekpyrotic scenario pointed out in , remain unresolved.

In this paper we would like to consider some other aspects of the ekpyrotic scenario. First of all, it was claimed in  that the action and the cosmological solution describing three static branes in the ekpyrotic scenario have been obtained in . However, the 5d bulk & brane action in  was found for the case of the two boundary branes only. We will show that it cannot be generalized to the case of boundary and bulk branes using the formalism of Ref. .

To add the bulk brane to this construction one must use the 4-form field, which was introduced in the context of 5d supergravity in singular spaces by Bergshoeff, Kallosh and Van Proeyen  and recently generalized in . We will show that the corresponding part of the 5d action and the part of the solution given in  are not quite correct as far as coefficients are concerned; we will present the corrected action and static BPS solution describing boundary and bulk brane in Section III (some technical details are given also in Appendix).

In Section IV we will discuss an effective 4d description of the 5d cosmology and give an argument that the 4d space in the ekpyrotic scenario can only collapse.Various versions of this argument were suggested to us independently by T. Banks, G. Dvali, and J. Maldacena . According to , this is indeed the case. The scale factor of the universe, which was obtained after integrating the 5d action over $`y`$ (the 5th dimension) and solving equations in the effective 4d theory, decreases. The authors of  argued that one can go back from the effective 4d theory to 5d theory and show that the scale factor of the visible brane grows. We believe, however, that in order to go back from the effective 4d description to the 5d theory one should perform an explicit investigation of the 5d geometry before the integrating over the 5th dimension.

This is not an easy task since the bulk brane potential $`V(Y)`$ driving expansion of the universe in the ekpyrotic scenario was added in  by hand to the 4d formulation of the theory, rather than to the 5d theory. One can only speculate how one should interpret this term from the point of view of the 5d theory.

In  this term was interpreted as a correction to the bulk brane action, which means that it represents an effective delta-functional addition to the total energy concentrated on the bulk brane. We will show in Sect. V that in this case the time-dependent ansatz for the metric and fields used in  does not provide a consistent solution to the 5d gravitational equations and to the equation for the dilaton field. In such a situation it is dangerous to take averages of the 5d action or of the 5d equations over $`y`$. If one does so with the 5d dilaton equation, one finds that the whole universe, including the visible brane, should exponentially collapse rather than expand.

In order to find a consistent solution of 5d equations one needs to make at least two important modifications discussed in Section VI. First of all, when one considers moving branes, or branes having energy momentum tensor different from the 4d cosmological constant, one should use a general ansatz for the metric and for the fields which respects the (planar) symmetry in the problem  rather than the factorized ansatz used in . The typical situation there when the junction conditions on the branes are taken into account is that the 5d metric has a non-factorizable dependence on time and fifth direction.

Moreover, we believe that the brane potential $`V(Y)`$ should be interpreted as the “radion potential.” It represents the total energy of the configuration involving the bulk brane standing at the distance $`Y`$ from the visible brane . In such a situation it may be incorrect to represent $`V(Y)`$ as a delta-functional contribution to the energy density localized on the bulk brane. Instead of that, one should find out the distribution of the fields responsible for the emergence of the long range interaction described by the potential $`V(Y)`$, and substitute their $`y`$-dependent energy-momentum tensor (rather than $`V(Y)`$) into the 5d equations. A simplest example of a similar situation is given by the energy of an electric capacitor. The energy of interaction of the charged plates of a capacitor $`V(R)`$ is proportional to the distance $`R`$ between the plates. However the electrostatic energy density $`E^2`$ is concentrated not at the plates but in the bulk.

We conclude that in addition to many other problems of the ekpyrotic scenario discussed in , there exists another one. The ansatz for the metric and fields used in  does not provide a correct solution to the dilaton equation and the 5d Einstein equations.

## II Brane tension and HW phenomenology

According to the Hor̆ava-Witten theory , the universe consists of two branes in 5d space, which appeared after 6 dimensions of the 11d space were compactified on Calabi-Yau (CY) space. The unification of weak, strong and electromagnetic interactions is achieved on the visible brane, with the gauge coupling $`\alpha _{GUT}0.04`$. The standard Hor̆ava-Witten phenomenology , including its versions with non-standard embedding , is based on the assumption that we live on the positive tension brane, and the volume of the CY space decreases towards the hidden brane, whereas the gauge coupling constant increases. This leads to the strong gauge coupling on the hidden brane, $`\alpha _{\mathrm{hid}}=O(1)`$. In the strong coupling regime one can obtain gaugino condensation, which plays an important role in the HW phenomenology.

The ratio of the gauge couplings is inversely proportional to the ratio of Calabi-Yau volumes at the positions of the branes ,

$$\alpha _{\mathrm{hid}}=\alpha _{\mathrm{GUT}}\frac{V_{\mathrm{vis}}}{V_{\mathrm{hid}}}.$$
(1)
The visible brane is located at $`y=0`$ and the hidden one at $`y=R`$ so that $`V_{\mathrm{vis}}=e^{\varphi (0)}`$ and $`V_{\mathrm{hid}}=e^{\varphi (R)}`$. This formula is valid for the versions of HW theory with the standard embedding without M5-branes, as well as for the versions with the non-standard embedding with M5-branes present in the bulk, see e.g. a discussion of the phenomenology of the theories with the non-standard embedding by Lukas, Ovrut and Waldram .

In  it was assumed that the tension of the visible brane $`\alpha `$ is negative. The volume of the CY space in  is proportional to $`D^3(y)`$, where $`D(y)=C+\alpha y`$. The function $`D(y)`$ grows at large $`y`$. As a result, the value of $`D^3(y)`$ near the hidden brane is 50 times greater than near the visible brane for the parameters used in ; it is 27 times greater for the parameters chosen in the replaced version of their paper. The authors of the ekpyrotic scenario did not calculate the value of $`\alpha _{\mathrm{hid}}`$ in their scenario , so we will estimate it now. If one has $`\alpha _{GUT}0.04`$ on the visible brane, one finds $`\alpha _{\mathrm{hid}}=O(10^3)`$. This is way too small to lead to the usual HW phenomenology, which requires $`\alpha _{\mathrm{hid}}1`$ and gaugino condensation . We are not saying that a consistent phenomenology with $`\alpha _{\mathrm{hid}}\alpha _{GUT}`$ is impossible, see e.g. , but this is a rather unconventional possibility.

Thus, we believe that the original version of the ekpyrotic scenario was at odds with the standard HW phenomenology as defined in . So why was the tension of the visible brane chosen to be negative in , despite all complications associated with such a choice? In the original version of Ref.  we read: As we will see in Section VB, it will be necessary for the visible brane to be in the small-volume region of space-time. This statement, as well as the related conclusion that the spectrum of perturbations in the ekpyrotic scenario is blue, was emphasized in many places of the text. The reason of this statement, as explained in Section VB of , was rooted in the idea that the decrease of $`D(Y)`$ is required for generation of density perturbations in the ekpyrotic scenario.

In ref.  a simple description of generation of density perturbations in the ekpyrotic scenario was presented§§§ Recently it was claimed  that if one takes into account gravitational backreaction, no perturbations of metric are produced in the ekpyrotic scenario. If this is the case, there is no need to continue discussion of this scenario. However, it is not obvious to us that the absence of fluctuations of the effective 4d metric discussed in  implies the absence of the 5d metric perturbations and the absence of the time delay of the brane collision. Perturbations investigated in  appeared not because of the fluctuations of the curvature of the 4d spacetime prior to the brane collision, but because of the “radion” perturbations $`\delta Y`$ (related to $`\delta g_{yy}`$). These perturbations, describing an inhomogeneous embedding of the bulk brane in the 5d spacetime, cannot be reduced to perturbations of the 4d geometry. In any case, before studying perturbations, one should first carefully examine the behavior of the nonperturbed metric. This is what we are going to do in our paper. using the methods developed in the theory of tachyonic preheating . It was shown in , in particular, that the requirement that the visible brane must be in the small-volume region of space-time is not necessary. Thus there is no reason to abandon the standard HW phenomenology and assume that we live on the negative tension brane. An improved version of the ekpyrotic scenario based on the assumption that we live on the positive tension brane was called “pyrotechnic universe” .

## III  M/string theory and 5d BPS domain walls

Compactification of M/string theory with extended objects down to 5d sometimes leads to appearance of supersymmetric domain walls. It has been explained in  that the supersymmetric domain walls in 5d must be charged. The relevant 4-form $`A`$ and 5-form field strength $`F=dA`$ play an important role in the supersymmetry of the bulk & brane construction of . It is the consequence of the M/string theory supersymmetry where RR-fields and M-theory form-field provide the balance of forces between BPS extended objects.

Ekpyrotic scenario is based on 11d theory with Hor̆ava-Witten domain walls and some M5-branes between them. A compactification of this system may lead to 5d theory with charged 3-branes. A complete M-theory derivation of the 5d theory with boundary HW domain walls and bulk branes between them was not actually worked out. Only the part of it with 2 HW walls was reduced to 5d in . The 4-form field describing HW domain walls and compactified M5-branes was introduced recently in  and it was argued in  that the 4-form formulation of the action in  is equivalent to the action presented in . It was claimed that this is easily seen by eliminating the 4-form using its equation of motion.

Since the purpose of these notes is to find the correct equations in 5d which are related to M/string theory and BPS construction, we have examined this claim and found that it is not quite correct. We use in this analysis the supersymmetric 5d theory with the 4-form field introduced in .

* In  there are only 2 HW walls, no branes in between. We will show that in case when in addition to the 2 boundary branes there are bulk branes present, the covariant action with the bulk dilaton potential (replacing the 5-form contribution) of the type  does not exist. Only the formulation with the 4-form field can describe the two boundary branes and a brane in between. We refer the reader to the Appendix A where all technical details are explained.
* In case of only 2 HW walls the action in  is not equivalent to the action in  unless the factor in front of the $`^2`$ kinetic term is changed.
* The expression for the 4-form for the static solution in  is not correct. The numerical factor has to be modified and the sign has to be changed to make it a BPS solution.

We start with the corrected version of the action in  (note the factor $`3/2`$ in front of $`^2`$)

$`S={\displaystyle \frac{M_5^3}{2}}{\displaystyle __5}d^5x\sqrt{g}\left(R{\displaystyle \frac{1}{2}}(\varphi )^2{\displaystyle \frac{3}{2}}{\displaystyle \frac{e^{2\varphi }_5^2}{5!}}\right)`$ (2)
$`3{\displaystyle \underset{i=1}{\overset{3}{}}}\alpha _iM_5^3{\displaystyle _{_4^{(i)}}}d^4\xi _{(i)}\left(\sqrt{h_{(i)}}e^\varphi {\displaystyle \frac{ϵ^{\mu \nu \kappa \lambda }}{4!}}𝒜_{\gamma \delta ϵ\zeta }_\mu X_{(i)}^\gamma _\nu X_{(i)}^\delta _\kappa X_{(i)}^ϵ_\lambda X_{(i)}^\zeta \right),`$ (3)
where $`\alpha _1=\alpha `$, $`\alpha _2=\alpha \beta `$, $`\alpha _3=\beta `$. The corrected form of the static solution isHere $`y`$ is a point of $`S^1/Z_2`$, i.e. $`R<yR`$, with $`0`$ and $`R`$ as fixed points, $`R`$ identified with $`R`$. This explains the factor of 2 for the fixed-point brane sources at $`0,R`$ (accounting for their images), and the presence of the brane in the bulk at $`y=Y`$ and its image at $`y=Y`$.

$`ds^2=D(y)(N^2dt^2+A^2d\stackrel{}{x}^2)+B^2D^4(y)dy^2,`$ (4)
$`e^\varphi =BD^3(y),`$ (5)
$`D(y)=\alpha y+C\mathrm{for}y<Y,`$ (6)
$`=(\alpha \beta )y+C+\beta Y\mathrm{for}y>Y,`$ (7)
where $`A,B,C,N`$ are constants, $`C>0`$,

$$𝒜_{0123}=+A^3NB^1D^1(y),_{0123y}=A^3NB^1D^2(y)D^{}(y),$$
(8)
and

$$[D(y)]^{^{\prime \prime }}=2[\alpha \delta (y)(\alpha \beta )\delta (yR)]\beta \delta (yY)\beta \delta (y+Y).$$
(9)
The factor $`A^3NB^1`$ in $``$ in eq. (8) was absent in . The sign of $`𝒜_{0123}`$ and $`_{0123y}`$ (which differs from the one in ), is easily checked by observing that the force on a static probe brane parallel to the source branes must vanish. The world-volume term cancels against the Wess-Zumino term in the static probe brane action. We are assuming the standard convention $`ϵ^{0123}=+1`$.

The static BPS solution is valid for branes that are not moving. It may serve as a starting point for finding time-dependent cosmological solutions.

## IV  4d view on ekpyrotic universe

Before investigating time-dependent cosmological solutions in 5d, let us see what could be expected from the point of view of the effective 4d theory. Indeed, if one considers the situation when the distance between the branes is very small and their motion is slow, one could expect that in the first approximation it should be possible to describe low energy theory from the point of view of 4d Einstein gravity (or Brans-Dicke theory) coupled to matter. Deviations from this description could occur if there were some processes with the energies comparable to the inverse distance between the branes. However, in the ekpyrotic scenario all energy scales (reheating temperature, Hubble constant, etc.) are several orders of magnitude smaller than the inverse distance between the branes $`1/R`$.

Therefore one may expect that the ekpyrotic scenario can be described entirely in terms of the 4d theory. But in such a case the universe, which was static in the beginning of the process, can only collapse. We present here, in a slightly modified form, the basic argument of Ref. ; see also .

Let us write down the Einstein equations for a homogeneous flat universe. The first equation is

$$H^2=\frac{8\pi G}{3}\rho .$$
(10)
The second equation, which follows from the first one and the energy conservation, can be represented in the following form:

$$\dot{H}=4\pi G(\rho +p).$$
(11)
Here $`\rho `$ and $`p`$ are the density and pressure in the effective 4d theory, and $`H=\dot{a}/a`$, where $`a`$ is the scale factor in the 4d space on the visible brane.

In the beginning, the branes do not expand, $`H=0`$, and the effective energy density and pressure vanish for the static brane configuration considered in the previous section.

This situation changes when one adds by hand the potential energy $`V(Y)`$ associated with the position of the bulk brane. According to , $`V(Y)<0`$. This, however, would be inconsistent with Eq. (10) unless one assumes that the bulk brane has nonzero velocity from the very beginning, so that the total energy density is non-negative. We do not want to speculate on how this configuration could emerge; see Ref.  for a discussion of the problem of initial conditions in the ekpyrotic scenario.

What is more important, Eq. (10) implies that $`\dot{H}0`$ because $`\rho +p0`$ in accordance with the null energy condition. Thus, if the universe begins in a static state, $`H=0`$, then it can only collapse, since $`H=\dot{a}/a0`$. Therefore one may argue that the ekpyrotic scenario cannot describe an expanding universe .

Let us compare these expectations with the results obtained in . To describe the motion of the bulk brane in the ekpyrotic scenario the authors started with the factorized ansatz based on (4)–(7) where it was assumed that some of the parameters of the static solution become functions of time but not of $`y`$ ($`A,N,YA(t),N(t),Y(t)`$), whereas some other parameters remain constant ($`\dot{B}=\dot{C}=0`$). They substituted this modified ansatz into the action (3), and integrated over $`y`$. In this way they obtained the 4d “moduli space” action with the Lagrangian density $`=_{bulk}^{4d}+_\beta ,`$ where

$$_{bulk}^{4d}=2\frac{3A^3BM_5^3}{N}_0^R𝑑yD^3(y,Y)\left[\left(\frac{\dot{A}}{A}\right)^2+3\frac{\dot{A}}{A}\frac{\dot{D}}{D}+\frac{1}{2}\left(\frac{\dot{D}}{D}\right)^2\right],$$
(12)
and

$$_\beta =\frac{3\beta M_5^3A^3B}{N}\left[\frac{1}{2}D^2(Y)\dot{Y}^2\right].$$
(13)
We gave here a form of the Lagrangian in which it is clear that it is integrated over the full range of the 5th coordinate.

In  the notation $`I_k(Y)_R^R𝑑yD^k(y,Y)=2_0^R𝑑yD^k(y,Y)`$ ($`k=1,2,\mathrm{}`$) and the explicit form of $`D`$ was used so that $`\dot{D}=0`$ at $`|y|<Y`$ and $`\dot{D}=\beta \dot{Y}`$ at $`|y|Y`$. After obtaining the effective 4d action by integration over $`y`$, the authors added by hand an important term, which plays a crucial role in their scenario:

$$\mathrm{\Delta }_\beta =3\beta M_5^3A^3BNV(Y).$$
(14)
The effective potential $`V(Y)`$ may appear, e.g., as a result of the nonperturbative effects associated with open M2-brane instantons . It was assumed in  that $`V(Y)e^{\alpha mY}`$ at large $`Y`$, and that it vanishes at $`Y=0`$. Here $`m`$ is some positive numerical constant specified (together with the tension $`\alpha `$) in .

The next step was to replace $`A`$ and $`N`$ by $`a=A(BI_3M_5)^{1/2}`$ and $`n=N(BI_3M_5)^{1/2}`$. This gives the effective 4d theory in the following form:

$`={\displaystyle \frac{3a^3M_5^2}{n}}\left\{\left({\displaystyle \frac{\dot{a}}{a}}\right)^2+J(Y)\beta ^2\dot{Y}^2+{\displaystyle \frac{\beta }{I_3}}\left[{\displaystyle \frac{1}{2}}D(Y)^2\dot{Y}^2n^2{\displaystyle \frac{V(Y)}{BI_3M_5}}\right]\right\},`$ (15)
where $`J(Y)\left(\frac{9I_{2b}^2}{4I_3^2}\frac{I_{1b}}{2I_3}\right)`$ and $`D(Y)=\alpha Y+C`$. Then the authors of  neglected the terms proportional to $`\beta ^2`$ (we will return to this point in the next section) and studied the 4d theory with the effective Lagrangian

$`_{\mathrm{eff}}={\displaystyle \frac{3a^3M_5^2}{n}}\left\{\left({\displaystyle \frac{\dot{a}}{a}}\right)^2+{\displaystyle \frac{\beta }{I_3}}\left[{\displaystyle \frac{1}{2}}D(Y)^2\dot{Y}^2n^2{\displaystyle \frac{V(Y)}{BI_3M_5}}\right]\right\}.`$ (16)
Since this Lagrangian leads to the equations which look exactly like the usual Friedmann equations for the scale factor $`a`$, the authors obtained the result which we expected on the basis of our general arguments: $`H=\frac{\dot{a}}{a}<0`$, i.e the 4d universe contracts.

At this stage it is very important to realize that we are not talking here about contraction of the bulk brane, or the visible brane, or the hidden brane. We are investigating an effective 4d geometry where the distinction between these branes completely disappeared. It was “washed away” by the integration over $`y`$. In a certain sense, one may imagine that integration over $`y`$ “glues” the three branes together. Thus we are talking about the contraction of the whole 4d space rather than of the one of the branes. All the difference between the branes in 4d must be encoded in the dynamics of the moduli fields, but not in the different rate of expansion or contraction of different branes.

Note that the contraction of the universe was crucial for ekpyrotic scenario in order to gain a residual kinetic energy of the moving brane and transfer it to the radiation upon the collision. This residual energy was suppressed by the small coefficient $`O(\beta /\alpha )`$ as compared to the maximal kinetic energy of the brane.

So where did the expansion of the universe come from in ? After solving the effective 4d equations, the authors decided to “unglue” the branes and go back to 5d. They studied the difference between the expansion of the universe as seen by different observers living on different branes and concluded that whereas the universe described by the overall scale factor $`a`$ collapses, the visible brane, described by the scale factor $`a_{}`$, may expand.

It is this point that is in an apparent contradiction with our 4d expectations suggesting that the scale factor of the universe in the ekpyrotic scenario cannot expand. However, the statement that $`a_{}`$ expands was based on the assumption that the function $`A`$ in the metric (4) depends only on $`t`$ but not on $`y`$, and that the effective 4d description with the metric ansatz of  provides a correct solution for the scale factor of the universe in the full 5d theory. As we will see in the next section, our analysis of the 5d solutions does not confirm this assumption.

## V 5d theory in the bulk

An assumption of the ekpyrotic cosmology in  was that the time dependent solution can be obtained from the static solution using the “moduli space” approximation, i.e. by replacing the constant moduli parameters of the static solution by time-dependent functions. More precisely, the ansatz used in  was to introduce only time but not $`y`$ dependence into $`A,N,Y`$ which were constants in the static BPS solution, so that they become $`A(t)`$, $`N(t)`$, $`Y(t)`$. All dependence on time in $`D(y)`$ then enters only via $`Y(t)`$, i.e. for the cosmological solution one takes $`D=D(y,Y(t))`$. The two extra parameters $`B,C`$ which appeared in the metric, in the dilaton and in the 4-form field were assumed to be time-independent.

As already discussed above, in  this ansatz was plugged into the 5-d action, and then integration over $`y`$ gave an action for time-dependent functions only. That action was taken as a starting point for a cosmological analysis. However, it is not clear a priori why this procedure is actually consistent, i.e. why the solutions of the resulting effective equations represent at the same time the solutions of the original 5-dimensional gravitational equations.

Indeed, it is well-known (see e.g. ) that the moduli space approximation, i.e. replacing moduli by time-dependent function may not always be a consistent. Unless the dependence on the internal coordinate ($`y`$ in the present case) is “homogeneous” so that it effectively “scales out” of the time-dependent higher-dimensional equations, one cannot simply replace these equations by their integrated (averaged over $`y`$) version – the moduli space approximation ansatz will not be consistent with the full set of the gravitational equations.

The question now arises whether it is possible to set up a 5d theory of ekpyrotic cosmology before integrating over $`y`$.

Since in  the interaction between the branes in 5d bulk was not specified but only the potential $`V(Y)`$ which lives on the bulk brane at $`y=Y`$ was introduced, we will also take this ansatz as part of the definition of the 5d theory. As we shall see, this will lead to a major problem with ekpyrotic cosmology, a contradiction with 5d solutions of equations of motion.

It could be expected that the 5d theory (3) does not describe the physics of the brane collision and does not address the moduli stabilization problem . The problem we are discussing now is completely different. Here we will analyse equations in the bulk during the roll of the bulk brane until it reaches the minimum of its potential $`V(Y)`$. If the 5d equations do not work even at this stage, then we do not have a consistent scenario not only during the brane collision but even before it.

To make the discussion as clear as possible, we shall start with the simplest equation – the one for the dilaton $`\varphi `$. The dilaton equation of motion in the bulk (away from the branes where the source terms vanish) is given by

$$\frac{2}{M_5^3}\frac{\delta S_{bulk}}{\delta \varphi (t,y)}=_\gamma (\sqrt{g}g^{\gamma \delta }_\delta \varphi )3\sqrt{g}\frac{e^{2\varphi }^2}{5!}=0.$$
(17)
It is satisfied for the above static solution (4)-(9). Plugging in the time-dependent ansatz of , we find that the new term in the dilaton equation, which was absent in static case, is simply the time-derivative one ($`\dot{D}_tD`$):

$$_t(\sqrt{g}g^{tt}_t\varphi )=3_t(A^3D^2BN^1\dot{D}).$$
(18)
The other terms in the dilaton equation, $`_y(\sqrt{g}g^{yy}_y\varphi )3\frac{e^{2\varphi }^2}{5!}`$, cancel not only for the static ansatz but for the time dependent ansatz as well, for the corrected action and solution given in (3)-(8). Following  in assuming that $`B`$ and $`C`$ are constant, we may study the dilaton equation which then reduces to the condition that (18) should vanish for all values of $`y`$ away from the branes. For $`y<Y`$ the term (18) vanishes since for $`D=C+\alpha y`$ we find that $`\dot{D}=0`$. However, for $`y`$ behind the moving brane, i.e. $`y>Y`$, we find that $`\dot{D}=\beta \dot{Y}0`$, and therefore the correction to the dilaton equation due to the time evolution does not vanish automatically. In this case Eq. (17) looks as follows:

$$\frac{2}{3M_5^3}\frac{\delta S_{bulk}}{\delta \varphi }=_t(A^3D^2BN^1\beta \dot{Y})=0,y>Y(t).$$
(19)
If we assume as in  that $`B`$ is constant and $`N(t)=A(t)`$ (which corresponds to the choice of $`t`$ as a conformal time), we get the condition

$$_t(A^2D^2\dot{Y})=0,y>Y(t).$$
(20)
The solution of this equation is

$$A^2(t)D^2(y,Y(t))\dot{Y}(t)=f(y),y>Y(t),$$
(21)
where $`f(y)`$ is an arbitrary function of $`y`$ that does not depend on $`t`$.

Let us first look at this equation in spirit of , replacing $`D^2(y,Y(t))`$ in the first approximation by its average value. When $`Y`$ changes from $`R`$ to $`0`$ in the ekpyrotic scenario, $`D(y,Y)`$ for any $`y>Y`$ changes just few times. Indeed, consider the maximum and minimum values of $`D(y,Y)`$. At $`y=0`$ in ekpyrotic cosmology it takes the minimum value $`D=C`$. At $`y=R`$ it takes the maximum value $`D(y,Y)=C+(\alpha \beta )R+\beta Y`$. Thus $`CD(y,Y)C+\alpha R`$. In the examples of  $`D_{min}=10^2`$ and $`D_{max}=310^2`$ or $`D_{min}=10^3`$ and $`D_{max}=210^3`$. Therefore, the function $`D`$ changes 3 and 2 times in these two examples, respectively, so one can in the first approximation replace $`D(y,Y)`$ by its average value.

On the other hand, in the ekpyrotic cosmology one had $`D\dot{Y}=\sqrt{2V(Y)}=\sqrt{2v}e^{\alpha mY/2}`$ . According to , when $`Y`$ changes from $`R`$ to $`0`$, the function $`e^{\alpha mY(t)/2}`$ changes from some value greater than $`e^{60}`$ to $`1`$. This means that during this time the scale factor $`A(t)\frac{f(y)}{D^2\dot{Y}(t)}`$ contracts at least $`e^{30}`$ times.

To compare this with the conclusions of  note that they define the scale factor of the visible brane as

$$a_{}(t)=A(t)\sqrt{C},$$
(22)
where $`C`$ is the constant in $`D`$. They conclude, using the 4d approximation, that $`a_{}(t)`$ grows in time and the visible universe expands. But as we have seen, this contradicts strongly to the solution of the 5d dilaton equation in the bulk, which shows that $`a_{}(t)`$ exponentially contracts.

However, the statement about expansion of the universe , as well as our statement about its exponential contraction, was based on averaging over the 5th dimension, assuming that the solutions of the 5d equations satisfying the metric ansatz of  do actually exist. Now let us look at our exact result, Eq. (21), more carefully. Since $`D(y,Y(t))=(\alpha \beta )y+C+\beta Y(t)`$ is a function of both $`y`$ and $`t`$, one can easily check that this equation is formally inconsistent, i.e. it does not have any solutions at all! This simply means that the ansatz for the metric and the fields used in  does not solve the 5d equations.

One can come to a similar conclusion using the gravitational equations in the bulk. First, note that general covariance leads to a relation between 5d Einstein equations, dilaton and 4-form equations. One can verify that the 4-form equation of motion in the bulk $`\frac{\delta S_{bulk}}{\delta A_{\alpha \beta \gamma \delta }}=0`$ is satisfied by the time dependent ansatz of Ref. . Thus we get an identity:

$$^\gamma \frac{\delta S_{bulk}}{\delta g_{\gamma \delta }}^\delta \varphi \frac{\delta S_{bulk}}{\delta \varphi }=0.$$
(23)
For a solution of bulk equations of motion each term in this identity should vanish. However, if the dilaton equation is not satisfied, this identity shows that a particular combination of the gravitational field equations cannot be satisfied by the ansatz of  (note that the derivatives $`^\delta \varphi `$ in the time and 5th direction do not vanish).

To see this more explicitly let us insert the ansatz of  directly into the action, with $`\dot{B}=\dot{C}=0`$. We shall change the variables so that

$$\stackrel{~}{A}^2(t,y,Y(t))=A^2(t)D^3(y,Y(t)),\stackrel{~}{N}^2(t,y,Y(t))=N^2(t)D^3(y,Y(t)).$$
(24)
Then, up to total derivatives,

$$_{\mathrm{bulk}}(t,y)=\frac{M_3^2}{2}\sqrt{g}\left(R\frac{1}{2}(\varphi )^2\frac{3}{2}\frac{e^{2\varphi }^2}{5!}\right)=\frac{3M_3^2\stackrel{~}{A}^3B}{\stackrel{~}{N}}\left[\left(\frac{\dot{\stackrel{~}{A}}}{\stackrel{~}{A}}\right)^2\frac{7}{4}\left(\frac{\dot{D}}{D}\right)^2\right].$$
(25)
Taking into account that $`\dot{\varphi }=3\dot{D}/D`$ we may also rewrite the bulk action as

$$_{\mathrm{bulk}}(t,y)=\frac{M_3^2\stackrel{~}{A}^3B}{\stackrel{~}{N}}\left[3\left(\frac{\dot{\stackrel{~}{A}}}{\stackrel{~}{A}}\right)^2\frac{7}{12}\dot{\varphi }^2\right].$$
(26)
This is a local action that should be varied to derive the local equations. If we perform the variation of this action over $`\varphi `$, we will get the same dilaton equation as in (20) (note that $`\stackrel{~}{A}`$ and $`\stackrel{~}{N}`$ are independent of $`\varphi `$), namely, $`_t(\frac{\stackrel{~}{A}^3B}{\stackrel{~}{N}}\dot{\varphi })`$=0 for $`y>Y`$.

The reason why this equation was not satisfied in  is that they took the action of the form (26), integrated it over $`y`$, added the brane action and neglected the term $`J(Y)\beta ^2\dot{Y}^2(t)𝑑y\left(\frac{M_3^2\stackrel{~}{A}^3B}{\stackrel{~}{N}}\dot{\varphi }^2\right)`$, as we have explained in the previous section (see eqs. (15), (16)). One may argue about whether this term is small or not as compared to the brane contribution in the integrated action. However, in 5 dimensions, i.e. before the integration over $`y`$, this term is given by $`\frac{\stackrel{~}{A}^3B\dot{\varphi }^2}{\stackrel{~}{N}}(t,y)`$. For $`y>Y`$ there is no other contribution to the dilaton equation. Therefore the variation of this term must vanish as Eq. (18) states.

## VI  Towards a consistent brane cosmology

As we have seen, different ways of taking averages in the situation where the 5d solutions do not exist can lead to dramatically different conclusions regarding expansion versus contraction of the universe. Thus, if one really wants to investigate cosmological consequences of the ekpyrotic scenario, one should find exact solutions of the corresponding 5d equations. This is especially important in the situation where the results of the averaging over the 5th dimension lead to the conclusions that are in an apparent contradiction with the Einstein equations in 4d, see Sect. IV.

The fact that the ansatz for the metric and the fields used in  does not solve the 5d equations is not very surprising. Indeed, it was shown in  that even in the simplest versions of brane cosmology describing one or two branes one should use a more general ansatz for the metric in order to satisfy the Israel junction conditions on the branes. A generic metric which respects the planar symmetry of the problem has the form

$$ds^2=n^2(t,y)dt^2+a^2(t,y)d\stackrel{}{x}^2+b^2(t,y)dy^2,$$
(27)
Note that here the functions $`a`$, $`b`$, and $`n`$ depend both on $`t`$ and $`y`$, and there still is a residual freedom of transformation of the coordinates $`(t,y)`$ Similarly, one may need to consider a more general ansatz for the fields as in .

Before one begins looking for exact solutions using a more general metric ansatz, one should re-examine other assumptions of the theory. Indeed, in  the potential $`V(Y)`$ was added by hand to the bulk brane action, whereas the bulk supergravity action remained unchanged. However, it is not obvious to us whether this is the proper way to introduce the inter-brane interactions in 5d.

As an illustrative example, consider two charged plates of a capacitor in ordinary electrodynamics, positioned at $`y=0`$ and $`y=R`$. If they have charges $`q`$ and $`q`$, and the electric field between the plates is $`E`$, then the potential energy of the interaction between the plates can be represented as the “brane potential” $`V(R)=qER`$. However, it would be incorrect to think that this energy is localized on the plates. Rather it is concentrated in the electric field between the plates. It is possible to use the potential $`V(R)`$ to study the motion of the branes. For example, if each brane has mass $`M`$, one can write $`m\ddot{R}=V^{}(R)`$, just as one does for the bulk brane acceleration in the ekpyrotic scenario. But if one studies gravitational backreaction of the electric field, it would be completely incorrect to replace the contribution of the electric field to the energy-momentum tensor in the bulk by the delta-functional term proportional to $`V(R)`$.

Similarly, if the potential $`V(Y)`$ appears due to the membrane instantons stretched between the branes, one should check whether the energy-momentum tensor in the bulk, as well as the dilaton and the 4-form field, changes due to these nonperturbative effects. Otherwise the appearance of the potential depending on the inter-brane separation would look as an example of the action at a distance.

In a certain sense, the potential $`V(Y)`$ is analogous to the effective potential $`V(r)`$ of the radion field introduced by Goldberger and Wise . It is a very useful concept if the only goal is to describe the forces acting on the bulk brane, ignoring the change of the metric produced by these forces. However, if one wants to study the corresponding changes in space-time geometry (and this was the main goal of ), one should perform a full investigation of the inter-brane interactions in 5d  and check whether one can add the fields responsible for the radion potential without additional fine-tuning and strong modification of the 4d geometry .

Thus one has a lot of things to do. First of all, one needs to find a theory with the potential $`V(Y)`$ which behaves as $`e^{\alpha mY}`$ at large $`Y`$ (the functional form is important). This potential should be smaller in absolute value than $`e^{120}`$ near the hidden brane, and should not have any positive contributions there with this accuracy. This fine-tuning is necessary to produce desirable density perturbations and avoid inflation. Also, this potential should not receive any contributions proportional to $`e^{\alpha m(RY)}`$ due to the interaction with the hidden brane . One must make sure that this potential vanishes at $`y=0`$, to avoid the cosmological constant problem. Then one must take into account that the visible brane and the hidden brane should be stabilized by some strong forces so that the effective potential of the corresponding moduli field could have mass on the TeV scale or even greater. One should also check that the strong forces leading to the brane stabilization do not interfere with the extremely weak interaction responsible for the potential $`e^{\alpha mY}`$. One cannot ignore the unresolved problem of brane stabilization (which was the position taken in ) and speculate about the inter-brane potentials suppressed by a factor of $`e^{120}`$.

When/if the theory with the desirable $`V(Y)`$ is found, one should investigate its $`5d`$ nature and make sure that the effects producing $`V(Y)`$ do not induce large curvature on the branes. Then one should write down equations in 5d taking into account the energy-momentum tensor in the bulk together with the junction conditions, and solve them.

And finally, one should find out what happens at the moment of the brane collision: whether the visible brane collapses, expands, stays at the same place or oscillates, etc. These issues have not been addressed in , and they cannot be fully analysed until the brane stabilization mechanism is understood.

Thus, if one wants to propose a consistent alternative to inflationary cosmology, one would need first to give a consistent formulation of the alternative theory, and then find a correct solution of the corresponding equations. As we have seen, this is a rather nontrivial task.

In this paper we studied a very limited part of this problem. We tried to check whether the basic assumptions of the ekpyrotic scenario (the ansatz for the metric and for the fields, and the modification of the bulk brane action proposed in ) can lead to a consistent 5d description of an expanding visible brane. We have found that this is not the case.

The authors are grateful to T. Banks, K. Benakli, M. Dine, N. Kaloper, J. Maldacena, V. Mukhanov, R. Myers, P. Nilles, A. Peet, K. Stelle, E. Witten for stimulating discussions and to A. Van Proeyen for confirming the calculations of the 4-form. We thank NATO Linkage Grant 975389 for support. L.K. was supported by NSERC and by CIAR; L.K. and A.L. were supported by NATO Linkage Grant 975389; R.K. and A.L. and were supported by NSF grant PHY-9870115; A.L. was also supported by the Templeton Foundation grant 938-COS273. A.T. was supported by the DOE grant DE-FG02-91ER40690, PPARC SPG grant 00613 and INTAS project 991590.

## VII Appendix: The 4-form story

There are two points about the 4-form dependence in  which must be changed to get the correct set up for charged BPS 5d domain walls. One has to correct the coefficient in the action and the coefficient in the solution for the 4-form. These corrections in the form sector are important in order to test the time-dependent ansatz of .

### A Action

The action in eq. (10) of  is not the one to which they refer as given in , . The one in  does not have a 4-form. We questioned the origin of their action in  and they replied in : The 4-form formulation of the action is equivalent to the action presented in . This is easily seen by eliminating the 4-form using its equation of motion.

In  there were indeed the relevant world-volume terms in the action (only 2 branes at the fixed points are present there, so there $`\beta =0`$), but there were no Wess-Zumino terms with the 4-form field:

$$S=\frac{1}{2\kappa _5^2}__5d^5x\sqrt{g}\left[R+\frac{1}{2V^2}_\alpha V^\alpha V+\frac{1}{3V^2}\alpha ^2\right]\underset{i=1}{\overset{2}{}}\frac{\sqrt{2}}{\kappa _5^2}\alpha _i_{_4^{(i)}}d^4\xi _{(i)}\sqrt{g}V^1,$$
(28)
We take the constants in the bulk and boundary actions as $`\alpha =\alpha ^{(1)}=\alpha ^{(2)}`$. Then we make the redefinitions $`Ve^\varphi ,\alpha \frac{3}{\sqrt{2}}\alpha ,\kappa _5^2M_5^3,RR`$ so that our notation agree with . That leads to

$$S=\frac{M_5^3}{2}__5d^5x\sqrt{g}\left[R\frac{1}{2}_\alpha \varphi ^\alpha \varphi \frac{3}{2}e^{2\varphi }\alpha ^2\right]3\underset{i=1}{\overset{2}{}}M_5^3\alpha _i_{_4^{(i)}}d^4\xi _{(i)}\sqrt{g}e^\varphi .$$
(29)
The kinetic terms match the ones in (10) in . Now we may perform the procedure suggested in . We promote the constant $`\alpha `$ to a function $`G(x)`$ and add a Lagrange multiplier of the form $`\epsilon ^{\alpha \beta \gamma \delta ϵ}(_ϵG)𝒜_{\alpha \beta \gamma \delta }`$ and assign the supersymmetry transformation to the 4-form so that its variation will compensate the variation of terms in the rest of the action with derivatives of $`G`$. If there are no sources, from the equation for the 4-form we learn that $`G(x)`$ is a constant, as before in usual gauged supergravity in d=5 without a 4-form where there is a constant gauge coupling. In the presence of sources we will find that $`G`$ is piecewise constant. If there are charged sources as in case of  we may also add the WZ term.

Thus we add a Lagrange multiplier term and a WZ term to the action of  given above in the form (29). Its normalization is arbitrary, we thus put a constant $`c`$ in front:

$$S_A=\frac{c}{4!}__5d^5x\left[\epsilon ^{\alpha \beta \gamma \delta ϵ}(_ϵG)𝒜_{\alpha \beta \gamma \delta }+\underset{i=1}{\overset{2}{}}2\alpha _i\epsilon ^{\mu \nu \rho \sigma }𝒜_{\mu \nu \rho \sigma }\delta (x^5x_i^5)\right].$$
(30)
The relative normalization between the two terms is arranged such that $`\alpha `$ jumps by 2$`\alpha _i`$ at brane $`i`$. If we differentiate over the 4-form we are back to (29). To have agreement with  we have to choose $`c=\frac{3}{2}M_5^3.`$ We may however do something else, namely, add and subtract a term quadratic in $``$. We find

$$S_G=__5d^5x\left[\frac{3M_5^3}{4}\sqrt{g}e^{2\varphi }\left(G\frac{2e^{2\varphi }}{3M_5^3\sqrt{g}}\frac{c}{5!}\epsilon ^{\alpha \beta \gamma \delta ϵ}_{\alpha \beta \gamma \delta ϵ}\right)^2\sqrt{g}\frac{e^{2\varphi }}{5!3M_5^3}c^2^2\right].$$
(31)
Now we vary over the field $`G`$ and from its equation we find that

$$G=\frac{2e^{2\varphi }}{3M_5^3\sqrt{g}}\frac{c}{5!}\epsilon ^{\alpha \beta \gamma \delta ϵ}_{\alpha \beta \gamma \delta ϵ}.$$
(32)
We thus find the action given in  but where the factor in front of the $`^2`$-term has to be changed by $`3/2`$. The correct action (i.e. the one which is equivalent to an action without 4-form in ) is thus

$`S={\displaystyle \frac{M_5^3}{2}}{\displaystyle __5}d^5x\sqrt{g}\left(R{\displaystyle \frac{1}{2}}(\varphi )^2{\displaystyle \frac{3}{2}}{\displaystyle \frac{e^{2\varphi }^2}{5!}}\right)`$ (33)
$`3{\displaystyle \underset{i=1}{\overset{3}{}}}\alpha _iM_5^3{\displaystyle _{_4^{(i)}}}d^4\xi _{(i)}\left(\sqrt{h_{(i)}}e^\varphi {\displaystyle \frac{ϵ^{\mu \nu \kappa \lambda }}{4!}}𝒜_{\gamma \delta ϵ\zeta }_\mu X_{(i)}^\gamma _\nu X_{(i)}^\delta _\kappa X_{(i)}^ϵ_\lambda X_{(i)}^\zeta \right).`$ (34)
This action for the two branes at the fixed points (i.e. without the third brane) taken in the static gauge with $`X^\mu =\xi ^\mu `$ and $`X^4=Y=const`$ does agree with the action in .

### B Bulk brane and 4-form

As was promised, we will show here that the formulation of the 5d supersymmetric theory with a bulk brane between the orbifold planes is impossible without the use of the 4-form or, equivalently, without the gauge coupling field $`G`$ dual to the field strength $`F=dA`$ .

In the case of two orbifold planes we may derive the action (28) by solving the equation of motion for the 4-form. This leads to

$$G^{}=2\alpha \left(\delta (y)\delta (yR)\right).$$
(35)
The solution for the gauge coupling field is $`G(y)=\alpha ϵ(y)`$. It simply means that the gauge coupling is positive for positive $`y`$ and negative for negative $`y`$. Therefore, $`(G(y))^2=\alpha ^2`$ for all $`y`$ and we find the bulk potential proportional to $`e^{2\varphi }\alpha ^2`$ as it was given in the original form of the action (29) where there was no 4-form, neither in the bulk, nor on the branes.

Now we may try to perform the same procedure of getting rid of the 4-form in case when the bulk brane is present. We find the gauge-coupling field dual to the 5-form field strength

$$G^{}=2\alpha \left(\delta (y)\delta (yR)\right)\beta \delta (yY)\beta \delta (y+Y).$$
(36)
In this case the solution for $`G`$ which is dual to the 5-form takes values

$`G^2`$ $`=`$ $`\alpha ^2\mathrm{at}0<|y|<Y,`$ (37)
$`G^2`$ $`=`$ $`(\alpha \beta )^2\mathrm{at}Y<|y|<R.`$ (38)
If we would try to plug this solution back into the action (into the term $`e^{2\varphi }G^2`$), we would find that the bulk potential is a piecewise function, and, therefore, there is no local general covariant theory. If, however, we keep the 4-form, there exists a local general covariant 5d bulk action. In the presence of sources, this will lead to piecewise values of the 5-form for the solutions, but the local general covariant action is available.

### C Solution

Let us suppose first that the action given in  were correct and let us try to check the solution for the 4-form. The corresponding equation is

$$_y\left(\sqrt{g}e^{2\varphi }^{0123y}\right)+3\left[\alpha \delta (y)+(\alpha \beta )\delta (yR)+\frac{\beta }{2}\delta (yY)+\frac{\beta }{2}\delta (y+Y)\right]=0.$$
(39)
For the ansatz in  we have $`\sqrt{g}=A^3BN[D(y)]^4`$ while $`\sqrt{g_{(4)}}=A^3N[D(y)]^2`$ and $`e^{2\varphi }=B[D(y)]^3`$. We find then that the solution is

$$_{0123y}=\frac{3}{2}A^3NB^1D^2(y)D^{}(y),$$
(40)
where

$$[D(y)]^{^{\prime \prime }}=2[\alpha \delta (y)(\alpha \beta )\delta (yR)]\beta \delta (yY)\beta \delta (y+Y).$$
(41)
This differs from solution given in  by a factor of $`\frac{3}{2}A^3NB^1`$ and by a sign. If we now start with the corrected form of the action in eq. (3) (with extra kinetic term factor $`\frac{3}{2}`$), we find the corrected expression for the 4-form

$$_{0123y}=A^3NB^1D^2(y)D^{}(y).$$
(42)
This solution differs from the one in  by the sign and by the factor $`A^3NB^1`$.


format: text

text:
# The SNR G106.3+2.7 and its Pulsar Wind Nebula: relics of triggered star formation in a complex environment

## 1 Introduction

The extended radio continuum source G106.3+2.7 was first classified as a supernova remnant (SNR) by Joncas & Higgs (1990). This SNR was later described by Pineault & Joncas (2000) in more detail. They separated the comet shaped structure into a small head with high radio surface brightness and a large extended tail of low surface brightness. A gap in the atomic hydrogen apparently correlated with the SNR suggested a distance of about 12 kpc. This would give the SNR a length of 200 pc and they concluded that G106.3+2.7 is a supernova remnant in the late stage of its isothermal evolution.

At the northern edge of the head there is a small shell-like structure with a diameter of about 3.5’. This source coincides with strong X-ray emission and the $`\gamma `$ray source 3EG J2227+6122. It was interpreted as the bow shock nebula of a fast moving pulsar or a pulsar wind blown bubble (Halpern et al., 2001a). These authors estimate the distance as $``$ 3 kpc from X-ray absorption values. Subsequently this pulsar, PSR J2229+6114, was detected in radio and X-rays (Halpern et al., 2001b). It has a period of 51.6 ms and a rotational energy loss rate of $`\dot{\mathrm{E}}=2.210^{37}`$ erg/s.

In this paper, by utilizing the data provided by the Canadian Galactic Plane Survey (CGPS), we analyze the possible relationship between the pulsar nebula and the supernova remnant. Since the number of known combined SNRs is rather low, such an analysis would improve our conceptual understanding on the interaction of a pulsar nebula with its environment.

## 2 Observations and Data Reduction

The observations discussed here form part of the CGPS, which is described in detail by Taylor et al. (2001). In carrying out the research presented here, we have used those parts of the CGPS database which derive from observations with the DRAO Synthesis Telescope (Landecker et al., 2000) and CO observations from the Five College Radio Astronomy Observatory (FCRAO) (Heyer et al., 1998). Angular resolution varies as cosec(declination) and therefore changes slowly across the final mosaics. Before assembly into a mosaic, the data for the individual fields are carefully processed to remove artifacts and to obtain the highest dynamic range, using the routines described by Willis (1999). Accurate representation of all structures to the largest scales is assured by incorporating data from large single antennas with data from the Synthesis Telescope, after suitable filtering in the Fourier domain. Continuum data are derived from the 408 MHz all-sky survey of Haslam et al. (1982), and from the 1.4 GHz Effelsberg survey (Reich et al., 1997). Single antenna H I data are obtained from a survey of the CGPS region made with the DRAO 26-m Telescope (Higgs & Tapping, 2000). To improve the resolution of the 1420 MHz continuum observations we also included data taken from the NVSS (Condon et al., 1998). The NVSS map was convolved to the resolution of our CGPS data. The resulting map was subtracted from the CGPS map and the difference added to the original NVSS map. The final image has a resolution of 45”.

## 3 The structure of G106.3+2.7

In Fig. 1 the radio image of G106.3+2.7 at 1420 MHz is shown. In the lower image the overall structure of the source is clearly visible. The major parts are a comet shaped structure with a diffuse extended tail of low radio surface brightness and a smaller more structured head with higher radio surface brightness. The upper image is an enlargement of the latter. The brightest feature is the “Boomerang” to the north, described by Halpern et al. (2001a). The point-like source just to its south is most likely an extragalactic source (Halpern et al., 2001a). The head itself shows two prominent features. There is a peak in the center, south of the extragalactic source and a shell-like structure on its southern boundary. Also visible is a diffuse bridge extending from the head to the pulsar nebula, suggesting an interaction of those structures. The SNR’s radio continuum emission has a spectral index of $`\alpha 0.57`$ (S$`\nu ^\alpha `$) (Pineault & Joncas, 2000). Neither the tail nor the head show significant spectral fluctuations. However, the pulsar nebula seems to have a flat radio spectrum (Halpern et al., 2001a). We obtained a flux density of $`90\pm 5`$ mJy by integration of our 1420 MHz data. Together with a flux density of 80 mJy at 4850 MHz taken from Becker et al. (1991) we derive a spectral index of $`\alpha =0.1\pm 0.1`$ for the pulsar nebula. This confirms the flat radio spectrum proposed by Halpern et al. (2001a) and supports its interpretation as a pulsar-powered synchrotron nebula.

Fig. 2 displays a map of polarized intensity at 1420 MHz in greyscale together with white contours representing the total power radio continuum. We find linear polarization in all parts of the object. The tail has only a small number of polarized patches but these show high percentage polarization of up to 70 %. The head is polarized except for an area between the pulsar nebula and the center of the head where we also find a peak in total power emission. The head shows peak polarization of up to 70 %. The most prominent polarization feature, however, is the pulsar nebula with an integrated percentage polarization of 40 %. Parts of the radio continuum bridge to the pulsar nebula appear to be highly polarized. These high percentage polarizations indicate low depolarization between the SNR and us. The mean rotation measure towards the “Boomerang” and the SNR G106.3+2.7 is around $`RM200`$ rad/m2. This high rotation measure results in Faraday rotation of $`500^{}`$ at a frequency of 1420 MHz. Since we find only low depolarization on the “Boomerang” it requires a very small distance which makes this object a local phenomenon, probably within 2 kpc. If it was further away, beam depolarization would occur because of superimposed regions of differing Faraday rotation within the beam.

The large OB association Cep OB2, which is located at a distance of about 950 pc (Garmany & Stencel, 1992), lies in the direction of G106.3+2.7 and partly overlaps with the supernova remnant. The members of this OB association ionize the material between them which leads to the depolarization of extended polarized emission from behind them at low radio frequencies. This gives us another upper limit for the distance to G106.3+2.7.

## 4 The structure of the cold environment

### 4.1 The head’s hydrogen envelope

In Fig. 3 we plot channel maps of our H I data cube representing an area around the head of G106.3+2.7 at low negative radial velocities. There is a shell-like H I structure around the eastern edge of the head best seen at -6.41 km/s but also quite prominent at -5.59 km/s and -7.23 km/s. A shell of H I is wrapped around the “Boomerang”, peaking at the same velocity as the shell around the head of the SNR. The coincidence in velocity clearly demonstrates that the SNR and the pulsar nebula are associated. This is the first time that an interaction between a pulsar nebula and cold interstellar material has been observed. It looks like the pulsar is either pushing the neutral hydrogen outwards or ionizing cold atomic hydrogen in its vicinity with its pulsar wind and thus creating the shell-like structure.

The center velocity of -6.4 km/s puts the SNR at a distance of 800 pc. At this distance the SNR has a length of 14 pc and a maximum width of 6 pc. The pulsar nebula is 0.8 pc wide. The small shell of atomic hydrogen contains about 0.4 M and the bigger shell around the head about 5 M. This leads to densities of about 200 cm-3 within both structures.

Halpern et al. (2001a) calculated the foreground H I column density to be $`N_H=6.310^{21}`$ cm-2 from the spectral fit to their X-ray data. Based on this value they estimated a distance of 3 kpc. We integrated the foreground H I column density from our data using $`N_H[\mathrm{cm}^2]=1.82310^{18}T_B(v)𝑑v`$ (Kerr, 1968) with $`T_B`$ representing the H I excess brightness temperature and $`v`$ its radial velocity. The result is $`N_H210^{21}`$ cm-2. This equation, however, is only valid if the neutral hydrogen in the foreground is optically thin. Comparison of the brightness temperature of two strong point like background sources in the vicinity of G106.3+2.7 ($`l=105.68^{}`$, $`b=2.42^{}`$ and $`l=107.84^{}`$, $`b=2.31^{}`$) with their H I absorption profiles indicates that most of the foreground material is optically thick with opacities higher than 3. This easily explains the factor of 3 between the calculated foreground $`N_H`$ and the one derived from the X-ray data.

### 4.2 The head’s molecular envelope

In Fig. 4 we plot channel maps of our CO data cube in the same area and at the same radial velocities as the H I data in Fig. 3. At the same velocities where we found the eastern hydrogen envelope there is a molecular envelope around the western edge of the head structure. But in this case we do not see the same structure in all channels. The location of the CO envelope is to the west at -5.6 km/s but rotates to the northwest at -7.2 km/s. Since around the eastern edge we find neutral hydrogen in all channels we can assume that the hydrogen completely surrounds this area. The CO, on the other hand, seems to be incomplete giving the expanding supernova remnant space to break through and form the smooth extended tail. The narrow shell of molecular material contains about 40 M of molecular hydrogen assuming $`N_{H_2}[\mathrm{cm}^2]=2.310^{20}T_{CO}\delta v`$. Here $`N_{H_2}`$ is the column density of molecular hydrogen, $`T_{CO}`$ \[K\] is the CO brightness temperature and $`\delta v`$ \[km/s\] is the velocity width of these structures.

### 4.3 The big picture

Fig. 5 shows an image of the entire region surrounding G106.3+2.7. The head of the SNR and its molecular and atomic envelope are visible in the lower left part of the image as is the small shell surrounding the pulsar wind. Apparently these structures represent only a small part of the entire cold environment. There is a large bubble centered at $`l106.1^{}`$ and $`b3.1^{}`$, with an inner diameter of about $`1^{}`$ which translates to about 14 pc at a distance of 800 pc. The outer diameter is about 20 pc. At the northern and eastern edge of the bubble we find patches of molecular material. The supernova remnant has apparently created a little bubble at the lower left edge of the large bubble with its shock wave.

## 5 Discussion

### 5.1 The creation of the cold environment

We suggest that a first generation of stars was formed in the center of the big H I bubble which at that time was a cloudy dense medium. With stellar wind effects and/or supernova explosions the large bubble was formed. However, some very dense cloudlets remained because they were too dense and probably also too far away from the location of the exciting stars. The fact that the remaining patches we observe are located at the outer edges of the bubble supports this, as does an H I structure surrounding the northeastern patch of molecular material, which is most likely the remains of dissociated molecular hydrogen at the surface of this cloud. Stellar winds and supernova explosions within the large bubble apparently also triggered the creation of the progenitor star of our SNR G106.3+2.7. At the time the progenitor star of G106.3+2.7 was formed, the environment most likely looked very much like the two northern patches of molecular material. The new born star was located inside this molecular environment. It was a massive star, probably an early B or an O-type star. The small bubble in the cold environment in which we now find the head of the SNR was either formed by the progenitor star’s stellar wind or by the expanding shockwave of the supernova explosion itself.

### 5.2 A bow shock or a pulsar wind nebula

Usually pulsar nebulae are expanding inside a bubble created by the shock wave of a supernova explosion. In this case there is no material left to interact with, because all the matter has been swept up and carried away already. The pulsar in G106.3+2.7 is apparently interacting with its dense cold environment. In the following we investigate two possible scenarios leading to the shell-like shape of the “Boomerang” pulsar nebula.

#### 5.2.1 The “Boomerang” as a bow shock nebula

Since the shape of the object bears a striking resemblance to a bow-shock structure, we first investigate the constraints which are imposed by the observations in order for this hypothesis to be viable. For simplicity, we assume that the explosion which caused the SNR took place near the center of the head component and resulted in the candidate pulsar getting a kick velocity which took it to its present position near the bow structure. In view of the observed morphology, we also suppose the motion to be nearly in the plane of the sky and neglect projection effects. The distance travelled, $`L`$, is comparable to the radius of the head component $`R_s`$, so we set $`LR_s`$. Letting $`t`$ represent the pulsar travel time and SNR age in years, $`v_p`$ the pulsar space velocity in km/s, and $`R_s`$ in pc, we have

$$R_s1.010^6v_pt$$
(1)
The equality between the pulsar relativistic wind pressure and ISM ram pressure implies that

$$\frac{\dot{E}}{4\pi l^2c}=n_o\mu m_Hv_p^2\dot{E}_{36}6.2\mu n_ol^2v_p$$
(2)
Here $`n_o`$ is the ambient number density in cm-3, $`l`$ is the distance between the pulsar and the apex of the bow shock in pc, and $`\dot{E}_{36}`$ is the pulsar rotational energy loss rate in $`10^{36}`$ erg/s. SNR dynamics give us:

$$R_s=\{\begin{array}{cc}0.34\left(\frac{E_{51}}{\mu n_o}\right)^{0.2}t^{0.4}\hfill & \text{for adiabatic expansion}\hfill \\ 1.4\left(\frac{ϵE_{51}}{\mu n_o}\right)^{5/21}t^{2/7}\hfill & \text{for isothermal expansion}\hfill \end{array}$$
(3)
Here the parameter $`ϵ`$ is the ratio of thermal energy in the SNR to the initial explosion energy $`E_{51}`$ \[$`10^{51}`$ erg\].

The pulsar distance from the center of the head is about 1.5 pc, its distance from the apex of the bowshock is about 0.4 pc and its rotational energy loss rate $`\dot{E}=2.210^{37}`$ erg/s (Halpern et al., 2001b). With an estimated value of 0.3 for $`ϵ`$ we get:

$$E_{51}\{\begin{array}{cc}1.710^8v_p\hfill & \text{adiabatic case}\hfill \\ 1.210^6v_p^{0.2}\hfill & \text{isothermal case}\hfill \end{array}$$
(4)
Any reasonable value for the pulsar velocity (100 km/s $`v_p`$ 1000 km/s) will give us by far the lowest ever recorded explosion energy for the supernova explosion. Allowing for inclination effects or substantial variations in the observed parameters still leaves us with highly unlikely values of the energy parameters.

#### 5.2.2 The “Boomerang” as a pulsar wind nebula

As an alternative, Halpern et al. (2001a) have suggested that the bow-shaped structure can be confined by the thermal pressure of the surrounding medium. This is indeed possible, but it implies an anomalously low velocity for an object belonging to a high-velocity class. The fact that the spectral index is flat makes it almost unavoidable that the emission originates from the relativistic pulsar wind itself. However we believe that the confinement is simply achieved because the relativistic wind is encountering the dense SNR shell in the northeast. In our picture the SN explosion took place at or very near the present position of the pulsar. In the northeast the shock has been quickly decelerated by its encounter with a relatively nearby and particularly high density medium, whereas towards the south it has expanded into a moderately dense medium (SNR head component) while breaking out in the southwest into a region of much lower density and giving rise to the tail component.

The morphology of the head and the tail structure supports this hypothesis. The head shows only the “Boomerang” and the southern shell as indications for the encounter of a shockwave with dense material. We do not find a shell structure towards the east of the head but we find an H I shell all around the east and the south of the head. We can explain this by moving the location of the explosion to the current position of the pulsar. To the north and to the east the shockwave is running into dense material. The first free path for the shockwave is to the southeast where the southern shell begins. Another point for this scenario is that the distance between the location of the supernova and the shell to the south is now larger that its distance to the beginning of the tail, which should be expected in our picture. If the explosion occured in the center of the head this would not be the case.

What seems peculiar however is the fact that the little bow-shaped shell in the north does not show enhanced non-thermal radio emission. This would appear to go against the expected view that an SNR should be brighter where the blast wave is moving against dense regions of the ISM. A similar situation has been encountered by Pineault et al. (1997) in their analysis of the CTA1 SNR. These authors argued that, in overdense regions where the shock velocity is considerably reduced, the Fermi acceleration mechanism (Bell, 1978) might either be too slow or simply break down. The time efficiency can be quantified by evaluating the timescale $`\tau _{\mathrm{acc}}`$ for acceleration of an electron to energy $`E`$, given by $`\tau _{\mathrm{acc}}E/Bv^2`$ (Ellison et al., 1990; Reynolds, 1996). Hence either a low shock velocity or low ambient magnetic field, or a combination of both may imply that the electrons radiating at a given observing frequency $`\nu `$, where $`E(\nu /B)^{1/2}`$, are not yet present in the energy distribution. An alternative approach is to consider the growth and decay rates of the Alfvén waves responsible for the repeated scattering of the relativistic particles back and forth across the shock (Bell, 1978; Blandford & Ostriker, 1978). On one hand, the growth rate is proportional to the shock velocity and thus decreases as the shock slows down and, on the other hand, the decay rate depends on the first power of the number density of neutral atoms. Again this means that the acceleration mechanism may be quenched whenever the shock wave encounters sufficiently dense clouds. It is also possible that a thin steep (spectral index of order -0.5) non-thermal layer is present but not detectable because of beam smearing.

The interpretation of the tail component as resulting from a breakout of the SNR shock wave is also compatible with theoretical predictions. The numerical studies of breakout by Tenorio-Tagle et al. (1985) have shown that the asymptotic ratio of the velocity of the material breaking out into a region of lower density to that of the shock wave still propagating into the original medium could be expressed as

$$B_{\mathrm{}}=4.221.89(\mathrm{log}_{10}\lambda +2),$$
(5)
where $`\lambda 1`$ is the ratio of the densities. The same ratio at any normalized time $`\tau =t/t_b`$, where $`t_b`$ is the time at which breakout occurred, was shown by Pineault et al. (1987) to be satisfactorily represented by the approximate analytical relation

$$B(\tau )=B_{\mathrm{}}(B_{\mathrm{}}1)/\tau .$$
(6)
It follows from this relation that the ratio $`r=L_t/R_sB_{\mathrm{}}`$ where $`L_t`$ is the length of the tail component, as long as $`\tau >>1`$. Using $`R_s3`$ pc and $`L_t12`$ pc gives us a value of 4 for $`r`$. The length of the tail is difficult to determine since there is no sharp outer boundary. The observed value implies $`\lambda 0.01`$. A value of 8 for $`r`$ would correspond to $`\lambda 10^4`$. It can also be shown that the breaking out material attains a maximum velocity

$$v_{bm}=v_b[0.35B_{\mathrm{}}(B_{\mathrm{}}/(B_{\mathrm{}}1))^{0.6}]$$
(7)
at $`\tau _m=(8/3)(11/B_{\mathrm{}})`$, where $`v_b`$ is the velocity just before breakout. For $`B_{\mathrm{}}=4`$ , the maximum velocity is about $`1.7v_b`$. Assuming that the age of the SNR is comparable to the characteristic age of the pulsar we take $`t10460`$ yr (Halpern et al., 2001b). For an adiabatic expanding SNR we derive $`E_{51}/\mu n_04.910^4`$. With a density of $`n_0100`$ cm-3 the explosion energy would be $`E710^{49}`$ erg. Usually the characteristic age of a pulsar is higher than its real age which means the calculated explosion energy is a lower limit. From the length of the tail, an average velocity of $`1125\mathrm{km}\mathrm{s}^1`$ is inferred. Equating this to the maximum velocity taken as $`1.7v_b`$ gives $`t_b540`$ yr, corresponding to $`R_{sb}1`$ pc. Although these numbers should not be taken at face value, they nevertheless indicate that the breakout hypothesis is consistent with the observed morphology.

## 6 Conclusion

Using the data from the CGPS we studied the morphology and the environment of the SNR G106.3+2.7 and the pulsar nebula associated with the pulsar J2229+6114 (the “Boomerang”). We conclude that these sources are the result of the same SN event. The SNR morphologically consists of two parts giving the object its cometary shape, in which the tail is due to an outbreak into the interior of a big H I bubble, whereas the head is created by the expanding shock wave interacting with dense ambient material. The kinematics of the associated neutral hydrogen and molecular material, $`v_{kin}6.4`$ km/s, suggests that both objects are located at the same distance of 800 pc, thereby the SNR is 14 pc long and 6 pc wide and the pulsar nebula has a diameter of 0.8 pc. This close distance is also supported by the presence of high polarized emission.

By studying the energetics of the SNR we found that the “Boomerang” nebula is created by the relativistic wind of the pulsar. For the “Boomerang” we found a spectral index of $`\alpha =0.1`$ and $``$ 40 % integrated percentage polarization at 1420 MHz. The analysis of the cold environment suggests that the creation of the progenitor star was triggered by stellar wind and/or SN explosions of a group of stars which also created the H I bubble. In this way this is the first detection of a pulsar nebula apparently interacting with its cold environment. It is also an example of a pulsar “displaced” from the center of its associated SNR. In this case, however, the displaced position of the pulsar is not due to high pulsar velocity.

The Dominion Radio Astrophysical Observatory is a National Facility operated by the National Research Council. The Canadian Galactic Plane Survey is a Canadian project with international partners, and is supported by the Natural Sciences and Engineering Research Council (NSERC). SP was supported by NSERC and the Fonds FCAR of Québec. We thank Lloyd Higgs and Tom Landecker for careful reading of the manuscript and discussions.


text:
# Glass transition in secondary structures formed by random RNA sequences

## Abstract

Formation of RNA secondary structures is an example of the sequence-structure problem omnipresent in biopolymers. A theoretical question of recent interest is whether a random RNA sequence undergoes a finite temperature glass transition. We answer this question affirmatively by first establishing the perturbative stability of the high temperature phase via a two replica calculation. Subsequently, we show that this phase cannot persist down to zero temperature by considering energetic contributions due to rare regions of complementary subsequences.

Introduction: RNA is an important biopolymer critical to all living systems rna . Like in DNA, there are four types of nucleotides (or bases) A, C, G, and U which, when polymerized can form double helical structures consisting of stacks of stable Watson-Crick (A–U or G–C) pairs. However unlike a long polymer of DNA, which is often accompanied by a complementary strand and forms otherwise featureless double helical structures, a polymer of RNA is usually single-stranded. It bends onto itself and forms elaborate spatial structures in order for bases located on different parts of the backbone to pair with each other. The structures encoded by the primary sequences often have important biological functions, much like how the primary sequence of a protein encodes its structure.

Understanding the encoding of structure from the primary sequence has been an outstanding problem of theoretical biophysics. Most work in the past decade has been focused on the problem of protein folding, which is very difficult analytically and numerically protein . Here, we study the problem of RNA folding, specifically the formation of RNA secondary structures which is more amenable to analytical and numerical approaches due to a separation of energy scales tino99 . Efficient algorithms mcca90 ; rnaalg  together with carefully measured free energy parameters frei86  describing the formation of various microscopic structures (e.g., stacks, loops, hairpins, etc.) allow the exact calculation of the ensemble of secondary structures formed by a given RNA molecule of up to a few thousand bases.

In this work, we are not concerned with the structure formed by a specific sequence. Instead, we will study the statistics of secondary structures formed by the ensemble of long random RNA sequences. It has been debated recently whether such an ensemble undergoes a finite temperature glass transition higg96 ; pagnani . However, the numerical results these studies are based on are not clear enough to allow unambiguous interpretation. In this letter, we provide analytical evidence supporting the existence of a finite temperature glass transition by studying some toy models of RNA folding. We characterize the behavior of RNA in the absence of disorder and show with the help of a two-replica calculation that disorder is perturbatively irrelevant. We then show that the assumption that the high-temperature phase exists at all temperatures leads to a contradiction below some finite temperature, thereby necessitating a finite temperature phase transition.

Model: A secondary structure $`S`$ of a polymer of RNA is the set of pairings formed between all of its monomers (or bases), with each base allowed in at most one pairing. We denote the pairing between the $`i^{\mathrm{th}}`$ and $`j^{\mathrm{th}}`$ monomer by $`(i,j)`$, with $`1i<jN`$. Each such structure can be represented by a diagram like the one shown in Fig. 1(a). In addition, it is common to exclude “pseudo-knots” like the one shown in Fig. 1(b) from the definition of secondary structures, so that any two base pairs $`(i,j)`$ and $`(k,l)`$ are either independent, i.e., $`i<j<k<l`$, or nested, i.e., $`i<k<l<j`$. This is permissible since long pseudo-knots are kinetically difficult to form and even the short ones occur infrequently due to energetic reasons tino99 . Experimentally, it is possible to “turn off” the pseudo-knots and other complicated tertiary contacts tino99  and study exclusively the class of secondary structures defined above.

In order to calculate Boltzmann factors for an ensemble of secondary structures, we need to assign an energy $`E[S]`$ to each structure $`S`$. For the purpose of secondary structure prediction, it is essential to model the energy as accurately as possible mcca90 ; rnaalg . However, for our interest in the universal statistical properties of long, random RNA sequences far below the denaturation temperature, it suffices to consider simplified models along the line of those used in earlier studies higg96 ; pagnani ; bund99 ; higg00 .

We associate an interaction energy $`\epsilon _{i,j}`$ with every pairing $`(i,j)`$ and assign $`E[S]=_{(i,j)S}\epsilon _{i,j}`$ as the total energy of the structure $`S`$. To retain the spirit of Watson-Crick pairing, we choose random sequences $`b_1\mathrm{}b_N`$ of the four bases $`A`$, $`U`$, $`C`$, and $`G`$ and then assign

$$\epsilon _{i,j}=\{\begin{array}{cc}\hfill u_m& (b_i,b_j)\text{ is a Watson-Crick base pair}\hfill \\ \hfill u_{mm}& \text{otherwise}\hfill \end{array}$$
(1)
with $`u_m,u_{mm}>0`$. Here, $`u_m`$ sets the energy scale. The value of $`u_{mm}`$ is not essential as long as it is repulsive, since there is always the option to not bind at all (with energy “0”) rather than to bind with a repulsive energy.

For analytical work, it is convenient to take the $`\epsilon _{i,j}`$ to be independent Gaussian random variables specified by the mean $`\overline{\epsilon }`$ and variance

$$\overline{(\epsilon _{i,j}\overline{\epsilon })(\epsilon _{k,l}\overline{\epsilon })}=D\delta _{i,k}\delta _{j,l}.$$
(2)
Throughout the text, we will use the overline to denote averages over the ensemble of random pairing energies. Here, the parameter $`D`$ provides a measure of the strength of the randomness. It is an approximation to the model (1) in two respects: First, it replaces the discrete distribution of energies by a Gaussian distribution. Moreover, it neglects the correlations between $`\epsilon _{i,j}`$ and $`\epsilon _{i,k}`$ induced by the shared base $`b_i`$. We do not anticipate universal quantities to depend on the subtle differences in the statistics of the $`\epsilon _{i,j}`$’s. This will be tested numerically by comparing the scaling behavior produced by the two models.

Given the energy of each secondary structure, we can study the partition function $`Z(N)=_{S𝒮(N)}e^{\beta E[S]}`$ where $`𝒮(N)`$ comprises all valid secondary structures of a molecule of length $`N`$. This function can be conveniently computed in terms of the restricted partition function $`Z_{i,j}`$ for the substrand $`b_i\mathrm{}b_j`$. For the simple energy model $`E[S]`$ described above, $`Z_{i,j}`$ can be split up according to the possible pairings of position $`j`$ to yield the exact recursion relation higg96 ; dege68 ; wate78

$$Z_{i,j}=Z_{i,j1}+\underset{k=i}{\overset{j1}{}}Z_{i,k1}e^{\beta \epsilon _{k,j}}Z_{k+1,j1}.$$
(3)
This equation can then be iterated to compute the full partition function $`Z(N)=Z_{1,N}`$ for arbitrary interaction energies $`\epsilon _{i,j}`$’s in $`O(N^3)`$ time mcca90 ; wate78 . It also forms the basis of analytical approaches to the problem.

The molten phase: If sequence disorder does not play an important role, we may describe the RNA molecule by replacing all the binding energies $`\epsilon _{i,j}`$ by some effective value $`\epsilon _0`$. As we will see later, this is an adequate description of random RNA at high enough temperatures (but below denaturation.) Below we briefly review the properties of RNA in this high temperature “molten phase”.

Since the $`Z_{i,j}`$’s become translationally invariant in the absence of disorder, it is straightforward to solve Eq. (3) for the molten phase partition function $`Z_0(N)`$ using the $`z`$-transform. For large $`N`$, it has the form bund99 ; dege68 ; wate78

$$Z_0(N)\stackrel{N1}{}N^{\theta _0}\mathrm{exp}[\beta Nf_0(q)]$$
(4)
where $`qe^{\beta \epsilon _0}`$ is the only parameter, $`\theta _0=3/2`$ is a universal exponent, and $`f_0(q)=k_BT\mathrm{ln}(1+2\sqrt{q})`$ is the free energy per length.

One useful observable characterizing the state of the RNA is the free energy cost $`\mathrm{\Delta }F(N)`$ of pinching together one end (say $`i=1`$) and the mid-point ($`i=N/2`$) of the polymer relative to the unperturbed state. The pinch effectively separates the polymer into two pieces of length $`N/2`$. In the molten phase, we simply have

$$\mathrm{\Delta }F=k_BT\mathrm{ln}[Z_0(N)/Z_0^2(N/2)]\frac{3}{2}k_BT\mathrm{ln}N.$$
(5)
It reflects the loss of configurational entropy in the allowed secondary structures due to the pinch.

Numerics: Before we delve into the analysis, we first present some numerical evidence for a phase transition between the high- and low- temperature behavior for the ensemble of random RNA. To this end, we generate many configurations of interaction energies $`\epsilon _{i,j}`$’s according to both models (1) and (2), with $`u_m=u_{mm}=1`$ and $`\overline{\epsilon }=1/2`$, $`D=3/4`$ respectively. Then for a wide range of temperatures we calculate the pinching free energy, which is given by $`\mathrm{\Delta }F=k_\mathrm{B}T\mathrm{ln}(e^{\beta \epsilon _{1,N/2}}Z_{2,N/21}Z_{N/2+1,N}/Z_{1,N})`$ in terms of the quantity $`Z_{i,j}`$ in (3) for sequences without translational invariance. The result is then averaged over $`1000`$ to $`10000`$ disorder configurations and illustrated in Fig. 2 for two representative temperatures.

At high temperature ($`k_BT=2`$), the pinch energy follows precisely the molten phase behavior expected according to Eq. (5) up to an additive constant. Thus, at high temperatures the molten phase description is applicable even if the interaction energies $`\epsilon _{i,j}`$ are not uniform. Moreover, we find no difference between the two models of disorder at $`k_BT=2`$. At low temperature ($`k_BT=0.025`$), the picture is different. The length dependence of $`\mathrm{\Delta }F`$ is still logarithmic. However, the prefactors differ between the two disorder choices and they are both by far larger than the prefactor $`0.0375`$ expected if the molten phase result is extrapolated to this temperature. This suggests that there is a distinct low temperature phase; this finding is reinforced by similar changes detected in other observables bund01 .

High temperature behavior: Now we will establish the stability of the molten phase against weak disorder by a perturbative analysis. Assuming that the specific choice of disorder does not matter at weak disorder (or high temperature) as supported by the numerical results above, we will use the uncorrelated Gaussian disorder characterized by Eq. (2) for the purpose of this analysis. The behavior of the system at weak disorder is determined by the lowest order terms in the perturbative expansion of the ensemble averaged free energy in the disorder strength $`D`$. This term is given by the two-replica partition function $`\overline{Z^2(N)}`$ of two RNA molecules sharing the same disorder. With the uncorrelated Gaussian energies (2), the ensemble average can be explicitly performed, yielding after some algebra

$$\overline{Z^2(N)}=_{\{S_1,S_2𝒮(N)\}}q^{|S_1|}q^{|S_2|}\stackrel{~}{q}^{|S_1S_2|}$$
(6)
where $`q\mathrm{exp}\left(\beta \overline{\epsilon }+\frac{1}{2}\beta ^2D\right)`$ and $`\stackrel{~}{q}\mathrm{exp}\left(\beta ^2D\right)`$ are the two relevant “Boltzmann factors”, and $`|S_i|`$ and $`|S_iS_j|`$ are the number of bases contained in structure $`S_i`$ or common to $`S_i`$ and $`S_j`$ respectively. This effective partition function has a simple physical interpretation: It describes two RNA molecules subject to a homogeneous attraction with effective interaction energy $`\epsilon _0\beta ^1\mathrm{ln}q=\overline{\epsilon }\frac{1}{2}\beta D`$ between any two bases of the same molecule. In addition, there is an inter-replica attraction characterized by the factor $`\stackrel{~}{q}`$ for each bond shared between the two replicas. The inter-replica attraction is induced by the same disorder shared by the replicas. It can potentially force the replicas to “lock” together, i.e., to become correlated.

It turns out that this two-replica problem can be solved exactly bund01 . The key idea leading to the solution is that the sum over the pairs of secondary structures in Eq. (6) can be reordered by first summing over the bonds which are common to the two structures, noting that the possible configurations of the common bonds are themselves the set $`𝒮(N)`$ of valid secondary structures. Within a given configuration of common bonds, all possibilities to place non-common bonds in the two individual structures can then be summed over, leading to an effective single RNA problem. However, the necessary algebra is quite involved bund01 ; here we just quote the results. The solution has the form $`\overline{Z^2(N)}N^\theta \mathrm{exp}[\beta Nf(q,\stackrel{~}{q})]`$, with two different expressions for $`\theta `$ and $`f`$ depending on whether $`\stackrel{~}{q}`$ is above or below a critical value $`\stackrel{~}{q}_c=1+1/[q^2_{N=1}^{\mathrm{}}Ng^2(N)]`$, where $`g(N)=Z_0(N)/(1+2\sqrt{q})^{N1}N^{3/2}`$ for large $`N`$.

For $`\stackrel{~}{q}<\stackrel{~}{q}_c`$, we have $`\theta =2\theta _0=3`$. In addition, $`f(q,\stackrel{~}{q})`$ is basically a modified version of $`2f_0(q)`$. In this regime, the two-replica partition function $`\overline{Z^2(N)}`$ is essentially a product of two single-replica partition functions $`Z_0(N)`$. Since there is no coupling between the two replicas, we conclude that the effect of disorder is irrelevant in this regime. For $`\stackrel{~}{q}>\stackrel{~}{q}_c`$, we have $`\theta =3/2`$ and $`f`$ given by a complicated function of $`Z_0(N)`$ bund01 . Here, the two-replica partition function has the same asymptotic form as that of the single-replica system in (4). This implies that the disorder coupling locks the two replicas together.

Of course, as already explained above, only the weak-disorder limit (i.e., $`\stackrel{~}{q}1`$) of the two-replica problem can be applied to the full random RNA problem. While $`\stackrel{~}{q}_c`$ itself depends on the disorder strength $`\beta ^2D`$, it converges in the weak disorder limit ($`\beta ^2D1`$) towards a constant $`\stackrel{~}{q}_c(\beta ^2D=0)>1`$. Therefore, for $`\beta ^2D1`$ we always have $`\stackrel{~}{q}<\stackrel{~}{q}_c`$. The molten phase is perturbatively stable upon the introduction of disorder, making it the appropriate description of random RNA at high temperatures.

Low temperature behavior: Next we determine whether the molten phase persists for all temperatures down to $`T=0^+`$. In the following, we will assume that long random RNA is in the molten phase for all temperatures, i.e., that the partition function for any substrand of large length is given by Eq. (4), with some effective value of $`q`$. Then, we will show that this assumption leads to a contradiction below some temperature $`T^{}>0`$. This contradiction implies that the molten phase description breaks down at some finite $`T_cT^{}`$. To be specific, we will consider the sequence disorder model (1) in this analysis.

We will again focus on the pinching free energy $`\mathrm{\Delta }F`$. Under the assumption that the random sequences are described by the molten phase, it is given by Eq. (5) for large $`N`$ and all $`T`$ independently of the effective value of $`q`$. On the other hand, we can study this pinching free energy for each given sequence of bases drawn from the ensemble of random sequences. For each such sequence, we can look for a continuous segment of $`\mathrm{}N`$ Watson-Crick pairs $`(b_i,b_j)(b_{i+1},b_{j1})\mathrm{}(b_{i+\mathrm{}1},b_{j\mathrm{}+1})`$ where the bases $`b_i\mathrm{}b_{i+\mathrm{}1}`$ are within the first half of the molecule and the bases $`b_{j\mathrm{}+1}\mathrm{}b_j`$ are in the second half. For random sequences, the probability of finding such complementary segments decreases exponentially with the length $`\mathrm{}`$ in a given region, with the largest $`\mathrm{}`$ in a sequence of length $`N`$ being typically $`\mathrm{ln}N/\mathrm{ln}2`$ arra90 .

Now we calculate the pinching free energy $`\mathrm{\Delta }F=F_{\mathrm{pinched}}F_{\mathrm{unpinched}}`$ by evaluating the two terms separately. The partition function corresponding to the unpinched free energy contains at least all the configurations in which the two complementary segments $`b_i\mathrm{}b_{i+\mathrm{}1}`$ and $`b_{j\mathrm{}+1}\mathrm{}b_j`$ are completely paired. Thus,

$$F_{\mathrm{unpinched}}F_{\mathrm{paired}}$$
(7)
where $`F_{\mathrm{paired}}`$ is the free energy of the ensemble of structures in which the two complementary segments are paired. The latter is the sum of the energy of the paired segments and those of the two remaining substrands of lengths $`L_1=ji2\mathrm{}+1`$ and $`L_2=N+ij1`$, i.e.,

$$F_{\mathrm{paired}}=\mathrm{}u_m+(N2\mathrm{})f_0+\frac{3}{2}k_BT[\mathrm{ln}(L_1)+\mathrm{ln}(L_2)].$$
(8)
The free energy $`F_{\mathrm{pinched}}`$ is, by the assumption of the molten phase, the interaction energy of the pinched base pair $`(b_1,b_{N/2})`$ plus the free energy of the two substrands $`b_2\mathrm{}b_{N/21}`$ and $`b_{N/2+1}\mathrm{}b_N`$. According to Eq. (4), this is $`F_{\mathrm{pinched}}=f_0N+2\times \frac{3}{2}k_BT\mathrm{ln}N`$ up to terms independent of $`N`$. Combining this with Eqs. (7) and (8), and noting that $`\mathrm{}`$ is typically of the order $`\mathrm{ln}N/\mathrm{ln}2`$ and $`L_1`$, $`L_2`$ are typically proportional to $`N`$, we finally obtain $`\overline{\mathrm{\Delta }F}[u_m+2f_0]\mathrm{ln}N/\mathrm{ln}2`$ for very large $`N`$. This is only consistent with Eq. (5) if

$$\frac{3}{2}k_BT[u_m+2f_0]/\mathrm{ln}2.$$
(9)
Since $`f_0`$ is a free energy per length, the right hand side of this equation is a decreasing function of temperature. Moreover, $`f_0(T=0)`$ is the average free energy per length of the minimum energy secondary structures of random sequences. The lowest possible energy any structure of an RNA molecule of $`N`$ bases can achieve is $`\frac{N}{2}u_m`$, i.e., if all bases form Watson-Crick pairs. However, if the sequence disorder leads to frustration, there is always a finite fraction of bases which cannot be incorporated into Watson-Crick base pairs 111This is not the case for the sequence disorder model at an alphabet size $`2`$ unless a minimum hairpin length constraint is introduced as in Ref. pagnani . and thus $`f_0(T=0)>u_m/2`$. Therefore, the right hand side of (9) is a decreasing function of $`T`$ starting at some positive value at $`T=0`$. It follows that there is some unique temperature $`T^{}`$ below which the consistency condition (9) breaks down, implying the inconsistency of the molten phase assumption in this regime. From this we conclude that there must be a phase transition away from the molten phase at some critical temperature $`T_cT^{}>0`$. Numerically, we find that $`f_0(T=0)0.46`$ for $`u_m=1`$, yielding $`k_BT^{}0.075`$. This is consistent with the numerically observed change between the low- and high- temperature behaviors at $`k_BT_c0.25`$ bund01 . Improved estimates of $`T_c`$ can be made by relaxing the condition of perfect complementarity of the two segments; this as well as a detailed characterization of the low temperature phase will be discussed elsewhere.

Conclusions: We studied the behavior of random RNA sequences in the regime below the denaturation transition. A two-replica calculation shows that disorder is perturbatively irrelevant, i.e., an RNA molecule with weak sequence disorder is in the molten phase where many secondary structures with comparable total energy coexist. By further considering the rare regions of strong sequence complementarity, we show that the molten phase cannot exist down to arbitrarily low temperatures, and thereby establish the existence of a distinct low temperature phase below some critical temperature $`T_c>0`$. Our analysis follows the approach used in Ref. tang00  to show the relevance of disorder in the denaturation of double-stranded DNA. It will be interesting to see whether a renormalization group theory along the line of Ref. tang00  may be constructed to elucidate the low-temperature phase of the random RNA problem.

Acknowledgments: The authors benefited from discussions with U. Gerland, D. Moroz, and especially L.-H. Tang regarding the role of rare regions. This work is supported by the NSF through grants no. DMR-9971456 and DBI-9970199.