text:
# Optimal Renormalization-Group Improvement of the Perturbative Series for the 𝑒⁺𝑒⁻-Annihilation Cross-Section

## Abstract

Using renormalization-group methods, we derive differential equations for the all-orders summation of logarithmic corrections to the QCD series for $`R(s)=\sigma (e^+e^{}\mathrm{hadrons})/\sigma (e^+e^{}\mu ^+\mu ^{})`$, as obtained from the imaginary part of the purely-perturbative vector-current correlation function. We present explicit solutions for the summation of leading and up to three subsequent subleading orders of logarithms. The summations accessible from the four-loop vector-correlator not only lead to a substantial reduction in sensitivity to the renormalization scale, but necessarily impose a common infrared bound on perturbative approximations to $`R(s)`$, regardless of the infrared behaviour of the true QCD couplant.

For center-of-mass squared-energy $`s`$, QCD corrections to $`R(s)\sigma (e^+e^{}\mathrm{hadrons})/\sigma (e^+e^{}\mu ^+\mu ^{})`$ are scaled by a perturbative QCD series ($`S`$):

$$R(s)=3\underset{f}{}Q_f^2S[\frac{\alpha _s(\mu )}{\pi },\mathrm{log}\left(\frac{\mu ^2}{s}\right)].$$
(1)
This series is extracted from the imaginary part of the $`\overline{\mathrm{MS}}`$ vector-current correlation function sgg ; ces ,

$$S[x,L]=1+\underset{n=1}{\overset{\mathrm{}}{}}x^n\underset{m=0}{\overset{n1}{}}T_{n,m}L^m,$$
(2)
with coefficients $`T_{n,m}`$ tabulated in Table 1 for 3–5 active flavors, as appropriate for the choice of the center-of-mass squared-energy $`s`$. Each order of this series depends upon the $`\overline{\mathrm{MS}}`$ renormalization scale parameter $`\mu `$, both through the couplant

$$x(\mu )\alpha _s(\mu )/\pi $$
(3)
and through powers of the logarithm

$$L(\mu )\mathrm{log}(\mu ^2/s).$$
(4)
Nevertheless, the all-orders series $`S`$ must ultimately be independent of renormalization scale. $`R(s)`$ is a measurable physical quantity necessarily independent of $`\mu `$, the artificial scale entering QCD calculations as a by-product of the regulation of Feynman-diagrammatic infinities. Hence,

$$0=\mu ^2\frac{\mathrm{d}S[x(\mu ),L(\mu )]}{\mathrm{d}\mu ^2}=\left(\frac{}{L}+\beta (x)\frac{}{x}\right)S[x,L].$$
(5)
The above renormalization group equation (RGE) is simply a chain-rule relation in which

$$\beta (x)\mu ^2\frac{\mathrm{d}x(\mu )}{\mathrm{d}\mu ^2}=x^2\underset{k=0}{\overset{\mathrm{}}{}}\beta _kx^k,$$
(6)
where known rvl  $`\overline{\mathrm{MS}}`$ $`\beta `$-function coefficients $`\beta _k`$ are also tabulated in Table 1. Thus, the RGE (5) is generally employed to provide scale dependence to the couplant $`x`$.

“Optimal” renormalization-group (RG) improvement is the inclusion of every term in a perturbative series of the form (2) that can be extracted by RG-methods from a perturbative computation to a given order utica . For example, a next-to-next-to-leading (NNL) order perturbative calculation determines only the coefficients $`T_{1,0}`$, $`T_{2,0}`$, and $`T_{2,1}`$ of the series (2). However, the RGE (5) can be utilized to determine all coefficients $`T_{k+1,k}`$ and $`T_{k+2,k}`$ within the series (2). The contributions of this infinite set of coefficients may then be summed analytically, as described below, thereby providing an “optimal RG-improvement” of the NNL expression.

The series $`S[x,L]`$, as defined in Eq. (2), may be rearranged in the following form:

$$S[x,L]=1+\underset{n=1}{\overset{\mathrm{}}{}}x^nS_n(xL),$$
(7)
where

$$S_k(u)\underset{n=k}{\overset{\mathrm{}}{}}T_{n,nk}u^{nk}.$$
(8)
Given knowledge of the $`k^{th}`$-order series coefficient $`T_{k,0}=S_k(0)`$, one can obtain $`S_k(x(\mu )L(\mu ))`$ explicitly, thereby summing over the entire set of $`k^{th}`$-order subleading logarithms contributing to the series (7). If we substitute the $`\beta `$-function series (6) into the RGE (5), we find that the aggregate coefficient of $`x^nL^{np}`$ vanishes $`(np)`$ provided the following recursion relation is upheld:

$$0=(np+1)T_{n,np+1}\underset{\mathrm{}=0}{\overset{p2}{}}(n\mathrm{}1)\beta _{\mathrm{}}T_{n\mathrm{}1,np}.$$
(9)
For example, if $`p=2`$, this recursion relation $`[T_{n,n1}=\beta _0T_{n1,n2}]`$ relates all leading-logarithm coefficients $`T_{n,n1}`$ within the series (2) to the known coefficient $`T_{1,0}=1`$, thereby enabling one to sum all orders of the leading-logarithm contributions

$$xS_1(xL)=x\underset{n=1}{\overset{\mathrm{}}{}}T_{n,n1}(xL)^{n1}=\frac{x}{1\beta _0xL}$$
(10)
to the series $`S[x,L]`$.

More generally, the recursion relation (9) may be utilized to obtain a succession of first-order inhomogeneous linear differential equations for the functions $`S_k(u)`$ within Eq. (7). If one multiplies Eq. (9) by $`u^{np}`$ and then sums over $`n`$ from $`n=p`$ to $`\mathrm{}`$, one finds from the definition (8) of $`S_k(u)`$ that

$$0=\frac{\mathrm{d}S_{p1}}{\mathrm{d}u}u\underset{\mathrm{}=0}{\overset{p2}{}}\beta _{\mathrm{}}\frac{\mathrm{d}S_{p\mathrm{}1}}{\mathrm{d}u}\underset{\mathrm{}=0}{\overset{p2}{}}(p\mathrm{}1)\beta _{\mathrm{}}S_{p\mathrm{}1},$$
(11)
which can be trivially rearranged $`(k=p1)`$ into a set of first-order linear differential equations

$$\frac{\mathrm{d}S_k}{\mathrm{d}u}\frac{k\beta _0}{(1\beta _0u)}S_k=\frac{1}{(1\beta _0u)}\underset{\mathrm{}=1}{\overset{k1}{}}\beta _{\mathrm{}}\left(u\frac{\mathrm{d}}{\mathrm{d}u}+k\mathrm{}\right)S_k\mathrm{},$$
(12)
with initial conditions $`S_k(0)=T_{k,0}`$.

Noting that $`S_0(u)=1`$ and that $`T_{1,0}=1`$ regardless of the number of active flavors, we find the first four solutions of (12) to be

$$S_1(xL)=\frac{1}{(1\beta _0xL)},$$
(13)
$$S_2(xL)=\frac{T_{2,0}\frac{\beta _1}{\beta _0}\mathrm{log}\left(1\beta _0xL\right)}{\left(1\beta _0xL\right)^2},$$
(14)
$$\begin{array}{cc}\hfill S_3(xL)=& \left(\frac{\beta _1^2}{\beta _0^2}\frac{\beta _2}{\beta _0}\right)/\left(1\beta _0xL\right)^2\hfill \\ & +\frac{\left\{T_{3,0}\left(\frac{\beta _1^2}{\beta _0^2}\frac{\beta _2}{\beta _0}\right)\frac{\beta _1}{\beta _0}\left(2T_{2,0}+\frac{\beta _1}{\beta _0}\right)\mathrm{log}(1\beta _0xL)+\frac{\beta _1^2}{\beta _0^2}\mathrm{log}^2(1\beta _0xL)\right\}}{\left(1\beta _0xL\right)^3}\hfill \end{array}$$
(15)
$$\begin{array}{cc}\hfill S_4(xL)=& \frac{1}{2}\left[\frac{\beta _1}{\beta _0}\left(\frac{\beta _1^2}{\beta _0^2}2\frac{\beta _2}{\beta _0}\right)+\frac{\beta _3}{\beta _0}\right](1\beta _0xL)^2+\left(2T_{2,0}+\frac{\beta _1}{\beta _0}\right)\left(\frac{\beta _1^2}{\beta _0^2}\frac{\beta _2}{\beta _0}\right)(1\beta _0xL)^3\hfill \\ & +2\frac{\beta _1}{\beta _0}\left(\frac{\beta _2}{\beta _0}\frac{\beta _1^2}{\beta _0^2}\right)(1\beta _0xL)^3\mathrm{log}(1\beta _0xL)\hfill \\ & +\left[T_{4,0}+\frac{\beta _3}{2\beta _0}\frac{1}{2}\frac{\beta _1^3}{\beta _0^3}2T_{2,0}\left(\frac{\beta _1^2}{\beta _0^2}\frac{\beta _2}{\beta _0}\right)\right](1\beta _0xL)^4\hfill \\ & +\frac{\beta _1}{\beta _0}\left[2\frac{\beta _1^2}{\beta _0^2}3\frac{\beta _2}{\beta _0}2T_{2,0}\frac{\beta _1}{\beta _0}3T_{3,0}\right](1\beta _0xL)^4\mathrm{log}(1\beta _0xL)\hfill \\ & +\frac{\beta _1}{\beta _0}\left[\frac{5\beta _1^2}{2\beta _0^2}+3T_{2,0}\frac{\beta _1}{\beta _0}\right](1\beta _0xL)^4\mathrm{log}^2(1\beta _0xL)\frac{\beta _1^3}{\beta _0^3}(1\beta _0xL)^4\mathrm{log}^3(1\beta _0xL).\hfill \end{array}$$
(16)
To explore the near-infrared regime of perturbative QCD, we specialize to the case of three active flavors. Using Table 1 $`n_f=3`$ values within Eqs. (13)–(16), we find that the version of the series (7) which incorporates full summation of leading and two subsequent subleading orders of logarithms is given by

$$S^{(\mathrm{\Sigma })}[x(\mu ),L(\mu )]=1+xS_1(xL)+x^2S_2(xL)+x^3S_3(xL),$$
(17)
where $`(u=xL)`$

$$S_1(u)=1/(19u/4),$$
(18)
$$S_2(u)=\frac{1.63982\frac{16}{9}\mathrm{log}(19u/4)}{(19u/4)^2},$$
(19)
$$S_3(u)=\frac{1.31057}{(1\frac{9}{4}u)^2}+\frac{\left\{8.973338.99096\mathrm{log}\left(1\frac{9u}{4}\right)+3.16049\mathrm{log}^2\left(1\frac{9u}{4}\right)\right\}}{(1\frac{9u}{4})^3}.$$
(20)
Moreover, we note from Eq. (16) that

$$\begin{array}{cc}\hfill S_4(xL)=& \frac{5.35589}{\left(1\frac{9}{4}xL\right)^2}+\frac{\left[6.62811+4.65981\mathrm{log}\left(1\frac{9}{4}xL\right)\right]}{\left(1\frac{9}{4}xL\right)^3}\hfill \\ & +\frac{\left[T_{4,0}+11.9840+31.8738\mathrm{log}\left(1\frac{9}{4}xL\right)+29.5946\mathrm{log}^2\left(1\frac{9}{4}xL\right)5.61866\mathrm{log}^3\left(1\frac{9}{4}xL\right)\right]}{\left(1\frac{9}{4}xL\right)^4},\hfill \end{array}$$
(21)
thereby providing for inclusion of the $`x^4S_4(xL)`$ contribution to the series (7) for three active flavors. The series coefficient $`T_{4,0}`$ appearing in Eq. (21) has not yet been calculated perturbatively, which is why we have not included the $`x^4S_4(xL)`$ contribution to $`S[x,L]`$ in Eq. (17). (An asymptotic Padé approximant estimate $`T_{4,0}1.90`$ for the $`n_f=3`$ case is presented in ref. ces .)

To examine whether the summation of leading and subsequent subleading logarithm factors decreases dependence on the unphysical renormalization-scale parameter $`\mu `$, we compare the $`\mu `$-dependence of Eq. (17) for a fixed value of $`s`$ to that of the $`n_f=3`$ version of the series (2) truncated after four-loop-order ($`4\mathrm{}`$) contributions to the vector-current correlation function:

$$S^{(4\mathrm{})}[x(\mu ),L(\mu )]=1+x+(1.63982+9L/4)x^2+(10.2839+11.3792L+81L^2/16)x^3.$$
(22)
Such $`\mu `$-dependence enters Eqs. (17) and (22) both through $`L=\mathrm{log}(\mu ^2/s)`$ and through $`x=x(\mu )`$, which is assumed to evolve via Eq. (6) (with $`n_f=3`$ choices for $`\beta _{03}`$) from an initial value choice $`x(m_\tau )=\alpha _s(m_\tau )/\pi =0.33/\pi `$ aleph ; ptg . Figure 1 displays a comparison of the $`\mu `$-dependence of Eqs. (17) and (22) at fixed $`s=1.5\mathrm{GeV}^2`$ . Although both expressions exhibit little variation with $`\mu `$ over the $`1.3\mathrm{GeV}\mu 3\mathrm{GeV}`$ range, we see that Eq. (17) exhibits much less variation with $`\mu `$ in the near-infrared regime below $`1.3\mathrm{GeV}`$. These results clearly indicate that renormalization-scale-invariance is more effectively upheld via the summations of leading and subsequent subleading orders of logarithms that occur within Eq. (17). We emphasize that Eqs. (17) and (22) both follow from “RG-improvement” of the same calculational information \[the coefficients $`T_{1,0},T_{2,0},`$ and $`T_{3,0}`$\]; simply put, such RG-improvement is more effectively implemented in Eq. (17) than in Eq. (22).

The usual prescription for obtaining the purely-perturbative (non-power-law) QCD contributions to $`R(s)`$ at four-loop order ptg ; fjy  is to set $`\mu =\sqrt{s}`$ within Eq. (22), and then to substitute the resulting series,

$$S^{(4\mathrm{})}[x(\sqrt{s}),L(\sqrt{s})]=1+x(\sqrt{s})+1.63982x^2(\sqrt{s})10.2839x^3(\sqrt{s}),$$
(23)
into Eq. (1). \[Note from Eq. (4) that $`L(\sqrt{s})=0`$.\] This prescription follows from the presumed renormalization-scale invariance of the truncated series (22), thereby leading to an expression that depends only on the physical scale $`s`$ for the electron-positron annihilation process. Note that all $`s`$-dependence of Eq. (17) resides entirely in the variable $`L`$. Let us first fix $`\mu =m_\tau `$ within Eq. (17) so as to incorporate the benchmark couplant value $`x(m_\tau )=0.33/\pi `$ everywhere $`x`$ appears in Eqs. (18)–(20). With this choice, the following summation-of-logarithms series can be substituted into Eq. (1) to obtain $`R(s)`$:

$$\begin{array}{cc}\hfill S^{(\mathrm{\Sigma })}[x(m_\tau ),L(m_\tau )]=& 1+\frac{0.33}{\pi }S_1\left(\frac{0.33}{\pi }\mathrm{log}\left(\frac{m_\tau ^2}{s}\right)\right)\hfill \\ & +\left(\frac{0.33}{\pi }\right)^2S_2\left(\frac{0.33}{\pi }\mathrm{log}\left(\frac{m_\tau ^2}{s}\right)\right)+\left(\frac{0.33}{\pi }\right)^3S_3\left(\frac{0.33}{\pi }\mathrm{log}\left(\frac{m_\tau ^2}{s}\right)\right).\hfill \end{array}$$
(24)
In Figure 2 we compare the $`s`$-dependence of this series to that of Eq. (23), for which all $`s`$ dependence resides in the evolution of $`x(\sqrt{s})`$. To make this comparison, such evolution is anchored to the initial value $`x(m_\tau )=0.33/\pi `$ via the differential equation (6) with $`n_f=3`$ values for $`\beta _{03}`$ \[Table 1\]. This initial value ensures that the series (23) and (24) coincide when $`\sqrt{s}=m_\tau `$. Figure 2 shows that both series continue to coincide over the range $`750\mathrm{MeV}\sqrt{s}m_\tau `$. For values of $`\sqrt{s}`$ less than 750 MeV, however, the truncated series (23) drops off quite suddenly at $`\sqrt{s}650\mathrm{MeV}`$, a consequence of the large negative coefficient of $`x^3\left(\sqrt{s}\right)`$, whereas the series (24) \[as obtained from the full summation-of-logarithms series (17)\] continues to probe the infrared domain of $`R(s)`$ even for values of $`\sqrt{s}400\mathrm{MeV}`$. In short, the summation of all leading and subsequent two subleading logarithms within the perturbative series (2) serves to extend the domain of the $`R(s)`$ series further into the infrared. This property, as well as the reduced renormalization-scale dependence evident in Figure 1, suggests that such summation is particularly appropriate for the near-infrared region characterizing sum-rule applications of purely-perturbative QCD corrections to current-correlation functions.

It is also evident from Figure 2 that the domain of the summation-of-logarithms series (24) manifests a singularity below $`\sqrt{s}=400\mathrm{MeV}`$, despite the fact that the $`s`$-dependence of the series (24) is decoupled entirely from any infrared behavior of the couplant $`x`$, which is held constant at $`x(m_\tau )`$. To understand this restriction on the domain of $`R(s)`$, we first note that each summation (13)–(16) becomes singular when $`1\beta _0xL0`$. Such resummation singularities have also been observed to occur in completely different contexts, including deep inelastic structure functions vogt .

The singularity property of eqs. (13)–(16) is upheld for all summations $`S_k(xL)`$. The solution to the differential equation (12) is necessarily of the form

$$S_k(xL)=\frac{T_{k,0}}{(1\beta _0xL)^k}+\left(\text{Particular solution depending on}\{S_{k1},S_{k2},\mathrm{},S_1\}\right).$$
(25)
Since the coefficients $`T_{k,0}`$ are results of kth-order Feynman diagram calculations, the kth-order pole in (25) at $`1\beta _0xL=0`$ is genuine and will not be canceled by particular-solution contributions that are sensitive to at most $`(k1)^{th}`$-order Feynman-diagrammatic coefficients $`\{T_{1,0},T_{2,0},\mathrm{},T_{k1,0}\}`$. For a given choice of renormalization scale $`\mu `$, this singularity implies \[via Eqs. (3) and (4)\] that each summation $`S_k(xL)`$ within the full series (7) becomes singular for a sufficiently small value of $`s`$:

$$1\beta _0\frac{\alpha _s(\mu )}{\pi }\mathrm{log}\left(\frac{\mu ^2}{s_{min}}\right)=0s_{min}=\mu ^2\mathrm{exp}\left(\frac{\pi }{\beta _0\alpha _s(\mu )}\right).$$
(26)
For example, if the renormalization scale $`\mu `$ is chosen (as in Figure 2) to be $`m_\tau `$, a choice for which $`\alpha _s(m_\tau )(=0.33\pm 0.02`$ ptg ) is phenomenologically accessible, then each term in the series (7) is seen to become progressively more singular as $`s`$ approaches $`m_\tau ^2\mathrm{exp}[4\pi /(9\times 0.33)]=(215\mathrm{MeV})^2`$ from above. Furthermore, Figure 3 shows that the low-energy behaviour of the resummed expression (17) is only weakly dependent on $`\mu `$ in the region $`\sqrt{s}>600\mathrm{MeV}`$ for choices of $`\mu `$ between $`0.6m_\tau `$ and $`1.6m_\tau `$, with the appearance of singular points corresponding to $`195\mathrm{MeV}<s_{min}<250\mathrm{MeV}`$ from (26).

It is to be emphasized that this infrared boundary on the physical scale $`s`$ entering term-by-term within the series (7) is not a manifestation of any infrared boundary cem  on the evolution of the QCD couplant $`x(\mu )`$. Even if the higher order contributions to the $`\beta `$-function (6) were somehow to conspire to allow the couplant to be well-behaved in the infrared region \[e.g., to have infrared-stable fixed point behavior\], the restriction (26) would still apply upon making a specific choice of the renormalization-scale parameter $`\mu `$ and its corresponding value of $`\alpha _s(\mu )`$. Curiously, though, this infrared restriction on $`s`$ can be easily shown to coincide with the “infrared-slavery” Landau singularity $`\mathrm{\Lambda }`$ associated with naive evolution of the QCD couplant via a one-loop $`\beta `$-function. The one loop version of Eq. (6),

$$\mu ^2\frac{\mathrm{d}x}{\mathrm{d}\mu ^2}=\beta _0x^2,x=\frac{\alpha _s}{\pi },$$
(27)
is satisfied by the relation

$$\frac{\alpha _s}{\pi }=\frac{1}{\beta _0\mathrm{log}(\mu ^2/\mathrm{\Lambda }^2)}$$
(28)
which is equivalent to Eq. (26) provided $`s_{min}`$ is identified with $`\mathrm{\Lambda }^2`$. Indeed, in a one-loop world \[$`x^1\mathrm{}(\mu )=1/\beta _0\mathrm{log}\left(\mu ^2/\mathrm{\Lambda }^2\right)`$\] where $`\mathrm{\Lambda }`$ serves as a universal infrared boundary, the one-loop analogues of the summation-of-logarithms series (17) and the truncated series (23) are necessarily equivalent:

$$\begin{array}{cc}\hfill S_1\mathrm{}^{(\mathrm{\Sigma })}& =1+x^{(1\mathrm{})}(\mu )S_1\left[x^{(1\mathrm{})}(\mu )\mathrm{log}\left(\mu ^2/s\right)\right]=1+\frac{1}{\beta _0\mathrm{log}\left(\frac{\mu ^2}{\mathrm{\Lambda }^2}\right)}\left\{\frac{1}{1\beta _0\left[\frac{1}{\beta _0\mathrm{log}\left(\frac{\mu ^2}{\mathrm{\Lambda }^2}\right)}\right]\mathrm{log}\left(\frac{\mu ^2}{s}\right)}\right\}\hfill \\ & =1+\frac{1}{\beta _0\mathrm{log}\left(\frac{s}{\mathrm{\Lambda }^2}\right)}=1+x^{(1\mathrm{})}\left(\sqrt{s}\right)\hfill \end{array}$$
(29)
We find it remarkable that $`\mathrm{\Lambda }`$, the one-loop couplant’s Landau pole, persists as an infrared boundary on the domain of each summation contributing to Eq. (7), the summation-of-logarithms formulation of the perturbative series within $`R(s)`$. Consequently, $`s=\mathrm{\Lambda }^2`$ serves as an infrared boundary for any approximation to $`R(s)`$ involving the truncation of the series (7), such as the expression (17) obtained via optimal RG-improvement of the four-loop vector-current correlation function’s imaginary part.

###### Acknowledgements.

We are grateful to V.A. Miransky for useful discussions, and to the Natural Sciences & Engineering Research Council of Canada (NSERC) for financial support. We also wish to thank SUNY Institute of Technology for hospitality in the preparation of this research. The work of AHF has been supported in part by the 2002 State of New York/UUP Professional Development Committee, and the 2002 Summer Grant from the School of Arts and Sciences, SUNY Institute of Technology. The work of MRA has been supported in part by an internal summer research grant from Mount Allison University. AS wishes to acknowledge support from an NSERC Summer Undergraduate Research Award.


format: text

text:
# Problems of G and Multidimensional Models

## 1 Introduction

Dirac’s Large Numbers Hypothesis (LNH) is the origin of many theoretical explorations of time-varying G. According to LNH, the value of $`\dot{G}/G`$ should be approximately the Hubble rate. Although it has become clear in recent decades that the Hubble rate is too high to be compatible with experiment, the enduring legacy of Dirac’s bold stroke is the acceptance by modern theories of non-zero values of $`\dot{G}/G`$ as being potentially consistent with physical reality. There are three problems related to $`G`$, which origin lies mainly in unified models predictions: 1) absolute $`G`$ measurements, 2) possible time variations of $`G`$, 3) possible range variations of $`G`$ – non-Newtonian, or new interactions. For 1) and 3) see .

After the original Dirac hypothesis some new ones appeared and also some generalized theories of gravitation admitting the variations of an effective gravitational coupling. We may single out three stages in the development of this field:

1. Study of theories and hypotheses with variations of FPC, their predictions and confrontation with experiments (1937-1977).
2. Creation of theories admitting variations of an effective gravitational constant in a particular system of units, analyses of experimental and observational data within these theories  (1977-present).
3. Analyses of FPC variations within unified models  (present).

Different theoretical schemes lead to temporal variations of the effective gravitational constant:

1. Empirical models and theories of Dirac type, where $`G`$ is replaced by $`G(t)`$.
2. Numerous scalar-tensor theories of Jordan-Brans-Dicke type where $`G`$ depending on the scalar field $`\sigma (t)`$ appears.
3. Gravitational theories with a conformal scalar field arising in different approaches .
4. Multidimensional unified theories in which there are dilaton fields and effective scalar fields appearing in our 4-dimensional spacetime from additional dimensions . They may help also in solving the problem of a variable cosmological constant from Planckian to present values.

A striking feature of most modern scalar-tensor and unification theories, is that they do not admit a unique and universal constant values of physical constants and of the Newtonian gravitational coupling constant $`G`$ in particular. In this paper we briefly set out the results of some calculations which have been carried out for various theories, and we discuss various bounds that may be suggested by multidimensional theories. Although the bounds on G-dot and $`G(r)`$ are in some classes of theories rather wide on purely theoretical grounds as a result of adjustable parameters, we note that observational data concerning other phenomena may place limits on the possible range of these adjustable parameters.

Here we limit ourselves to the problem of G-dot (for $`G(r)`$ see ). We show that the various theories predict the value of $`\dot{G}/G`$ to be $`10^{12}/yr`$ or less. The significance of this fact for experimental and observational determinations of the value of or upper bound on G-dot is the following: any determination with error bounds significantly below $`10^{12}`$ will typically be compatible with only a small portion of existing theoretical models and will therefore cast serious doubt on the viability of all other models. In short, a tight bound on G-dot, in conjunction with other astrophysical observations, will be a very effective ”theory killer.”

Some estimations for G-dot were done long ago in the frames of general scalar tensor theories using values of cosmological parameters ($`\mathrm{\Omega }`$, $`H`$, $`q`$ etc) . With modern values they predict $`\dot{G}/G`$ at the level of $`10^{12}/yr`$ and less (see recent estimations of A. Miyazaki , predicting time variations of $`G`$ at the level of $`10^{13}yr^1`$) for the Machian-type cosmological solution in the Brans-Dicke theory).

The most reliable experimental bounds on $`\dot{G}/G`$ (radar ranging of spacecraft dynamics ) and laser lunar ranging  give the limit of $`10^{12}/yr`$).

## 2 G-dot in $`(4+N)`$-dimensional cosmology with multicomponent anisotropic fluid

We consider here a $`(4+N)`$-dimensional cosmology with an isotropic 3-space and an arbitrary Ricci-flat internal space. The Einstein equations provide a relation between $`\dot{G}/G`$ and other cosmological parameters.

### 2.1 The model

Let us consider $`(4+N)`$-dimensional theory described by the action

$$S_g=\frac{1}{2\kappa ^2}d^{4+N}x\sqrt{g}R,$$
(1)
where $`\kappa ^2`$ is the fundamental gravitational constant. Then the gravitational field equations are

$$R_P^M=\kappa ^2(T_P^M\delta _P^M\frac{T}{N+2}),$$
(2)
where $`T_P^M`$ is a $`(4+N)`$-dimensional energy-momentum tensor, $`T=T_M^M`$, and $`M,P=0,\mathrm{},N+3`$.

For the $`(4+N)`$-dimensional manifold we assume the structure

$$M^{4+N}=R\times M_k^3\times K^N,$$
(3)
where $`M_k^3`$ is a 3-dimensional space of constant curvature,
$`M_k^3=S^3,R^3,L^3`$ for $`k=+1,0,1`$, respectively, and $`K^N`$ is a $`N`$-dimensional compact Ricci-flat Riemann manifold.

The metric is taken in the form

$$g_{MN}dx^Mdx^N=dt^2+a^2(t)g_{ij}^{(3)}(x^k)dx^idx^j+b^2(t)g_{mn}^{(N)}(y^p)dy^mdy^n,$$
(4)
where $`i,j,k=1,2,3;m,n,p=4,\mathrm{},N+3;g_{ij}^{(3)},g_{mn}^{(N)},a(t)`$ and $`b(t)`$ are, respectively, the metrics and scale factors for $`M_k^3`$ and $`K^N`$. For $`T_P^M`$ we adopt the expression of the multicomponent (anisotropic) fluid form

$$(T_P^M)=\underset{\alpha =1}{\overset{m}{}}diag(\rho ^\alpha (t),p_3^\alpha (t)\delta _j^i,p_N^\alpha (t)\delta _n^m).$$
(5)
Under these assumptions the Einstein equations take the form

$`{\displaystyle \frac{3\ddot{a}}{a}}+{\displaystyle \frac{N\ddot{b}}{b}}={\displaystyle \frac{\kappa ^2}{N+2}}{\displaystyle \underset{\alpha =1}{\overset{m}{}}}[(N+1)\rho ^\alpha 3p_3^\alpha Np_N^\alpha ],`$ (6)
$`{\displaystyle \frac{2k}{a^2}}+{\displaystyle \frac{\ddot{a}}{a}}+{\displaystyle \frac{2\stackrel{2}{\stackrel{.}{a}}}{a^2}}+{\displaystyle \frac{N\stackrel{.}{a}\stackrel{.}{b}}{ab}}={\displaystyle \frac{\kappa ^2}{N+2}}{\displaystyle \underset{\alpha =1}{\overset{m}{}}}[\rho ^\alpha +(N1)p_3^\alpha Np_N^\alpha ],`$ (7)
$`{\displaystyle \frac{\ddot{b}}{b}}+(N1){\displaystyle \frac{\stackrel{2}{\stackrel{.}{b}}}{b^2}}+{\displaystyle \frac{3\stackrel{.}{a}\stackrel{.}{b}}{ab}}={\displaystyle \frac{\kappa ^2}{N+2}}{\displaystyle \underset{\alpha =1}{\overset{m}{}}}[\rho ^\alpha 3p_3^\alpha +2p_N^\alpha ].`$ (8)
The 4-dimensional density is

$$\rho ^{\alpha ,(4)}(t)=_Kd^Ny\sqrt{g^{(N)}}b^N(t)\rho ^\alpha (t)=\rho ^\alpha (t)b^(t),$$
(9)
where we have normalized the factor $`b(t)`$ by putting

$$_Kd^Ny\sqrt{g^{(N)}}=1.$$
(10)
On the other hand, to get the 4-dimensional gravity equations one should put $`8\pi G(t)\rho ^{\alpha (4)}(t)=\kappa ^2\rho ^\alpha (t)`$. Consequently, the effective 4-dimensional gravitational “constant” $`G(t)`$ is defined by

$$8\pi G(t)=\kappa ^2b^N(t)$$
(11)
whence its time variation is expressed as

$$\dot{G}/G=N\dot{b}/b.$$
(12)
### 2.2 Cosmological parameters

Some inferences concerning the observational cosmological parameters can be extracted just from the equations without solving them . Indeed, let us define the Hubble parameter $`H`$, the density parameters $`\mathrm{\Omega }^\alpha `$ and the ”deceleration” parameter $`q`$ referring to a fixed instant $`t_0`$ in the usual way

$$H=\dot{a}/a,\mathrm{\Omega }^\alpha =8\pi G\rho ^{\alpha ,(4)}/3H^2=\kappa ^2\rho ^\alpha /3H^2,q=a\ddot{a}/\dot{a}^2.$$
(13)
Besides, instead of $`G`$ let us introduce the dimensionless parameter

$$g=\dot{G}/GH=Na\dot{b}/\dot{a}b.$$
(14)
Then, excluding $`b`$ from (6) and (8), we get

$$\frac{N1}{3N}g^2g+q\underset{\alpha =1}{\overset{m}{}}A^\alpha \mathrm{\Omega }^\alpha =0$$
(15)
with

$$A^\alpha =\frac{1}{N+2}[2N+1+3(1N)\nu _3^\alpha +3N\nu _N^\alpha ],$$
(16)
where

$$\nu _3^\alpha =p_3^\alpha /\rho ^\alpha ,\nu _N^\alpha =p_N^\alpha /\rho ^\alpha ,\rho ^\alpha >0.$$
(17)
When $`g`$ is small we get from (15)

$$gq\underset{\alpha =1}{\overset{m}{}}A^\alpha \mathrm{\Omega }^\alpha .$$
(18)
Note that (18) for $`N=6`$, $`m=1`$, $`\nu _3^1=\nu _6^1=0`$ (so that $`A^1=13/8`$) coincides with the corresponding relation of Wu and Wang  obtained for large times in case $`k=1`$ (see also ).

If $`k=0`$, then in addition to (18), one can obtain a separate relation between $`g`$ and $`\mathrm{\Omega }^\alpha `$, namely,

$$\frac{N1}{6N}g^2g+1\underset{\alpha =1}{\overset{m}{}}\mathrm{\Omega }^\alpha =0$$
(19)
(this follows from the Einstein equation $`R_0^0\frac{1}{2}R=\kappa ^2T_0^0`$, which is certainly a linear combination of (6)-(8).

The present observational upper bound on $`g`$ is

$$g\begin{array}{c}<\hfill \\ \hfill \end{array}0.1$$
(20)
if we take in accord with

$$\stackrel{.}{G}/G\begin{array}{c}<\hfill \\ \hfill \end{array}0.6\times 10^{11}(y^1)$$
(21)
and $`H=(0.7\pm 0.1)\times 10^{11}(y^1)70\pm 10(km/s.Mpc)`$.

### 2.3 Two-component example: dust + $`(N1)`$-brane

Let us consider two component case: $`m=2`$. Let the first component (called ”matter”) be a dust, i.e.

$$\nu _3^1=\nu _N^1=0,$$
(22)
and the second one (called ”quintessence”) be a $`(N1)`$-brane, i.e.

$$\nu _3^2=1,\nu _N^2=1.$$
(23)
We remind that as it was mentioned in  the multidimensional cosmological model on product manifold $`𝐑\times M_1\times \mathrm{}\times M_n`$ with fields of forms (for review see ) may be described in terms of multicomponent ”perfect” fluid  with the following equations of state for $`\alpha `$-s component: $`p_i^\alpha =\rho ^\alpha `$ if $`p`$-brane worldvolume contains $`M_i`$ and $`p_i^\alpha =\rho ^\alpha `$ in opposite case. Thus, the field of form matter leads us either to $`\mathrm{\Lambda }`$-term, or to stiff matter equations of state in internal spaces.

In this case we get from (18) for small $`g`$

$$gq\frac{2N+1}{N+2}\mathrm{\Omega }^1+4\frac{N1}{N+2}\mathrm{\Omega }^2,$$
(24)
and for $`k=0`$ and small $`g`$ we obtain from (19)

$$1g\mathrm{\Omega }^1+\mathrm{\Omega }^2.$$
(25)
Now we illustrate the formulas by the following example when $`N=6`$ ($`K^6`$ may be a Calabi-Yau manifold) and

$$q=\mathrm{\Omega }^1=\mathrm{\Omega }^2=0.5.$$
(26)
We get from (24)

$$g\frac{1}{16}0.06$$
(27)
in agreement with (20).

In this case the second fluid component corresponds to magnetic (Euclidean) $`NS5`$-brane (in $`D=10`$ type I, Het or II A string models).


text:
# On the stability of vortex-plane solitons: The solution of the problem of Josephson-vortex structure in layered superconductors and stacked junctions

## I Introduction

In this paper, we present the solution of the problem of equilibrium vortex structure in layered superconductors and stacked Josephson junctions in the presence of a parallel, static, homogeneous external magnetic field $`H`$ and provide a unified interpretation of the available experimental data. Our approach consists in a rigorous mathematical analysis of the stability of all types of flux configurations, proposed in the literature, by means of exact variational methods for microscopic free-energy functionals.K99 ; K01 ; K1

Within the framework of these methods, we have previously obtained a complete classification of all possible static soliton solutions to coupled sine-Gordon (SG) equations for phase differences both in infinite \[$`N=\mathrm{}`$, $`N`$ is the number of superconducting(S) layers\] layered superconductorsK99 ; K01  and finite ($`N<\mathrm{}`$) Josephson-junction stacksK1  for $`H>0`$. Based on the fundamental argument of soliton physicsBP75 ; Bo78 ; R82 ; DEGM82 ; GM86 ; S89  that topological solitons in nonlinear field theories are minimizers of the energy functionals (free-energy functionals in our case),r2  we have identified these solutions with equilibrium Josephson-vortex configurations. Their magnetic field has symmetry typical of plane defects, hence the term ”vortex planes”. Physically, a vortex plane can be regarded as a bound state of interlayer vortices (one vortex per each insulating layer in the plane). In contrast to a deep-rooted beliefB73  in an ”analogy” with Abrikosov vortices in continuum type-II superconductors, the SG equations for $`H>0`$ do not admit static soliton solutions that can be identified with an ”isolated Josephson vortex” or a ”triangular Josephson-vortex lattice”.

Unfortunately, a wide-spread misunderstanding of the fact that equilibrium Josephson vortices are nothing but static soliton solutions to the SG equations incurred a misunderstanding of the exact mathematical results of Refs. K99, ; K01, ; K1, . In the critical comment, Ref. Kr02, , V. M. Krasnov and L. N. Bulaevskii ”disprove” the conclusions of Refs. K99, ; K01,  by claiming (in contradiction to the fundamentals of soliton physics) that vortex-plane solitons ”maximize the free energy”.r1  According to Ref. Kr02, , the ”instability” of vortex planes is ”similar to the instability of the laminar solutiondG  for type-II superconductors”. \[Alternating superconducting and normal layers, envisaged by the laminar model,dG  have nothing to do with soliton physics and do not possess the property of topological stability: see a proof in subsection V.B of the present paper.\] V. M. Krasnov and L. N. Bulaevskii insist on hypothetical ”isolated fluxons”, allegedly, having ”lower energy” for $`H>0`$ and characterized by a ”much smaller length scale”. No exact mathematical definition of the ”isolated fluxons” is given. We have not found any definition in the original papers by the critics of vortex planes, either: For example, it is claimed in Refs. BCG91, ; Kr01,  that Josephson vortices ”do not exist” in single Josephson junctions with $`W2\lambda _J`$ ($`W`$ is the junction width, $`\lambda _J`$ is the Josephson length). However, an exact, closed-form analytical solutionK1 ; K2  to the single static SG equation, appropriate for this case, clearly demonstrates the existence of phase-difference solitons for arbitrarily small $`W`$, provided the external field $`H`$ is sufficiently high. \[As shown in Refs. K99, ; K01, , exactly these solitons account for the well-known Fraunhofer pattern of the critical Josephson current $`I_c\left(H\right)`$ for $`W2\lambda _J`$.\]

Furthermore, the manuscript Ref. K1, , although submitted twice to Physical Review B (November 2000, August 2001), is still not published. In particular, one of the referees disputed the conclusions of Ref. K1, , because, in his opinion, the soliton boundary conditionsBo78 ; R82 ; DEGM82 ; GM86 ; S89  employed therein ”overdetermined” the problem of the classification of equilibrium Josephson-vortex configurations. He argued that certain numerical simulations for the SG equations had demonstrated, aside from vortex-plane solitons, the existence of ”single-vortex” solutions. According to the referee, these solutions had ”lower free energy” than the vortex planes for given $`H`$. As in Ref. Kr02, , no exact definition of such solutions was given.

It should be emphasized that the idea of an ”analogy” between the Josephson-vortex structure in layered superconductors and the Abrikosov-vortex structure in continuum type-II superconductors was not supported in Ref. B73,  and subsequent publicationsCC90 ; BC91 ; BLK92 ; Ko93  by any serious mathematical arguments.F98 ; K99 ; K01  Neither was it confirmed by direct experimental observations of the equilibrium Josephson-vortex structure in artificial stacked junctions at $`H>0`$.N96 ; Yu98  Unfortunately, most theoretical efforts were constrained by the idea of an ”analogy”, hence the use of mathematically ill-formulated methods, such as, e.g., a ”continuum-limit approximation”.B73 ; CC90 ; BC91 ; BLK92 ; Ko93  For instance, the exact SG equations for the phase differences were not even derived in Refs. B73, ; CC90, , concerned with a ”single Josephson vortex” at $`H>0`$. (As shown by Farid,F98  equations of Ref. CC90,  have no physical solution.) Since the problem of the stability of the proposed ”vortex configurations” (i.e., whether they are actual points of minima of the free-energy functionals) required the use of rigorous mathematical methods, it was not even posed in Refs. B73, ; CC90, ; BC91, ; BLK92, ; Ko93, .

Concerning numerical simulations for the static SG equations,Kr00 ; Kr01  there is an unjustified tendency to identify any kink-type feature of the phase difference with a ”Josephson vortex”, without any analysis of its stability. In contrast to the exact analytical methods of Refs. K99, ; K01, ; K1, , the numerical approach does not provide any means to establish a full set of necessary and sufficient conditions of the minimum of the free-energy functionals. Typically,Kr00  numerical simulations start with an incorrectly formulated (both mathematically and physically) boundary value problem that does not meet the criterion of uniqueness.CH

To close the issue of the equilibrium Josephson-vortex structure in layered superconductors and stacked junctions, we determine analytically (by means of exact methods of the calculus of variationsL62  and soliton physicsZK74 ; BP75 ; R82 ; DEGM82 ; GM86 ; S89 ) the type of all stationary points of the exact microscopic Gibbs free-energy functional,K1  generating the static SG equations. Our consideration applies to an arbitrary number of superconducting layers $`N`$, including the cases $`N=2`$ (a single junction) and $`N\mathrm{}`$ (an infinite layered superconductor). As a result, we obtain a complete classification of all nontrivial solutions to the SG equations, considered in the literature (both analytically and numerically), with respect to their stability. As could be expected from the general arguments of soliton physics, the only minimizers of the free-energy functional are the Meissner solution (the ”vacuum” state) and soliton vortex-plane solutions.K99 ; K01 ; K1  The latter represent the actual equilibrium Josephson-vortex configurations for $`H>0`$. In contrast, non-soliton configurations (e.g., ”isolated fluxons”, ”triangular Josephson-vortex lattices”, etc.) are absolutely unstable: They are nothing but saddle points of the Gibbs free-energy functional and are not even stationary points of the Helmholtz free-energy functional (obtained from the former by a Legendre transformation). \[Physically, non-soliton solutions violate conservation laws for the current and the flux, as was first noticed in the case of an infinite layered superconductor ($`N=\mathrm{}`$) in Refs. K99, ; K01, .\]

In section II, we clarify a relationship between the correct formulation of the boundary value problem to the SG equations and a full set of necessary and sufficient conditions of the minimum of the Gibbs free-energy functional. The proof of the stability of the Meissner solution and vortex-plane solitons is given in sections III, IV. In section III, we establish that the sufficient condition of the minimum of the Gibbs free-energy functional consists in the vanishing of the surface variation of the corresponding Helmholtz free-energy functional, which, in turn, yields conservation laws for the flux and the intralayer current. In subsection IV.A, we derive the soliton boundary conditions directly from the conservation law for the flux, which provides the sought proof of the stability of the Meissner solution and vortex-plane solitons. The main physical and mathematical properties of these solutions are discussed in subsection IV.B. In subsection IV.C, we analyze the obtained results from the point of view of general theory of topological defectsBP75 ; Bo78 ; M79 ; R82 ; DEGM82 ; GM86 ; S89  and explain dynamic stability of vortex planes, established in numerical simulations.SBP93 ; Kl94 ; Ma00

Some important physical and mathematical issues, related to the main results of the paper, are discussed in section V. In subsection V.A, we present a rigorous analytical description of unstable solutions (”isolated fluxons”, ”triangular Josephson-vortex lattices”, etc.), proposed in previous theoretical publications and numerical simulations. In subsection V.B, we draw a comparison with the Abrikosov vortices in type-II superconductors. In subsection V.C, we analyze the available experimental data from the point of view of the stable vortex-plane configurations.

In section VI, we summarize the obtained results, systematize our criticism of previous approaches and make some concluding remarks. In Appendix A, we give a list of mathematical formulas, relevant to the main text. In Appendix B, we establish a relationship to the variational principle of Refs. K99, ; K01,  for infinite layered superconductors.

Throughout the paper, we adhere to the dimensionless notation of Ref. K1, . The geometry of the problem is that of figures 1, 2 in Ref. K1, : The superconductor occupies the region $`\left[0xN1\right]\times \left[LyL\right]\times \left(\mathrm{}<z<\mathrm{}\right)`$, where $`2L=W`$; the layering axis (the $`c`$-axis) is $`x`$; the axis $`y`$ is along the layers; the external magnetic field is along the axis $`z`$: $`𝐇=\left(0,0,H0\right)`$. The phase difference between two successive S-layers is denoted as $`\varphi _n\phi _n\phi _{n1}`$ ($`\varphi _0=\varphi _N0`$), with $`n=0,1,\mathrm{},N1`$ being the S-layer number. The $`c`$-axis external current is not considered: $`I=0`$.

## II The formulation of the problem

In Ref. K1, , we have derived the coupled static SG equations for the phase differences $`\varphi _n`$ by minimizing the exact microscopic Gibbs free-energy functionalK01  $`\mathrm{\Omega }[f_n,\varphi _n,𝐀;H]`$ with respect to $`f_n`$ and $`𝐀`$($`f_n`$ is the reduced modulus of the order parameter in the $`n`$-th S-layer, $`𝐀`$ is the vector potential). In the limit $`r\left(T\right)1`$, $`HH_{c2}`$ \[$`r\left(T\right)`$ is the parameter of the interlayer coupling, $`H_{c2}`$ is the upper critical field\], when $`f_n=1`$, the SG equations appear as solubility conditions for the Maxwell equations in the gauge

$$𝐀=[0,A(x,y),0].$$
(1)
The SG equations read:K1

$$\frac{d^2\varphi _n\left(y\right)}{dy^2}=\frac{1}{ϵ^2}\underset{m=1}{\overset{N1}{}}G^1(n,m)\mathrm{sin}\varphi _m\left(y\right),n=1,\mathrm{},N1;$$
(2)
where $`G^1(n,m)`$ is a Jacobian matrix with elements $`G^1(n,n)=2+ϵ^2`$ ($`n=1,\mathrm{},N1`$), $`G^1(n+1,n)=G^1(n,n+1)=1`$ ($`n=1,\mathrm{},N2`$), and $`G^1(n,m)=0`$ for $`\left|nm\right|>1`$.

The requirement that the local field $`H_n\left(y\right)`$ ($`n1x<n`$, $`n=1,\mathrm{},N1`$) be equal to the applied one at $`y=\pm L`$ has led to the conditions

$$\frac{d\varphi _n}{dy}\left(L\right)=\frac{d\varphi _n}{dy}\left(L\right),n=1,\mathrm{},N1;$$
(3)
$$\frac{d\varphi _n}{dy}\left(L\right)=\frac{d\varphi _{n+1}}{dy}\left(L\right)\frac{d\varphi }{dy}\left(L\right)0,n=1,\mathrm{},N2.$$
(4)
\[The condition $`\frac{d\varphi }{dy}\left(L\right)0`$ merely reflects the fact that the local field is parallel to the applied one $`H0`$. \] From the requirement that the local field be equal to the applied one at $`x=0`$, $`x=N1`$, we have obtained

$$\varphi _n\left(y\right)=\varphi _{Nn}\left(y\right),n=1,\mathrm{},N1.$$
(5)
Equations (2) and boundary conditions (3) are satisfied by functions of the type

$$\varphi _n\left(y\right)=\varphi _n\left(y\right)+2\pi Z_n,$$
(6)
where the constants $`Z_n`$ can be arbitrarily chosen from the set $`0,\pm 1,\pm 2,\mathrm{}`$ As is pointed out in Ref. K1, , the fixation of the constants $`Z_n`$ requires imposition of boundary conditions on $`\varphi _n`$ at $`y=\pm L`$. Based on general arguments of soliton physics that soliton solutions are minimizers of corresponding energy functionals, we have imposed in Ref. K1,  standard soliton boundary conditions on $`\varphi _n`$.

As it has turned out, the only possible solutions, compatible with the requirement (4), are the Meissner solution and the soliton vortex-plane solutions, for which (5) is satisfied automatically and

$$Z_n=Z_{n+1}N_v,n=1,\mathrm{},N2;N_v=0,1,2,\mathrm{}$$
(7)
in (6), with $`N_v=0`$ representing the topologically trivial Meissner solution. Given thatK1

$$H_n\left(y\right)=H\left[G(n,1)+G(n,N1)\right]+\frac{ϵ^2}{2}\underset{m=1}{\overset{N1}{}}G(n,m)\frac{d\varphi _m\left(y\right)}{dy},$$
(8)
the field $`H=\overline{H}`$, corresponding to a concrete configuration $`\left\{\overline{\varphi }_n\right\}`$ with $`N_v=\overline{N}_v`$, is determined by

$$\frac{d\overline{\varphi }}{dy}\left(L\right)=2\overline{H}.$$
(9)
Note that the first, phase-independent, term in (8) is a contribution of the field penetrating through the boundaries $`x=0`$, $`x=N1`$, and the second term is a contribution of the field penetrating through the boundaries $`y=\pm L`$. The matrix $`G(n,m)`$is the inverse of $`G^1(n,m)`$: for its properties, see Appendix A.

In contrast to the above exact variational method, numerical simulationsKr00 ; Kr01  for (2) start with the imposition of the boundary conditions

$$\frac{d\varphi _n}{dy}\left(\pm L\right)=2H,$$
(10)
without any regard to sufficient conditions of the minimum of the Gibbs free-energy functional. Such an approach is based on an erroneous belief that all solutions to (2), (10) minimize the free-energy functional. (Manifestations of this belief are the naive ”energy arguments” of Ref. Kr02, , appealing to ”a difference in the length scales”, and calculations by means of combinatoricsKr00 ; Kr01  of the ”number of quasi-equilibrium fluxon modes”.) However, conditions (10) do not specify any boundary value problem for (2): By virtue of the symmetry relations (6), the imposition of the boundary condition on $`\frac{d\varphi _n}{dy}`$ at $`y=L`$ automatically ensures the fulfillment of the same boundary condition at $`y=L`$, whereas the constants $`Z_n`$ remain undetermined. Thus, the ”boundary value problem” (10) does not satisfy the criterion of uniqueness,CH  which is a sign of the presence of unphysical (i.e., unobservable) solutions.

The existence of redundant solutions to (2), (10) is already clear for physical reasons: This ”boundary value problem” does not take any account of the necessity to ensure the continuity of the local field at the boundaries $`x=0`$, $`x=N1`$. To understand at a rigorous mathematical level where the unphysical solutions come from, we have to consider all stationary points of the generating Gibbs free-energy functional, rewritten via $`\varphi _n`$ and $`\frac{d\varphi _n}{dy}`$.K1  In this way, we will derive a full set of necessary and sufficient conditions of the minimum directly from the variational principle, obtain an independent proof of the fact that the Meissner solution and the soliton vortex-plane solutionsK99 ; K01 ; K1  are the unique minimizers of the problem and establish the character of the instability of unphysical non-soliton solutions. In our consideration, we will employ the first integral of (2) that, taking account of (10), has the formK1

$$\underset{n=1}{\overset{N1}{}}\mathrm{cos}\varphi _n\left(y\right)+\frac{ϵ^2}{2}\underset{n=1}{\overset{N1}{}}\underset{m=1}{\overset{N1}{}}G(n,m)\frac{d\varphi _n\left(y\right)}{dy}\frac{d\varphi _m\left(y\right)}{dy}$$

$$=\frac{2H^2}{H_s^2}\left(N1\right)+\underset{n=1}{\overset{N1}{}}\mathrm{cos}\varphi _n\left(L\right),$$
(11)
where $`H_s`$ is the superheating (penetration) field of a semiinfinite ($`0y<+\mathrm{}`$) Josephson-junction stack, given by Eq. (66).

## III The necessary and sufficient conditions of the minimum of the Gibbs free-energy functional

The generating Gibbs free-energy functional for the SG equations (2) has the formK1

$$\mathrm{\Omega }[\varphi _n,\frac{d\varphi _n}{dy};H]=F[\varphi _n,\frac{d\varphi _n}{dy};H]4Hr\left(T\right)\underset{n=1}{\overset{N1}{}}\mathrm{\Phi }_n\frac{\varphi _n\left(L\right)\varphi _n\left(L\right)}{2\pi },$$
(12)
$$F[\varphi _n,\frac{d\varphi _n}{dy};H]$$

$$=r\left(T\right)[\frac{2H^2}{H_s^2}W(N1)+\frac{ϵ^2}{2}\underset{n=1}{\overset{N1}{}}\underset{m=1}{\overset{N1}{}}G(n,m)\underset{L}{\overset{L}{}}dy\frac{d\varphi _n\left(y\right)}{dy}\frac{d\varphi _m\left(y\right)}{dy}$$

$$+\underset{n=1}{\overset{N1}{}}\underset{L}{\overset{L}{}}dy[1\mathrm{cos}\varphi _n\left(y\right)]],$$
(13)
$$\mathrm{\Phi }_n=\pi \left[1G(n,1)G(n,N1)\right].$$
(14)
Note that the functional (12) is measured from the condensation energy $`\mathrm{\Omega }_0=\frac{NW}{2}`$. Moreover, the energy of the external field in the absence of the sample is subtracted.LL

Our task is to establish the necessary and sufficient conditions of the minimum of (12):

$$\mathrm{\Delta }\mathrm{\Omega }[\varphi _n,\frac{d\varphi _n}{dy};H]$$

$$=\mathrm{\Omega }[\varphi _n+\delta \varphi _n,\frac{d\varphi _n}{dy}+\frac{d\delta \varphi _n}{dy};H]\mathrm{\Omega }[\varphi _n,\frac{d\varphi _n}{dy};H]0$$
(15)
First, we observe that, in contrast to (12), the functional (13) is positive definite

$$F[\varphi _n,\frac{d\varphi _n}{dy};H]0,$$
(16)
since the matrix $`G(n,m)`$ is positive definite: see Appendix A. The absolute minimum of (13) is achieved for $`H=0`$, $`\varphi _n0`$ ($`n=1,2,\mathrm{},N1`$). Hence, the functional (13) necessarily has minima for any $`H0`$. The first variations of (12), (13) are

$$\delta \mathrm{\Omega }[\varphi _n,\frac{d\varphi _n}{dy};H]=\delta F[\varphi _n,\frac{d\varphi _n}{dy};H]4Hr\left(T\right)\underset{n=1}{\overset{N1}{}}\mathrm{\Phi }_n\frac{\delta \varphi _n\left(L\right)\delta \varphi _n\left(L\right)}{2\pi },$$
(17)
$$\delta F[\varphi _n,\frac{d\varphi _n}{dy};H]$$

$$=r\left(T\right)\underset{n=1}{\overset{N1}{}}\underset{L}{\overset{L}{}}𝑑y\left[\mathrm{sin}\varphi _n\left(y\right)ϵ^2\underset{m=1}{\overset{N1}{}}G(n,m)\frac{d^2\varphi _m\left(y\right)}{dy^2}\right]\delta \varphi _n\left(y\right)$$

$$+2\frac{d\varphi }{dy}\left(L\right)r\left(T\right)\underset{n=1}{\overset{N1}{}}\mathrm{\Phi }_n\frac{\delta \varphi _n\left(L\right)\delta \varphi _n\left(L\right)}{2\pi }.$$
(18)
The requirement that the volume variation in (18) vanish yields the SG equations (2), as expected. Of special interest to us are surface variations, i.e., the last terms in (17), (18): The requirement that these variations vanish determines boundary conditions to (1). (For a very clear discussion of the relationship between the surface variation and boundary conditions, see Ref. L62, , section II.15.) In the derivation of the surface variation in (18), we have used conditions (3), (4) that the local field be continuous at the boundaries $`y=\pm L`$. The requirement of the continuity of the local field at the boundaries $`x=0`$, $`x=N1`$ has not been so far employed. \[Recall our remark in section II that the disregard of this requirement is the reason for unphysical solutions to (1), (10).\]

If we simply enforce the conditions (10), the surface variations in (17) and (18) cancel out: Thus, all solutions to (2), (10) are stationary points of the Gibbs free-energy functional (12). However, under (10), the surface variation in (18) does not vanish. Therefore, not all solutions to (2), (10) are stationary points of (13). We have to examine conditions of the stationarity of (13) in more detail.

The requirement that the local field is fixed at the boundaries $`x=0`$, $`x=N1`$ is equivalent to the requirement that the vector potential $`𝐀`$ is fixed at $`x=0`$, $`x=N1`$. Consider now the total flux $`\mathrm{\Phi }`$. In the gauge (1), we have

$$\mathrm{\Phi }=\underset{L}{\overset{L}{}}𝑑y\left[A(N1,y)A(0,y)\right].$$
(19)
On the other hand,

$$\mathrm{\Phi }=\underset{n=1}{\overset{N1}{}}\underset{L}{\overset{L}{}}𝑑yH_n\left(y\right)=HW\left(N1\right)\frac{H_s^21}{H_s^2}+\underset{n=1}{\overset{N1}{}}\mathrm{\Phi }_n\frac{\varphi _n\left(L\right)\varphi _n\left(L\right)}{2\pi },$$
(20)
where the first term is the flux penetrating through the boundaries $`x=0`$, $`x=N1`$, and the second term is the flux penetrating through the boundaries $`y=\pm L`$. (We will call it the ”Josephson flux”, $`\mathrm{\Phi }_J`$.) Given that

$$\delta A(0,y)=\delta A(N1,y)=0,$$
(21)
we have

$$\delta \mathrm{\Phi }=\delta \mathrm{\Phi }_J=\underset{n=1}{\overset{N1}{}}\mathrm{\Phi }_n\frac{1}{2\pi }\underset{L}{\overset{L}{}}𝑑y\frac{d\delta \varphi _n\left(y\right)}{dy}=\underset{n=1}{\overset{N1}{}}\mathrm{\Phi }_n\frac{\delta \varphi _n\left(L\right)\delta \varphi _n\left(L\right)}{2\pi }=0.$$
(22)
Thus, the continuity of the field at $`x=0`$, $`x=N1`$ imposes a constraint on the variations:

$$\mathrm{\Phi }_J=\underset{n=1}{\overset{N1}{}}\mathrm{\Phi }_n\frac{1}{2\pi }\underset{L}{\overset{L}{}}𝑑y\frac{d\varphi _n\left(y\right)}{dy}=const0.$$
(23)
The result (23) is exactly what had to be expected:LL  By virtue of the Meissner effect, the flux $`\mathrm{\Phi }_J`$ (and, of course, $`\mathrm{\Phi }`$) is stable against any small perturbations, represented by variations $`\varphi _n\left(y\right)\varphi _n\left(y\right)+\delta \varphi _n\left(y\right)`$. Equivalent forms of (22), (23) are

$$\delta \varphi _n\left(L\right)=\delta \varphi _n\left(L\right),$$
(24)
$$\varphi _n\left(L\right)\varphi _n\left(L\right)=c_n=const0.$$
(25)
Note that the existence of conserved physical quantities of the type of $`\mathrm{\Phi }_J`$ is a precursor to the existence of soliton solutions in nonlinear field theories.ZK74 ; BP75 ; Bo78 ; R82 ; DEGM82 ; GM86 ; S89  In Appendix B, we establish a relationship between the conservation of $`\mathrm{\Phi }_J`$ and the conservation of the intralayer current, which, in turn, establishes a relationship to the variational principle for infinite ($`N=\mathrm{}`$) layered superconductors.K99 ; K01

What will be shown now is that all the stationary points of (13) are the unique minimizers of both (12) and (13). First, we notice that high-order variations of (12) and (13) coincide: $`\delta ^nF=\delta ^n\mathrm{\Omega }`$, $`n2`$. Thus, all the minimizers of (13) are minimizers of (12). On the other hand, the minimizers of (12) obeying (23) are minimizers of (13): From the condition of the minimum (15), we get

$$F[\varphi _n+\delta \varphi _n,\frac{d\varphi _n}{dy}+\frac{d\delta \varphi _n}{dy};H]F[\varphi _n,\frac{d\varphi _n}{dy};H]\delta \mathrm{\Phi }_J$$

$$=F[\varphi _n+\delta \varphi _n,\frac{d\varphi _n}{dy}+\frac{d\delta \varphi _n}{dy};H]F[\varphi _n,\frac{d\varphi _n}{dy};H]0.$$
(26)
\[Physically, this fact means the equivalence of the description in terms of the Gibbs free energy and the Helmholtz free energy: Because of (23), the functional $`F[\varphi _n,\frac{d\varphi _n}{dy};\mathrm{\Phi }_J]F[\varphi _n,\frac{d\varphi _n}{dy};0]`$ can be regarded as the Helmholtz free-energy functional.\] Using the standard technique,ZK74  it is straightforward to prove that all stationary points of (12), obeying (23), are minimizers of (12) \[and, thus, of (13)\]. Indeed, let $`\left\{\overline{\varphi }_n\right\}`$ be the stationary point for $`\mathrm{\Phi }_J=\overline{\mathrm{\Phi }}_J`$ and corresponding $`H=\overline{H}`$. In the vicinity of $`\left\{\overline{\varphi }_n\right\}`$, i.e., for $`\varphi _n=\overline{\varphi }_n+\delta \overline{\varphi }_n`$, we have the following estimate:

$$\mathrm{\Omega }[\varphi _n,\frac{d\varphi _n}{dy};\overline{H}]$$

$$r\left(T\right)\left[\frac{ϵ^2}{2}\underset{n=1}{\overset{N1}{}}\underset{m=1}{\overset{N1}{}}G(n,m)\underset{L}{\overset{L}{}}𝑑y\frac{d\varphi _n\left(y\right)}{dy}\frac{d\varphi _m\left(y\right)}{dy}4\overline{H}\underset{n=1}{\overset{N1}{}}\mathrm{\Phi }_n\frac{1}{2\pi }\underset{L}{\overset{L}{}}𝑑y\frac{d\varphi _n\left(y\right)}{dy}\right]$$

$$4\overline{H}r\left(T\right)\overline{\mathrm{\Phi }}_J.$$
(27)
Inequalities (27) show that $`\mathrm{\Omega }`$ has a lower bound in the vicinity of any stationary point $`\left\{\overline{\varphi }_n\right\}`$, obeying (23); hence, $`\left\{\overline{\varphi }_n\right\}`$ is a minimizer of $`\mathrm{\Omega }`$ and $`F`$ , Q.E.D.

To strengthen (27), we minimize the right-hand side of the first inequality with respect to $`\frac{d\varphi _n}{dy}`$, obtaining

$$ϵ^2\underset{m=1}{\overset{N1}{}}G(n,m)\left[\frac{d\varphi _m}{dy}\right]_{\mathrm{min}}=\frac{2\overline{H}}{\pi }\mathrm{\Phi }_n,$$
(28)
$$\mathrm{\Omega }[\varphi _n,\frac{d\varphi _n}{dy};\overline{H}]2\overline{H}r\left(T\right)\overline{\mathrm{\Phi }}_J.$$
(29)
Taking into account that $`\left\{\overline{\varphi }_n\right\}`$ is a solution of (2), making use of (11) and (28), we get:

$$\mathrm{\Omega }[\varphi _n,\frac{d\varphi _n}{dy};\overline{H}]\mathrm{\Omega }[\overline{\varphi }_n,\frac{d\overline{\varphi }_n}{dy};\overline{H}]$$

$$=r\left(T\right)\left[W\underset{n=1}{\overset{N1}{}}\left[1\mathrm{cos}\overline{\varphi }_n\left(L\right)\right]+ϵ^2\underset{n=1}{\overset{N1}{}}\underset{m=1}{\overset{N1}{}}G(n,m)\underset{L}{\overset{L}{}}𝑑y\frac{d\overline{\varphi }_n\left(y\right)}{dy}\frac{d\overline{\varphi }_m\left(y\right)}{dy}4\overline{H}\overline{\mathrm{\Phi }}_J\right]$$

$$r\left(T\right)W\underset{n=1}{\overset{N1}{}}\left[1\mathrm{cos}\overline{\varphi }_n\left(L\right)\right]0.$$
(30)
The inequality $`\mathrm{\Omega }(\overline{\varphi }_n,\frac{d\overline{\varphi }_n}{dy};\overline{H})0`$ is a manifestation of the Meissner effect and had to be expected from general thermodynamic arguments.LL  Relations (30) immediately yield

$$F[\varphi _n,\frac{d\varphi _n}{dy};\overline{H}]F[\overline{\varphi }_n,\frac{d\overline{\varphi }_n}{dy};\overline{H}]$$

$$r\left(T\right)W\underset{n=1}{\overset{N1}{}}\left[1\mathrm{cos}\overline{\varphi }_n\left(L\right)\right]+4\overline{H}r\left(T\right)\overline{\mathrm{\Phi }}_J4\overline{H}r\left(T\right)\overline{\mathrm{\Phi }}_J.$$
(31)
We want to emphasize that inequalities of the type (29), (31) are typical of soliton physics: They are used to establish the existence and stability of soliton solutions.ZK74 ; BP75 ; R82 ; DEGM82 ; GM86 ; S89

Note an alternative interpretation of the variational principle

$$\delta \mathrm{\Omega }[\varphi _n,\frac{d\varphi _n}{dy};H]=\delta F[\varphi _n,\frac{d\varphi _n}{dy};0]4Hr\left(T\right)\underset{n=1}{\overset{N1}{}}\mathrm{\Phi }_n\frac{1}{2\pi }\delta \underset{L}{\overset{L}{}}𝑑y\frac{d\varphi _n\left(y\right)}{dy}=0.$$
(32)
The field $`H`$ in (32) can be considered as a Lagrange multiplier, implying that variation can be performed without any restrictions on $`\delta \varphi _n\left(L\right)`$, $`\delta \varphi _n\left(L\right)`$. In this case, the requirement that the surface variation vanish yields conditions (3), (4). Boundary conditions on $`\varphi _n`$ are uniquely determined by (23): see the next section. The value $`H=\overline{H}`$ for a concrete minimizer $`\{\overline{\varphi }_n,\overline{\mathrm{\Phi }}_J\}`$ should be found from the condition of thermodynamic equilibrium

$$\frac{\mathrm{\Omega }(\overline{\varphi }_n,\frac{d\overline{\varphi }_n}{dy};\overline{H})}{\overline{\mathrm{\Phi }}_J}=0.$$
(33)
Indeed, $`\mathrm{\Omega }(\overline{\varphi }_n,\frac{d\overline{\varphi }_n}{dy};\overline{H})`$ can be written as

$$\mathrm{\Omega }(\overline{\varphi }_n,\frac{d\overline{\varphi }_n}{dy};\overline{H})=r\left(T\right)[W\underset{n=1}{\overset{N1}{}}[1\mathrm{cos}\overline{\varphi }_n\left(L\right)\frac{1}{W}\underset{L}{\overset{L}{}}dy\overline{\varphi }_n\left(y\right)\mathrm{sin}\overline{\varphi }_n\left(y\right)]$$

$$+2[\frac{d\overline{\varphi }}{dy}\left(L\right)2\overline{H}]\overline{\mathrm{\Phi }}_J],$$
(34)
which by virtue of (33) immediately yields (9).

In summary, we have proved the following: The SG equations (1) and relations (10) ensure only the stationarity of the Gibbs free-energy functional (12). The necessary and sufficient conditions of the minimum of both (12) and (13) (which is the Helmholtz free-energy functional for $`H=0`$) are given by (3), (4) and the constraint (23). Solutions to (2), (10) that do not obey this constraint are absolutely unstable. The character of this instability can be easily established. Indeed, such solutions are not even stationary points of the Helmholtz free-energy functional, therefore $`\delta ^2\mathrm{\Omega }=\delta ^2F`$ need not have definite sign. Moreover, the functional (12) is unbounded in the vicinity of these solutions. Thus, they are nothing but saddle points of (12).

## IV The proof of the stability of the Meissner solution and vortex-plane solitons

### IV.1 Boundary conditions on $`\varphi _n`$

Our aim now is to establish boundary conditions on $`\varphi _n`$ directly from the constraint (23). Given that the SG equations (2), boundary conditions (3), (4) and the constraint (23) represent a full set of necessary and sufficient conditions of the minimum of the Gibbs and Helmholtz free-energy functionals, we will obtain, in this manner, the sought proof of the stability of the Meissner solution and soliton vortex-plane solutions.

First, we observe that since a minimizer of (12), (13) must ensure the vanishing of both the surface and volume variations in (17) and (18), it should necessarily belong to the class of functions that satisfy (2), (3), (4) and the symmetry relations (6). Thus, the variation of the surface terms in (17) and (18) is performed with respect to trial functions that take only discrete values at $`y=0`$:

$$\varphi _n\left(0\right)=\pi Z_n,$$
(35)
where $`Z_n`$ can be arbitrarily chosen from the set $`0,\pm 1,\pm 2,\mathrm{}`$ .These functions can be subdivided into classes parameterized by an $`\left(N1\right)`$-dimensional ”vector”

$$𝐐=(Z_1,Z_2,\mathrm{},Z_{N1}).$$
(36)
In view of (35), the requirement of the continuity of variations can be met if and only if

$$\delta \varphi _n\left(0\right)=0,$$
(37)
which means that all the minima of (12), (13) are parameterized by the vector $`𝐐`$, and the variation of the surface terms in (17) and (18) is performed with respect to trial functions that belong to a certain class (36). Moreover, the symmetry relations (6) imply $`\delta \varphi _n\left(L\right)=\delta \varphi _n\left(L\right)`$. Combined with (14), this yields

$$\delta \varphi _n\left(L\right)=\delta \varphi _n\left(L\right)=0.$$
(38)
Now we combine (15) with (6) to obtain

$$\varphi _n\left(L\right)\pi Z_n=\frac{c_n}{2}0,$$
(39)
$$\varphi _n\left(L\right)+\pi Z_n=\frac{c_n}{2}0.$$
(40)
Since the inequalities in (39), (40) should hold for any $`Z_n`$, including $`Z_n=0`$, we get

$$\varphi _n\left(L\right)0,\varphi _n\left(L\right)0.$$

Moreover, since for any fixed set $`\left\{c_n\right\}`$ the set $`\left\{\varphi _n\left(\pm L\right)\right\}`$ must belong to a certain unique class (36),

$$\pi <\varphi _n\left(L\right)0,$$
(41)
$$2\pi Z_nc_n<2\pi \left(Z_n+1\right),Z_n=\left[\frac{1}{2\pi }\underset{L}{\overset{L}{}}𝑑y\frac{d\varphi _n\left(y\right)}{dy}\right]=0,1,2,\mathrm{},$$
(42)
where $`\left[u\right]`$ is the integer part of $`u`$.

Given (37), (38), the boundary conditions (41) and

$$\varphi _n\left(0\right)=\pi Z_n,Z_n=0,1,2,\mathrm{}$$
(43)
together with (3), (4) determine, in principle, a complete set of conditions for the minimizer of (12), (13). The solution belonging to a certain class (36) \[with $`Z_n`$ as in (43)\] first appears when all $`\varphi _n\left(L\right)=0`$, and $`c_n=2\pi Z_n`$ \[see the estimates (30), (31): under these conditions, the density of the Josephson energy at the boundary is a minimum\]. Thus, we are confronted with the standardBo78 ; R82 ; DEGM82 ; GM86 ; S89  soliton boundary value problem

$$\varphi _n\left(L\right)=0,\varphi _n\left(0\right)=\pi Z_n,Z_n=0,1,2,\mathrm{},$$
(44)
plus the boundary conditions (4). To find out what type of minimizers can be realized in reality, we have to solve Eqs. (2).K1 ; K2  Note that $`N2`$ relations (4) implicitly impose $`N2`$ conditions on $`N1`$ constants $`Z_n`$. Therefore, as shown in Refs. K1, ; K2, , the only solutions to (2), (44), compatible with (4), are those that satisfy (7), i.e., the Meissner solution and the soliton vortex-plane solutions:

$$𝐐_v=(N_v,N_v,\mathrm{},N_v),N_v=\left[\frac{1}{2\pi }\underset{L}{\overset{L}{}}𝑑y\frac{d\varphi _n\left(y\right)}{dy}\right]=0,1,2,\mathrm{}$$
(45)
The properties of these solutions will be discussed in more detail in what follows. Here we want to brief on the results for $`H=0`$.K1 ; K2  For $`H=0`$, $`L<\mathrm{}`$, there are no soliton solutions at all (including the vortex planes), and the only stable solution is the trivial Meissner solution $`\varphi _1=\varphi _2=\mathrm{}\varphi _{N1}0`$. The situation changes drastically for $`H=0`$, $`L=\mathrm{}`$: The imposition of soliton boundary conditions on $`\varphi _n\left(\pm \mathrm{}\right)`$ automatically ensures the fulfillment of the boundary conditions $`\frac{d\varphi _n}{dy}\left(\pm \mathrm{}\right)=0`$ by virtue of Eqs. (2) themselves and some elementary theorems of mathematical analysis. Aside from the vortex-plane solution with $`N_v=1`$ in (45) (where $`L=\mathrm{}`$), we have a variety of soliton solutions (36) with $`Z_n`$ arbitrarily chosen from the set $`0,\pm 1`$. The fact that for $`H=0`$, $`L=\mathrm{}`$ each $`\varphi _n`$ can ”accommodate” no more than one vortex or anivortex is a generalization of the well-knownBo78 ; R82  result for the single static SG equation and can be easily proved by the use of the first integral (11) with the right-hand side equal to $`N1`$.

### IV.2 The Meissner solution and soliton vortex-plane solutions

The range of the existence of the solutions, parameterized by (45), is determined by the boundary value problem (41), (43) (with all $`Z_n=N_v`$) and relation (45):K1 ; K2

$$0H<H_0H_{sL},forN_v=0,$$
(46)
$$\sqrt{H_{N_v1}^2H_s^2}H<H_{N_v},forN_v=1,2,\mathrm{},$$
(47)
where $`H_{sL}`$ has the meaning of the superheating field of the Meissner state ($`N_v=0`$) for $`L<\mathrm{}`$ ($`H_{sL}>H_s`$ for $`L<\mathrm{}`$, and $`H_s\mathrm{}H_s`$). The lower bound in (46), (47) is determined by the exact upper bound in (41). At $`H=supH=H_{N_v}`$, when all $`\varphi _n\left(L\right)=inf\varphi _n\left(L\right)=\pi `$, there is instability of the saddle-point type (see the end of section III). Note that both the Meissner solution ($`N_v=0`$) and the vortex-plane solutions ($`N_v1`$) automatically satisfy the symmetry relations (5).K1 ; K2

It is instructive to verify the general inequalities (30), (31). Mathematically, it is sufficient to do this only for $`H`$ equal to the lower bounds in (46), (47): By continuity arguments, the result will be valid in the whole field range. For the Meissner solution the verification is trivial. For $`N_v1`$, we employ the exact expression

$$F[\overline{\varphi }_n,\frac{d\overline{\varphi }_n}{dy};\overline{H}]=r\left(T\right)ϵ^2\underset{n=1}{\overset{N1}{}}\underset{m=1}{\overset{N1}{}}G(n,m)\underset{L}{\overset{L}{}}𝑑y\frac{d\overline{\varphi }_n\left(y\right)}{dy}\frac{d\overline{\varphi }_m\left(y\right)}{dy},$$
(48)
where $`\overline{H}\sqrt{H_{N_v1}^2H_s^2}`$. As shown in Ref. K1, , in this case

$$\left[\frac{d\overline{\varphi }_n\left(y\right)}{dy}\right]_{\mathrm{min}}=\frac{d\overline{\varphi }_n}{dy}\left(\pm L\right)=2\sqrt{H_{N_v1}^2H_s^2},foralln.$$
(49)
Combining Eqs. (48), (49), we get exactly the lower bound in (31).

As is emphasized in Ref. K1, , the obtained solutions are valid for any $`N`$, including the cases $`N=2`$ (a single junction) and $`N=3`$ (a double-junction stack). (For $`N=2,3`$, we have derived in Ref. K1,  exact, closed-form analytical expressions.) The solutions with $`N_v1`$ are pure solitons only at $`H=\sqrt{H_{N_v1}^2H_s^2}`$, when $`j_n\left(\pm L\right)=\mathrm{sin}\varphi _n\left(\pm L\right)=0`$ \[$`j_n\left(y\right)`$ is the density of the Josephson currentK1  for $`f_n=1`$\]. In the rest of the regions (47), we have solitons ”dressed” by the Meissner field. In the case of $`N=2`$ (the single junction), Owen and ScalapinoOS67  called these regions the ”$`N_v`$ to $`N_v+1`$ vortex mode”. \[Because the principle of superposition does not apply to the nonlinear Eqs. (2), the Meissner and the vortex-plane fields cannot be separated from each other.\] It is clear that the vortex-plane solutions for $`N>2`$ are a direct generalization of ordinary vortices in single junctions.

Of special interest is the overlap of the regions (46), (47) for $`N_v=\overline{N}_v`$ and $`N_v=\overline{N}_v+1`$. As a result, the obtained solutions cover the whole field range $`0H<\mathrm{}`$, as they should. Mathematically, the overlap is related to the fact that the solution with $`N_v=\overline{N}_v`$ cannot be continuously transformed into the solution with $`N_v=\overline{N}_v+1`$ by changing $`H`$, as is always the case for solitons. For the single junction, the overlap was first established numerically in Ref. OS67,  and discussed qualitatively in Ref. K70, . The overlap practically vanishes for $`H_{N_v}H_s`$ . Given that all $`H_{N_v}`$ decrease when $`W=2L`$ increases,K1  the overlap is stronger for large $`W`$ and can involve several neighboring states. Physically, the actual equilibrium state is the one that corresponds to the absolute minimum of the Gibbs free energy for given $`H`$. The rest of the allowed states are metastable. In view of the above-mentioned discontinuity, a transition from the state with $`N_v=\overline{N}_v`$ to the state with lower Gibbs free energy $`N_v=\overline{N}_v+1`$ will necessarily be a phase transition of the first-order type.K70 ; K99 ; K01  It particular, the lower critical field $`H_{c1}`$ is determined from the requirement that the Gibbs free energy of the state $`N_v=1`$ be equal to that of the Meissner state ($`N_v=0`$) and satisfies the relation $`\sqrt{H_{sL}^2H_s^2}<H_{c1}<H_{sL}`$.

### IV.3 Topological considerations and stability in the dynamic regime

The stability of the Meissner solution and vortex-plane solitons can be better understood, if we analyze the obtained results from the general point of view of the stability of topological defects in continuum media.BP75 ; Bo78 ; M79 ; R82 ; DEGM82 ; GM86 ; S89  To this end, we consider the density of the Gibbs free energy (12) at the boundaries $`y=\pm L`$.

Because of the general symmetry relations (6), valid for any solution to (6), (10), the density of the Josephson energy is equal at $`y=L`$ and $`y=+L`$:

$$1\mathrm{cos}\varphi _n\left(L\right)=1\mathrm{cos}\varphi _n\left(+L\right),n1x<n,n=1,2,\mathrm{},N1.$$
(50)
Taking into account (11), we conclude that also the density of the total free energy is equal at the boundaries $`y=L`$ and $`y=+L`$ and thus corresponds to the degenerate equilibrium (”vacuum”) state, unperturbed by topological defects (solitons). Mathematically, the boundary of the interval $`LyL`$ can be considered as a $`0`$-dimensional sphere: $`S^0=\{L,+L\}`$. Given that configurations $`\varphi _n`$ and $`\varphi _n+2\pi Z_n`$ ($`Z_n=0,\pm 1,\pm 2,\mathrm{}`$) are physically indistinguishable, we can fix the values $`\varphi _n\left(L\right)`$ as in (41) and regard the functions

$$\psi _n\left(+L\right)\frac{\varphi _n\left(+L\right)+\varphi _n\left(L\right)}{2\pi }=Z_n$$
(51)
as continuous maps of the boundary into the additive group of the integers, $`𝐙`$: $`S^0\stackrel{\psi _n}{}𝐙`$. ($`𝐙`$ is the group of the degeneracy of the equilibrium state, or the order-parameter space.) The fact of the existence of topologically nontrivial maps of this type, realized by soliton solutions, is often expressed in terms of the ”zeroth homotopy group”BP75 ; Bo78 ; M79 ; R82 ; DEGM82 ; GM86 ; S89  $`\pi _0\left(M\right)`$, where the index ”$`0`$” stands for the boundary $`S^0`$ and $`M`$ is the order-parameter space:

$$\pi _0\left(𝐙\right)=𝐙.$$
(52)
\[Note that $`\pi _0\left(M\right)`$ is merely the set of disconnected components of the space $`M`$.\] Because of the boundary conditions (4), all $`\psi _n`$ at $`H>0`$ realize the same mapping : $`Z_1=Z_2=\mathrm{}=Z_{N1}Z`$. The external field $`H>0`$ breaks the symmetry $`\varphi _n\varphi _n`$ \[see the second term in (12)\]. Therefore, only the values $`ZN_v=0,1,2,\mathrm{}`$ are allowed, with $`N_v=0`$ being the ”vacuum”, Meissner state. In this way, we arrive at the natural topological classification (45) of the minimizers of (12). Owing to the continuity conditions (36), (37), variation in (17), (18) is allowed only with respect to trial functions $`\left\{\varphi _n\right\}`$ that have the same end points $`\varphi _n\left(\pm L\right)`$ and the middle point $`\varphi _n\left(0\right)`$ as the minimizer $`\left\{\overline{\varphi }_n\right\}`$, i.e., $`\left\{\varphi _n\right\}`$ are homotopic to $`\left\{\overline{\varphi }_n\right\}`$ and belong to the same class (45): hence the stability of $`\left\{\overline{\varphi }_n\right\}`$ against continuous perturbations.

Numerical simulations for time-dependent coupled SG equations have revealed exceptional stability of vortex planes in the dynamic regime as well: see figure 7 in Ref. SBP93, , figure 8 in Ref. Kl94, , and figure 3 in Ref. Ma00, . (The authors of Refs. SBP93, ; Kl94, ; Ma00,  employ the terms ”coherent”, ”in-phase” or ”phase-locked modes” instead of our term ”vortex planes” that we prefer for physical reasons.) Although a detailed analysis of the dynamics of vortex planes is beyond the scope of this paper and will be done elsewhere, the results of Refs. SBP93, ; Kl94, ; Ma00,  can be explained already at this stage.

In the absence of dissipation, the dynamic SG equations, describing an evolution of the system in the time interval $`t_itt_f`$, can be derived from a corresponding Lagrangian by use of a variational principle. The requirement that the surface variation vanish on the whole perimeter of the space-time boundary leads to the conditions $`\delta \varphi _n(y,t_i)=\delta \varphi _n(y,t_f)=0`$ and a generalization of the conservation law for the Josephson flux $`\mathrm{\Phi }_J`$, Eq. (23):

$$\mathrm{\Phi }_J=\underset{n=1}{\overset{N1}{}}\mathrm{\Phi }_n\frac{1}{2\pi }\underset{L}{\overset{L}{}}𝑑y\frac{\varphi _n(y,t)}{y}=\underset{n=1}{\overset{N1}{}}\mathrm{\Phi }_n\frac{\varphi _n(L,t)\varphi _n(L,t)}{2\pi }=const,$$
(53)
which means that the differences $`\varphi _n(L,t)\varphi _n(L,t)`$ do not depend on $`t`$. Thus, by fixing the boundary conditions $`\varphi _n(\pm L,t)`$ at $`t=t_i`$ as in subsection IV.A, we fix the initial value of the flux $`\mathrm{\Phi }_J=\overline{\mathrm{\Phi }}_J`$ that will not change in the course of the evolution of the system from $`t=t_i`$ to $`t=t_f`$. The topological type of the solution \[see (46)\] will not change, either:

$$N_v=\left[\frac{1}{2\pi }\underset{L}{\overset{L}{}}𝑑y\frac{\varphi _n(y,t)}{y}\right]=\left[\frac{\varphi _n(L,t)\varphi _n(L,t)}{2\pi }\right]=const.$$
(54)
As usual,Bo78 ; R82 ; DEGM82 ; GM86 ; S89  this situation can be formalized in terms of the conserved topological current

$$j_\mu =\underset{n=1}{\overset{N1}{}}\frac{\mathrm{\Phi }_n}{2\pi }ϵ_{\mu \tau }_\gamma \varphi _n,_\mu j_\mu =0,$$
(55)
where $`\mu ,\gamma =0,1`$; $`_\mu =(_t,_x)`$; and $`ϵ_{\mu \tau }`$ is the antisymmetric symbol on two indices, $`ϵ_{01}=ϵ_{10}=1`$; with $`\mathrm{\Phi }_J=\underset{L}{\overset{L}{}}𝑑yj_0`$ being the topological charge.

As should be clear from these results, time-dependent SG equations alone cannot describe Josephson-vortex penetration, i.e., an evolution of the system from the topologically trivial Meissner state, $`N_v=0`$, to a state with $`N_v0`$. Unfortunately, this important issue has not been realized in Ref. Kr02,  that claims to have ”demonstrated a dynamic process of vortex penetration” by means of numerical simulations for time-dependent SG equations.

## V Discussion

### V.1 Unstable solutions to the SG equations

As is proved in section III, all non-topological, non-soliton solutions to (2), (10) that do not meet the requirement of the flux conservation (23) are absolutely unstable: They are nothing but saddle points of the Gibbs free-energy functional (12), cannot be assigned any ”free energy” and are therefore unobservable. Since the requirement of the continuity of variations in (17) does not impose on such solutions any constraints of the type (36), (37), they can be continuously transformed into the stable Meissner solution or a vortex-plane solution, representing the actual minimum of (12) at a given $`H`$, by a series of infinitesimal deformations of $`\varphi _n`$ without a violation of the boundary conditions (10). \[A clear illustration of such a transformation for non-topological defects in a system of planar spins see in Ref. M79, , section II.B.\]

Analytically, all unstable configurations for $`H0`$ can be obtained using the symmetry relations (6), as solutions to the boundary value problem

$$\frac{d\varphi _n}{dy}\left(L\right)=2H,\varphi _n\left(0\right)=\pi Z_n,Z_n=0,1,2,\mathrm{}$$
(56)
that violate topological boundary conditions, derived in subsection IV.A. For example, one can set in (56) $`Z_1=Z_2=\mathrm{}=Z_{N1}Z`$ and increase $`H`$ beyond the upper bound $`H_{N_v}`$ of the stability regions (46), (47) for a given $`N_v=Z`$. By continuously increasing $`H`$ beyond the stability region of the Meissner state, unstable configurations with $`Z=0`$, interpreted as ”Josephson-vortex penetration”, were obtained in numerical simulations for static SG equations (Ref. Kl94, , figure 7) and time-dependent SG equations (Ref. Kr02, , figure 2). Analogous instability for $`Z1`$ is demonstrated by numerical simulations in the dynamic regime \[$`H=`$const $`>0`$, and an increasing transport current $`I>I_c\left(H\right)`$\] in Ref. Ma00, : see region IV in figure 3 therein; region III corresponds to dynamically stable vortex planes.

Unstable solutions appear also when not all $`Z_n`$ in (56) are equal to each other. Thus, an unstable ”single Josephson vortex” corresponds to the choice $`Z_l=1`$, $`Z_{nl}=0`$. Solutions of this type were obtained in several numerical simulations.SBP93 ; Kr00 ; Kr01  By way of illustration, we consider here only the case $`H=0`$. As is explained at the end of subsection IV.A, the only stable configuration for $`H=0`$, $`L<\mathrm{}`$ is the trivial Meissner state $`\varphi _1=\varphi _2=\mathrm{}=\varphi _{N1}0`$. Figure 5 in Ref. SBP93, , and figures 1, 2 in Ref. Kr01,  clearly show that the solutions presented therein, in reality, are characterized by all $`Z_n=\left[\frac{1}{2\pi }\underset{L}{\overset{L}{}}𝑑y\frac{d\varphi _n\left(y\right)}{dy}\right]=0`$ and, thus, belong to the class $`𝐐_v=(0,0,\mathrm{},0)`$ of the general topological classification (45). By means of continuous deformations, they can be transformed into the trivial Meissner solution.

Other unstable solutions, available in the literature, can be analyzed along the same lines. In particular, the ”triangular Josephson-vortex lattice with the period $`x_p=1`$”, proposed in Ref. BC91, , corresponds to the case $`Z_{odd}=Z`$, $`Z_{even}=Z+1`$.

### V.2 A comparison with Abrikosov vortices in type-II superconductors

As the formation of a vortex-plane soliton involves only phase differences between successive S-layers, it does not affect the topology of the layered superconductor. In contrast, the appearance of a linear ($`𝐑^1`$) singularity of the order parameter $`\mathrm{\Delta }\left(𝐫\right)=\left|\mathrm{\Delta }\left(𝐫\right)\right|\mathrm{exp}\left[i\phi \left(𝐫\right)\right]`$ is necessary for the formation of an Abrikosov vortex in continuum type-II superconductors. Thus, in the presence of a single Abrikosov vortex, the topology of the continuum type-II superconductor changes from $`𝐑^3`$ (the three-dimensional Euclidean space) to $`𝐑^3/𝐑^1=S^1`$. Therefore, the notion of the ”vortex core”dG  is inherent (both physically and mathematically) to the Abrikosov vortex and is meaningless in the case of the vortex plane.

An isolated Abrikosov vortex is itself a stable object, both topologically and energetically. (The latter can be proved by the use of the same mathematical methods as those employed in our section III: see, e.g., Refs. DEGM82, ; GM86, ; S89, .) Therefore, an equilibrium state of $`N_v`$ Abrikosov vortices is determined by comparing the values of the Ginzburg-Landau free-energy functional for different spatial configurations, which yields the well-known triangular lattice as the most favorable one.dG  In contrast, the notion of the ”Josephson-vortex lattice” is senseless for layered superconductors: One can only speak of $`N_v`$-soliton (vortex-plane) states, with $`N_v=0`$ representing the Meissner state, and each vortex plane being a ”Josephson vortex” itself.

In the case of extreme type-II superconductors, the linear singularities, associated with Abrikosov vortices, can be easily incorporated into the framework of the simple London model.dG  The resulting equation is a linear inhomogeneous partial differential equation for the local field. Owing to linearity, the local field is a superposition of the Meissner and vortex fields. As is emphasized in subsection IV.B, this is not the case for layered superconductors because of the nonlinearity of the SG equations (2). Unfortunately, this important issue was not understood in some publications concerned with Josephson-vortex penetration.BF92

As we can see, there is no ”analogy” between the Abrikosov-vortex structure in continuum type-II superconductors and the Josephson-vortex structure in layered superconductors in the naive sense.B73  Instead, there is a much subtler mathematical analogy: The topological classification of vortex configurations in type-II superconductors is isomorphic to that in layered superconductors. A proof is straightforward. For the reasons explained above, the boundary of a type-II superconductor is, in general, topologically equivalent to a one-dimensional sphere (a circle) $`S^1`$. The order parameter space is $`M=U\left(1\right)`$ (the symmetry group of quantum electrodynamics). Topologically, $`U\left(1\right)=S^1`$. Thus, soliton solutions, in this case, realize nontrivial maps $`S^1S^1`$. All the continuous maps $`S^1S^1`$ have a group structure of the fundamental (or first homotopy) group $`\pi _1\left(S^1\right)`$.BP75 ; Bo78 ; M79 ; R82 ; DEGM82 ; GM86 ; S89  Given that $`S^1=𝐑/𝐙`$ ($`𝐑`$ is the additive group of the real numbers), we can write

$$\pi _1\left(S^1\right)=\pi _1\left(𝐑/𝐙\right)=\pi _0\left(𝐙\right)=𝐙,$$
(57)
which should be compared with (52). As in the case of layered superconductors, the external magnetic field $`H>0`$ breaks the symmetry $`\phi \phi `$. Thus, only the states parameterized by $`N_v=0,1,\mathrm{}`$ are possible, with $`N_v=0`$ being the ”vacuum”, Meissner state.

To conclude this discussion, we have to clarify a typical misunderstandingKr02  concerning the role of the laminar modeldG  in type-II superconductivity. In reality, the order-parameter space of a continuum type-II superconductor, $`M=S^1`$, precludes the existence of topologically stable plane defects, envisaged by the laminar model. Indeed, consider two points $`P_1=(x_0,a,z_0)`$, $`P_2=(x_0,b,z_0)`$ on the opposite sides of such a defect, in unperturbed regions of the superconductor. Join these points by a continuous path, parameterized by $`ayb`$. The boundary of the interval $`[a,b]`$ is a $`0`$-dimensional sphere $`S^0=\{a,b\}`$, which leads us to a consideration of the maps $`S^0S^1`$. However, the pertinent homotopy group $`\pi _0\left(S^1\right)`$ is trivial,BP75 ; Bo78 ; M79 ; R82 ; GM86 ; S89  i.e.,

$$\pi _0\left(S^1\right)=0,$$
(58)
in contrast to (52) and (57). Topological instability of the ”laminar solution” hardly allows one to expect that this solution corresponds to any local minimum of the Ginzburg-Landau free-energy functional. Therefore, a comparison with unstable ”isolated fluxons” or ”triangular Josephson-vortex lattices” in layered superconductors is much more appropriate than the far-fetchedKr02  ”similarity” to the vortex plane.

### V.3 The interpretation of experimental data

Experimental observations of the vortex structure in layered superconductors can be roughly subdivided into two groups: (i) direct observations, allowing one to”visualize” the flux distribution;N96 ; Yu98 ; Mo98  (ii) indirect observations (i.e., measurements of $`c`$-axis transport properties,So96 ; La96 ; Kr00  magnetization,De99  and the upper critical fieldRu80 ). Here, we present an overview of these observations, showing that all the experimental data available up to now can be explained in terms of the stable vortex-plane configurations. A detailed quantitative analysis can be done with the use of the results of Refs. K99, ; K01, ; K1, .

Josephson-flux distribution, characteristic of vortex planes for $`H>0`$, was directly observed on artificial low-$`T_c`$ stacked junctions in Ref. N96,  (by low-temperature scanning electron microscopy) and in Ref. Yu98,  (by polarized neutron reflection). In particular, the double-junction stackN96  ($`N=3`$) has revealed the phase-difference symmetry $`\varphi _1=\varphi _2`$, exactly as could be expected from the general relations (5) for vortex planes. In both the experiments,N96 ; Yu98  penetration of the flux occurred simultaneously and coherently into all the junctions, in full agreement with the scenario for the vortex planes.K99 ; K01  Moreover, accompanying measurements of magnetization in Ref. Yu98,  have shown typical oscillations and hysteresis. The oscillations should be viewed as a manifestation of a series of first-order phase transitions, discussed in section IV.2 and Refs. K99, ; K01, , whereas the hysteresis is implied by the overlap of the regions (46), (47).

We draw attention to a possible application to high-$`T_c`$ superconductivity. Oscillations of magnetization in parallel fields, interpreted as evidence of Josephson nature of the flux, have been reported for YBCO in Ref. De99, . According to Ref. De99, , ”the temperature dependence of the magnetization contradicts the present theoretical expectations”.

As shown by our self-consistent calculations,K99 ; K01  the oscillating behavior of the critical Josephson current $`I_c\left(H\right)`$ (the Fraunhofer pattern) is a result of successive penetration of vortex planes and their pinning by the edges of the superconductor. Oscillating $`I_c\left(H\right)`$ dependencies have been observed both on artificial low-$`T_c`$ stacked junctionsSo96 ; Kr00  and high-$`T_c`$ layered superconductors BSCCO.La96 ; Kr00  ”Irregularities” of the dependence $`I_c\left(H\right)`$,So96 ; Kr00  such as, e.g., multivaluedness and aperiodicity, can be easily explained by the overlap of the regions (46), (47). Behavior of this type was observed a long time ago on the single Josephson junction,Sch70  which confirmed the theoretical prediction of the overlap for ordinary Josephson vortices.OS67

The most ”ancient” experimental confirmation of the stability of vortex planes is provided by observations of the ”crossover” behavior of $`H_{c2}\left(T\right)`$ in artificial low-$`T_c`$ stacked junctions.Ru80  For $`HH_{c2}`$, the condition $`f_n=1`$, employed in the derivation of Eqs. (2), in no longer valid. However, periodic modulations of $`f_n\left(y\right)`$, caused by the presence of vortex planes and therefore identical in all the S-layers, account for the observed behavior of $`H_{c2}\left(T\right)`$.K99 ; K01

Finally, we want to comment on direct observations of non-equilibrium isolated vortices in layered high-$`T_c`$ superconductors at $`H=0`$.Mo98  As is explained at the end of subsection IV.A, nontrivial flux configurations cannot exist in a layered superconductor with ideal periodicity at $`H=0`$, $`L<\mathrm{}`$. However, the presence of structural defects (e.g., stacking faults, as is hinted in Ref. Mo98, ) violates the condition of ideal periodicity and should stabilize energetically an otherwise unstable configuration. In this situation, we indeed expect to obtain non-equilibrium isolated vortices, because their self-energy is lower and the $`c`$-axis extent is smaller than those of vortex planes at $`H=0`$, $`L<\mathrm{}`$.K01  A detailed mathematical analysis of this case can be done on the basis of the results of Ref. K1, .

## VI Conclusions

In brief, we have solved the problem of the classification of all solutions to (2), (10) with respect to their stability. In our consideration, we have employed exact methods of the calculus of variations, soliton physics and the exact results of Ref. K1,  \[the expression for the Gibbs free-energy functional (12), the first integral (11) and the solution to the soliton boundary value problem (2), (4), (44)\]. In view of obvious mathematical complexity, the problem of the stability of vortex configurations could not be solved by use of inadequate methods, employed, e.g., in Refs. B73, ; CC90, ; BC91, ; BLK92, ; Ko93, ; Kr02, , and was not even posed in any of these publications.

In full agreement with the fundamentals of soliton physics, we have proved that the only minimizers of both the Gibbs and Helmholtz free energy functionals are the Meissner solution (the ”vacuum” state) and soliton vortex-plane solutions. They represent the actual equilibrium field configurations. The obtained results allowed us to explain exceptional stability of vortex planes, established in numerical simulations, and to provide a unified interpretation of the experimental data available up to now.

In contrast, non-soliton configurations (”isolated fluxons”, triangular Josephson-vortex lattices”, etc.), proposed in previous publications, turned out to be absolutely unstable and unobservable: They are nothing but saddle points of the Gibbs free-energy functional and are not even stationary points of the Helmholtz free-energy functional. Physically, these configurations violate conservation laws for the flux and the current.

One may ask a natural question why exactly the unstable configurations were previously proposed as the ”equilibrium state”, whereas the actual minimizers of the free energy (vortex-plane solitons) were neglected. We think that the answer lies in the following:

i) the hypothesis of an ”analogy” with the Abrikosov-vortex structure in type-II superconductors, accepted without any mathematical justification;B73

ii) the absence of an exact mathematical definition of the ”Josephson vortex”. In Refs. B73, ; CC90, , ”isolated Josephson vortices” were discussed without any consideration of the SG equations. The fact that Josephson vortices are nothing but static solitons of the SG equations was not realized in subsequent publications, either. For example, the existence of Josephson vortices in the case $`W2\lambda _J`$ was denied in Refs. BCG91, ; Kr01, , which is refuted by our Eqs. (47), valid for any $`W`$;

iii) obvious mathematical mistakes in the treatment of the Lawrence-Doniach modelLD  for infinite ($`N=\mathrm{}`$) layered superconductors. The neglectBC91 ; BCG91 ; BLK92 ; Ko93  of the surface variation in the variational principle for the Lawrence-Doniach functional resulted in a loss of the conservation laws for the current and the flux, as was first pointed out in Ref. K99, ; K01, : Soliton solutions are a corollary of these conservation laws;

iv) the absence of any investigation of analytical properties of the coupled static SG equations for $`H>0`$, $`W<\mathrm{}`$. Pertinent soliton solutions were obtained in our papers: in the exactly solvable cases $`N=\mathrm{}`$ (Refs. K99, ; K01, ), $`N=2,3`$ (Refs. K1, ; K2, ), and in the general case $`2N<\mathrm{}`$ (Ref. K1, ). Standard methods of soliton physicsZK74 ; BP75 ; Bo78 ; R82 ; DEGM82 ; GM86 ; S89  as well as advanced methods of the calculus of variations and of the theory of differential equations, employed in our analysis, were completely disregarded in previous theoretical publications;BC91 ; BCG91 ; BLK92 ; Ko93

v) the absence of any attempts to analyze the stability of the proposed ”vortex” configurations, both in theoretical publicationsB73 ; CC90 ; BC91 ; BLK92 ; Ko93  and numerical simulations.Kr01  Dynamic stability of vortex planes, noticed in numerical simulations,SBP93 ; Kl94 ; Ma00  was not understood and neglected;

vi) the neglect of direct experimental observationsN96 ; Yu98  of the Josephson-vortex structure at $`H>0`$: These observations have clearly revealed that exactly the vortex planes (not ”isolated fluxons” or ”triangular lattices”) are the actual equilibrium field configurations;

vii) long-term domination of the subjective point of viewB73  and the absence of any pluralism of opinion. As a result, the critical remarksF98  are neglected, whereas the attempts to clarify the situation within the framework of a rigorous mathematical approachK99 ; K01  are immediately attackedKr02  with the use of inappropriate methods.r1

We hope that this paper will finally convince both theorists and experimentalists, specializing in the field of weak superconductivity, of the necessity to give up the old, unsound theoretical prejudices: The wealth of magnetic properties of layered superconductors (both low- and high-$`T_c`$) cannot be understood without the solitons. One should also think of possible practical applications of the vortex-plane solitons, e.g., in submillimeter-wave generators, as is proposed in Ref. Ma00, . Given the role of the single SG equation in different fields of physics (quantum optics;DEGM82  the Skyrme and the Thirring models in elementary particle physics;Bo78 ; R82 ; DEGM82  the theory of dislocations and magnetism, let alone the Josephson effect, in condensed matter physics,DEGM82  etc.), we expect that our exact results for the coupled SG equations (including the single one as a particular case) may find applications in these fields as well.

The coupled SG equations for $`H>0`$, $`W<\mathrm{}`$ have not been studied in mathematical literature, either. Our exact analytical results for the static case constitute only the first step in this direction. The next stage should be analytical properties of time-dependent equations. Our paper may stimulate interest in this problem of specialists in applied mathematics as well.

## Appendix A The properties of the matrix $`G(n,m)`$

The explicit form of $`G(n,m)`$ is

$$G(n,m)=\frac{1}{2ϵ\sqrt{1+\frac{ϵ^2}{4}}}\left[\mu ^{\left|nm\right|}\frac{\mu ^n\left(\mu ^{mN}\mu ^{Nm}\right)+\mu ^{Nn}\left(\mu ^m\mu ^m\right)}{\mu ^N\mu ^N}\right],$$
(59)
where

$$\mu =1+\frac{ϵ^2}{2}ϵ\sqrt{1+\frac{ϵ^2}{4}}.$$
(60)
The following properties of $`G(n,m)`$ are obvious:

$$G(n,m)=G(m,n),$$
(61)
$$G(n,Nm)=G(Nn,m).$$
(62)
The matrix $`G(n,m)`$ is positive definite, since all its eigenvalues $`e_j`$ are positive:

$$e_j=\frac{\lambda _j^2}{ϵ^2},\lambda _j=\frac{ϵ}{\sqrt{2+ϵ^22\mathrm{cos}\frac{\pi j}{N}}},j=1,2,\mathrm{},N1.$$
(63)
Of importance are the summation rules:

$$\underset{m=1}{\overset{N1}{}}G(n,m)=\frac{1}{ϵ^2}\left[1G(n,1)G(n,N1)\right]$$

$$=\frac{1}{ϵ^2}\left[1\frac{\mu ^n+\mu ^{N+n}\mu ^n\mu ^{Nn}}{\mu ^N\mu ^N}\right],1nN1;$$
(64)
$$\underset{n=1}{\overset{N1}{}}\underset{m=1}{\overset{N1}{}}G(n,m)=\frac{1}{ϵ^2}\left[N1\frac{2\sqrt{1+\frac{ϵ^2}{4}}ϵ}{ϵ}\frac{1\mu ^{N1}}{1+\mu ^N}\right]\frac{N1}{ϵ^2H_s^2},$$
(65)
where

$$H_s=\left[1\frac{2\sqrt{1+\frac{ϵ^2}{4}}ϵ}{ϵ\left(N1\right)}\frac{1\mu ^{N1}}{1+\mu ^N}\right]^{\frac{1}{2}}$$

$$\sqrt{\frac{\left(N1\right)N}{2}}\left[\underset{k=0}{\overset{\left[\frac{N}{2}\right]1}{}}\lambda _{2k+1}^2\mathrm{cos}{}_{}{}^{2}\frac{\pi \left(2k+1\right)}{2N}\right]^{\frac{1}{2}}$$
(66)
is the superheating (penetration) field of a semiinfinite ($`0y<\mathrm{}`$) layered superconductor.K1

## Appendix B A relationship to the variational principle for infinite layered superconductors ($`N=\mathrm{}`$)

The intralayer currents for $`f_n=1`$ are given byK1

$$J_n\left(y\right)=\frac{d\phi _n\left(y\right)}{dy}2A(n,y)$$

$$=\frac{1}{4\pi }\left[H_n\left(y\right)H_{n+1}\left(y\right)\right],n=0,1,\mathrm{},N1,$$
(67)
where $`H_0\left(y\right)=H_N\left(y\right)H`$. Using the second relation in (67), we get

$$\underset{n=0}{\overset{N1}{}}J_n\left(y\right)=0,$$
(68)
which is the conservation law for the total intralayer current. Moreover, in view of (5) and (61), (62),

$$H_n\left(y\right)=H_{Nn}\left(y\right).$$
(69)
Hence,

$$J_n\left(y\right)=J_{Nn1}\left(y\right).$$
(70)
Using the first relation in (67), we write

$$\underset{L}{\overset{L}{}}𝑑y\left[J_n\left(y\right)J_{n1}\left(y\right)\right]$$

$$=\varphi _n\left(L\right)\varphi _n\left(L\right)\underset{L}{\overset{L}{}}𝑑y\left[A(n,y)A(n1,y)\right],n=1,2,\mathrm{},N1.$$
(71)
The second term on the right-hand side of (71) is the flux between the S-layers $`n`$ and $`n1`$. We can therefore rewrite (71) using (8):

$$\underset{L}{\overset{L}{}}𝑑y\left[J_n\left(y\right)J_{n1}\left(y\right)\right]=HW\left[G(n,1)+G(n,N1)\right]$$

$$+\varphi _n\left(L\right)\varphi _n\left(L\right)+\frac{ϵ^2}{2}\underset{m=1}{\overset{N1}{}}G(n,m)\left[\varphi _m\left(L\right)\varphi _m\left(L\right)\right],n=1,2,\mathrm{},N1.$$
(72)
In view of the flux-conservation conditions (24), the variation of the right-hand side of (72) vanishes, hence

$$\delta J_n\left(y\right)=\delta J_{n1}\left(y\right),n=1,2,\mathrm{},N1.$$
(73)
Combined with the current-conservation law (68), relations (73) yield:

$$\delta J_n\left(y\right)=0,n=0,1,\mathrm{},N1,$$
(74)
which means that partial intralayer currents are also conserved.

Consider the case $`N1`$. For $`n`$ satisfying the condition $`\left[ϵ^1\right]nN1\left[ϵ^1\right]`$, we can proceed to the limit $`N\mathrm{}`$ in the second relation (67), obtaining

$$J_n\left(y\right)=\frac{d\phi _n\left(y\right)}{dy}2A(n,y)=0.$$
(75)
This is exactly the result derived for the infinite ($`N=\mathrm{}`$) layered superconductor in Refs. K99, ; K01,  by means of an exact variational principle, based on the use of the conservation law for the total intralayer current.