text:
# Comment on “Density perturbations in the ekpyrotic scenario”

## I Introduction

Recently, a model called the ekpyrotic scenario ekp  was proposed as a possible competitor to the inflationary paradigm inflation  on the basis that it could solve the problems of the hot big bang model and also produce an almost scale invariant spectrum of primordial scalar perturbations. Although the foundations of the model are still under debate pyro , we only shall comment on the controversial Lyth ; BF ; HT  claim, crucial for the ekpyrotic scenario, that the power spectrum is scale-invariant perturbekp .

In Ref. perturbekp , the complicated five-dimensional evolution of the background is modeled by an effective four-dimensional spacetime which experiences a singular bounce at $`\eta =0`$ ($`\eta `$ being the conformal time) as the effective scale factor vanishes. The calculation of the power spectrum is a controversial issue as different authors would use different matching conditions matching ; MS ; Ruth  at the singular bounce. There, Bardeen’s gravitational potential Bardeen  $`\mathrm{\Phi }`$ diverges. The authors of Ref. perturbekp  suggested to work with the comoving density contrast $`ϵ_\mathrm{m}`$ which is regular at $`\eta =0`$ and has two linearly independent modes $`D`$ and $`E`$, such that $`ϵ_\mathrm{m}=ϵ_0D(\eta )+ϵ_2E(\eta )`$. At the singular point the equality of the coefficients before and after the bounce is assumed perturbekp : $`ϵ_0(0^+)=ϵ_0(0^{})`$ and $`ϵ_2(0^+)=ϵ_2(0^{})`$. Although this rule is given without justification (recall that the standard junction conditions follow from the Einstein equations), we shall refer to such equalities as matching conditions hereafter.

In Ref. us  we criticized the matching conditions suggested in Ref. perturbekp . We showed that they lead to an ambiguity, and thus an arbitrary choice needs to be made that present-day physics cannot make. Ref. perturbekp  was subsequently modified, see Ref. perturbekp2 , to address some critical comments that have been made on the ekpyrotic model, one of them being the above mentioned ambiguity exhibited in Ref. us . In the new version perturbekp2 , more details about the matching conditions are given and one claim has been added (claim 1): “One situation of special interest is … where the potential is irrelevant at $`\varphi \mathrm{}`$ and there is no radiation in the incoming state. In this case, $`ϵ_2(0^{})=0`$. In this case, we would obtain the same final result from any matching rule which set $`ϵ_2(0^+)=Aϵ_2(0^{})`$, with any constant $`A`$.” It concerns a particular case that we had not considered and for which our result supposedly does not apply. It should be noted however that the rest of Ref. perturbekp2  inconsistently relies on the general case instead of this particular one, as it was the case in the preceding version perturbekp . Claim 1 clearly utilizes, in an essential manner, the ideas developed in Ref. us  and represents a tentative a priori response to the points Ref. us  raised.

References perturbekp ; perturbekp2  also state (claim 2): “The prescription is invariant under redefining the independent solutions, e.g. by adding an arbitrary amount of the solution $`E(\eta )`$ to $`D(\eta )`$. Matching any other non-singular perturbation variable, defined to be an arbitrary linear combination of $`ϵ_\mathrm{m}`$ and $`ϵ_\mathrm{m}^{}`$ with coefficients which are non-singular background variables (defined to possess power series expansions in $`\tau `$, as above) will, with the same prescription of matching the amplitudes of both linearly independent solutions, also yield precisely the same result”. This is an important point in the ekpyrotic scenario because, if true, it would support the use of the variable $`ϵ_\mathrm{m}`$ as a tool to define matching conditions.

Below we show that both quoted claims are incorrect.

## II Is the Ekpyrotic spectrum unique?

Let us start with claim 1, i.e. the question concerning the invariance of the spectrum under rescaling of the parameter $`ϵ_2`$ and let us first define our notation. The authors of Ref. perturbekp2  consider the equation of state close to the singular bounce to be $`\omega p/\rho =\omega _0+\omega _1\eta +\omega _2\eta ^2+\mathrm{}`$. Before the bounce, the kinetic energy of the scalar field dominates and the dynamics can be approximated by that of a free scalar field for which $`\omega _0=1`$ and $`\omega _1=\omega _2=0`$ (note that one can show that the values of $`\omega _i`$, $`i>2`$, are irrelevant for this discussion). This should apply for $`\eta <0`$ close to the singularity ($`|\eta |1`$) and corresponds precisely to the situation that is assumed in claim 1 and to the case where the equations of Ref. us  cannot be applied (as mentioned in that article), hence the argument of Ref. perturbekp2 . Below, we complete the arguments of Ref. us  to include this particular case and show that the conclusions of Ref. us  remain unchanged.

At leading order, the behavior of the scale factor is a simple power-law, $`a(\eta )=\mathrm{}_0(\eta )^{1/2}/(2\sqrt{2})`$, in agreement with the notation of Eq. (24) in Ref. us . The comoving Hubble rate and the sound speed are $`a^{}/a=1/(2\eta )`$ and $`c__\mathrm{S}^2=1`$. The solution for the Bardeen potential Bardeen ; perturb  in the long wavelength limit reads

$$\mathrm{\Phi }=\frac{3}{4}B_1(k)\frac{}{a^2}+\frac{3}{4}B_2(k)\frac{}{a^2}^\eta \frac{\mathrm{d}\tau }{\theta ^2},$$
(1)
with $`\theta 1/(a\sqrt{1+\omega })`$ and the integral sign stands for the primitive of the integration kernel. The coefficients $`B_1`$ and $`B_2`$ are functions of the comoving wavenumber $`k`$, to be determined by means of a matching with the initial vacuum condition perturbekp ; us . In the case at hand,

$$\mathrm{\Phi }=\frac{3}{\mathrm{}_0^2\eta ^2}B_1+\frac{3}{8}B_2.$$
(2)
Note that the equation of motion for the Bardeen potential for a free scalar field reads $`\mathrm{\Phi }^{\prime \prime }+(3/\eta )\mathrm{\Phi }^{}+k^2\mathrm{\Phi }=0`$, which can be solved exactly in terms of Bessel functions,

$$\mathrm{\Phi }=\frac{3}{(k\eta )}\left[\frac{1}{4}B_2J_1(k\eta )+\frac{\pi }{2}\frac{k^2}{\mathrm{}_0^2}B_1N_1(k\eta )\right],$$
(3)
whose expansion in powers of $`\eta `$ permits to recover the special solution (43-44) of Ref. perturbekp  as well as Eq. (2).

Before the bounce the density contrast on comoving hypersurfaces, $`ϵ_\mathrm{m}=2k^2\mathrm{\Phi }/(3^2)`$, is given by

$$ϵ_\mathrm{m}=\frac{8k^2}{\mathrm{}_0^2}B_1k^2\eta ^2B_2.$$
(4)
Comparison of Eq. (4) with Eqs. (43-44) of Ref. perturbekp  yields

$$ϵ_0^<(k)=\frac{8k^2}{\mathrm{}_0^2}B_1^<,ϵ_2^<(k)=k^2B_2^<,$$
(5)
showing that the two “modes” decouple. It is interesting to compare with the case $`\omega _1^<0`$, for which one has us  (here we assume $`\omega _2^>=0`$ for the sake of simplicity; this does not change in any way the conclusion)

$`ϵ_0^<(k)`$ $`=`$ $`{\displaystyle \frac{8k^2}{\mathrm{}_0^2}}B_1^<{\displaystyle \frac{128\mathrm{ln}2k^2}{9\omega _1^{<2}}}B_2^<,`$ (6)
$`ϵ_2^<(k)`$ $`=`$ $`{\displaystyle \frac{9k^2\omega _1^{<2}}{8\mathrm{}_0^2}}B_1^<\left(1+2\mathrm{ln}2\right)k^2B_2^<.`$ (7)
Let us note at this point that, with the choice we made of the normalization of the integral in Eq. (1), the limit $`\omega _10`$ is singular, as can be seen from Eqs. (6) – (7) and as already mentioned in Ref. us . Another choice is possible as one is free to add an arbitrary constant amount of $`B_2`$ into $`B_1`$, simply by making the integral a definite one. In particular, one can arrange that the constant factor $`(\mathrm{ln}2)`$ in the above equations be made to vanish from the outset, and it is easy to convince oneself that such a choice does not modify anything if the correct equations are used.

In the ekpyrotic scenario, and in the long wavelength limit, the coefficient $`B_1^<k^{3/2}`$ is dominant with respect to $`B_2^<k^{1/2}`$. This means that the relation $`ϵ_29\omega _1^2ϵ_0/64`$ holds true if $`\omega _10`$, whereas it does not if $`\omega _1=0`$. In the latter case, one has $`ϵ_2=\mathrm{}_0^2B_2ϵ_0/(8B_1)kϵ_00`$, see Eq. (5). The authors of Ref. perturbekp2  missed that the limit $`\omega _10`$ does in no way imply $`ϵ_2=0`$ exactly, as implied by claim 1, in which this value is said to be used to perform the matching. Moreover, as demonstrated below, the fact that $`ϵ_2`$ becomes subdominant does not remove the ambiguity noticed in Ref. us .

This is also related to the excessive claim in perturbekp2  that “there is no long wavelength contribution to $`\zeta `$ in the collapsing phase”. Inserting (3) into the gauge-invariant variable $`\zeta (2/3)(^1\mathrm{\Phi }^{}+\mathrm{\Phi })/(1+\omega )+\mathrm{\Phi }`$ (see e.g. MS ; Bardeen ; perturb ), gives

$$\zeta =\frac{1}{2}B_2J_0(k\eta )+\frac{\pi }{4}\frac{k^2}{\mathrm{}_0^2}B_1N_0(k\eta )$$
(8)
which, in the limit $`\eta 0`$, yields

$$\zeta \frac{1}{2}B_2\frac{1}{2}\frac{k^2}{\mathrm{}_0^2}B_1[\mathrm{ln}(k\eta )+\gamma __\mathrm{E}],$$
(9)
where $`\gamma __\mathrm{E}`$ is Euler constant coming from the expansion of the Bessel function. Since the second term is singular the first term is subdominant, but that does not imply that it vanishes altogether.

After the singular bounce ($`\eta >0`$), $`\omega _1^>0`$ and one recovers the relations us

$`ϵ_0^>(k)`$ $`=`$ $`{\displaystyle \frac{8k^2}{\mathrm{}_0^2}}B_1^>{\displaystyle \frac{128\mathrm{ln}2k^2}{9\omega _1^{>2}}}B_2^>,`$ (10)
$`ϵ_2^>(k)`$ $`=`$ $`{\displaystyle \frac{9k^2\omega _1^{>2}}{8\mathrm{}_0^2}}B_1^>\left(1+2\mathrm{ln}2\right)k^2B_2^>.`$ (11)
The authors of Ref. perturbekp2  propose to glue the epoch after the bounce to the epoch before the bounce by imposing that $`ϵ_0`$ and $`ϵ_2`$ are the same before and after the singularity. Contrary to standard junction conditions, and until some more microphysics is specified, this recipe does not rest on any physical principle and one may wonder what is the rationale behind it. The reason why one is pushed to adopt a new rule in this case is linked to the fact that one is dealing with a perturbative approach around a singular background whose meaning is questionable Lyth . Moreover, we argued in Ref. us  that this prescription is ambiguous since an arbitrary rescaling factor $`f(\omega ,k)`$ shows up in the final result. Therefore, the scenario contains an arbitrary function in a physically measurable quantity: until this function can be calculated unambiguously by first principles, the model cannot be falsified.

As in Ref. us , we rescale $`ϵ_2`$ by the completely arbitrary factor $`1/f`$. This gives the constant part of the gravitational potential

$$B_2^>=B_2^<\frac{f^>}{f^<}+\frac{9\omega _1^{>2}}{8\mathrm{}_0^2}B_1^<.$$
(12)
This is the main equation of Sec. II. It represents the equivalent of Eq. (54) of Ref. us  for the case $`\omega _1^<=0`$. On the other hand, the decaying mode amplitude reads

$$B_1^>=(1+2\mathrm{ln}2)B_1^<\frac{16\mathrm{ln}2\mathrm{}_0^2}{9\omega _1^{>2}}B_2^<\frac{f^>}{f^<},$$
(13)
given here for the sake of completeness.

Note that if the expansion of the equation of state parameter is done to higher order in $`\eta `$, all the relations of this section are only changed by inclusion of $`\omega _2`$ through the replacement $`(3/8)\omega _1^2w^{(2)}\omega _2+(3/8)\omega _1^2`$, in agreement with Ref. perturbekp ; perturbekp2 , while none of the $`\omega _i`$ for $`i>2`$ does contribute.

Eq. (12) shows that the spectrum effectively acquires a scale invariant piece $`B_1^<`$, together with an arbitrary piece $`B_2^<f^>/f^<`$. Indeed, as discussed in Ref. us , the functions $`f^<`$ and $`f^>`$ can depend on the background constants (e.g., $`\omega `$) and on the considered scale $`k`$ in a way which is not yet given from first principles. The choice $`f^>=f^<`$ fixes the $`k`$dependence of the spectrum, but, without physical justification, this remains an arbitrary choice, equivalent to assume a scale invariant spectrum from the outset.

## III An example

Let us now consider claim 2 according to which the choice of $`ϵ_\mathrm{m}`$ is essentially unique. The ambiguity inherent to the method proposed in Refs. perturbekp ; perturbekp2  can be most clearly exhibited by means of an example suggested at the Euro-conference in Annecy in December 2001 by G. Veneziano Ve1 . This section relies completely upon his idea.

In Ref. perturbekp2  it is argued that the density contrast $`ϵ_\mathrm{m}`$ is the quantity of interest because it is finite at $`\eta =0`$ and we have shown above how this quantity is used to propagate the spectrum through the singularity. However, $`ϵ_\mathrm{m}`$ is not the only finite, physically relevant, quantity. Based on the method developed in Ref. Ve , we could equally well consider and use the conjugate momentum $`\mathrm{\Pi }`$ to the conserved quantity $`\zeta `$ MS  which obeys the Hamilton like equations, valid for isentropic perturbations,

$$\mathrm{\Pi }=z^2\zeta ^{},\mathrm{\Pi }^{}=k^2c__\mathrm{S}^2z^2\zeta ,$$
(14)
where $`z^2a^2(1+\omega )/c__\mathrm{S}^2`$. It is worth pointing out that since the quantity $`\zeta `$ diverges at the singularity, it may be more useful to combine Eqs. (14) into the single one

$$\mathrm{\Pi }^{\prime \prime }+\left[3\left(c__\mathrm{S}^2\omega \right)2\right]\mathrm{\Pi }^{}+k^2c__\mathrm{S}^2\mathrm{\Pi }=0.$$
(15)
In the neighborhood of the singularity, $`c__\mathrm{S}^2=\omega =1`$, and the solutions of Eq. (15) are $`k\eta J_1(k\eta )`$ and $`k\eta N_1(k\eta )`$, where $`J_1`$ and $`N_1`$ are Bessel function of the first and second kind respectively. These solutions are completely regular at $`\eta =0`$. Note that Eqs. (14) and (15) apply precisely in the two cases of interest here, namely that of purely hydrodynamical perturbations and that of a free scalar field, by replacing $`c__\mathrm{S}^21`$ MS  in this last case. The quantity $`\mathrm{\Pi }`$ is related to $`ϵ_\mathrm{m}`$ as

$$\mathrm{\Pi }=a^2ϵ_\mathrm{m}.$$
(16)
Note also that the use of $`\mathrm{\Pi }`$ may be argued to be more appropriate in view of the fact that in the regular bounce case bounce , even though admittedly a different case as the ekpyrotic model, the variable $`ϵ_\mathrm{m}`$ ends up being an odd function of $`\eta `$, whereas $`\mathrm{\Pi }`$ is even. In the short-duration limit, that would imply that $`ϵ_\mathrm{m}`$ experiences a jump whereas $`\mathrm{\Pi }`$ is continuous.

Since $`|a^2|\mathrm{}_0^2/16`$ is finite at the singularity $`\eta =0`$, this relation clearly shows that the variable $`\mathrm{\Pi }`$, being a linear combination of $`ϵ_\mathrm{m}`$ and $`ϵ_\mathrm{m}^{}`$ (with vanishing coefficient for the latter), satisfies all the requirements of claim 2. Accordingly, it seems that there is no convincing reason to use $`ϵ_\mathrm{m}`$ rather than $`\mathrm{\Pi }`$. In the vicinity of the bounce, $`\mathrm{\Pi }`$ reads

$$\mathrm{\Pi }=\frac{\mathrm{}_0^2s}{16}[ϵ_0\frac{3}{4}\omega _1ϵ_0\eta +(ϵ_2\frac{3w^{(2)}}{8}ϵ_0)\eta ^2+\mathrm{}],$$
(17)
where $`s=1`$ before the bounce and $`s=+1`$ after the bounce. The proposal of Ref. perturbekp  then would consist in assuming that the coefficients in front of the constant term on the one hand and in front of the $`\eta ^2`$ term on the other hand are the same before and after the bounce. In the present context, this reduces to

$$ϵ_0^<=ϵ_0^>,ϵ_2^<=ϵ_2^>\frac{3w^{(2)>}}{8}ϵ_0^>,$$
(18)
since $`w^{(2)<}=0`$ by definition. From these relations, it is easy to establish that

$$ϵ_0^>=\frac{8k^2}{\mathrm{}_0^2}B_1^<,ϵ_2^>=k^2B_2^<\frac{3w^{(2)>}k^2}{\mathrm{}_0^2}B_1^<.$$
(19)
Finally, the computation of $`B_2^>`$ shows that the spectrum no longer contains a scale invariant piece,

$$B_2^>=B_2^<k^{1/2},$$
(20)
but gives a spectral index equal to $`3`$. It should be noticed that the same calculation goes through, with the same resulting spectral index if, in Eq. (16), one replaces the factor $`a^2`$ by its magnitude $`|a^2|`$ which is not only regular but also symmetric across the bounce. Therefore, the choice of the variable used to apply the rule proposed in Ref. perturbekp  plays a crucial role. In the absence of an underlying reason to choose a variable rather than another at the present stage, $`\mathrm{\Pi }`$ seems to be as relevant as $`ϵ_\mathrm{m}`$, and we are led to the conclusion that claim 2 is definitely erroneous. Moreover, $`\mathrm{\Pi }`$ is not the only well-defined variable that leads to such a conclusion: there is an infinite number of such possible combinations.

## IV Conclusions

To conclude, let us mention yet another possibility. In Ref. Ruth , it was suggested that the condition used in Ref. perturbekp  may be cast in the form of a matching of the energy perturbation as seen by an observer comoving with the fluid, and that such a matching is “at least as natural as matching the energy in the longitudinal gauge”. Firstly, one should notice that no rigorous proof exists that the conditions of Ref. Ruth  are equivalent to the prescription utilized in Ref. perturbekp . Secondly, we have shown in Ref. us  that, in the well known case of radiation to matter transition for which one can compare with the exact solution, this matching condition yields an incorrect result. A similar conclusion has been reached in Ref. BF  in the case of the reheating transition in which the equation of state also jumps. However, it can be argued that a bouncing situation is in no way comparable to a jump in the equation of state and that the usual matching conditions may then not apply. This is indeed what was found in Ref. bounce  for the case of a non singular bounce. The case for the ekpyrotic scenario is still more involved, as it was also argued that the bounce must be singular if the brane collision is to produce radiation Syksy . In any case, matching through an unspecified bounce is still an open question as there is no geometrical or physical argument imposing some particular choice, and no conclusion can be made about the power spectrum of perturbations in general, unless some non-singular bounce is specified, as it is the case of Ref. bounce .

Finally, we would like to insist on the fact that the overall gauge-invariant perturbation theory may be completely meaningless when a singularity is reached, since some gauge transformations that are admissible in any other situation could turn out to be singular. In such a situation, a “gauge-invariant” variable will, at the singular point, become “gauge-dependent”, and its use rather than the use of any other variable become an arbitrary, physically meaningless, choice. In fact, the transformation between the zero shear (conformal Newtonian) slicing and the comoving slicing is indeed singular at the bounce, which actually means that one or both slices are unphysical. Lyth Lyth  has shown that there is no slicing on which the density contrast and the intrinsic curvature perturbation are finite simultaneously. In such a situation, linear perturbation theory becomes meaningless.

Acknowledgments

We would like to thank Robert Brandenberger, Ruth Durrer, Justin Khoury, David Lyth and Raymond Schutz for many enlightening discussions. We are especially indebted with Gabriele Veneziano for proposing the variable $`\mathrm{\Pi }`$ of Sec. III and numerous other discussions. NPN would like to acknowledge IAP for warm hospitality during the time this work was being done, and CNRS and CNPq for financial support. DJS acknowledges financial support from the Austrian Academy of Sciences.


format: text

text:
# Reply to Comment on: ”Exact solutions of the Lawrence-Doniach model for layered superconductors”

## I Some background information on analytical properties of coupled static SG equations

The difference in mathematical description of the finite and the infinite LD models, emphasized in Ref.  , is not a ”speculation”: It results from translational invariance of the barrier potential and the absence of external boundaries in the layering direction in the latter case. The additional symmetry of the infinite model manifests itself in the appearance of an unphysical degree of freedom. The elimination of this unphysical degree of freedom requires minimization of the Gibbs free energy with respect to the phases $`\phi _n`$. In contrast, minimization with respect to $`\phi _n`$ is prohibited for the finite model by the physical requirement that the local field be equal to $`H`$ at $`x=0`$, $`x=\left(N1\right)p`$.

As an illustration, we compare coupled static SG equations for $`\varphi _n`$ in both the cases. \[For $`r(T)1`$ and $`HH_{c2}`$, these equations constitute solvability conditions for the Maxwell equations.\] For $`N=\mathrm{}`$, they read

$$\lambda _J^2\frac{d^2\varphi _n}{dy^2}=\frac{1}{ϵ^2}\left[\left(2+ϵ^2\right)\mathrm{sin}\varphi _n\mathrm{sin}\varphi _{n+1}\mathrm{sin}\varphi _{n1}\right],$$
(1)
$$n=0,\pm 1,\pm 2,\mathrm{}$$
(2)
Physical solutions to (1), (2) are subject to the requirement that the local field be equal to $`H`$ at $`y=\pm L`$, which yields

$$\frac{d\varphi _n}{dy}\left(\pm L\right)=2epH.$$
(3)
Conditions (3) allow only for solutions of the type

$$\varphi _n(y)=\varphi _n(y)+0\text{mod}\left(2\pi \right).$$
(4)
Note that due to the symmetry relations (4) conditions (3) alone do not specify any boundary value problem for (1), (2): The additive constants on the right-hand side of (4) remain undetermined. The determination of the constants $`0`$mod$`\left(2\pi \right)`$ that define a topological type of the solution requires imposition of boundary conditions on $`\varphi _n`$. A misunderstanding of this issue misled the author of the Comment into thinking that the solution to the boundary value problem for static SG equations could be ”nonunique”. However, the actual mathematical arbitrariness of (1), (2), related to the unphysical degree of freedom, lies in an infinite number of these equations: For an infinite number of ordinary differential equations one cannot formulate any existence and uniqueness theorems. The exact minimization procedure of Refs.  (or, alternatively, of Appendix A) leads to relations

$$\varphi _n(y)=\varphi _{n+1}(y)\varphi (y)$$
(5)
that reduce (1), (2) to the well-defined single static SG equation, with the single Josephson length $`\lambda _J`$. Unfortunately, in spite of ”extensive theoretical studies” of this equation, the fact that it admits an exact analytical solution, valid for any $`H`$ and $`L`$, has not been realized in any previous publications:

$$\varphi (y)=\pi \left(N_v1\right)+2\text{am}(\frac{y}{k\lambda _J}+K\left(k^2\right),k),$$
(6)
$$\text{dn}(\frac{L}{k\lambda _J},k)=\frac{\sqrt{1k^2}}{k}\frac{H_s}{H},N_v=2m,m=0,1,\mathrm{};$$
(7)
$$\varphi (y)=\pi N_v+2\text{am}(\frac{y}{k\lambda _J},k),$$
(8)
$$\text{dn}(\frac{L}{k\lambda _J},k)=k\frac{H}{H_s},N_v=2m+1,m=0,1,\mathrm{},$$
(9)
where am$`\left(u\right)`$ and dn$`\left(u\right)=\frac{d}{du}`$am$`\left(u\right)`$ are the Jakobi elliptic functions, and $`K\left(k^2\right)`$ is the elliptic integral of the first kind. The superheating (penetration) field $`H_sH_s\mathrm{}=\left(ep\lambda _J\right)^1`$ determines the upper bound of the existence of the Meissner state in the semiinfinite (along the layers) system $`0y<+\mathrm{}`$. The topological number $`N_v=0,1,\mathrm{}`$ specifies the number of vortex planes, with $`N_v=0`$ for the Meissner state. The range of the existence of the soliton solution with $`N_v1`$ is given by

$$\sqrt{H_{N_v1}^2H_s^2}HH_{N_v},$$
(10)
where $`H_{N_v}`$ is determined by the implicit equation

$$\frac{L}{\lambda _J}=\left(N_v+1\right)\frac{H_s}{H_{N_v}}K\left(\frac{H_s^2}{H_{N_v}^2}\right),$$
(11)
with $`H_0H_{sL}>H_s`$ being the superheating (penetration) field for $`L<\mathrm{}`$. The lower bound in (10) (the appearance of the soliton) is determined by applying (3) to the solution of the boundary value problem with topological conditions

$$\varphi (L)=0,\varphi (0)=\pi N_v,$$
(12)
which corresponds to the requirement that the density of the Josephson energy \[and, essentially, the free energy: see (A2)\] be a minimum at the boundaries $`y=\pm L`$. \[Owing to the symmetry relations (4), it is convenient to impose a topological condition at $`y=0`$ instead of $`y=+L`$. The role of topological boundary conditions is exhaustively discussed in the literature on soliton physics.\] The topological solution becomes unstable, when the density of the Josephson energy is a maximum at $`y=\pm L`$, i.e.,

$$\varphi (L)=\pi ,\varphi (0)=\pi N_v,$$
(13)
which determines, by (3), the upper bound in (10). \[For the Meissner solution, an analogous condition, determining $`H_{sL}`$, is obtained from (13) by setting $`N_v=0`$.\] In the intermediate field range, the appropriate boundary conditions follow from (3) and continuous dependence on $`H`$:

$$\frac{d\varphi }{dy}\left(L\right)=2epH,\varphi (0)=\pi N_v.$$
(14)
\[For the Meissner solution, we should set $`N_v=0`$ in (14).\]

Under corresponding redefinition of the parameters $`\lambda _J`$, $`H_s`$, the solution (6)-(11) describes the single junction ($`N=2`$) and the double-junction stack ($`N=3`$, see section II), which shows that vortex planes are a direct generalization of ordinary Josephson vortices. (In the single-junction case, $`N_v`$ should be interpreted as the number of ordinary vortices.) Equations (6)-(11) cover all well-known limiting cases.

An important feature of the exact solution is an overlap of the regions (10) with different topological numbers $`N_v`$. For the single junction, this overlap (not the alleged ”nonuniqueness of the solution”) was first established by a numerical analysis. It is also qualitatively discussed in Ref.  . Mathematically, the overlap is related to the fact that the solution with $`N_v`$ cannot be continuously transformed into the solution with $`N_v+1`$ by changing $`H`$. As shown in section II, the overlap is typical of vortex-plane solutions.

For $`N<\mathrm{}`$, equations (1) still hold, however, condition (2) should be replaced by the conditions

$$n=1,\mathrm{},N1,\varphi _0=\varphi _N0\left(2N<\mathrm{}\right).$$
(15)
Relations (3), (4) also hold, however, relations (5) should be replace by

$$\varphi _n=\varphi _{Nn},$$
(16)
obtained from the boundary conditions on the local field at $`x=0,`$ $`x=\left(N1\right)p`$. In contrast to the case $`N=\mathrm{}`$, we can now take advantage of powerful analytical methods of the theory of ordinary differential equations. Using Picard’s uniqueness and existence theorem, we establish the central analytical property of (1), (15):

Lemma. The initial value problem for (1), (15) with arbitrary initial conditions $`\varphi _n(y_0)=\alpha _n`$, $`\frac{d\varphi _n}{dy}(y_0)=\beta _n`$ has a unique solution in the whole interval $`\mathrm{}<y<+\mathrm{}`$. This solution has continuous derivatives with respect to $`y`$ of arbitrary order and continuously depends on the initial data.

This fundamental existence and uniqueness Lemma will allow us to obtain a complete classification of all soliton (vortex) solutions to (1), (15).

## II The classification of soliton (vortex) solutions to static coupled SG equations for $`N<\mathrm{}`$

By analogy with (6)-(11), we should find solutions to (1), (15) subject to topological boundary conditions

$$\varphi _n(L)=0,\varphi _n(0)=0\text{mod}\left(\pi \right),$$
(17)
with all possible nonzero sets of the constants $`0`$mod$`\left(\pi \right)`$, compatible with the requirement (3) that all $`\frac{d\varphi _n}{dy}\left(L\right)`$ be equal to each other. This task can be reformulated in terms of the initial value problem at $`y=L`$:

$$\varphi _n(L)=0,\frac{d\varphi _n}{dy}\left(L\right)=2ep\stackrel{~}{H}>0,$$
(18)
where $`\stackrel{~}{H}`$ is so far unknown. According to the Lemma of section I, the solution to (18), satisfying topological conditions at $`y=0`$, will represent in $`[L,L]`$ the unique type of the solution to (17), compatible with the boundary conditions on $`\frac{d\varphi _n}{dy}`$. Using continuous dependence of this solution on $`H`$, we can establish a full range of its existence. This program is carried out to a full extent in Ref.  . Here we present the main results. As could be expected, the only possible soliton solutions for $`H>0`$, $`L<\mathrm{}`$ turn out to be of the vortex-plane type. In addition to the symmetry relations (4), (16), they satisfy the conditions

$$\pi \varphi _n\left(L\right)0,\varphi _n(0)=\pi N_v,\frac{d\varphi _n}{dy}\left(y\right)>0,n=1,\mathrm{},N1.$$
(19)
The range of the existence of the solution with $`N_v1`$ is given by (10), where the superheating (penetration) field $`H_s`$ is now an explicit analytical function of $`N`$ \[$`H_s\left(N<\mathrm{}\right)>H_s\left(N=\mathrm{}\right)\left(ep\lambda _J\right)^1`$\], and $`H_{N_v}`$ is determined by the solution to the boundary value problem

$$\varphi _n\left(L\right)=\pi ,\varphi _n(0)=\pi N_v,n=1,\mathrm{},N1.$$
(20)
Using analogous considerations, it is straightforward to establish the absence of any soliton solutions (including the vortex planes ) for $`H=0`$, $`L<\mathrm{}`$. (The latter result is indirectly confirmed by the numerical simulations of Refs.  for $`H=0`$, performed in a finite interval on the $`y`$ axis: Appropriate topological boundary conditions are not fulfilled in the figures presented therein.)

As a simple illustration of the above mathematical techniques, we present an explicit proof of the absence of single-vortex solutions to a pair of coupled static SG equations ($`N=3`$) for $`H>0`$. In this case, the boundary value problem for a single vortex, positioned at $`n=1`$, is formulated as follows:

$$\varphi _{1,2}\left(L\right)=0,\varphi _1\left(0\right)=\pi ,\varphi _2\left(0\right)=0.$$
(21)
However, this boundary value problem has no solution, meeting the requirement $`\frac{d\varphi _1}{dy}\left(L\right)=\frac{d\varphi _2}{dy}\left(L\right)`$, for any $`\stackrel{~}{H}>0`$. Indeed, the initial value problem (18) for $`N=3`$ has an exact analytical solution on the whole axis $`\mathrm{}<y<+\mathrm{}`$, valid for any $`\stackrel{~}{H}>0`$:

$$\varphi _1(y)=\varphi _2(y)\varphi (y)=\pi +2\text{am}(\frac{y+L}{k\lambda _{J1}}+K\left(k^2\right),k),$$
(22)
$`\lambda _{J1}={\displaystyle \frac{ϵ\lambda _J}{\sqrt{1+ϵ^2}}},k={\displaystyle \frac{1}{\sqrt{1+\left(ep\lambda _{J1}\stackrel{~}{H}\right)^2}}},`$

where $`\lambda _{J1}`$ is the unique Josephson length for $`N=3`$. (Contrary to a remark in the Comment, both the Meissner solution and the vortex planes for $`3N<\mathrm{}`$ are characterized by a distribution of Josephson lengths, whose number is equal to the integer part of $`\frac{N}{2}`$.) According to the Lemma, expression (22) is the unique solution, compatible with the boundary conditions at $`y=L`$ for any $`\stackrel{~}{H}>0`$. As anticipated by (19), it satisfies the condition $`\frac{d\varphi }{dy}>0`$ on the whole axis $`y`$, whereas (3) and (21) imply $`\frac{d\varphi _2}{dy}\left(0\right)<0`$. Note that although the equality $`\varphi _1=\varphi _2`$ could be expected by (16), these symmetry relations were not employed in our proof. (The fact that $`\varphi _1=\varphi _2`$ for $`H>0`$ was verified experimentally on artificial double-junction stacks.)

Although vortex planes are unique soliton solutions for $`H>0`$, $`L<\mathrm{}`$, the situation changes drastically at $`H=0`$, $`L=\mathrm{}`$, when there is no need to bother about boundary conditions on $`\frac{d\varphi _n}{dy}`$. The imposition of topological boundary conditions

$$\varphi _n(y)\underset{y\mathrm{}}{}0,\varphi _n(y)=0\text{mod}(\pi )$$
(23)
automatically ensures the fulfillment of the asymptotic boundary conditions

$$\frac{d\varphi _n}{dy}\left(y\right)\underset{y\mathrm{}}{}0,$$
(24)
by virtue of Eqs. (1), (15) themselves and some elementary theorems of mathematical analysis. As shown in Refs.  , apart from the vortex-plane solution, equations (1), (15) admit for $`H=0`$, $`L=\mathrm{}`$ a variety of topological configurations: single-vortex, vortex-vortex and vortex-antivortex solutions. However, all these solutions are subject to the requirement that the constants $`0`$mod$`\left(\pi \right)`$ be chosen from the set $`0,\pm \pi `$, i.e., each $`\varphi _n`$ can ”accommodate” for $`H=0`$, $`L=\mathrm{}`$ no more than one vortex or antivortex. (This is a typical result of the theory of nonlinear ordinary differential equations that could not be obtained by means of numerical simulations, restricted to a finite interval on the $`y`$ axis.) In particular, the solution representing a single vortex, positioned at $`n=l`$, is specified by the conditions

$`\varphi _l(y)\underset{y\mathrm{}}{}0,{\displaystyle \frac{d\varphi _l}{dy}}\left(y\right)\underset{y\mathrm{}}{}0,\varphi _l(0)=\pi \text{,}{\displaystyle \frac{d\varphi _l}{dy}}\left(0\right)>0\text{,}`$

$$\varphi _n(y)\underset{y\mathrm{}}{}0,\frac{d\varphi _n}{dy}\left(y\right)\underset{y\mathrm{}}{}0,\varphi _n(0)=0\text{,}\frac{d\varphi _n}{dy}\left(0\right)<0\text{, for }nl,$$
(25)
plus a certain existence condition on all $`\frac{d\varphi _n}{dy}\left(0\right)`$. The main distinctive feature of this solution, the negative sign of $`\frac{d\varphi _{nl}}{dy}\left(0\right)`$ in (25), is clearly reproduced by figures 1, 2 in Ref.  . In this context, the remark that the numerical solution ”agrees well with Bulaevskii and Clem-Coffey solutions” is inappropriate: Instead of exact Eqs. (1), references  employ a mathematically ill-defined ”continuum-limit approximation” that completely neglects this feature. (As shown by Farid,  ”solutions” of this type are, in fact, a mathematical fiction.)

Summarizing, the actual domain of the existence of the single-vortex solution is restricted to $`H=0`$, $`L=\mathrm{}`$. Concerning figure 1 in the Comment, based on the numerical calculations for $`H=0`$, the ratio of the self-energy of the vortex plane to that of the single vortex at $`H=0`$ is unimportant, in light of the absence of the latter at $`H>0`$. For the same reason, the expression$`H_{c1}=4e\pi E_{\text{sing}}/\mathrm{\Phi }_0`$” is physically senseless. The actual lower critical field $`H_{c1}`$ is determined by the condition that the Gibbs free energy of the state with a single vortex plane be equal to that of the Meissner state. In the two exactly solvable cases, $`N=\mathrm{}`$ and $`N=3`$ with $`(epH_s)^1L<\mathrm{}`$, it is given by $`H_{c1}=\frac{2}{\pi }H_s`$, as for the single junction. (In general, $`\sqrt{H_{sL}^2H_s^2}<H_{c1}<H_{sL}`$.) Calculations of the number of ”quasiequilibrium fluxon modes” by means of combinatorics are also invalid, because they do not take into account the necessity to satisfy topological boundary conditions together with the boundary conditions on $`\frac{d\varphi _n}{dy}`$. (Figure 6 in Ref.  does not provide any evidence that such conditions for isolated fluxons are fulfilled in the numerical simulations, either.) The comparison  of the soliton vortex-plane solutions with the laminar model for type-II superconductors is a misunderstanding: Alternating superconducting and normal layers, envisaged by this model, have nothing to do with soliton physics and do not possess the property of topological stability.

The statement that ”free energy of any isolated solution is twice the Josephson energy” is incorrect: See, e.g., (A2). However, the self-energy of soliton solutions for $`H=0`$, $`L=\mathrm{}`$ is indeed twice the Josephson energy, which has been established by the use of the first integral of (1) in Ref.  (not Ref.  , as is claimed in the Comment). (From a field-theoretical point of view, this fact is a manifestation of the virial theorem.)

## III Josephson-vortex penetration

Isolated Josephson vortices cannot penetrate layered superconductors and stacked junctions at any $`H>0`$, because they do not form any equilibrium configurations. This problem cannot be circumvented by any consideration of penetration as a ”dynamic process” whose final stage should be an equilibrium flux configuration. The experiments on artificial stacks have clearly shown that Josephson vortices penetrate all the junctions simultaneously and coherently (i.e., in the form of vortex planes), in complete agreement with the scenario of Refs.  .

A description in terms of static equations, accepted in Refs.  and previous publications on single Josephson junctions, implies that a flux relaxation time is much smaller than a characteristic time of the change of the external field, which corresponds to experimental conditions of thermodynamic equilibrium. As in the case of Abrikosov vortices in type-II superconductors, penetration occurs when vortex planes appear at the surface $`y=\pm L`$ (i.e., at $`H=H_{N_v}`$), which signifies the vanishing of the surface barrier. (In contrast, the static SG equations with $`N3`$ do not admit solutions with a single Josephson vortex at the surface.) Mathematically, penetration constitutes a change of the topological type of the solution and should be interpreted as a series of first-order phase transition. Discontinuities related to these phase transitions manifest themselves, in particular, in jumps of magnetization. Such jumps of magnetization were indeed observed experimentally in stacked junctions.

What figure 2 in the Comment actually shows is a violation of the conditions of thermodynamic equilibrium in the numerical simulations. The hypothesis  of ’decomposition of a ”breather”’ is also irrelevant. Certainly, the single ($`N=2`$) time-dependent SG equation for $`H=0`$, $`L=\mathrm{}`$ has an exact analytical solution, known as the ”breather”. However, apart from the fact that time-dependent SG equations have nothing to do with thermodynamic equilibrium, the author of the Comment does not demonstrate the existence of such a solution for $`N3`$, $`H>0`$, $`L<\mathrm{}`$ and does not explain why it should ”decompose”.

## IV Experimental observations of Josephson-vortex configurations: Some concluding remarks

There is no contradiction between our conclusion about the absence of single-vortex configurations at $`H>0`$ and the observations of non-equilibrium isolated vortices in layered high-$`T_c`$ superconductors at $`H=0`$. In the case of Ref.  , isolated vortices can be pinned by structural defects, because their energy is lower and the $`c`$-axis extent is smaller than those of vortex planes at $`H=0`$.

Measurements of $`c`$-axis transport properties do not provide any direct information on flux distribution in stacked junctions. Both the multivaluedness and the aperiodicity of $`I_c(H)`$ dependence, reported in Ref.  , can be explained by the overlap of states with different numbers of vortex planes, discussed in section I.

The Josephson-vortex structure in artificial stacked junctions for $`H>0`$, $`I=0`$ was directly observed in Ref.  (by low-temperature scanning electron microscopy) and in Ref.  (by polarized neutron reflection). As is emphasized in the introduction and sections II, III, all our theoretical conclusions fully agree with these observations. In contrast, the concept of ”isolated fluxons” stands in direct contradiction to the experiments of Refs.  .

In conclusion, a revision of some old prejudices in the theory of Josephson systems should begin with the case of the single junction ($`N=2`$). In spite of ”extensive theoretical studies” for almost four decades, there is still an erroneous belief that Josephson vortices ”do not exist” in single junctions with $`L\lambda _J`$, which is refuted by the exact, closed-form solution (6)-(11). \[As shown in Refs.  , exactly the ordinary Josephson vortices and the vortex planes account for the occurrence of the well-known Fraunhofer pattern of $`I_c(H)`$ for small $`L.`$\] The lack of the understanding of analytical properties of the single SG equation is one of the reasons why the solution to the more complex coupled SG equations was not obtained in previous publications. Adequate mathematical techniques to tackle this problem are proposed in Refs.  . The author of the Comment himself could profit from the application of these techniques to his own research.

## A A new method of exact minimization of the Gibbs free energy of the infinite LD model

Equations (1), (2) have the first integral

$$C\left(H\right)\underset{n}{}\mathrm{cos}\varphi _n(y)=\frac{ϵ^2\lambda _J^2}{2}\underset{n}{}G_{\mathrm{}}(n,m)\frac{d\varphi _n(y)}{dy}\frac{d\varphi _n(y)}{dy},$$
(A1)
$`C\left(H\right)=2\left(ep\lambda _JH\right)^2{\displaystyle \frac{V}{pL}}+{\displaystyle \underset{n}{}}\mathrm{cos}\varphi _n(L),`$

$`G_{\mathrm{}}(n,m)={\displaystyle \frac{\mu ^{\left|nm\right|}}{2ϵ\sqrt{1+\frac{ϵ^2}{4}}}},\mu =1+{\displaystyle \frac{ϵ^2}{2}}ϵ\sqrt{1+{\displaystyle \frac{ϵ^2}{4}}}.`$

The term on the right-hand side of (A1) is the density of electromagnetic energy. For $`r(T)1`$ and $`HH_{c2}`$, one can eliminate the electromagnetic-energy term from the LD free-energy functional by the use of the Maxwell equations and (A1), obtaining

$`\mathrm{\Omega }_{LD}[\phi _n;H]={\displaystyle \frac{H_c^2(T)V}{4\pi }}[{\displaystyle \frac{1}{2}}+r(T)[1+4\left(ep\lambda _JH\right)^2+{\displaystyle \frac{2pL}{V}}{\displaystyle \underset{n}{}}\mathrm{cos}\varphi _n\left(L\right)`$

$$\frac{2p}{V}\underset{n}{}\underset{L}{\overset{L}{}}dy\mathrm{cos}\varphi _n(y)\frac{4ep^2\lambda _JLH}{V}\underset{n}{}[\varphi _n(L)\varphi _n(L)]]].$$
(A2)
Taking the variation of (A2) with respect to $`\phi _n(y)`$, one immediately arrives at the relations (5). (Mathematically, the elimination of the electromagnetic-energy term is equivalent to the elimination of the intrinsic constraint in Ref. .) Unfortunately, the existence of the first integral (A1) was not noticed in any previous publications on the infinite LD model, which partly explains the difficulties with the minimization of the Gibbs free energy with respect to $`\phi _n`$.