text:
# The bifurcation diagrams for the Ginzburg-Landau system for superconductivity

## 1 Introduction

In this paper, we want to give a detailed description of the bifurcation diagrams of an infinite superconducting cylinder of cross section $`\mathrm{\Omega }`$, submitted to an exterior magnetic field $`h_0`$. The response of the material varies greatly according to the value of $`h_0`$, the size of the cross section and the Ginzburg-Landau parameter $`\kappa `$ that characterizes the material: superconductivity appears in the volume of the sample for low fields and small samples, under the form of vortices for higher fields, bigger samples and larger values of $`\kappa `$, and is destroyed for high fields. The type of response of a superconducting material has been studied numerically and theoretically by various authors in various asymptotic regimes . Here, we want to give a complete numerical picture of the bifurcation diagrams for all values of the parameters. This type of computation has already been done in dimension 1 by one of the authors using auto . Here, the computation is made in a 2 dimensional domain using numerical solutions of the well-known Ginzburg-Landau model , based on a code first developed in . We examine the behavior of the energy, the magnitude of the order parameter and the magnetization versus the magnetic field for various solution branches. We also provide some analysis on the detailed findings such as the stability of solutions.

The paper is organized as follows: the Ginzburg-Landau model is briefly stated in Section 2 and the main features of the numerical codes used in the computation are described Section 3. The complete phase diagrams are given in Section 4, along with detailed analysis. A conclusion is given in Section 5.

## 2 The Ginzburg-Landau model

Let $`\mathrm{\Omega }`$ denote a two dimensional bounded domain which represents the cross section of a three dimensional cylinder occupied by the superconducting sample. Assume that the cylinder is homogeneous along its axis and a constant applied field $`H_0`$ is placed along the axis direction as well. Then, the Gibbs free energy $`𝒢`$ may be written in the following form :

$$\begin{array}{c}𝒢(\psi ,A)=_\mathrm{\Omega }\left(f_n+\alpha |\psi |^2+\frac{\beta }{2}|\psi |^4\right)𝑑\mathrm{\Omega }\hfill \\ +_\mathrm{\Omega }\left[\frac{1}{2m_s}\left|\left(i\mathrm{}+\frac{e_sA}{c}\right)\psi \right|^2+\frac{|\text{curl }A|^2}{8\pi }\frac{\text{curl }AH_0}{4\pi }\right]𝑑\mathrm{\Omega }.\hfill \end{array}$$

Here, $`\psi `$ is the (complex-valued) order parameter, $`A`$ is the magnetic potential, $`\text{curl }A`$ is the magnetic field, $`\alpha `$ and $`\beta `$ are constants (with respect to the space variable $`x`$) whose values depend on the temperature, $`c`$ is the speed of light, $`e_s`$ and $`m_s`$ are the charge and mass, respectively, of the superconducting charge-carriers, and $`2\pi \mathrm{}`$ is Planck’s constant.

After proper nondimensionalization, we can reformulate the free energy functional as:

$$𝒢(\psi ,A)=_\mathrm{\Omega }\left|\left(iA\right)\psi \right|^2+\frac{\kappa ^2}{2}\left(1|\psi |^2\right)^2+|\text{curl }Ah_0|^2d\mathrm{\Omega },$$

where $`\kappa `$ is the Ginzburg-Landau parameter representing the ratio of the penetration depth and the coherence length, $`h_0`$ the applied magnetic field and $`d`$ the characteristic size of the domain $`\mathrm{\Omega }`$, that is $`\mathrm{\Omega }=dD`$ where $`D`$ is a fixed domain.

The system that we are going to study is the following Ginzburg-Landau equations derived as the Euler-Lagrange equations for the minimizers of the functional $`𝒢`$ :

$$\{\begin{array}{c}(iA)^2\psi =\kappa ^2\psi (1|\psi |^2)\text{ in }\mathrm{\Omega },\hfill \\ \text{curl curl }A=|\psi |^2A+\frac{i}{2}(\psi ^{}\psi \psi \psi ^{})\text{ in }\mathrm{\Omega },\hfill \end{array}$$
(1)
which are supplemented by the boundary conditions

$$\{\begin{array}{c}(\psi iA\psi )\text{n}=0\text{ on }\mathrm{\Omega },\hfill \\ \text{curl }A=h_0\text{ on }\mathrm{\Omega },\hfill \end{array}$$

and gauge constraints

$$\{\begin{array}{c}\text{div }A=0\text{ in }\mathrm{\Omega },\hfill \\ A\text{n}=0\text{ on }\mathrm{\Omega }.\hfill \end{array}$$

Here, $`\mathrm{\Omega }`$ is the boundary of $`\mathrm{\Omega }`$ and n is its unit outer normal. With the above nondimensionalization, $`|\psi |`$ takes values between $`0`$ and $`1`$: the normal state corresponds to $`|\psi |=0`$ while the Messiner state corresponds to $`|\psi |=1`$.

This Ginzburg-Landau model has a special family of solutions called the normal solutions: $`\psi =0`$ and $`\text{curl }A=h_0`$, which correspond to the situation where superconductivity is destroyed. According to the values of the different parameters $`\kappa `$, $`d`$ and $`h_0`$, the system may have other solutions: superconducting solutions, for which $`\psi `$ is never 0 and vortex solutions for which $`\psi `$ has isolated zeroes (see Figure 1 and 2). For a complete introduction to the topic, one may refer to .

## 3 Numerical codes

The numerical code that we are using is based on the finite element approximations of the Ginzburg-Landau model, first proposed in  and later used in many settings, see  for instance. The codes have many different variants that can be used to simulate the solutions of the Ginzburg-Landau models in the high-$`\kappa `$ setting , thin films with variable thickness , samples with normal inclusions  as well as layered materials. Here, we choose the standard version that solves the Ginzburg-Landau equations on a rectangular domain $`\mathrm{\Omega }`$ and $`d`$ will denote the characteristic size of the rectangle. We note that other numerical methods, such as the gauge invariant difference methods , i.e, the so-called bond-and-link variable methods or the method of eigen-functions, have also been used to compute the phase diagrams for the Ginzburg-Landau equations .

In our implementation of the finite element approximation, we use a uniform triangular grid with piecewise quadratic polynomials for both $`\psi `$ and $`A`$. It is shown that (see  for instance), if $`(\psi _\delta ,A_\delta )`$ is the finite element solution on a given mesh with mesh size $`\delta `$, the convergence of the approximation is assured, and asymptotically, we have

$$\psi \psi _\delta _2+AA_\delta _2=O(\delta ^3)$$

where $`_2`$ denote the standard mean square $`L^2`$ norm. For each set of calculation, we refine the mesh size until numerical convergence is evident.

In Figures 1 and 2, we present a few typical plots for the numerical solutions of the equation (1). For each solution, the plots include a surface plot of the magnitude of the order parameter, a surface plot of the magnetic field given by $`\mathrm{curl}A`$ and a vector plot of the superconducting current. In Figure 1, we have a solution with a single vortex at the center of the domain which corresponds to the parameter values $`\kappa =0.23`$, $`d=16.8`$ and $`h_0=0.563`$.

In Figure 2, we present the plots for a solution with two vortices corresponding to $`\kappa =0.8`$, $`d=4`$ and $`h_0=1.2`$.

For fixed $`\kappa `$, $`d`$ and $`h_0`$, we are interested in finding the number of solutions of (1) and their stability. A continuation in the parameter spaces is used for getting solutions with different parameter values. With $`\kappa `$, $`d`$ given but $`h_0`$ allowed to vary, for a computed solution branch, we plot $`\psi _{\mathrm{}}`$, the maximum magnitude of the order parameter $`\psi `$ in the domain, the free energy $`𝒢`$ and the magnetization versus the applied field $`h_0`$. These phase diagrams or bifurcation diagrams will give us information on the solutions (number and stability) for each $`h_0`$. They were drawn by Ginzburg  in some limiting cases of the parameters. Here we want to give a more complete description of these diagrams for all values of the parameters.

The results of our numerical computations allow us to separate the $`\kappa `$-$`d`$ plane into different regions depending on the shape of the bifurcation diagram. A solution branch may exist in one region but may cease to exist in another one, it may also develop hysteresis in some regions.

## 4 The bifurcation diagrams

It is well known that for large fields ($`h>h_{}`$), the only solution is the normal solution . For smaller fields ($`h<h_{}`$), the normal solution always exists but there are other solutions which display four different types of behaviors. These behaviors depend on the values of $`\kappa `$ and $`d`$. In Figure 3, we have plotted four curve segments $`\{\kappa _i(d)\}_1^4`$ separating the $`\kappa d`$ plane into four regions $`\{R_i\}_1^4`$.

The four curve segments in Figure 3 are the results of the computation described in the earlier section. All four curves meet close to $`\kappa =1/\sqrt{2}`$, $`\kappa _2(d)`$ is of the form $`2.112/d`$, $`\kappa _3(d)`$ is tending to 0.4 at infinity. The four regions correspond to the four types of behaviors for the bifurcation diagrams. For convenience, for each $`i=1,2,3,4`$, we use $`d=d_i(\kappa )`$ to denote the inverse function of the function $`\kappa =\kappa _i(d)`$ wherever the inverse is well-defined. What distinguishes the different regions are features like the existence (or the lack of existence) of vortex solutions, the global and local stability of solutions, and the hysteresis phenomena.

We now provide detailed descriptions of the solution behavior for each region in Figure 3.

Region 1: $`d<d_1(\kappa )`$ and $`d<d_2(\kappa )`$.

This corresponds to the situation where the cross section of the superconducting sample is small enough. The bifurcation diagram is illustrated in Figure 4. The corresponding plot of the energy is given in Figure 5 and the magnetization curve in Figure 6.

Throughout this region, there is a unique non-normal solution for $`h<h_{}`$. This solution is a superconducting solution which is the global minimizer of the free energy $`𝒢`$. The curve $`|\psi |_{\mathrm{}}`$ against $`h`$ is monotonically decreasing. When increasing the field, the magnitude of the superconducting solution decreases until it turns normal at $`h_0=h_{}`$. Conversely, when decreasing the field, the normal solution turns superconducting also for $`h_0=h_{}`$. The transition to the normal solution is of second order, that is the energy of the superconducting solution tends to the energy of the normal solution at the transition and there is no hysteresis phenomenon. This can also be seen on the magnetization curve of Figure 6.

There is no vortex solution for the parameters $`(d,\kappa )`$ in this region. This reflects the fact that $`d`$ is too small to allow enough room for a vortex to exist since a vortex core is of typical size $`C/\kappa `$.

In , it is proved that for $`d<\mathrm{min}\{d_0,d_1/\kappa \}`$, a solution of the Ginzburg-Landau equations has no vortex, and for such a solution, $`\psi `$ is almost constant in the domain and is unique up to multiplication by a constant of modulus 1. Also it is proved that the bifurcation curve for $`|\psi |_{\mathrm{}}`$ is decreasing and that $`h_{}`$ is asymptotically $`C/d`$ when $`d`$ tends to 0. More precisely, for a disc, when $`d`$ tends to 0, the limiting curve is $`\psi _{\mathrm{}}^2=1d^2h_0^2/2\kappa ^2`$.

Region 2: $`\kappa >\kappa _2(d)`$ and $`\kappa >\kappa _3(d)`$ .

In contrast with Region 1, this region corresponds to the situation where the typical size of vortices ($`C/\kappa `$) is small enough compared to the size of the domain. This region displays the typical type II behavior of superconductors. The bifurcation diagram is illustrated in Figure 7, the energy in Figure 8 and the magnetization curve is given in Figure 9.

It has been well understood both physically and in more recent years mathematically that, for sufficiently large $`\kappa `$, there are vortex solutions which are the global minimizers of the free energy for a certain range of fields. The number of vortices depends on the strength of the applied field. The maximum number of vortices increases with $`d`$ and $`\kappa `$, a type-II behavior.

For very small fields, the global minimizer is the superconducting solution (solid line). As the field is increased, the superconducting solution loses its global stability ($`h=h_{c_1}`$ and for even larger fields loses its local stability, $`h=h_0^{}`$). Then the global minimizer starts to nucleate vortices. In Figures 7, 8, 9, the solution branch corresponding to a one vortex solution is illustrated by a dotted line and to a two vortex solution by a dashed line.

Let $`n`$ be a non-negative integer (the flux quanta). Hysteresis phenomena occur in our experiments when going from a solution with $`n`$ vortices to a solution with $`n+1`$ vortices: when increasing the field, the solution with $`n`$ vortices remains locally stable, though it is no longer a global minimizer until $`h_n^{}`$ where it jumps to a solution with $`n+1`$ vortices. Similarly when decreasing the field, the solution with $`n+1`$ vortices remains locally stable down to $`h_{}^{n+1}<h_n^{}`$. Such a hysteresis pattern is typical for parameters in this region. The fields $`h_n^{}`$ and $`h_{}^{n+1}`$ correspond to the critical field where the vortex solution loses its local stability and their values can be identified from Figures 7 and 8 (transition from superconducting solution to 1 vortex solution and from one vortex to two vortices).

Nevertheless, the transition from normal to the last vortex solution is of second order since the vortex solution is stable, as can be seen from the energy and magnetization curves. This transition occurs at a field that is usually called $`h_{c_3}`$ in the literature.

There have been many numerical works on these vortex solutions: see for instance . Recently, a rigorous analysis of vortices, their number and the critical thermodynamic fields in the high kappa limit has been done by Sandier and Serfaty : in particular they obtain in the high kappa limit that $`h_{c_1}`$ (the field for which the energy of the superconducting solution is equal to the energy of the one vortex solution) is of order $`C\mathrm{log}\kappa `$ and the field $`h_n^{}`$ is asymptotic to $`h_{c_1}+n\mathrm{log}\mathrm{log}\kappa `$.

The onset of superconductivity in decreasing fields (instability of normal solutions and computation of the fields of nucleation) has been analyzed by Bernoff and Sternberg  and del Pino, Felmer, Sternberg . Their works provide, as $`d`$ and $`\kappa `$ tends to $`\mathrm{}`$, an asymptotic development of $`h_{c3}`$, the field at which the normal solution bifurcates to a vortex solution. This is what is called surface superconductivity. A linearizing the Ginzburg-Landau equation has been done near the normal solution. In the high kappa limit, their computation yields

$$h_{c_3}\frac{\kappa ^2}{\lambda _1}+C\kappa \kappa _{max}+o(\kappa )$$
(2)
where $`\kappa _{max}`$ is the maximal curvature of the domain and $`\lambda _1`$ is the first eigenvalue of the linearized problem and is approximately equal to $`0.59`$. In the high $`d`$ limit, it yields

$$h_{c_3}\frac{\kappa ^2}{\lambda _1}+\frac{C\kappa \kappa _{max}}{d^2}+o(\frac{1}{d^2}).$$
(3)
This expansion is consistent with the work of Saint James and de Gennes  who got the first term of this expansion in the case of an infinite plane in one dimension. In two dimensions, one has to take into account the curvature of the cross section. In the case of the disc, the equivalent of expansions (2) and (3) have been carried out by  in the limit $`\kappa d`$ large.

The hysteresis phenomenon has been rigorously analyzed in  for the superconducting solution ($`n=0`$).

Some interesting questions to be addressed here are: what are the critical fields for the vortex solutions losing their local stability ($`h_n^{}`$ and $`h_{}^n`$)? can we find them in some asymptotic limit such as $`d`$ large or $`\kappa `$ large or $`\kappa d`$ large? and can we prove the existence of these hysteresis phenomena for the general case?

Curve $`\kappa _2(d)`$:

For fixed $`\kappa `$ above 0.7, when $`d`$ is increased from 0, the point $`(d,\kappa )`$ is first in Region 1. Then it reaches the critical value $`d_2(\kappa )`$. For $`d=d_2(\kappa )`$, the bifurcation diagram $`\psi _{\mathrm{}}`$ vs $`h`$ is decreasing and the superconducting solution bifurcates from the normal solution at $`h=h_{}`$. For $`d`$ a little bigger than $`d_2(\kappa )`$, there is a vortex solution bifurcating from the normal solution close to $`h_{}`$. Hence for $`d=d_2(\kappa )`$, at $`h=h_{}`$, the linearized problem near the normal state has two eigenfunctions: one without vortices and one with a vortex. This is how uniqueness of solution is lost when increasing $`d`$, though the assertion needs to be proved mathematically.

Let $`D`$ be the fixed domain such that $`\mathrm{\Omega }=dD`$. Then a bifurcated solution near the normal state $`(0,h_0a_0)`$ (where $`a_0`$ is such that $`\text{curl }a_0=1`$ in $`D`$ and $`a_0n=0`$ on $`\mathrm{\Omega }`$) is of the form $`(\epsilon \varphi ,h_0a_0+\epsilon B)`$. Let $`\omega =h_0d`$. The second variation of $`𝒢`$ near the normal state is

$$\frac{𝒢^2}{\epsilon ^2}(\epsilon \varphi ,h_0a_0+\epsilon B)=\frac{1}{d^2}_D|(i+\omega ^2a_0)\varphi |^2\kappa ^2d^2|\varphi |^2+|\text{curl }B|^2.$$

Let

$$\lambda (\omega )=inf(_D|(i+\omega ^2a_0)\varphi |^2,\varphi _{L^2}=1,\varphi H^1(D,\text{ }\mathrm{C})).$$

Hence if $`\lambda (\omega )>\kappa ^2d^2`$, the normal solution is stable, if $`\lambda (\omega )<\kappa ^2d^2`$, the normal solution is unstable and if $`\lambda (\omega )=\kappa ^2d^2`$, degenerate stability occurs. For the eigenvalue $`\lambda (\omega )=\kappa ^2d^2`$, bifurcation of non normal solutions occurs. Thus, one has to study

$$\{\begin{array}{c}(i\omega ^2a_0)^2\varphi =\lambda (\omega )\varphi \text{in}D,\hfill \\ \frac{\varphi }{n}=0\text{on}D,\hfill \end{array}$$
(4)
with $`\lambda (\omega )=\kappa ^2d^2`$. For most values of $`\kappa `$ and $`d`$, the field $`h`$ such that the first eigenvalue $`\lambda (\omega )`$ is equal to $`\kappa ^2d^2`$ yields a single eigenfunction. In region 1, the eigenfunction has no vortex while in region 2, the eigenfunction has a vortex. Thus, the curve $`\kappa _2(d)`$ corresponds to those values of $`\kappa `$ and $`d`$ for which the eigenvalue has two different eigenfunctions, one without vortices and one with a vortex. That is, on $`\kappa _2(d)`$, the vortex state starts to exist.

This situation has been studied in the case of a ball in . See also  for more recent developments. In this case, the solutions of (4) with $`n`$ vortices are of the form $`\xi _n(r)exp(in\theta )`$ and have eigenvalue $`\lambda (\omega ,n)`$. In particular in  they draw the function $`\lambda (\omega ,n)`$ versus $`\omega `$. The curves $`\lambda (\omega ,0)`$ and $`\lambda (\omega ,1)`$ intersect exactly once for $`\omega =\omega _0`$ and $`\lambda =\lambda _0`$. Because of the bifurcation condition $`\lambda (\omega )=\kappa ^2d^2`$, it implies that $`\kappa ^2d^2=\lambda _0`$, hence the curve $`\kappa _2(d)`$ is of the form $`\kappa d=constant`$. It would be interesting to give a rigorous proof that the curves $`\lambda (\omega ,0)`$ and $`\lambda (\omega ,1)`$ intersect only once for the case of the disc and for the case of a more general domain. In region 1, that is below $`\kappa _2(d)`$, the first eigenfunction is simple and leads to a solution without vortex. In region 2, that is above $`\kappa _2(d)`$, the first eigenfunction is simple and leads to a solution with one vortex, but we expect that there is also an eigenfunction with no vortex for a lower field $`h`$.

Similarly, the curves $`\lambda (\omega ,n)`$ and $`\lambda (\omega ,n+1)`$ also intersect only once on the numerics of  which means in our setting that there are curves $`\kappa d=C_n`$ at which the eigenvalue has two eigenfunctions with $`n`$ and $`n+1`$ vortices. Above $`\kappa d=C_n`$, a solution with $`n+1`$ vortices starts to bifurcate from the normal state and below it, a solution with $`n`$ vortices starts to bifurcate from the normal state, so that the curve $`\kappa d=C_n`$ are the critical curves for the existence of $`n+1`$ vortices.

In the general case that we are studying, it is totally open to prove that there is a unique value of $`\omega `$ such that $`\lambda (\omega )`$ has two eigenvalues, one with a vortex and one without. This would yield to $`\kappa _2(d)=C/d`$, which is what we have found numerically. Moreover, we observe that the field of bifurcation $`h_{}`$ satisfies $`h_{}d=\omega `$ hence is constant along $`\kappa _2(d)`$.

Taking this analysis of bifurcation a little further allows us to define

$$(\kappa ,d)=\{h,s.t.\lambda (\sqrt{hd})=\kappa ^2d^2\}.$$

In region 1, we expect that $``$ has a single element while in Region 2, we expect this set to have several elements corresponding to the various branches of solutions with several vortices bifurcating from the normal state. But this analysis is open even in the case of the disc.

Region 3: $`\kappa _4(d)<\kappa <\kappa _3(d)`$.

For parameters in this region, that is large domains and intermediate $`\kappa `$ (in a relative sense), a typical phase diagram is illustrated in Figure 10 with the energy in Figure 11 and the magnetization in Figure 12. Three solution branches are shown which represent the normal solution, the superconducting solution (solid line) and a solution with a single vortex (dashed line). A profile for one of the vortex solutions is given in Figure 1.

The superconducting solution displays a hysteresis behavior as before, but when increasing the field, it turns normal instead of going on the vortex branch as in Region 2. More precisely, as the field is increased, the superconducting solution loses its global stability, then its local stability and its drops to the normal branch when the super-heating field is reached. Conversely when decreasing the field from the normal state, the solution gets on the vortex branch, though it is only a local minimizer. So the transition when going down the field is of second order as can be seen on the magnetization curve (Figure 12), but when going up, there is a hysteresis. When decreasing the field further, the solution jumps to the superconducting branch. The vortex solution is never a global minimizer. In fact, our numerical experiments indicate that there are only locally stable vortex solutions.

We note that, for values of the Ginzburg-Landau parameter $`\kappa `$ in this region, the vortex state has not been frequently studied in the literature, except for superconducting sample with extreme geometrical conditions such as thin films, disks and rings. In the latter cases, the material displays typical type-II behavior for all ranges of $`\kappa `$ as the Ginzburg-Landau models can be simplified to allow an almost uniform penetration of the magnetic field . However, the current study is done for three dimensional infinite cylinders and the simplified models are not directly applicable. In fact, from the plot of the magnetic field given in Figure 1, we see that there is considerable variation in the field strength over the cross section.

Curve $`\kappa _3(d)`$:

Let us call $`H_c`$ the thermodynamic critical field introduced by Ginzburg : the energy of the superconducting solution is equal to the energy of the normal solution at this field (in our nondimensionalization, it means that the energy of the superconducting solution is zero). The curve $`\kappa _3(d)`$ corresponds to the situation where there is a small amplitude vortex solution bifurcating from the normal solution exactly at $`H_c`$. One could hope to determine this curve mathematically.

We notice that as $`d`$ tends to infinity, $`\kappa _3(d)`$ tends to a finite limit close to 0.4 which is less than $`1/\sqrt{2}`$. Using (3), one can get that in the high $`d`$ limit, $`\kappa `$ is close to $`\lambda _1H_c`$, where $`\lambda _1`$ is the first eigenvalue of the linearized problem (4) in an infinite domain and is close to 0.59. In Ginzburg’s computations  (see also ), $`H_c`$ is exactly $`1/\sqrt{2}`$ in the limit $`d=\mathrm{}`$. This gives $`\lambda _1/\sqrt{2}`$ as a limit for $`\kappa _3(d)`$ at infinity, which is close to $`0.42`$. A rigorous mathematical justification of this asymptotic behavior remains to be provided.

Let us recall that the critical value of $`\kappa `$ equal to $`1/\sqrt{2}`$ that separates type I and type II superconductors was discovered by Ginzburg studying the bifurcation of the superconducting  solution from the normal state for an infinite domain. Then Saint James and de Gennes discovered surface superconductivity, that is superconductivity is nucleated first in the surface and thus appears for higher fields than was calculated by Ginzburg. Here, we see that for $`\kappa `$ less than $`\kappa _3(d)`$, the vortex solution is no longer stable but is locally stable near the bifurcation.

Region 4: $`\kappa <\kappa _4(d)`$ and $`d>d_1(\kappa )`$.

For parameters in this region, that is $`\kappa `$ small but domains large enough, the typical bifurcation diagram is illustrated in Figure 13 with the energy in Figure 14 and the magnetization in Figure 15.

There are superconducting solutions displaying a hysteresis phenomenon and no locally stable vortices. The superconducting solution is not always the global minimizer, but when increasing the field, the sample remains superconducting until reaching a super-heating field $`h^{}`$, where the solution becomes normal with a discontinuous transition. Similarly, when decreasing the field, the sample stays normal until the field $`h_{}`$ which is less than $`h^{}`$, where it turns superconducting by a discontinuous transition. Mathematically, we believe that there is a range of fields, between $`h_{}`$ and $`h^{}`$ where there are multiple superconducting solutions. This is in analogy with what happens for the one dimensional case. For more rigorous analysis of the one dimensional models, we refer to .

Next, we also note that there is no locally stable vortex solution in the region. It is well-known that asymptotically for small $`\kappa `$, the vortex solution is not energetically favorable, and the material belongs to the typical type-I regime  where the phase transition is characterized by the superconducting/normal interface rather than the vortex state. When varying the field, the superconducting or normal solutions will not turn into a one vortex solution. However, if we do continuation from a vortex solution for a bigger value of $`\kappa `$ (in Region 3), and decrease $`\kappa `$, we can still find existence of solutions with vortices when we reduce $`\kappa `$ to values in Region 4 despite of the instability of vortex solutions. For even smaller $`\kappa `$, continuations in $`\kappa `$ or in other parameters from the vortex solutions fails to produce any new vortex solutions. We believe that when decreasing $`\kappa `$, the vortex solution first loses its stability near the normal solution (on the curve $`\kappa _4(d)`$, but it remains locally stable for $`\psi _{\mathrm{}}`$ a little higher in the branch. For very small $`\kappa `$, (especially less than $`C/d`$), we believe that there is no vortex solution at all. This has been proved in  in the case of a disc. It is an open problem to prove that for fixed $`d`$ and for $`\kappa `$ small enough, vortex solutions do not exist.

Curve $`\kappa _1(d)`$:

If $`\kappa `$ is fixed below 0.7 and $`d`$ is increased, then the curve $`d_1(\kappa )`$ is crossed. It remains establish the mathematical existence of this curve. In the particular case where $`\kappa `$ is very small, the order parameter $`\psi `$ is almost constant, there are no vortices in the domain so that system (1) simplifies to

$$\mathrm{\Delta }A=|\psi |^2A,\text{in}\mathrm{\Omega }\text{curl }A=1\text{on}\mathrm{\Omega }$$

where $`|\psi |`$ is a constant that depends on $`h_0`$. The boundary condition for $`\psi `$ yields $`|\psi |(|\psi |^2+h_0^2A^21)=0`$, which is the equation of the bifurcation curve. One has to find the critical $`d`$ for which the curve $`|\psi |(h_0)`$ changes direction of bifurcation near $`|\psi |=0`$, so that the bifurcation goes from stable for small $`d`$ to unstable for larger $`d`$. Another way to study this curve is to make the bifurcation analysis near the normal state, described in the analysis of $`\kappa _2(d)`$, which yields to (4). Then one would need to take this development to higher order to get the sign of energy of the bifurcated branch. This sign changes on $`\kappa _1(d)`$. The fact that the bifurcation from the normal state is unstable for large $`d`$ has not been studied.

Curve $`\kappa _4(d)`$:

Another open problem is to determine the behavior of $`\kappa _4(d)`$ as $`d`$ tends to infinity. We expect it to be of the order $`C/d`$ for some constant $`C`$. We believe that the analysis that we have explained for the curve $`\kappa _2(d)`$ is the same here, that is on $`\kappa _4(d)`$ as well the eigenvalue has two eigenfunctions. The same analysis of the linearized problem needs to be performed. The difference with $`\kappa _2(d)`$ is that the eigenfunction with no vortex is stable on $`\kappa _2(d)`$ and unstable on $`\kappa _4(d)`$.

The point of intersection

Note that all curves $`\kappa _i(d)`$ intersect at the same point. Indeed the point of intersection of $`\kappa _2(d)`$ and $`\kappa _4(d)`$ has an eigenvalue with two eigenfunctions, one of which (the one without vortices) changes stability. Hence this point also belongs to $`\kappa _1(d)`$ since on $`\kappa _1(d)`$ the stability of the solution without vortex changes. Finally this point belongs to $`\kappa _3(d)`$ since the energy of the bifurcated solution is zero for both eigenfunctions, in particular for the vortex solution. We want to point out that a similar analysis for the intersection of these curves has been performed in  in the one dimensional setting and it yields to a solvability condition of fourth order at the point of intersection.

## 5 Conclusion

We have obtained very detailed bifurcation diagrams for the Ginzburg-Landau model of a two dimensional cross section of a three dimensional superconducting cylinder when the applied field is along the direction of the axis. Detailed analysis are provided for the solution branches in the different regions of the $`(\kappa ,d)`$ plane. What distinguishes the different regions are features like the existence (or the lack of existence) of vortex solutions, the global and local stability of solutions, and the hysteresis phenomena. We note that the analysis includes regions with small and intermediate values of $`\kappa `$ which is not often featured in the existing studies, as most of the works in the literature focus on the vortex state which appears for larger values of $`\kappa `$ or for thin films. Though we have computed the bifurcation diagrams with rectangular cross sections, they are very representative of the general cases. Naturally, for small samples, the geometric conditions may have a stronger effect on leading to detailed alternations to the bifurcation curves. A change of topology, such as rings or shells, may present other complications  but we expect similar bifurcation diagrams remain valid. In addition, comparisons with existing theories have been made in the paper. Some remaining questions have also been raised.

Acknowledgments: This work was made possible thanks to the joint research project France HongKong (Procore) grant.


format: text

text:
# 1 Introduction

## 1 Introduction

Motivated by recent developments in string theory, the idea of confining matter in a three-brane embedded in a higher dimensional spacetime has stimulated lately a lot of activity in several fields of research. In cosmology, it has been shown that, quite generically, the matter confinement in a brane entails a non-standard evolution of the brane geometry . In other words, the Friedmann equations governing the expansion of our ‘apparent’ universe are modified. This could lead to trouble because nucleosynthesis would give very different light element abundances with a modified scale factor evolution. Therefore, a strong constraint on brane models is the recovery of standard cosmological evolution at least since nucleosynthesis. A nice solution to this problem has been provided by applying the idea of the Randall-Sundrum model  that a tension in the brane can be compensated by a negative cosmological constant in the bulk, as far as the internal brane expansion is concerned. Applied to cosmology, this idea leads to the recovery of the usual Friedmann law when the ordinary cosmological energy density is much less than the brane tension . By contrast, at earlier times, when the energy density is higher than the tension one finds a non-standard evolution in the brane .

Another solution that has been suggested in the case of finite size extra-dimension (see e.g.  for the case of infinite extra-dimension) is to consider a stabilizing potential for the so-called radion . The precise definition of the radion is the source of some confusion in the literature. It is usually described as the metric component along the extra-dimension(s), but as emphasized in  it can be reinterpreted as the relative distance between the two branes. Stabilization of the radion can be explicitly realized  by adding some matter in the bulk, a scalar field for simplicity, so that the relative distance between the two branes will be stabilized. This mechanism has been applied  to the two-brane model of Randall and Sundrum , in order to recover gravity consistent with the observations on the negative tension brane (without stabilization, one gets a Brans-Dicke gravity with a negative parameter ). In cosmology, this mechanism has been discussed in  and  but the claim that this would lead to the usual Friedmann equations has been challenged in . Let us mention that the radion has also been discussed in the context of holography .

The purpose of the present work will not be to study the cosmological consequences of a stabilizing mechanism for the radion, task which will be left for future work, but more modestly to characterize the radion in cosmology (see also a brief discussion in ), to write down the equations of motion that govern its dynamics, which will be analysed in some simple examples.

Previous treatments of the radion in brane systems  start from the beginning with a perturbative approach: the radion is thus hidden among the metric fluctuations, which have been the subject of active research lately in cosmology (see  and references therein). Our approach is quite different in the sense that the radion dynamics appears directly from the junction conditions. It also gives us the full non-linear dynamics of the homogeneous radion (i.e. ignoring the ordinary spatial dependence of the radion). We are thus able to show that the dynamics of the radion is more subtle than that of a scalar field coupled to the trace of the matter energy-momentum tensor.

Our plan is the following. We begin in section 2 by defining the coordinate system in which we describe the two-brane configuration. We then focus on the internal geometries inside the two branes in section 3. Section 4 is the central part of this paper, where the equations of motion of the radion are derived. Starting in section 5, we apply our general results to the case of an empty anti-de Sitter ($`AdS_5`$) bulk spacetime. The last two sections are devoted to the analysis of the radion evolution in this geometry, first using a perturbative approach (section 6), then in the full non-linear case (section 7).

## 2 The two-brane cosmological system

The purpose of this section will be to present in a geometrical way the space-time configuration we wish to study. The spacetime will be supposed to be five-dimensional and to satisfy three-dimensional spatial homogeneity and isotropy in order to obtain (unperturbed) cosmological configurations. Two spatially homogeneous and isotropic three-branes live in this spacetime. The fifth dimension is periodic and we assume a mirror (orbifold) symmetry across each of the branes. The space between the two branes represents half of the periodic space along the fifth dimension, the other half being the mirror image with respect to any of the two branes.

In any coordinate system, even adapted to the three-dimensional homogeneity and isotropy, the two branes would be in general in motion along the fifth dimension. It turns out that it is possible, and convenient, to define the fifth dimensional coordinate $`y`$ as part of a Gaussian normal coordinate system with respect to one of the branes, chosen arbitrarily and called the reference brane. The reference brane is at rest in the fifth dimension, in $`y=0`$ say, and the spacetime metric in this coordinate system (which is not unique) is of the form

$$ds^2=g_{AB}dx^Adx^B=n(t,y)^2dt^2+a(t,y)^2\gamma _{ij}dx^idx^j+dy^2,$$
(1)
where the $`x^i`$ are coordinates for the maximally symmetric three-dimensional submanifolds, $`t`$ is a time coordinate. The metric components are necessarily of this form because of the symmetries, but their explicit dependence of $`t`$ and $`y`$ can be obtained only by solving the five-dimensional Einstein’s equations (with the appropriate boundary conditions),

$$G_{AB}R_{AB}\frac{1}{2}Rg_{AB}=\kappa ^2T_{AB},$$
(2)
where $`R_{AB}`$ is the Ricci tensor, $`R`$ its trace, and $`T_{AB}`$ the energy-momentum tensor. The coupling of the matter to the geometry, $`\kappa ^2`$, can be written as the inverse of the third power of the five dimensional reduced Planck mass.

If the reference brane is now at rest, the second brane will be in general moving in the same coordinate system. Its position at any time $`t`$ will be denoted

$$y=(t),$$
(3)
and we call this function $`(t)`$, which represents the relative distance between the two branes, the cosmological radion. Note that the rôle of the reference brane and of the second brane can of course be exchanged. Finally we wish to stress that here, unlike most of the literature, our radion is defined in a non-perturbative way. We will show later that our definition coincides, in the linearized approach, with the traditional definition of the radion.

## 3 Induced geometries in the two branes

The induced metric in the reference brane immediately follows from the spacetime metric (1) and can be written in the usual Robertson-Walker form,

$$ds_{induced}^2=dt^2+a_0^2d𝐱^2,$$
(4)
where the time coordinate has been redefined such that $`n(t,0)=1`$, and $`a_0(t)a(t,y=0)`$. $`t`$ is thus the proper time associated with the reference brane.

To get the induced geometry for the second brane, one must take into account its displacement in the fifth dimension, and one finds from (1) the induced line element

$$ds_{induced}^2=\left[n^2(t,(t))\dot{}^2\right]dt^2+a^2(t,(t))d𝐱^2,$$
(5)
where a dot will always stand for a (possibly partial) derivative with respect to $`t`$. This can be rewritten in the form

$$ds_{induced}^2=d\tau ^2+a_2^2(\tau )d𝐱^2,$$
(6)
where $`\tau `$ is the second brane cosmic time (i.e. the proper time for a comoving observer) which is related to the reference brane cosmic time $`t`$ by

$$d\tau =n(t,(t))\sqrt{1\frac{\dot{}^2}{n^2(t,(t))}}dtn_2\gamma ^1dt.$$
(7)
And the scale factor in the second brane is simply

$$a_2(t)=a(t,(t)),$$
(8)
or $`a_2(\tau )`$ if one expresses $`t`$ in terms of $`\tau `$.

It is also useful to introduce the expansion rate for the second brane, which can be defined as

$$H_2(t)\frac{da_2/dt}{a_2}=\left(\frac{\dot{a}}{a}+\frac{a^{}}{a}\dot{}\right)_2.$$
(9)
However, $`H_2`$ does not coincide with the standard definition of the Hubble parameter for a cosmological observer in the second brane because it is defined with respect to the time $`t`$, which is the cosmic time in the reference brane but not in the second brane. The traditional Hubble parameter for the second brane corresponds to the definition

$$_2(\tau )\frac{da_2/d\tau }{a_2}=\frac{\gamma }{n_2}H_2$$
(10)
Finally, one can introduce the unit velocity vector corresponding to the second brane, whose components read

$$u^A=\{\frac{dt}{d\tau },0,0,0,\frac{dy}{d\tau }\}=n_2^1\gamma \{1,0,0,0,\frac{d}{dt}\}.$$
(11)
At this point, we can easily infer the cosmological evolution in the second brane from the knowledge of the cosmological evolution in the first brane and the position of the second brane with respect to the first. What will be studied next is the connection with the matter content of the second brane, which has been ignored up to this point.

## 4 Junction conditions

The matter content of the branes is directly related to the jump of the extrinsic curvature tensor across the brane. This relation has been established explicitly in the case of a brane at rest with respect to the coordinate system (1) in our previous works . We generalize our result to include the case of a brane moving with respect to the coordinate system.

We will now consider explicitly the second brane, but all the following expressions apply as well to the reference brane by taking $`=0`$. The extrinsic curvature tensor on the brane is defined by the expression

$$K_{AB}=h_A^C_Cn_B,$$
(12)
where $`n^A`$ is a unit vector field normal to the brane worlsheet and

$$h_{AB}=g_{AB}n_An_B$$
(13)
is the induced metric on the brane.

It is easy to compute the components of the unit normal vector from the components of the brane velocity (11), by using the fact that $`n^A`$ must satisfy the two following conditions:

$$g_{AB}n^An^B=1,g_{AB}n^Au^B=0,$$
(14)
and we get

$$n^A=\{\frac{\dot{}}{n^2\sqrt{1\frac{\dot{}^2}{n^2}}},0,0,0,\frac{1}{\sqrt{1\frac{\dot{}^2}{n^2}}}\}.$$
(15)
Then, a straightforward calculation consisting in the substitution of the components (15) into the definition (12) yields the following (non zero) components for the extrinsic curvature tensor:

$`K_0^0`$ $`=`$ $`{\displaystyle \frac{\ddot{}+nn^{}2\frac{n^{}}{n}\dot{}^2\frac{\dot{n}}{n}\dot{}}{n^2\left(1\frac{\dot{}^2}{n^2}\right)^{5/2}}},`$ (16)
$`K_0^5`$ $`=`$ $`\dot{}K_0^0,K_5^0=K_0^5/n^2,`$ (17)
$`K_j^i`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{1\frac{\dot{}^2}{n^2}}}}\left({\displaystyle \frac{a^{}}{a}}+{\displaystyle \frac{\dot{a}\dot{}}{an^2}}\right)\delta _j^i,`$ (18)
$`K_5^5`$ $`=`$ $`{\displaystyle \frac{\dot{}^2\left(\ddot{}+nn^{}2\frac{n^{}}{n}\dot{}^2\frac{\dot{n}}{n}\dot{}\right)}{n^4\left(1\frac{\dot{}^2}{n^2}\right)^{5/2}}},`$ (19)
where all coefficients take their value on the brane, i.e. at $`t`$ and $`y=(t)`$.

Let us now introduce the energy-momentum tensor on the second brane, which, because of the symmetries of the setup, is necessarily of the perfect fluid form

$$S_{AB}^{(2)}=(\rho _2+P_2)u_Au_B+P_2h_{AB}.$$
(20)
We also define $`\widehat{S}_{AB}`$ as

$$\widehat{S}_{AB}S_{AB}\frac{1}{3}Sh_{AB}.$$
(21)
The junction conditions  for a hypersurface in a five-dimensional world read

$$\left[K_{\mu \nu }\right]=\kappa ^2\widehat{S}_{\mu \nu },$$
(22)
where the brackets stand for the jump across the brane, and $`K_{\mu \nu }=e_\mu ^Ae_\nu ^BK_{AB}`$ (and a similar expression for $`\widehat{S}_{\mu \nu }`$) where the four vectors $`e_\mu ^A`$ ($`\mu =0,1,2,3`$) form a basis of the vector space tangent to the brane worldvolume. For the second brane, this gives explicitly

$$\left[K_{\mu \nu }\right]_2=K_{\mu \nu }(t,(t)^+)K_{\mu \nu }(t,(t)^{})=2K_{\mu \nu }(t,(t)),$$
(23)
where the last equality follows from the identification of the point $`y=(t)`$ with the point $`y=(t)`$ and from the $`Z_2`$ symmetry across the second brane. Note that, whereas all expressions (15-22) apply also to the reference brane (by taking $`=0`$), the last expression differs by a sign for the two branes: one would find $`\left[K_{\mu \nu }\right]_0=2K_{\mu \nu }(t,y=0)`$ for the reference brane.

It is easy to compute the components of $`\widehat{S}_A^B`$, using the velocity components (11) and one finds

$`\widehat{S}_0^0`$ $`=`$ $`{\displaystyle \frac{2\rho _2+3P_2}{3}}\gamma ^2`$ (24)
$`\widehat{S}_5^0`$ $`=`$ $`{\displaystyle \frac{2\rho _2+3P_2}{3}}\gamma ^2{\displaystyle \frac{\dot{}}{n^2}},\widehat{S}_0^5={\displaystyle \frac{2\rho _2+3P_2}{3}}\gamma ^2\dot{}`$ (25)
$`\widehat{S}_j^i`$ $`=`$ $`{\displaystyle \frac{1}{3}}\rho _2\delta _j^i`$ (26)
$`\widehat{S}_5^5`$ $`=`$ $`{\displaystyle \frac{2\rho _2+3P_2}{3}}\gamma ^2{\displaystyle \frac{\dot{}^2}{n^2}}.`$ (27)
Equating this with the extrinsic curvature tensor, according to (22), all equations reduce to only two equations, which read

$`{\displaystyle \frac{\ddot{}}{n^2}}+{\displaystyle \frac{n^{}}{n}}\left(12{\displaystyle \frac{\dot{}^2}{n^2}}\right){\displaystyle \frac{\dot{n}}{n}}{\displaystyle \frac{\dot{}}{n^2}}`$ $`=`$ $`{\displaystyle \frac{\kappa ^2}{6}}\left(2\rho _2+3P_2\right)\left(1{\displaystyle \frac{\dot{}^2}{n^2}}\right)^{3/2}`$ (28)
$`{\displaystyle \frac{a^{}}{a}}+{\displaystyle \frac{\dot{a}}{a}}{\displaystyle \frac{\dot{}}{n^2}}`$ $`=`$ $`{\displaystyle \frac{\kappa ^2}{6}}\rho _2\left(1{\displaystyle \frac{\dot{}^2}{n^2}}\right)^{1/2},`$ (29)
where, once more, all the metric coefficients take their value on the brane. The above system describes the full non-linear dynamics of the homogeneous cosmological radion. The non-linearity is manifest in the terms quadratic in $`\dot{}`$ but is also hidden in the fact that all the metric coefficients depend on $``$.

It is interesting to note that the left hand side of equation (28) is simply the (covariant) acceleration in the fifth direction of a particle with the velocity (11). This can be understood by noticing that $`u^AK_{AB}u^B=n_Bu^C_Cu^B`$, as a consequence of (12) and (14). $`u^AK_{AB}u^B`$, which is also $`K_0^0`$, therefore vanishes for a geodesic motion. The right hand side of equation (28) will thus induce a deviation from geodesic motion for a “comoving” observer on the second brane (at fixed coordinate $`x^i`$).

The second order equation (28) can in fact be derived from the first order equation (29). Indeed, if one takes the time derivative of (29) and adds it to $`H_2`$ times (29), one gets, after an overal multiplication by $`\gamma ^2`$, the expression

$`H_2\left[{\displaystyle \frac{\ddot{}}{n^2}}+{\displaystyle \frac{n^{}}{n}}\left(12{\displaystyle \frac{\dot{}^2}{n^2}}\right){\displaystyle \frac{\dot{n}}{n}}{\displaystyle \frac{\dot{}}{n^2}}\right]`$ (31)
$`{\displaystyle \frac{1}{3}}G_{05}\left(1{\displaystyle \frac{\dot{}^4}{n^4}}\right){\displaystyle \frac{1}{3}}\left(n^2G_{00}+G_{55}\right)\left(1{\displaystyle \frac{\dot{}^2}{n^2}}\right)\dot{}`$
$`=`$ $`{\displaystyle \frac{\kappa ^2}{6}}H_2\left(2\rho _2+3P_2\right)\left(1{\displaystyle \frac{\dot{}^2}{n^2}}\right)^{3/2},`$ (32)
where we have used

$$\dot{\rho _2}=3H_2\left(\rho _2+P_2\right),$$
(33)
which expresses the energy conservation on the second brane (assuming there is no inflow/outflow of energy from/to the bulk). A bulk energy-momentum tensor, satisfying the three-dimensional symmetries of homogeneity and isotropy, can be written in the form

$$T_{AB}=\rho _Bv_Av_B+P_T\left(h_{AB}^{(2)}+v_Av_B\right)+P_Bh_{AB}^{(3)},$$
(34)
where $`v^A`$ is a unit time-like vector (whose ordinary spatial components are zero) describing the fluid velocity, $`P_T`$ represents the transverse pressure (along the fifth dimension) and $`P_B`$ represents the pressure in the three ordinary dimensions; $`h_{AB}^{(2)}`$ is the projection tensor on the two-dimensional $`(t,y)`$ subspace and $`h_{AB}^{(3)}`$ is the projection tensor on the three-dimensional ordinary space. Because of the $`Z_2`$ symmetry, the bulk fluid cannot go “across” the branes, and we must thus impose the further condition on the bulk fluid that the vector field $`v^A`$ coincides, on the branes, with the brane velocity $`u^A`$ 333 Simultaneously $`h_{AB}^{(2)}+v_Av_B=n_An_B`$. . It is then not difficult to check explicitly that the second and third terms on the left hand side of (32), expressed in terms of $`T_{AB}`$ by use of Einstein’s equations (2), just cancel. Equation (28) is thus derived, provided $`H_2`$ is not zero.

Equation (29) can be solved algebraically to express $`\dot{}`$ in terms of the other quantities. One finds

$`\dot{}`$ $`=`$ $`n\left({\displaystyle \frac{a^{}\dot{a}}{a^2n}}\pm {\displaystyle \frac{\kappa ^2\rho _2}{6}}\sqrt{{\displaystyle \frac{\dot{a}^2}{a^2n^2}}{\displaystyle \frac{a^2}{a^2}}+{\displaystyle \frac{\kappa ^4\rho _2^2}{36}}}\right)\left({\displaystyle \frac{\dot{a}^2}{n^2a^2}}+{\displaystyle \frac{\kappa ^4\rho _2^2}{36}}\right)^1.`$ (35)
If $`a^{}/a`$ and $`\rho _2`$ are of the opposite sign, then (29) has one solution if $`|\dot{a}/(an)|>|a^{}/a|`$, no solution otherwise. If $`a^{}/a`$ and $`\rho _2`$ are of the same sign, then (29) has two solutions when $`|\dot{a}/(an)|<|a^{}/a|`$, and only one solution otherwise. In both cases, when there is a unique solution, it corresponds to the root in (35) with the sign being the same as that of $`\dot{a}/a`$.

To conclude this section, let us mention that the equations of motion for the radion can be rewritten in terms of the cosmic time on the second brane, $`\tau `$, instead of the cosmic time on the reference brane $`t`$ :

$`_{,\tau \tau }+{\displaystyle \frac{n^{}}{n}}\left(1+_{,\tau }^2\right)`$ $`=`$ $`{\displaystyle \frac{\kappa ^2}{6}}\left(2\rho _2+3P_2\right)\sqrt{1+_{,\tau }^2},`$ (36)
$`{\displaystyle \frac{a^{}}{a}}+{\displaystyle \frac{\dot{a}}{an}}{\displaystyle \frac{_{,\tau }}{\sqrt{1+_{,\tau }^2}}}`$ $`=`$ $`{\displaystyle \frac{\kappa ^2}{6}}\rho _2\left(1+_{,\tau }^2\right)^{1/2}.`$ (37)
The second equation has been obtained by  in the particular case where the reference brane is Minkowski (hence $`\dot{a}/a=0`$ and $`\rho _2`$ is a constant.

## 5 Radion dynamics in an $`AdS_5`$ bulk spacetime

In order to explore the dynamical evolution of the radion, we now need to be more specific and to resort to explicit expressions for the spacetime metric (1). We will thus choose to work with an $`AdS_5`$ bulk spacetime, for which it is possible to recover (approximately) , at late times, standard cosmology, and which corresponds to the cosmological, i.e. non static, generalization of the Randall-Sundrum model .

Keeping aside the two branes, we assume that the bulk spacetime is empty, except for a negative cosmological constant $`\mathrm{\Lambda }`$. It is useful to define a mass scale $`\mu `$ associated with $`\mathrm{\Lambda }`$, according to the expression

$$\mu =\sqrt{\mathrm{\Lambda }/6},$$
(38)
($`\mu ^1`$ is the $`AdS_5`$ curvature radius) and an energy density scale $`\sigma `$, defined by

$$\sigma =6\frac{\mu }{\kappa ^2}.$$
(39)
Note that, in the Randall-Sundrum model , the energy density in the brane is supposed to be exactly this $`\sigma `$. In the present work, we consider any value for the energy densities in the two branes.

We will also assume that the five-dimensional Weyl tensor is zero. Our global spacetime is therefore made of two pieces of $`AdS_5`$ separated by the two branes. Assuming the usual $`Z_2`$ symmetry with respect to each brane, the two $`AdS_5`$ pieces are then simply mirror symmetric.

In , we have solved explicitly the five-dimensional Einstein’s equations

$$G_{AB}\mathrm{\Lambda }g_{AB}=\kappa ^2T_{AB},$$
(40)
in the case where the only matter in spacetime is the one confined in a brane. We obtained the following form for the component of (1):

$`a(t,y)`$ $`=`$ $`a_0(t)\left(\mathrm{cosh}\mu y\eta _0\mathrm{sinh}\mu |y|\right)`$ (41)
$`n(t,y)`$ $`=`$ $`\mathrm{cosh}\mu y\stackrel{~}{\eta }_0\mathrm{sinh}\mu |y|,`$ (42)
where $`\eta _0`$ and $`\stackrel{~}{\eta }_0`$ are two functions of time, related to the energy density in the reference brane by

$$\eta _0(t)=\frac{\rho _0(t)}{\sigma },\stackrel{~}{\eta }_0(t)=\eta _0+\frac{\dot{\eta }_0}{H_0},$$
(43)
where $`H_0\dot{a}_0/a_0`$ is the Hubble parameter. For future use we also define $`\eta _2`$ as

$$\eta _2(t)=\frac{\rho _2(t)}{\sigma },$$
(44)
where $`\rho _2(t)`$ is the energy density of the second brane. It must be stressed that, whereas the functional form of the metric components (41,42) is obtained by solving the Einstein equations in the bulk, the relationship between the coefficients and the reference brane matter given in (43) comes from the junction conditions  across the reference brane.

Finally, the evolution of the scale factor $`a_0(t)`$ and of the energy density $`\rho _0(t)`$ is determined by solving the coupled system consisting of the generalized Friedmann equation

$$H_0^2=\frac{\mathrm{\Lambda }}{6}+\frac{\kappa ^4}{36}\rho _0^2\frac{k}{a_0^2}=\mu ^2\left(\eta _0^21\right)\frac{k}{a_0^2}$$
(45)
where $`k=1,0,1`$ according to the three-dimensional spatial curvature (the term $`𝒞/a_0^4`$ does not appear because we have assumed a vanishing Weyl tensor) and the usual energy conservation law

$$\dot{\rho }_0=3H_0\left(\rho _0+P_0\right).$$
(46)
A useful identity is an integral of Einstein’s equations

$$\left(\frac{\dot{a}}{na}\right)^2=\frac{1}{6}\kappa ^2\rho _B+\left(\frac{a^{}}{a}\right)^2\frac{k}{a^2}$$
(47)
which was obtained in . The generalized Friedmann equation (45) on the reference brane is simply the application of this identity (with the junction conditions) in $`y=0`$.

The fully non-linear evolution of the radion, coupled to the evolution of the matter on the two branes, is thus given by the first-order system of differential equations consisting of (45), (46), (33) and (35), which can be rewritten in the simpler form

$`\dot{}`$ $`=`$ $`n\left({\displaystyle \frac{a^{}\dot{a}}{a^2n}}\pm {\displaystyle \frac{\kappa ^2\rho _2}{6}}|_2|\right)\left({\displaystyle \frac{\dot{a}^2}{n^2a^2}}+{\displaystyle \frac{\kappa ^4\rho _2^2}{36}}\right)^1,`$ (48)
where $`_2`$ is the Hubble parameter of the second brane defined in (10).

One can also check, starting from the expression (9) in terms of the metric component (41) and using (29) and (45), all terms in $`\dot{}`$ cancel and one recovers the generalized Friedmann equation on the second brane, that is

$$_2^2=\frac{\mathrm{\Lambda }}{6}+\frac{\kappa ^4}{36}\rho _{2}^{}{}_{}{}^{2}\frac{k}{a_2^2}=\mu ^2\left(\eta _2^21\right)\frac{k}{a_2^2}.$$
(49)
It has to be so because the generalized Friedmann equations comes from a purely local analysis, and the rôle of the two branes can be exchanged.

## 6 Radion fluctuations

In this section, we study perturbatively the dynamics of the radion fluctuations. We first define the background solution as the case where the radion is frozen, i.e. does not depend on time. This means that the second brane is at rest with respect to the first brane, i.e.

$$y_2(t)=\overline{},$$
(50)
where $`\overline{}`$ is a constant. One can then find explicit solutions of the Einstein’s equations for the whole two-brane system. However, to ensure the continuity of the metric along the compact fifth dimension, one finds that the matter content in the second brane is necessarily related to the matter content of the first brane. This result was already noticed for a vanishing bulk cosmological constant . We generalize it here to the case of a non-zero bulk cosmological constant.

If $``$ is constant, then the system (28-29) simplifies considerably to yield

$$\left(\frac{n^{}}{n}\right)_2=\frac{\kappa ^2}{6}\left(2\overline{\rho _2}+3\overline{P_2}\right)$$
(51)
and

$$\left(\frac{a^{}}{a}\right)_2=\frac{\kappa ^2}{6}\overline{\rho _2}.$$
(52)
Substituting the explicit expressions (41) and (42) for $`a`$ and $`n`$, which depend on the energy content on the first brane, one can express the energy density and pressure of the second brane as functions of the energy density and pressure on the first brane and of the position $``$. One finds

$$\kappa ^2\overline{\rho _2}(t)=6\mu \frac{\mathrm{sinh}\mu \overline{}\eta _0\mathrm{cosh}\mu \overline{}}{\mathrm{cosh}\mu \overline{}\eta _0\mathrm{sinh}\mu \overline{}},$$
(53)
and

$$\kappa ^2\overline{P_2}(t)=2\mu \frac{\mathrm{sinh}\mu \overline{}\stackrel{~}{\eta }_0\mathrm{cosh}\mu \overline{}}{\mathrm{cosh}\mu \overline{}\stackrel{~}{\eta }_0\mathrm{sinh}\mu \overline{}}\frac{2}{3}\kappa ^2\overline{\rho _2}.$$
(54)
Note that if one considers an equation of state of the form $`P_0=w_0\rho _0`$ with $`w_0`$ constant for the matter on the reference brane (which includes the standard cases of radiation, non-relativistic matter and a cosmological constant), the equation of state of the matter on the second brane, which is completely determined here as a function of the relative distance between the two branes, will correspond to an equation of state $`\overline{P}_2=w_2\overline{\rho }_2`$ where, in general, $`w_2`$ is time dependent. Indeed $`\eta _0`$ and

$$\stackrel{~}{\eta }_0=(2+3w_0)\eta _0$$
(55)
are time dependent in general.

One exception is when the reference brane contains only a cosmological constant, i.e. $`w_0=1`$, then one also gets $`w_2=1`$ on the second brane. There is in this case an interesting correspondence between the cosmological constants on the two branes. Introducing the effective cosmological constant from the brane point of view as

$$\lambda =\frac{\mathrm{\Lambda }}{6}+\frac{\kappa ^4}{36}\rho =\mu ^2(\eta ^21),$$
(56)
we have

$$\lambda _2=\frac{\lambda _0}{\left(\mathrm{cosh}\mu \overline{}\eta _0\mathrm{sinh}\mu \overline{}\right)^2}.$$
(57)
Thus, the two branes are simultaneously $`dS_4`$ ($`\lambda >0`$, $`|\eta |>1`$), $`M_4`$ ($`\lambda =0`$, $`|\eta |=1`$) or $`AdS_4`$ ($`\lambda <0`$, $`|\eta |<1`$) .

Another interesting case is when the reference brane undergoes ordinary cosmology, which occurs when $`\eta _0`$ is very small. One can then neglect $`\eta _0`$ as well as $`\stackrel{~}{\eta }_0`$ in the expressions (53-54) and the matter on the second brane reduces effectively to a cosmological constant.

In order to obtain the equation of motion for the radion fluctuations, we linearize the exact equations of motion (28-29) about the equilibrium configuration we have just defined. The linearized system reads

$`{\displaystyle \frac{\ddot{\delta }}{n^2}}+\delta \left({\displaystyle \frac{n^{}}{n}}\right){\displaystyle \frac{\dot{n}}{n}}{\displaystyle \frac{\dot{\delta }}{n^2}}`$ $`=`$ $`{\displaystyle \frac{\kappa ^2}{6}}\left(2\delta \rho _2+3\delta P_2\right)`$ (58)
$`\delta \left({\displaystyle \frac{a^{}}{a}}\right)+{\displaystyle \frac{\dot{a}}{a}}{\displaystyle \frac{\dot{\delta }}{n^2}}`$ $`=`$ $`{\displaystyle \frac{\kappa ^2}{6}}\delta \rho _2,`$ (59)
where we recall that all metric coefficients take their value on the second brane. Furthermore, one can use the explicit expressions (41-42) to find

$`\delta \left({\displaystyle \frac{n^{}}{n}}\right)`$ $`=`$ $`m_n^2\delta \left[\mu ^2{\displaystyle \frac{\kappa ^4}{36}}\left(2\overline{\rho _2}+3\overline{P_2}\right)^2\right]\delta `$ (60)
$`\delta \left({\displaystyle \frac{a^{}}{a}}\right)`$ $`=`$ $`m_a^2\delta \left[\mu ^2{\displaystyle \frac{\kappa ^4}{36}}\overline{\rho _2}^2\right]\delta ,`$ (61)
where we have introduced two effective squared masses $`m_n^2`$ and $`m_a^2`$. Note that both $`m_n`$ and $`m_a`$ vanish in the particular case where the Randall-Sundrum condition between the tension and the bulk cosmological constant is satisfied on the second brane.

If one adds (58) with three times (59), one obtains the equation of motion

$$\frac{\ddot{\delta }}{n^2}+\left(3\frac{\dot{a}}{a}\frac{\dot{n}}{n}\right)\frac{\dot{\delta }}{n^2}+\left(m_n^2+3m_a^2\right)\delta =\frac{\kappa ^2}{6}\delta T_2,$$
(62)
where $`\delta T_2\delta \rho _2+3\delta P_2`$ is the linear perturbation of the trace of the matter energy-momentum tensor. One thus recognizes the equation of motion for a scalar field coupled to the (perturbed) trace of the matter energy-momentum tensor. The friction coefficient is the usual one (rewriting in terms of $`\tau `$ would give the natural $`3_2`$ friction coefficient, as can be seen when extracting this equation directly from (36), (37)), but the mass term is an effective mass that depends both on the bulk cosmological constant and on the background matter content of the second brane.

This effective mass is time-dependent except in the case where the second brane matter has for equation of state $`P_2=\rho _2`$. The two squared masses are then identical and equal to

$$m_a^2=m_n^2=\lambda _2,$$
(63)
where $`\lambda _2`$ is the effective cosmological constant from the brane point of view defined in (56). One thus recovers exactly previous results . An immediate consequence of this is that when the first brane is $`dS_4`$, the small perturbations around the equilibrium position of the second brane (then also $`dS_4`$) are unstable, whereas when the first brane is $`AdS_4`$ (the second brane is then also $`AdS_4`$) they are stable. We will see in section 7 that this perturbative analysis matches with the exact non linear evolution.

The evolution of the radion fluctuations is entangled with the matter fluctuations. It is however possible when the matter fluctuations are of the adiabatic type, i.e.

$$\delta P_2=c_{s2}^2\delta \rho _2,$$
(64)
to derive an equation for only the radion fluctuations by considering the appropriate linear combination of (58) and (59),

$$\frac{\ddot{\delta }}{n^2}+\left[\left(2+3c_{s2}^2\right)\frac{\dot{a}}{an^2}\frac{\dot{n}}{n^3}\right]\dot{\delta }+\left[m_n^2+\left(2+3c_{s2}^2\right)m_a^2\right]\delta =0.$$
(65)
The evolution of the matter fluctuations is then related to the radion fluctuation by (59), i.e.

$$\frac{\kappa ^2}{6}\delta \rho _2=\frac{\dot{a}}{a}\frac{\dot{\delta }}{n^2}+m_a^2\delta .$$
(66)
Assuming an equation of state $`\overline{P_2}=w_2\overline{\rho _2}`$ (and adiabatic fluctuations, i.e. $`c_{s2}^2=w_2`$), the equation of motion (65) takes the form, after substitution of $`m_a^2`$ and $`m_n^2`$ by their explicit expressions,

$$\delta _{,\tau \tau }+(2+3w_2)_2\delta _{,\tau }+3(1+w_2)\left(\mu ^2\frac{\kappa ^4}{36}\overline{\rho _2}^2(2+3w_2)\right)\delta =0,$$
(67)
where we have used here the second brane cosmic time $`\tau `$. One can solve this equation for de Sitter ($`w_2=1`$) on the two branes and one gets

$$\delta =\delta _ie^{_2(\tau \tau _i)}\frac{\kappa ^2}{6_2^2}\delta \rho _{2i}\left(1e^{_2(\tau \tau _i)}\right),$$
(68)
where $`\delta _i`$ is the initial radion fluctuation and $`\delta \rho _{2i}`$ is the initial density fluctuation (at time $`\tau _i`$). For de Sitter brane cosmology, we recover the instability mentioned above.

Note that in the case of the radion fluctuations in the strict Randall-Sundrum model, studied in , the friction coefficient in the above equation vanishes because the background is static, and both the effective squared masses $`m_a^2`$ and $`m_n^2`$ also vanish because of the RS condition relating the bulk cosmological constant and the brane tensions. One thus ends up with the equation of motion for a free scalar field.

## 7 Non linear evolution of the radion

In this section, we explore the non-perturbative aspects of the radion evolution in some simple cases, leaving a more systematic analysis for future work. We are going to consider the cases where the matter on the reference brane behaves like a cosmological constant. In these simple cases, the reference brane matter does not evolve, i.e. $`\dot{\eta }_0=0`$, and the metric components are separable and can be written in the form

$$n(t,y)=n(y),a(t,y)=a_0(t)n(y).$$
(69)
The explicit form of $`n(y)`$ depends on the sign of $`\eta _01`$. If $`\eta _0>1`$, the brane is $`dS_4`$; if $`\eta _0<1`$, the brane is $`AdS_4`$; and the case $`\eta _0=1`$ corresponds to the Randall-Sundrum model . We will then further specialize to the cases where the brane matter of the second also behaves as a cosmological constant, implying $`\dot{\eta }_2=0`$.

### 7.1 $`dS_4`$ brane

For a $`dS_4`$ reference brane, the dependence of the metric components on the fifth coordinate is simply (see e.g. )

$$n(y)=\sqrt{\eta _0^21}\mathrm{sinh}\mu (y_\mathrm{h}|y|),$$
(70)
where

$$y_\mathrm{h}=\frac{1}{\mu }\mathrm{arg}\mathrm{coth}\eta _0$$
(71)
Here we suppose that the second brane is always located at $`<y_\mathrm{h}`$. Otherwise, there would be a horizon at $`=y_\mathrm{h}`$.

Equation (35) can then be rewritten in terms of dimensionless quantities in the form

$$x_{,s}=𝒦_0^2\mathrm{sinh}x\frac{h_0\mathrm{cosh}x\pm \eta _2\mathrm{sinh}x\sqrt{\eta _2^2𝒦_0^2\mathrm{sinh}^2x𝒦_0^2\mathrm{cosh}^2x+h_0^2}}{h_0^2+\eta _2^2𝒦_0^2\mathrm{sinh}^2x}$$
(72)
where we have introduced the dimensionless distance of the second brane to the horizon $`x`$ and the dimensionless time $`s`$

$$x\mu (y_\mathrm{h}),s=\mu t.$$
(73)
The dimensionless Hubble parameter $`h_0`$ is defined by

$$h_0\frac{1}{a_0}\frac{da_0}{ds}=\frac{H_0}{\mu },$$
(74)
and $`𝒦_0`$ is defined by

$$𝒦_0=\sqrt{|\eta _0^21|}$$
(75)
In general equation (72) is non separable since both $`\eta _2`$ and $`h_0`$ are function of time. It gets however simplified when one assumes that the energy momentum tensor of the second brane is that of a cosmological constant (implying $`\dot{\eta }_2=0`$), and if one assumes moreover that $`k=0`$ (in order to account for a slicing of $`dS_4`$ by flat 3-spaces) so that $`h_0=𝒦_0^2`$ (as can be seen from (45)). In this case equation (72) reads

$$x_{,s}=\sqrt{\eta _0^21}\mathrm{sinh}x\frac{\mathrm{cosh}x\pm \eta _2\sqrt{\eta _2^21}\mathrm{sinh}^2x}{1+\eta _2^2\mathrm{sinh}^2x},$$
(76)
it becomes separable and can be solved analytically. Since the matter energy densities are frozen on both branes, the evolution of the radion is fully embodied in this unique first order equation. The evolution of the radion is thus completely determined by the knowledge of the initial position of the second brane, which we call $`_i`$. If the initial position $`_i`$ is located at the equilibrium position

$$\overline{x}=\mu (y_\mathrm{h}\overline{})=\mathrm{arg}\mathrm{coth}\eta _2,$$
(77)
then the brane does not move. However, if the brane is initially away from the equilibrium position, it moves. According to the perturbative analysis, a brane slightly off the equilibrium position moves further away since we found that the radion is unstable.

Here, we can determine the non-linear evolution. We have solved numerically the evolution of the radion in the cases where the brane is initially on the left ($`_i<\overline{}`$) and on the right ($`_i>\overline{}`$) of the equilibrium position. We observe the expected unstable behaviour. In the first case, illustrated in Fig. 1, the brane moves towards the reference brane until it collides with it (numerically it goes across it). In the second case, illustrated in Fig. 2, the brane moves further away but its motion freezes asymptotically at the horizon, as seen in the reference brane proper time. However, as can be seen from (7) and (76), this takes only a finite amount of the second brane proper time, as is typical for the crossing of a (Rindler or black hole) horizon.

One may obtain a simple analytical expression which illustrates this behaviour in the case of a empty second brane : $`\eta _2=0`$. The solution is

$$x=\mathrm{arg}\mathrm{tanh}\left[\mathrm{tanh}x_i\mathrm{exp}\left(\sqrt{\eta _0^21}(ss_i)\right)\right],$$
(78)
where $`x_i`$ is the initial position at time $`s_i`$. This is also the solution of (76) for non zero $`\eta _2`$ when the second brane approaches the horizon (i.e. for $`x`$ small).

### 7.2 $`AdS_4`$ brane

For an $`AdS_4`$ reference brane one finds for the metric component

$$n(y)=\sqrt{1\eta _0^2}\mathrm{cosh}\mu (|y|y_c),$$
(79)
where

$$y_c=\frac{1}{\mu }\mathrm{arg}\mathrm{tanh}\eta _0.$$
(80)
In this case the warp factor $`n`$ has a minimum in $`y=y_c`$ enabling to have a configuration with two static branes of positive tensions  (see also ). In any case, equation (35) reads

$$x_{,s}=𝒦_0^2\mathrm{cosh}x\frac{h_0\mathrm{sinh}x\pm \eta _2\mathrm{cosh}x\sqrt{\eta _2^2𝒦_0^2\mathrm{cosh}^2x𝒦_0^2\mathrm{sinh}^2x+h_0^2}}{h_0^2+\eta _2^2𝒦_0^2\mathrm{cosh}^2x},$$
(81)
where $`x`$ is defined by

$$x\mu (y_c),$$
(82)
and $`s`$ and $`𝒦_0`$ are defined respectively as in (73) and (75). The scale factor is given straighforwardly by solving (45); one finds $`a_0=\mathrm{sin}(𝒦_0s)/\mu 𝒦_0`$, so that in this FLRW parametrization of $`AdS_4`$ there are coordinate singularities in $`𝒦_0s=0`$ and $`𝒦_0s=\pi `$. The dimensionless Hubble parameter is then $`h_0=𝒦_0\text{cotan}(𝒦_0s)`$. The equation of motion (81) remains non separable even in the simple case where $`\dot{\eta _2}=0`$, with $`\eta _20`$. We have solved numerically for the evolution of the radion in the cases where the brane (with $`\eta _2`$ being constant and positive) is initially (at time $`s_i`$ such that $`𝒦_0s_i<\pi /2`$) on the left ($`x_i<\overline{x}`$) and on the right ($`x_i>\overline{x}`$) of the equilibrium position $`\overline{x}`$ defined by

$$\overline{x}=\mu (\overline{}y_c)=\mathrm{arg}\mathrm{tanh}\eta _2.$$
(83)
For small enough $`s_i`$, the only acceptable solution is given by taking the $`+`$ sign in front of the square root in equation (81), as can be seen using the arguments given at the end of section 4.

In both cases ($`x_i<\overline{x}`$ and $`x_i>\overline{x}`$), illustrated in figures 3 and 4, the brane first moves back toward its equilibrium position, and then reaches a point (at time $`s_c`$ close to $`\pi /2𝒦_0`$) where the square root444which is also the determinant of the quadratic equation obtained from equation (29) in equation (81) vanishes. The solution obtained from (81) with the $`+`$ sign in front of the square root must then be matched (after time $`s_c`$) to a solution of (81) with a $``$ sign is taken in front of the square root (the two solutions have the same tangent in $`s_c`$). Figures 3 and 4 are obtained with this matching.

Equation (81) becomes separable, and can be easily integrated, in the limiting case where the second brane is empty i.e. $`\eta _2=0`$. In this case, one finds:

$$x=\mathrm{arg}\mathrm{tanh}\left[\mathrm{tanh}x_i\frac{\mathrm{cos}\left(𝒦_0s\right)}{\mathrm{cos}\left(𝒦_0s_i\right)}\right].$$
(84)
### 7.3 $`M_4`$ brane

When the reference brane is Minkowskian ($`\eta _0`$ being 1, and $`a=n=\mathrm{exp}(\mu |y|)`$) , equation (35) is given by

$$x_{,s}=\pm \frac{e^x}{\eta _2}\sqrt{\eta _2^21},$$
(85)
where $`x`$ is given here by

$$x\mu $$
(86)
and $`s`$ is defined as above. The second brane equilibrium condition is simply given by the Randall-Sundrum  matching condition expressed here as $`\eta _2=1`$. For this particular $`\eta _2`$, one sees easily from equation (85) that the second brane must be at rest with respect to the first one, but that the equilibrium position $`\overline{x}`$ can be chosen arbitrarily. However a small departure from the equilibrium condition $`\eta _2=1`$ results in either a colliding of the two branes or a runaway behaviour.

## Acknowledgments

We thank Arthur Lue, who helped us to understand some subtleties associated with the numerical integration of equation (81). The work of C.D. is sponsored in part by NSF Award PHY 9803174, and by David and Lucile Packard Foundation Fellowship 99-1462.