text:
# The Meissner effect and the vortex structure in stacked junctions and layered superconductors: Exact analytical results

## I Introduction

We present a rigorous mathematical examination of the problem of the Meissner effect and of the vortex structure in thin-layer Josephson-junction stacks and layered superconductors, with a static, homogeneous external magnetic field $`H>0`$ applied parallel to the layers (along the $`z`$ axis, see Fig. 1.) We consider periodic systems composed of an arbitrary number $`N1`$ of identical superconductor/insulator (S/I) junctions ($`2N<\mathrm{}`$, where $`N`$ is the number of S-layers, with $`x`$ being the layering axis), occupying an arbitrary region $`[L,L]`$ along the $`y`$ axis.

We begin with the investigation of analytical properties of a finite set of coupled static sine-Gordon (SG) equations for phase differences $`\varphi _n`$($`1nN1`$), obtained by the minimization of a microscopic Gibbs free-energy functional. Equations of this type were first derived within the framework of a phenomenological approach. Under corresponding redefinition of the parameters, they also apply to the phenomenological Lawrence-Doniach (LD) model with $`N<\mathrm{}`$. Although particular numerical solutions to these equations were obtained in several publications, any analytical analysis of their general properties has not been undertaken up until now.

Making use of standard methods of the theory of differential equations, we arrive at a complete classification of the solutions to the SG equations subject to appropriate physical boundary conditions at $`y=\pm L`$. In particular, we prove that for $`H>0`$ these equations admit only Meissner solutions and topological ”vortex-plane” solutions. Both the types of solutions are characterized, in general, by $`\left[\frac{N}{2}\right]`$ Josephson lengths $`\lambda _{Jk}`$ ($`\left[\frac{N}{2}\right]`$ is the integer part of $`\frac{N}{2}`$, $`k=0,1,\mathrm{},\left[\frac{N}{2}\right]1`$). The Meissner solutions persist up to a certain superheating field of the Meissner state, $`H_{sL}`$. We show that the field $`H_{sL}`$ simultaneously determines the penetration field for the vortex planes. For $`L\lambda _{J\mathrm{max}}\lambda _{J0}`$, we derive an explicit analytical expression for the superheating (penetration) field $`H_sH_s\mathrm{}`$ as a function of $`N2`$. For $`HH_s`$, we obtain a closed-form analytical expression for the Meissner local field $`H_n(y)`$, valid for any $`N2`$, and investigate a transition to the infinite ($`N=\mathrm{}`$) layered-superconductor limit.

Thermodynamically stable vortex-plane solutions, physically, represent coherent chains of Josephson vortices (one vortex per each I-layer in a chain), positioned in planes parallel to the coordinate axes $`x`$, $`z`$. (See Fig. 2.) Such solutions were thoroughly studied for infinite layered superconductors in Refs.  . For $`H_{c1}HH_{c2}`$, they correspond to the ”transparent state” considered by Theodorakis within the framework of the infinite LD model. In the case of a double-junction stack ($`N=3`$), they coincide with the well-known ”in-phase mode”. Significantly, coherent Josephson vortex configurations that can be identified with the vortex-plane solutions have been directly observed in experiments on artificial double-junction stacks and weakly-coupled multilayers. Besides giving a proof of the existence and stability of the vortex-plane solutions for $`3`$$`N<\mathrm{}`$, we discuss their main physical properties, such as periodicity in the $`y`$ direction and the overlapping of states with different numbers of vortex planes $`N_v`$. In the limit $`L\lambda _{J\mathrm{max}}\lambda _{J0}`$, we derive exact analytical expressions for the self-energy $`E_v`$ of the state with $`N_v=1`$ and for the lower critical field $`H_{c1}`$. We also investigate a physical and mathematical relationship between vortex-plane solutions for $`N3`$ and the ordinary Josephson vortices in a single junction ($`N=2`$).

In contrast to a prevailing view, the coupled SG equations do not possess any single-vortex solutions for $`H>0`$ and $`L<\mathrm{}`$, as well as such vortex solutions as the ”triangular lattice” and the ”row mode”. Thus, single-vortex solutions appear only in the limit $`H=0`$, $`L=\mathrm{}`$, which, physically, signifies their absolute thermodynamic instability. In this regard, it should be emphasized that single-vortex configurations were originally proposed for infinite layered superconductors without an appropriate analysis of the actual conditions of their existence. (See Refs.  for the criticism.) An excessive preoccupation with single-vortex configurations in previous publications can be explained by the confusion of the problem of finding the lowest-energy topological solutions at $`H=0`$ with the problem of the minimization of the Gibbs free energy at $`H>0`$: Although at $`H=0`$, $`L=\mathrm{}`$ the energy of single-vortex configurations is lower than that of the vortex planes, the former ones are irrelevant to the problem of minimization because of their absolute instability at $`H>0`$, $`L<\mathrm{}`$. For the same reason, single-vortex solutions for $`H=0`$, $`L=\mathrm{}`$ cannot be used for estimates of the lower critical field $`H_{c1}`$. The absence of single-vortex configurations at $`H>0`$ was established for infinite layered superconductors by the use of exact variational methods in Refs.  . In order to clarify this issue for $`N<\mathrm{}`$, we analyze those features of single-vortex solutions at $`H=0`$, $`L=\mathrm{}`$ which preclude their existence at $`H>0`$, $`L<\mathrm{}`$.

Section II of the paper is devoted to exact mathematical formulation of the problem. In section III, we derive all major mathematical and physical results sketched above. Mathematical Lemmas 1, 2 and the Theorem of subsection IIIA form the basis for a subsequent physical analysis in subsections IIIB-IIID. The general consideration of section III is illustrated in section IV by two exactly-solvable examples, namely of a single thin-layer junction ($`N=2`$) and of a double-junction stack ($`N=3`$). (In section IV, we present a new exact solution to the single SG equation, valid for arbitrary $`H0`$ and $`L>0`$.) The obtained results are discussed in section V. Appendices A-C contain mathematical derivations and proofs omitted in the main text.

## II Formulation of the problem

We consider a periodic structure consisting of alternating $`N`$ superconducting(S) and $`N1`$ insulating (I) layers ($`2N<\mathrm{}`$) in a parallel, static, homogeneous external magnetic field $`0<HH_{c2}`$, where $`H_{c2}`$ is the upper critical field. (See Fig. 1.) The S-layer thickness is $`a`$, and $`p`$ is the period; the length of the structure in the $`y`$ direction is $`W=2L`$, and $`W_z`$ is the length of the structure in the $`z`$ direction ($`W_z\mathrm{}`$). We set $`\mathrm{}=c=1`$ and assume that

$$\frac{T_{c0}T}{T_{c0}}1,$$
(1)
$$\xi _0a,$$
(2)
$$a\mathrm{min}\{\zeta (T),\lambda (T),\alpha ^1\xi _0\},\alpha \frac{3\pi ^2}{7\zeta (3)}\underset{0}{\overset{1}{}}𝑑ttD(t)1,$$
(3)
$$ap.$$
(4)
Here $`T_{c0}`$ is the critical temperature of an isolated S-layer, $`\xi _0`$ is the BCS coherence length; $`\zeta (T)`$ and $`\lambda (T)`$ are, respectively, the Ginzburg-Landau (GL) coherence length and the penetration depth; $`D(\mathrm{cos}\theta )`$ is the incidence-angle-dependent tunneling probability of the I-layer between two successive S-layers.

Conditions (1) and (2) ensure the applicability of the GL-type expansion within each S-layer. Condition (3) corresponds to the thin S-layer limit, whereas condition (4) is employed here for the sake of mathematical simplicity only. To simplify further the mathematical description, we introduce dimensionless units by

$`{\displaystyle \frac{x}{p}}x,`$

$`{\displaystyle \frac{y}{\overline{\lambda }_J}}y,`$

$`{\displaystyle \frac{H}{\overline{H}_s}}H,`$

where the quantities on the left-hand side are dimensional, with $`\overline{\lambda }_J=\left(8\pi ej_0p\right)^{1/2}`$ being the Josephson penetration depth ($`j_0`$ is the density of the Josephson current in a single junction with thick electrodes) and $`\overline{H}_s=\left(ep\overline{\lambda }_J\right)^1`$ being the superheating (penetration) field of the infinite layered superconductor. In our dimensionless units, for example, the flux quantum is $`\mathrm{\Phi }_0=\pi `$, and the lower critical field of the infinite layered superconductor is $`\overline{H}_{c1}=\frac{2}{\pi }`$.

Under the above conditions, the structure is completely described by a closed set of self-consistent microscopic equations for the reduced modulus of the superconducting order parameter in the $`n`$th S-layer $`f_n(y)`$ ($`0f_n1`$, $`n=0,1,\mathrm{},N1`$) and the phase difference between the $`n`$th and $`(n1)`$th S-layers $`\varphi _n(y)`$ ($`\varphi _n\phi _n\phi _{n1}`$, $`n=1,2,\mathrm{},N1`$):

$`f_0(y)f_0^3(y)=r(T)\left[{\displaystyle \frac{ϵ^2}{2}}{\displaystyle \frac{d^2f_0(y)}{dy^2}}+{\displaystyle \frac{2\left[HH_1(y)\right]^2}{ϵ^2f_0^3(y)}}+{\displaystyle \frac{1}{2}}\left[f_0(y)f_1(y)\mathrm{cos}\varphi _1(y)\right]\right],`$

$`f_n(y)f_n^3(y)=r(T)[{\displaystyle \frac{ϵ^2}{2}}{\displaystyle \frac{d^2f_n(y)}{dy^2}}+{\displaystyle \frac{2\left[H_n(y)H_{n+1}(y)\right]^2}{ϵ^2f_n^3(y)}}`$

$$+\frac{1}{2}[2f_n(y)f_{n+1}(y)\mathrm{cos}\varphi _{n+1}(y)f_{n1}(y)\mathrm{cos}\varphi _n(y)]],1nN2,$$
(5)
$`f_{N1}(y)f_{N1}^3(y)=r(T)[{\displaystyle \frac{ϵ^2}{2}}{\displaystyle \frac{d^2f_{N1}(y)}{dy^2}}+{\displaystyle \frac{2\left[HH_{N1}(y)\right]^2}{ϵ^2f_{N1}^3(y)}}`$

$`+{\displaystyle \frac{1}{2}}[f_{N1}(y)f_{N2}(y)\mathrm{cos}\varphi _{N1}(y)]];`$

$$\frac{df_n}{dy}\left(\pm L\right)=0,0nN1;$$
(6)
$`{\displaystyle \frac{1}{f_n^2(y)}}\left[H_{n+1}(y)H_n(y)\right]{\displaystyle \frac{1}{f_{n1}^2(y)}}\left[H_n(y)H_{n1}(y)\right]ϵ^2H_n(y)`$

$$=\frac{ϵ^2}{2}\frac{d\varphi _n(y)}{dy},1nN1,$$
(7)
$$H_0(y)=H_N(y)=H,$$
(8)
$$\frac{d\varphi _n}{dy}\left(\pm L\right)=2H,1nN1,$$
(9)
where

$`r(T){\displaystyle \frac{\zeta ^2(T)\alpha }{a\xi _0}},ϵ{\displaystyle \frac{\sqrt{ap}}{\lambda }},`$

and the local magnetic field in the $`n`$th I-layer ($`n1<x<n`$) is given by

$`H_n(y)={\displaystyle \frac{1}{2}}{\displaystyle \underset{L}{\overset{y}{}}}𝑑uf_n(u)f_{n1}(u)\mathrm{sin}\varphi _n(u)+H`$

$$=\frac{1}{2}\underset{L}{\overset{y}{}}𝑑uf_n(u)f_{n1}(u)\mathrm{sin}\varphi _n(u)+H.$$
(10)
Note that although in most physical situations $`ϵ1`$, we will not need the smallness of $`ϵ`$ in our mathematical consideration. Because of the obvious property $`f_n(y)=f_n(y)`$, equation (10) implies that the phase differences $`\varphi _n`$ meet the condition

$$\varphi _n(y)=\varphi _n(y)+0\text{mod}2\pi .$$
(11)
The dimensionless Gibbs free energy $`\mathrm{\Omega }(H)`$, normalized via the relation

$`{\displaystyle \frac{4\pi \mathrm{\Omega }(H)}{H_c^2(T)a\lambda _J\mathrm{}W_z}}\mathrm{\Omega }(H),`$

in terms of the mean-field quantities $`f_n`$, $`\varphi _n`$ and $`H_n(y)`$ has the form

$`\mathrm{\Omega }(H)={\displaystyle \underset{n=0}{\overset{N1}{}}}{\displaystyle \underset{L}{\overset{L}{}}}dy[f_n^2(y)+{\displaystyle \frac{1}{2}}f_n^4(y)+{\displaystyle \frac{r(T)ϵ^2}{2}}\left({\displaystyle \frac{df_n(y)}{dy}}\right)^2`$

$`+{\displaystyle \frac{2r(T)}{ϵ^2f_n^2(y)}}[H_{n+1}(y)H_n(y)]^2+2r(T)[H_n(y)H]^2]`$

$$+\frac{r(T)}{2}\underset{n=1}{\overset{N1}{}}\underset{L}{\overset{L}{}}𝑑y\left[f_{n1}^2(y)+f_n^2(y)2f_n(y)f_{n1}(y)\mathrm{cos}\varphi _n(y)\right].$$
(12)
Here, the sum of the three terms in the first line on the right-hand side is the condensation energy, the sum of the two terms in the second line is the electromagnetic energy $`E_{em}`$, and the last term is the Josephson energy $`E_J`$. The intralayer current in the $`n`$th S-layer $`J_n(y)`$ (normalized to $`\overline{H}_s`$) and the density of the Josephson current between the $`n`$th and the $`(n1)`$th S-layers $`j_{n,n1}(y)`$ (normalized to $`j_0`$) are given by

$$J_n(y)=\frac{1}{4\pi }\left[H_n(y)H_{n+1}(y)\right],0nN1,$$
(13)
and

$$j_{n,n1}(y)=2\frac{dH_n(y)}{dy}=f_n(y)f_{n1}(y)\mathrm{sin}\varphi _n(y),1nN1,$$
(14)
respectively.

Assuming that the temperature range satisfies the condition of the weak-coupling limit

$$r(T)1,$$
(15)
one can obtain a perturbative solution for $`f_n`$ and $`\varphi _n`$ up to any desired order in $`r(T)`$, starting from the zero-order solution to (5), (6),

$$f_n=1,$$
(16)
and the zero-order equations for $`\varphi _n`$,

$$H_{n+1}(y)\left(2+ϵ^2\right)H_n(y)+H_{n1}(y)=\frac{ϵ^2}{2}\frac{d\varphi _n(y)}{dy},1nN1,$$
(17)
where $`H_n(y)`$ are given by (10) with $`f_n=1`$ and satisfy the boundary conditions (8). For most applications, it is sufficient to consider expressions for physical quantities only in leading order in $`r(T)`$. Thus, for example, substituting (16) and the solution of (17) into (12) immediately yields a first-order expansion for the Gibbs free energy, because first-order corrections to the condensation-energy term cancel out.

A detailed mathematical analysis of Eqs. (17) is the subject of section III. Here we point out that these equations can be transformed into a very useful for application form by solving for $`H_n(y)`$ (see Appendix A for mathematical details):

$$H_n(y)=h_n(y)+H_n,$$
(18)
$$h_n(y)=\frac{ϵ^2}{2}\underset{m=1}{\overset{N1}{}}G(n,m)\frac{d\varphi _m(y)}{dy},$$
(19)
$$H_n=\frac{H\left(\mu ^n+\mu ^{N+n}\mu ^n\mu ^{Nn}\right)}{\mu ^N\mu ^N},$$
(20)
where $`G(n,m)`$ are given by (A9), and $`\mu `$ is given by (A5). By (A8), and (9), (A13), expression (18) explicitly satisfies boundary conditions (8) and $`H_n(\pm L)=H`$. Moreover, the $`y`$-independent quantities $`H_n`$ in (18) have clear physical meaning: Being solutions of (17) with $`\frac{d\varphi _m}{dy}0`$, they describe distribution of the local magnetic field within the I-layers in the homogeneous Meissner state. (See section III of Ref. .) \[Note that in an infinite layered superconductor $`H_n0`$, and $`\underset{m=1}{\overset{N1}{}}G(n,m)\mathrm{}\underset{m=\mathrm{}}{\overset{+\mathrm{}}{}}G_{\mathrm{}}(n,m)\mathrm{},`$ where $`G_{\mathrm{}}(n,m)`$ are defined by (A17).\] By comparing Eqs. (18) with Eqs. (10) (where $`f_n=1`$), using the property $`H_n=H_{Nn}`$ and (A11), we establish the symmetry relations

$$\varphi _n(y)=\varphi _{Nn}(y),h_n(y)=h_{Nn}(y),H_n(y)=H_{Nn}(y).$$
(21)
Relations (21) replace the relations $`\varphi _n(y)=\varphi _{n+1}(y)\varphi (y)`$ and $`h_n(y)=h_{n+1}(y)h(y)`$ of the infinite layered superconductor.

By the use of Eqs. (18), the Gibbs free energy up to first order in $`r\left(T\right)1`$ can be given the form

$`\mathrm{\Omega }(H)={\displaystyle \frac{NW}{2}}+r(T)[{\displaystyle \frac{2H^2W}{H_s^2}}(N1)`$

$`+{\displaystyle \frac{ϵ^2}{2}}{\displaystyle \underset{n=1}{\overset{N1}{}}}{\displaystyle \underset{m=1}{\overset{N1}{}}}G(n,m){\displaystyle \underset{L}{\overset{L}{}}}𝑑y{\displaystyle \frac{d\varphi _n(y)}{dy}}{\displaystyle \frac{d\varphi _m(y)}{dy}}`$

$$+2\underset{n=1}{\overset{N1}{}}\underset{L}{\overset{L}{}}dy\mathrm{sin}{}_{}{}^{2}\frac{\varphi _n(y)}{2}4H\underset{n=1}{\overset{N1}{}}\mathrm{\Phi }_n\frac{\varphi _n(L)\varphi _n(L)}{2\pi }],$$
(22)
$$\mathrm{\Phi }_n=\pi \left[1\frac{\mu ^n+\mu ^{N+n}\mu ^n\mu ^{Nn}}{\mu ^N\mu ^N}\right],$$
(23)
where $`H_s`$ is the superheating (penetration) field defined in Eq. (A14), and $`\mathrm{\Phi }_n`$ is the total flux carried by a Josephson vortex positioned in the $`n`$th I-layer. (The quantities $`H_s`$ and $`\mathrm{\Phi }_n`$ are thoroughly discussed in section IIIB.) The quadratic form of the electromagnetic energy \[the second line on the right-hand side of Eq. (22)\] can be diagonalized with the help of (A15), (A16):

$$E_{em}=\frac{1}{N}\underset{k=0}{\overset{\left[\frac{N}{2}\right]1}{}}\lambda _{2k+1}^2\underset{n=1}{\overset{N1}{}}\underset{m=1}{\overset{N1}{}}\mathrm{sin}\frac{\pi n\left(2k+1\right)}{N}\mathrm{sin}\frac{\pi m\left(2k+1\right)}{N}\underset{L}{\overset{L}{}}𝑑y\frac{d\varphi _n(y)}{dy}\frac{d\varphi _m(y)}{dy},$$
(24)
where $`\left[u\right]`$ is the integer part of $`u`$. \[In deriving (24), we employed the symmetry relations (21).\] Equation (24) explicitly shows that physical solutions in a stack with $`N`$ S-layers are characterized by $`\left[\frac{N}{2}\right]`$ Josephson length scales $`\lambda _{Jk}\lambda _{2k+1}`$ ($`k=0,1,\mathrm{},\left[\frac{N}{2}\right]1`$) from the set (A15).

## III Major results

### A Analysis of the equations for the phase differences

By differentiation with respect to $`y`$, integrodifferential equations (17) reduce to a system of $`N1`$ ordinary nonlinear second-order differential equations

$`{\displaystyle \frac{d^2\varphi _1(y)}{dy^2}}={\displaystyle \frac{1}{ϵ^2}}\left[\left(2+ϵ^2\right)\mathrm{sin}\varphi _1(y)\mathrm{sin}\varphi _2(y)\right],`$

$`{\displaystyle \frac{d^2\varphi _n(y)}{dy^2}}={\displaystyle \frac{1}{ϵ^2}}\left[\left(2+ϵ^2\right)\mathrm{sin}\varphi _n(y)\mathrm{sin}\varphi _{n+1}(y)\mathrm{sin}\varphi _{n1}(y)\right],2nN2,`$

$$\frac{d^2\varphi _{N1}(y)}{dy^2}=\frac{1}{ϵ^2}\left[\left(2+ϵ^2\right)\mathrm{sin}\varphi _{N1}(y)\mathrm{sin}\varphi _{N2}(y)\right]$$
(25)
with boundary conditions (9). As shown in section II, physical solutions to (25) necessarily obey the symmetry relations (11) and (21).

A detailed analysis of the analytical structure of the coupled static SG equations (25) is the main task of this subsection. Such differential equations were first derived within the framework of a less rigorous, phenomenological approach in Ref.  . Particular numerical solutions (for $`H=0`$) were obtained in several publications. However, any analytical investigation of these equations has not been carried out up until now.

#### 1 Existence and uniqueness of the solution to the initial value problem

Consider Eqs. (25) on the whole axis $`\mathrm{}<y<+\mathrm{}`$. Two simple properties of (25) are quite obvious: If $`\varphi _n(y)`$ ($`1nN1`$) is a solution, the functions $`\overline{\varphi }_n(y)`$ given by

$$\overline{\varphi }_n(y)=\varphi _n(y)+2\pi k\text{(}k\text{ is an integer),}$$
(26)
and

$$\overline{\varphi }_n(y)=\varphi _n(y+c)\text{(}c\text{ is an arbitrary constant)}$$
(27)
are also solutions. \[The latter is a result of the fact that $`y`$ does not enter explicitly the right-hand side of (25).\] Our conclusions about the solution to (25) will be substantially based on another key property, which we formulate as a lemma:

Lemma 1. Consider an arbitrary interval $`I=[L_1,L_2]`$ and $`y_0I`$. The initial value problem for Eqs. (25) with arbitrary initial conditions $`\varphi _n(y_0)=\alpha _n`$, $`\frac{d\varphi _n}{dy}(y_0)=\beta _n`$ has a unique solution in the whole interval $`I`$. This solution has continuous derivatives with respect to $`y`$ of arbitrary order and continuously depends on the initial data. (For the proof of Lemma 1, see Appendix B.)

It is worth noting that the existence and uniqueness of a smooth solution to the initial value problem in the whole interval $`I`$ is rather nontrivial for nonlinear differential equations: For such equations, theorems of existence and uniqueness are usually valid only locally, in the neighborhood of initial data. In our case, global character of the solution and its infinite differentiability are ensured by the fact that $`\varphi _n`$ enter the right-hand side of Eqs. (25) only as arguments of the sine. Note that because of the arbitrariness of the interval $`I`$, the solution can be uniquely continued onto the whole axis $`\mathrm{}<`$$`y<+\mathrm{}`$.

Differentiating (18) with respect to $`y`$ yields

$$\mathrm{sin}\varphi _n(y)=ϵ^2\underset{m=1}{\overset{N1}{}}G(n,m)\frac{d^2\varphi _m(y)}{dy^2}.$$
(28)
Multiplying (28) by $`\frac{d\varphi _n(y)}{dy}`$, summing over the layer index $`n`$ with the use of (A10) and performing integration, we arrive at the first integral of Eqs. (25):

$$C\left(H\right)\underset{n=1}{\overset{N1}{}}\mathrm{cos}\varphi _n(y)=\frac{ϵ^2}{2}\underset{n=1}{\overset{N1}{}}\underset{m=1}{\overset{N1}{}}G(n,m)\frac{d\varphi _n(y)}{dy}\frac{d\varphi _m(y)}{dy},$$
(29)
$$C\left(H\right)=\frac{2H^2}{H_s^2}\left(N1\right)+\underset{n=1}{\overset{N1}{}}\mathrm{cos}\varphi _n(L),$$
(30)
where $`C\left(H\right)`$ is the constant of integration. Using (29), one can eliminate the electromagnetic-energy term from the Gibbs free energy (22).

The existence of the first integral of Eqs. (25) was first established in Ref. . However, the explicit form of the matrix elements $`G(n,m)`$ and the actual value of $`C\left(H\right)`$ was not determined by the authors of Ref. . Unfortunately, the existence of the first integral for the infinite ($`N=\mathrm{}`$) layered superconductor was not noticed in any previous publications, which partly explains the difficulties with the minimization of the Gibbs free energy.  \[Equation (29) holds for $`N=\mathrm{}`$ too, taking account of the substitution $`\underset{m,n=1}{\overset{N1}{}}G(n,m)\mathrm{}\underset{m,n=\mathrm{}}{\overset{+\mathrm{}}{}}G_{\mathrm{}}(n,m)\mathrm{}`$, where $`G_{\mathrm{}}(n,m)`$ are defined by (A17)\]. In the case $`N=\mathrm{}`$, the Gibbs free energy $`\mathrm{\Omega }\left(H\right)`$ should be additionally minimized with respect to the phases $`\phi _n`$, which by the use of the first integral immediately yields the exact solution with $`\varphi _n(y)=\varphi _{n+1}(y)\varphi (y)`$.

For our further consideration, we will need an auxiliary lemma:

Lemma 2. Consider Eqs. (25) on the whole axis $`\mathrm{}<`$$`y<+\mathrm{}`$. Conditions

$$\varphi _n(y_0)=0,\frac{d\varphi _n}{dy}(y_0)=0,n=1,2,\mathrm{},N1,$$
(31)
where $`y_0\pm \mathrm{}`$, cannot by satisfied simultaneously for all $`n`$ by any nontrivial solution of (25).

The proof of Lemma 2 is straightforward: Indeed, conditions (31) specify the initial value problem for (25). This initial value problem has the trivial solution $`\varphi _1=\varphi _2=\mathrm{}=\varphi _{N1}0`$. By Lemma 1, this solution is unique and there are no nontrivial solutions that satisfy (31).

The importance of Lemma 2 lies in the fact that it prohibits the existence of any topological (vortex) solutions in any finite interval $`[L_1,L_2]`$ in the absence of an external field ($`H=0`$). Indeed, any such solution must necessarily satisfy the boundary conditions $`\varphi _n(L_1)=0`$ and $`\varphi _n(L_2)=0`$mod$`2\pi `$ for all $`n`$, with $`\varphi _n(L_2)0`$ for at least one $`n=m`$. For $`H=0`$, we also have $`\frac{d\varphi _n}{dy}(L_1)=\frac{d\varphi _n}{dy}(L_2)=0`$ for all $`n`$. However, the former and the latter conditions are incompatible by virtue of Lemma 2.

Another important consequence of Lemma 2 is that for any finite interval $`[L,L]`$ the constant of integration in (29) satisfies the inequality

$$C(H)>N1,$$
(32)
for any $`H>0`$. Indeed, for the Meissner solution we have $`\varphi _n(y)=\varphi _n(y)`$ and $`\varphi _n(0)=0`$ \[see (11)\]. By setting $`y=0`$ in (29), applying Lemma 2 and using the fact that the quadratic form on the right-hand side of (29) is positively definite, we get (32). For topological solutions, inequality (32) is quite obvious.

#### 2 Localized solutions. The criterion of existence

Alongside solutions to (25) in a finite interval $`[L,L]`$, we will consider physical solutions to (25) in the infinite ($`\mathrm{}<y<+\mathrm{}`$) and a semiinfinite ($`0y<+\mathrm{}`$) intervals. To ensure the convergence of the integrals over $`y`$ in the expression for the free energy \[see (22)\], such solutions must necessarily satisfy the asymptotic boundary conditions

$$\varphi _n(y)\underset{y\pm \mathrm{}}{}0\text{mod }2\pi ,$$
(33)
$$\frac{d\varphi _n(y)}{dy}\underset{y\pm \mathrm{}}{}0,$$
(34)
for any $`1nN1`$. (Note that Lemma 2 does not obtain for $`y_0=\pm \mathrm{}`$.) These solutions with a square-integrable derivative will be called localized. Their consideration can be reduced to the consideration of the standard initial value problem at $`y=0`$, with the initial data subject to a certain existence condition.

By inserting (33) and (34) into (29), we find the value of the constant of integration: $`C\left(0\right)=N1`$. \[Compare with inequality (32) for $`L<\mathrm{}`$, $`H>0`$.\] Thus Eq. (29) becomes

$$\underset{n=1}{\overset{N1}{}}\mathrm{sin}{}_{}{}^{2}\frac{\varphi _n(y)}{2}=\frac{ϵ^2}{4}\underset{n=1}{\overset{N1}{}}\underset{m=1}{\overset{N1}{}}G(n,m)\frac{d\varphi _n(y)}{dy}\frac{d\varphi _m(y)}{dy}.$$
(35)
Substituting initial values $`\varphi _n(0)=\alpha _n`$, $`\frac{d\varphi _n}{dy}(0)=\beta _n`$ into (35), we obtain the desired condition on $`\alpha _n`$ and $`\beta _n`$ that we call the criterion of the existence of localized solutions:

$$\underset{n=1}{\overset{N1}{}}\mathrm{sin}{}_{}{}^{2}\frac{\alpha _n}{2}=\frac{ϵ^2}{4}\underset{n=1}{\overset{N1}{}}\underset{m=1}{\overset{N1}{}}G(n,m)\beta _n\beta _m.$$
(36)
Indeed, Lemma 1 guarantees the existence, uniqueness and differentiability of a solution for arbitrary $`\alpha _n`$ and $`\beta _n`$ in the whole interval $`(\mathrm{},+\mathrm{})`$ \[or $`[0,+\mathrm{})`$\]. Owing to our special choice of the constant of integration in (35), the solution determined by $`\alpha _n`$ and $`\beta _n`$ obeying (36) will necessarily satisfy asymptotic conditions (33), (34). Moreover, this solution will automatically satisfy the conditions

$$\frac{d^k\varphi _n(y)}{dy^k}\underset{y\pm \mathrm{}}{}0,$$
(37)
for any $`1nN1`$ and $`2k`$, which can be verified by repeated differentiation of (28) and application of (33), (34).

Note that localized solutions are characterized by the absence of any electromagnetic effects at $`y\pm \mathrm{}`$. Furthermore, owing to (35), the Josephson energy of localized solutions $`E_J`$ is exactly equal to their electromagnetic energy $`E_{em}`$. (From a field-theoretical point of view, this fact is a manifestation of the virial theorem.)

#### 3 Solutions with periodic derivatives. The major Theorem

Our conclusions about the type and the properties of physical solutions for $`H>0`$ will be based on the major Theorem:

Theorem. Consider Eqs. (25) on the whole axis $`\mathrm{}<`$$`y<+\mathrm{}`$. The initial value problem

$$\varphi _n(y_0)=0,\frac{d\varphi _n}{dy}(y_0)=2H>0,n=1,2,\mathrm{},N1,$$
(38)
for (25) has a unique solution on the whole axis $`\mathrm{}<`$$`y<+\mathrm{}`$. This solution is characterized by the properties

$$\varphi _n\left(y+P\right)=\varphi _n\left(y\right)+2\pi ,$$
(39)
$$\frac{d\varphi _n}{dy}\left(y+P\right)=\varphi _n\left(y\right)>0,$$
(40)
$$\varphi _n\left(y_0+\frac{1}{2}P\right)=\pi ,$$
(41)
$$\frac{d\varphi _n}{dy}\left(y_0+\frac{1}{2}P\right)=2\sqrt{H^2+H_s^2},\text{for all }n=1,2,\mathrm{},N1,$$
(42)
where $`H_s`$ is defined via (A14), and the period $`P`$ is given by

$$P(H)=\frac{2}{\sqrt{H^2+H_s^2}}\underset{(0,0,\mathrm{},0)}{\overset{(\frac{\pi }{2},\frac{\pi }{2},\mathrm{},\frac{\pi }{2})}{}}\frac{ds}{\sqrt{1\frac{1}{N1}\frac{H_s^2}{H^2+H_s^2}\underset{n=1}{\overset{N1}{}}\mathrm{sin}{}_{}{}^{2}\theta _{n}^{}}},$$
(43)
$`ds=ϵH_s\sqrt{{\displaystyle \frac{1}{N1}}{\displaystyle \underset{n=1}{\overset{N1}{}}}{\displaystyle \underset{m=1}{\overset{N1}{}}}G(n,m)d\theta _nd\theta _m}.`$

The integration in (43) is performed along the integral curve:

$`\theta _n={\displaystyle \frac{\varphi _n(y)}{2}},y_0yy_0+{\displaystyle \frac{1}{2}}P,n=1,2,\mathrm{},N1,`$

with $`ds`$ being the length increment. (See Appendix C for the proof.)

An immediate physical consequence of the Theorem is the absence of single-vortex solutions to (25) in a finite interval $`[L,L]`$ for any $`H>0`$. (For $`H=0`$, the absence of single-vortex solutions in a finite interval follows from Lemma 2.) Indeed, if such solutions existed, they would satisfy at a certain $`H=H_1>0`$ the boundary conditions

$`\varphi _n(L)=0,{\displaystyle \frac{d\varphi _n}{dy}}(L)=2H_1>0,n=1,2,\mathrm{},N1,`$

$`\varphi _l(L)=2\pi ,{\displaystyle \frac{d\varphi _l}{dy}}(L)=2H_1>0;\varphi _n(L)=0,{\displaystyle \frac{d\varphi _n}{dy}}(L)=2H_1>0,nl.`$

However, these boundary conditions are mutually incompatible: The boundary conditions at $`y=L`$ specify the initial value problem for (25). By the Theorem, this initial value problem admits a unique solution with $`\frac{d\varphi _n}{dy}>0`$ on the whole axis $`\mathrm{}<`$$`y<+\mathrm{}`$ for all $`n`$, whereas the boundary conditions at $`y=L`$ imply that $`\frac{d\varphi _n}{dy}`$ with $`nl`$ change the sign twice in the interval $`[L,L]`$.

Analogously, one can prove the absence for $`H>0`$ and $`L<\mathrm{}`$ of any other incoherent vortex solutions, as well as such vortex configurations as the ”triangular lattice” and the ”row mode”. In contrast, the properties (39)-(43) are inherent to the vortex-plane solutions.

Note that $`P(H)`$ $`\underset{H0}{}\mathrm{}`$, in agreement with Lemma 2. On the other hand, for $`HH_s`$, the path of integration in (43) is a straight line: $`\theta _n=\frac{\pi }{P}\left(yy_0\right)`$, $`y_0yy_0+\frac{1}{2}P`$, $`n=1,2,\mathrm{},N1`$, which yields $`P(H)`$ $`\underset{HH_s}{}`$ $`\frac{\pi }{H}`$.

For a single thin-layer junction ($`N=2`$), equations (25) reduce to a single SG equation for $`\varphi _1\varphi `$, with $`\lambda _{J0}=H_s^1\frac{ϵ}{\sqrt{2+ϵ^2}}`$. In agreement with Refs.  , the period (43) becomes

$$P(H)=\frac{2}{\sqrt{H^2+H_s^2}}K\left(\frac{H_s^2}{H^2+H_s^2}\right),$$
(44)
where $`K\left(k^2\right)`$ is the complete elliptic integral of the first kind. Relations (39)-(43), applied to $`\varphi `$, determine an exact vortex solution in terms of the Jakobi elliptic functions am$`(u)`$ and dn$`(u)=\frac{d}{du}`$am$`(u)`$.

The true vortex-plane solutions appear in a double-junction stack ($`N=3`$). Owing to the symmetry relations (21), equations (25) again reduce to a single SG equation for $`\varphi _1=\varphi _2\varphi `$, with $`\lambda _{J0}=H_s^1\frac{ϵ}{\sqrt{1+ϵ^2}}`$. The period $`P`$ is again given by (44) (with redefined $`H_s`$), and the exact vortex-plane solution is again expressed via the functions am$`(u)`$ and dn$`(u)`$. These two exactly-solvable examples ($`N=2`$ and $`N=3`$) not only demonstrate the power and generality of the Theorem, but also show that the vortex-plane solutions are a natural generalization of the well-known vortex solutions in a single junction. (See section IV for a more detailed discussion of the cases $`N=2,3`$.)

### B Meissner solutions. The superheating (penetration) fields $`H_{sL}`$ and $`H_sH_s\mathrm{}`$

In a finite interval $`[L,L]`$, the Meissner solution is characterized by the relations

$$\varphi _n(y)=\varphi _n(y),1nN1,$$
(45)
resulting from the general symmetry (11), and obeys the conditions (21). Thus, the Meissner boundary value problem is completely specified by the boundary conditions

$$\frac{d\varphi _n}{dy}\left(L\right)=2H>0,\varphi _n\left(0\right)=0,1nN1.$$
(46)
Up to a certain field $`H=H_{sL}`$, the boundary value problem (46) has a unique solution with $`\frac{d\varphi _n}{dy}>0`$ in the whole interval $`[L,L]`$ and $`\pi \varphi _n\left(L\right)<0`$ for all $`n`$. The field $`H_{sL}`$ is determined from the conditions $`\varphi _n\left(\pm L\right)=\pm \pi `$ for all $`n`$. \[As follows from (30), (32), $`H_{sL}>H_s`$.\] Indeed, for $`HH_s`$ the Meissner solution can be obtained in a closed analytical form as a linear combination of $`\left[\frac{N}{2}\right]`$ exponentials with $`\lambda _{Jk}`$ \[see (24)\]: According to (30), (32), this case allows of linearization. The existence of a unique Meissner solution for higher values of $`H`$ follows from continuous dependence of the solution on the initial data $`\varphi _n\left(L\right)=\alpha _n`$ and $`\frac{d\varphi _n}{dy}\left(L\right)=2H\beta _n`$. (See Lemma 1.) The Meissner solution ceases to exist at $`H=H_{sL}`$: For this field, both $`\frac{d\varphi _n}{dy}(y)`$ and $`h_n(y)`$ are local maxima at $`y=\pm L`$. \[See (19), (28).\]

Physically, the conditions $`\varphi _n\left(\pm L\right)=\pm \pi `$ correspond to the vanishing of the surface barrier induced by Josephson currents \[$`j_{n,n1}\left(\pm L\right)=0`$\] and the formation of a vortex plane at the side boundaries of the stack. (Thus, at $`H=H_{sL}`$, the total flux due to the field penetration through the $`y=\pm L`$ interfaces is exactly equal to the flux carried by a single vortex plane.) For $`H>H_{sL}`$, only topological vortex-plane solutions are possible. However, the vortex-plane solutions can exist and be favorable energetically already at fields $`H<H_{sL}`$. \[See the next subsection.\] For these reasons, $`H_{sL}`$ should be identified both with the superheating field of the Meissner state and the penetration field for the vortex planes. According to the Theorem, at $`H=H_{sL}`$, $`\frac{d\varphi _n}{dy}(0)=2\sqrt{H_{sL}^2H_s^2}`$ for all $`n`$, and hence $`H_n(0)=\sqrt{H_{sL}^2H_s^2}+\left(H_{sL}\sqrt{H_{sL}^2H_s^2}\right)\left[G(n,1)+G(n,N1)\right]`$. The field $`H_{sL}`$ itself is determined by the implicit equation

$$2L=P\left(\sqrt{H_{sL}^2H_s^2}\right),$$
(47)
where $`P`$ is given by (43), with the path of integration being determined by the solution to the boundary value problem $`\varphi _n\left(L\right)=\pi `$, $`\varphi _n\left(0\right)=0`$. Note that $`H_{sL}H_s`$ for $`LH_s^1`$, which will be explained below.

The consideration of the Meissner effect becomes simpler, when $`L\lambda _{J\mathrm{max}}\lambda _{J0}`$. In this case, the interval $`[L,L]`$ can be transformed into $`[0,+\mathrm{})`$ by changing the variable $`yyL`$ and proceeding to the limit $`L\mathrm{}`$. In the semiinfinite interval $`[0,+\mathrm{})`$, the Meissner solution necessarily satisfies the conditions $`\varphi _n(y)\underset{y+\mathrm{}}{}0`$, $`\frac{d\varphi _n(y)}{dy}\underset{y+\mathrm{}}{}0`$, and we can use the criterion (36) that now takes the form

$$\frac{1}{N1}\underset{n=1}{\overset{N1}{}}\mathrm{sin}{}_{}{}^{2}\frac{\varphi _n(0)}{2}=\frac{H^2}{H_s^2},$$
(48)
where, by (A14) and (A15), (A16),

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

$$\sqrt{\frac{\left(N1\right)N}{2}}\left[\underset{k=0}{\overset{\left[\frac{N}{2}\right]1}{}}\lambda _{Jk}^2\mathrm{cot}{}_{}{}^{2}\frac{\pi \left(2k+1\right)}{2N}\right]^{\frac{1}{2}}.$$
(49)
Physical interpretation of the quantity $`H_s`$ is straightforward. The maximum value of the left-hand side of (48) is achieved when all $`\varphi _n(0)=\pi `$ and is equal to unity, which corresponds to $`H=H_s`$ on the right-hand side. Hence $`H_s`$ should be identified with the superheating (penetration) field $`H_s\mathrm{}`$ ($`H_sH_s\mathrm{}`$). Note that $`H_s`$ for $`N<\mathrm{}`$ is always higher than the corresponding field $`\overline{H}_s=1`$ of the infinite ($`N=\mathrm{}`$) layered superconductor. For $`N\mathrm{}`$, $`H_s1`$.

According to (48), equations (25) can be linearized for $`HH_s`$. In this limit, the explicit Meissner solution in the interval $`[0,+\mathrm{})`$ is

$$\varphi _n(y)=\frac{4H}{N}\underset{k=0}{\overset{\left[\frac{N}{2}\right]1}{}}\lambda _{Jk}\mathrm{cot}\frac{\left(2k+1\right)\pi }{2N}\mathrm{sin}\frac{\left(2k+1\right)n\pi }{N}\mathrm{exp}\left[\frac{y}{\lambda _{Jk}}\right],$$
(50)
$$H_n(y)=\frac{2H}{N}\underset{k=0}{\overset{\left[\frac{N}{2}\right]1}{}}\lambda _{Jk}^2\mathrm{cot}\frac{\left(2k+1\right)\pi }{2N}\mathrm{sin}\frac{\left(2k+1\right)n\pi }{N}\mathrm{exp}\left[\frac{y}{\lambda _{Jk}}\right]+H_n,$$
(51)
where $`H_n`$ is given by (20). The distribution of the currents $`J_n`$ and $`j_{n,n1}`$ can be easily obtained from (13), (14) using (50), (51).

It is instructive to investigate a transition to the infinite layered-superconductor limit $`N12\left[ϵ^1\right]`$ ($`ϵ<1`$) in Eqs. (50), (51). Thus, for I-layers whose layer index $`n`$ satisfies the condition

$$\left[ϵ^1\right]nN1\left[ϵ^1\right],$$
(52)
we have $`H_n=H\left(\mu ^n+\mu ^{Nn}\right)H`$. Under the condition (52), the main contribution to the sums over $`k`$ in (50), (51) comes from $`k\left[\frac{N}{2}\right]1`$. For such $`k`$ and $`N1`$, $`\lambda _{Jk}\lambda _{J0}\overline{\lambda }_J=1`$. Therefore, for $`n`$ in the interval (52) we get

$`\varphi _n(y){\displaystyle \frac{8H}{\pi }}\mathrm{exp}\left(y\right){\displaystyle \underset{k=0}{\overset{+\mathrm{}}{}}}{\displaystyle \frac{\left(1\right)^k}{2k+1}}=2H\mathrm{exp}\left(y\right),`$

$`H_n(y)H\mathrm{exp}\left(y\right),`$

in complete agreement with the results of Refs.  .

Two final remarks would be in order here. The fact that the Meissner solution in a stack with $`N`$ S-layers is characterized by $`\left[\frac{N}{2}\right]`$ Josephson lengths $`\lambda _{Jk}`$ was not noticed in any previous publications. Moreover, in contrast to continuum type-II superconductors,  the local field $`H_n(y)`$ at $`H=H_s`$ cannot be represented as a sum of a purely Meissner field and a vortex field, because the principle of superposition does not hold for the nonlinear Eqs. (25). Unfortunately, the latter point was not realized in Ref.  concerned with infinite layered superconductors.

### C Vortex-plane solutions. The lower critical field $`H_{c1}`$

As explained in the Introduction, by a vortex plane we understand a coherent chain of $`N1`$ Josephson vortices (one vortex per each I-layer) positioned in a plane parallel to the coordinate axes $`x`$, $`z`$: see Fig. 2. (In this plane, all $`\frac{d\varphi _n}{dy}`$ and $`h_n`$ are local maxima.) Such solutions in an interval $`[L,L]`$ obey the general relations (11) with the same constant $`0`$mod$`2\pi 0`$ for all $`n`$ and the symmetry conditions (21): hence the same set of $`\left[\frac{N}{2}\right]`$ Josephson lengths $`\lambda _{Jk}`$ as in the Meissner case. The Theorem envisages the existence of vortex-plane solutions for any $`H>0`$. Using the Theorem and Lemma 1, we can easily establish all basic properties of the vortex-plane solutions. They can be summarized as follows.

The boundary value problem for the vortex-plane solutions is specified by the boundary conditions

$$\frac{d\varphi _n}{dy}\left(L\right)=2H>0,\varphi _n\left(0\right)=\pi N_v,1nN1,N_v=1,2,\mathrm{},$$
(53)
where $`N_v`$ is the number of vortex planes. As in the case of the Meissner boundary value problem (46), the boundary value problem (53) for a fixed $`N_v`$ and any $`L>0`$ has a unique solution, with $`\frac{d\varphi _n}{dy}>0`$ in the whole interval $`[L,L]`$ and $`\pi \varphi _n\left(L\right)0`$ for all $`n`$, in a certain field range

$$\sqrt{H_{N_v1}^2H_s^2}HH_{N_v},$$
(54)
where $`H_0H_{sL}`$. The lower bound of the existence of the $`N_v`$-vortex-plane solution is determined by the boundary value problem

$$\varphi _n\left(L\right)=0,\varphi _n\left(0\right)=\pi N_v,1nN1,$$
(55)
whereas the upper bound is determined by the boundary value problem

$$\varphi _n\left(L\right)=\pi ,\varphi _n\left(0\right)=\pi N_v,1nN1.$$
(56)
The field $`H_{N_v}`$ is determined by the implicit equation

$$2L=\left(N_v+1\right)P\left(\sqrt{H_{N_v}^2H_s^2}\right),$$
(57)
where $`P`$ is given by (43), with the path of integration being determined by the solution to the boundary value problem (56). Note that for $`N_v=0`$ (the Meissner state) Eq. (57) coincides with Eq. (47), as it should. Physically, $`P`$ determines the period of the distribution of $`H_n(y)`$ and $`j_{n,n1}(y)`$.

As follows from (54), the range of the existence of the state with $`N_v`$ vortex planes overlaps with that of the state with $`N_v1`$ vortex planes. Thus, the state with $`N_v=1`$ overlaps with the Meissner state. For $`LH_s^1`$, the overlapping can occur for several neighboring states. The overlapping is negligibly small only for $`LH_s^1`$ (with arbitrary $`N_v`$) and for $`N_v1`$ (with arbitrary $`L`$). (For ordinary Josephson vortices in a single junction, the overlapping is well known.) However, the actual number of vortex planes, $`N_v`$, for given $`H`$ should be determined from the requirement that the free energy be an absolute minimum. Note that for $`L\mathrm{}`$, $`\sqrt{H_0^2H_s^2}\sqrt{H_{sL}^2H_s^2}0`$ and the state with $`N_v=1`$ appears already in zero external field. This case is discussed in more detail below.

For $`LH_s^1`$, we can proceed to the limit $`L\mathrm{}`$ and consider the state with $`N_v=1`$ in the interval $`(\mathrm{},+\mathrm{})`$. With the asymptotic boundary conditions

$$\varphi _n(y)\underset{y\mathrm{}}{}0,\varphi _n(y)\underset{y+\mathrm{}}{}2\pi ,\frac{d\varphi _n(y)}{dy}\underset{y\pm \mathrm{}}{}0$$
(58)
for all $`n`$, this state satisfies the criterion (36) and can be obtained as the solution to an initial value problem at $`y=0`$. Indeed, it can be constructed from the Meissner solution in the interval $`[0,+\mathrm{})`$ at $`H=H_s`$: Using the property (26), by means of the transformation $`\varphi _n\varphi _n+2\pi `$ we obtain in the interval $`[0,+\mathrm{})`$ a solution that satisfies the initial conditions

$$\alpha _n\varphi _n(0)=\pi ,\beta _n\frac{d\varphi _n}{dy}(0)=2H_s,1nN1,$$
(59)
and the conditions (58) for $`y+\mathrm{}`$. By Lemma 1, this solutions can be uniquely continued into the whole interval $`(\mathrm{},+\mathrm{})`$. The solution thus obtained satisfies the conditions (58) for $`y\pm \mathrm{}`$ and hence represents the desired singe-vortex-plane solution. By the construction, $`h_n(0)=H_s\left[1G(n,1)G(n,N1)\right]`$. The total flux carried by the vortex plane is

$$\mathrm{\Phi }=\underset{n=1}{\overset{N1}{}}\mathrm{\Phi }_n=\pi \left(N1\right)\left[1\frac{2\sqrt{1+\frac{ϵ^2}{4}}ϵ}{ϵ\left(N1\right)}\frac{1\mu ^{N1}}{1+\mu ^N}\right],$$
(60)
where $`\mathrm{\Phi }_n`$ is given by (23). Note that in contrast to Josephson junctions with thick electrodes and infinite layered superconductors, the flux carried by a Josephson vortex in a finite thin-layer S/I structure is not quantized and is always smaller than the flux quantum $`\mathrm{\Phi }_0=\pi `$. For $`N12\left[ϵ^1\right]`$, $`\mathrm{\Phi }\pi \left(N1\right)`$, as it should.

To determine the thermodynamic lower critical field $`H_{c1}`$ at which the vortex-plane solutions for $`LH_s^1`$ become energetically favorable, we must calculate the difference between the Gibbs free energy in the presence of a single vortex plane, $`\mathrm{\Omega }_v(H)`$, and the Gibbs free energy of the homogeneous Meissner state, $`\mathrm{\Omega }_M(H)`$ \[the sum of the phase-independent terms in (22)\]:

$`\mathrm{\Omega }_v(H)\mathrm{\Omega }_M(H)`$

$$=r(T)\left[ϵ^2\underset{n=1}{\overset{N1}{}}\underset{m=1}{\overset{N1}{}}G(n,m)\underset{\mathrm{}}{\overset{+\mathrm{}}{}}𝑑y\frac{d\varphi _n(y)}{dy}\frac{d\varphi _m(y)}{dy}4\mathrm{\Phi }H\right],$$
(61)
where the total flux $`\mathrm{\Phi }`$ is given by (60). The first term on the right-hand side of (61) should be interpreted as the self-energy of the vortex plane: $`E_v=2E_{em}=2E_J`$. From (61), we get:

$$H_{c1}=\frac{E_v}{4r(T)\mathrm{\Phi }}.$$
(62)
For thermodynamically stable solutions, we must necessarily have $`H_{c1}<H_s`$. It is straightforward to verify that this condition is met by the vortex-plane solutions. Using (28), (58), the initial values $`\beta _n=2H_s`$ and integrating by parts, we convert $`E_v`$ into the form

$$E_v=2r(T)\left[\underset{n=1}{\overset{N1}{}}\underset{+\mathrm{}}{\overset{0}{}}𝑑y\varphi _n(y)\mathrm{sin}\varphi _n(y)2H_s\mathrm{\Phi }\right].$$
(63)
The first term on the right-hand side of (63) is positive, because in the region $`0y<+\mathrm{}`$ all $`\varphi _n`$ satisfy the relation $`\pi \varphi _n<2\pi `$. By the use of (28) and (A14), we obtain the following strict inequalities:

$`2H_s\mathrm{\Phi }<{\displaystyle \underset{n=1}{\overset{N1}{}}}{\displaystyle \underset{+\mathrm{}}{\overset{0}{}}}𝑑y\varphi _n(y)\mathrm{sin}\varphi _n(y)<4H_s\mathrm{\Phi },`$

$`0<E_v<4r(T)H_s\mathrm{\Phi }.`$

Hence,

$`0<H_{c1}<H_s,`$

as anticipated. Note that in all special cases admitting exact analytical solutions ($`N=\mathrm{}`$, and $`N=2,3`$, see section IV), $`H_{c1}=\frac{2}{\pi }H_s`$.

### D Single-vortex solutions for $`H=0`$, $`L=\mathrm{}`$, and other localized incoherent vortex solutions

As shown in section IIIA, the only topological (vortex) solutions admitted by Eqs. (25) for $`H>0`$ are the vortex-plane solutions. For $`H=0`$, equations (25) do not possess any topological solutions in a finite interval $`[L,L]`$. Therefore, incoherent vortex solutions to (25) can exist only for $`H=0`$, in the infinite interval $`(\mathrm{},+\mathrm{})`$, and must necessarily meet the requirements imposed on localized solutions. Below we consider in detail the most important type of such solutions, namely single-vortex solutions.

A single Josephson vortex positioned in the $`l`$th I-layer at $`y=0`$ obeys the symmetry relations

$$\varphi _l(y)=2\pi \varphi _l(y);\varphi _n(y)=\varphi _n(y),nl,$$
(64)
\[see (11)\] and the asymptotic boundary conditions

$$\varphi _l(y)\underset{y\mathrm{}}{}0,\varphi _l(y)\underset{y+\mathrm{}}{}2\pi ;\varphi _n(y)\underset{y\pm \mathrm{}}{}0,nl,$$
(65)
$$\frac{d\varphi _n(y)}{dy}\underset{y\pm \mathrm{}}{}0,\text{for all }1nN1.$$
(66)
Moreover, $`\frac{dh_n(y)}{dy}>0`$ in the region $`\mathrm{}<y<0`$, and $`\frac{dh_n(y)}{dy}<0`$ in the region $`0<y<+\mathrm{}`$. Hence, $`\varphi _n`$ must satisfy the relations

$$0<\varphi _n(y)<\pi ,\text{ }y(\mathrm{},0);\pi <\varphi _n(y)<0,\text{ }y(0,+\mathrm{}),\text{for }nl,$$
(67)
and the initial conditions

$$\alpha _l\varphi _l(0)=\pi ;\alpha _n\varphi _n(0)=0,nl,$$
(68)
$$\beta _l\frac{d\varphi _l(0)}{dy}>0;\beta _n\frac{d\varphi _n(0)}{dy}<0,nl.$$
(69)
A necessary condition of the existence of such a solution is provided by the general criterion (36) and has the form

$$\frac{ϵ^2}{4}\underset{n=1}{\overset{N1}{}}\underset{m=1}{\overset{N1}{}}G(n,m)\beta _n\beta _m=1.$$
(70)
The flux through the $`n`$th I-layer due to the vortex in the $`l`$th I-layer can be found using (19) and (65):

$$\mathrm{\Phi }_{nl}=\underset{\mathrm{}}{\overset{+\mathrm{}}{}}𝑑yh_n(y)=\frac{\pi ϵ}{2\sqrt{1+\frac{ϵ^2}{4}}}\left[\mu ^{\left|nl\right|}\frac{\mu ^n\left(\mu ^{lN}\mu ^{Nl}\right)+\mu ^{Nn}\left(\mu ^l\mu ^l\right)}{\mu ^N\mu ^N}\right].$$
(71)
The total flux carried by the vortex in the $`l`$th layer is $`\mathrm{\Phi }_l=_{n=1}^{N1}\mathrm{\Phi }_{nl}`$, where $`\mathrm{\Phi }_l`$ is given by (23) with $`n=l`$. For $`N<\mathrm{}`$, the total flux $`\mathrm{\Phi }_l`$ is not quantized and is less than the flux quantum $`\mathrm{\Phi }_0=\pi `$. (See the previous subsection.)

In contrast to the Meissner and the vortex-plane solutions, single-vortex solutions, in general, do not obey the symmetry (21), inherent to the original integrodifferential equations (7), (17) minimizing the Gibbs free-energy. Therefore, they are characterized by the full set (A15) of $`N1`$ length scales $`\lambda _i`$. (The only exclusion is a stack with an odd number of junctions $`N1`$ and $`l=\frac{N}{2}`$.)

Although the energy of single-vortex solutions is lower than the self-energy of the vortex-plane solutions at $`H=0`$, the former solutions cannot be used, e.g., for estimates of $`H_{c1}`$ because of their absolute instability at $`H>0`$. Unfortunately, this crucial issue was not understood in any previous publications. Thus, the main distinctive feature of single-vortex solutions that precludes their existence at $`H>0`$, i.e., $`\frac{d\varphi _n(0)}{dy}<0`$ for $`nl`$, was overlooked in Refs.  concerned with infinite layered superconductors. (The approach of Refs.  was criticized from different points of view in Refs.  .) The negative sign of $`\frac{d\varphi _n(0)}{dy}`$ for $`nl`$ was noticed in Ref.  . However, the effect of this property on the stability of single-vortex solutions at $`H>0`$ was not realized therein.

The feature (69) is clearly reproduced by the numerical studies of finite stacks. Unfortunately, in the absence of any analytical investigations of Eqs. (25), the authors of Refs.  failed to establish the actual domain of validity of their numerical results: $`H=0`$, $`L=\mathrm{}`$. The condition (70) was not found either. In this regard, it would be appropriate to point out an intrinsic limitation of the numerical calculations that does not allow one to consider them as a proof of the existence of single-vortex solutions even in the case $`H=0.`$ By necessity, these calculations are performed in a finite interval $`[L_1,L_2]`$. However, as shown in section IIIA, equations (25) do not admit any topological solutions (in a strict mathematical sense) in any finite interval at $`H=0`$. It should be also noted that the single-vortex solution cannot be represented by a linear combination of the terms $`4\mathrm{arctan}\mathrm{exp}\left[\frac{y}{\lambda _i}\right]`$, because the principle of superposition does not apply to the nonlinear Eqs. (25).

The consideration of other localized incoherent vortex solutions to (25) (i.e., solutions with $`2k<N1`$ vortices in the plane $`y=0`$, and vortex-antivortex configurations) can be done along the same lines. At $`H>0`$, all these solutions are unstable too.

## IV The two exactly-solvable cases

Below, we present an exact solution for the cases $`N=2,3`$. Valid for any $`H0`$ and $`L>0`$, this solution provides a clear illustration of the general results of sections IIIA-IIIC. The single-junction case ($`N=2`$) illuminates some typical specific features of the thin-S-layer limit and establishes a profound physical and mathematical relationship between the ordinary Josephson vortices and the vortex-plane solutions. In the double-junction case ($`N=3`$), our exact solution describes the true vortex planes, which is of considerable experimental interest: For a double-junction stack, we already have direct observations of Josephson-vortex configurations that can be unambiguously identified with the vortex-plane solutions.

### A A single thin-layer junction ($`N=2`$)

In this simplest case, a single phase difference $`\varphi _1(y)\varphi (y)`$ satisfies the usual static SG equation

$$\frac{d^2\varphi (y)}{dy^2}=\frac{1}{\lambda _{J0}^2}\mathrm{sin}\varphi (y),$$
(72)
with the Josephson length

$$\lambda _{J0}=\frac{ϵ}{\sqrt{2+ϵ^2}}.$$
(73)
Note that $`\lambda _{J0}`$, given by (73), for $`ϵ1`$ is much smaller than the Josephson length of a single junction with thick electrodes, which in our dimensionless units is $`\lambda _J=\sqrt{\frac{p}{2\lambda }}`$. From (20) and (49), we get the local field in the homogeneous Meissner state

$$H_1=\frac{2H}{2+ϵ^2},$$
(74)
and the superheating (penetration) field

$$H_s=\lambda _{J0}^1=\frac{\sqrt{2+ϵ^2}}{ϵ},$$
(75)
respectively. For $`ϵ1`$, the superheating (penetration) field $`H_s`$, given by (75), is much higher than the corresponding field of a single junction with thick electrodes $`H_s=\lambda _J`$.

Using the Theorem, we immediately find an exact solution to the boundary value problems (46) and (53):

$$\varphi (y)=\pi \left(N_v1\right)+2\text{am}(\frac{y}{k\lambda _{J0}}+K\left(k^2\right),k),$$
(76)
$$\text{dn}(\frac{L}{k\lambda _{J0}},k)=\frac{\sqrt{1k^2}}{k\lambda _{J0}H},N_v=2m,m=0,1,\mathrm{};$$
(77)
$$\varphi (y)=\pi N_v+2\text{am}(\frac{y}{k\lambda _{J0}},k),$$
(78)
$$\text{dn}(\frac{L}{k\lambda _{J0}},k)=k\lambda _{J0}H,N_v=2m+1,m=0,1,\mathrm{},$$
(79)
where am$`\left(u\right)`$ and dn$`\left(u\right)=\frac{d}{du}`$am$`\left(u\right)`$ are the Jakobi elliptic functions, and $`K\left(k^2\right)`$ is the elliptic integral of the first kind. The range of the existence of the solution with $`N_v=0,1,\mathrm{}`$ Josephson vortices is given by (54), where $`H_{N_v}`$ is determined by the implicit equation

$$\frac{L}{\lambda _{J0}}=\frac{\left(N_v+1\right)}{H_{N_v}\lambda _{J0}}K\left(\frac{1}{H_{N_v}^2\lambda _{J0}^2}\right).$$
(80)
Note that although the Meissner and the vortex solutions to (72) were studied a long time ago, a closed-form analytical solution of the type (76)-(80), valid for any $`0<L`$ in the whole field range $`0HH_{c2}`$, has not been found up until now. (Apparently for this reason, there exists an erroneous belief that Josephson vortices ”do not form” in single junctions with small $`L`$.) Using asymptotics of am$`\left(u\right)`$, dn$`\left(u\right)`$ and $`K\left(k^2\right)`$, we can obtain from (76)-(80) all physically interesting limiting cases.

#### 1 The Meissner solution for $`[0y<+\mathrm{})`$

For the fields $`0HH_s=\frac{\sqrt{2+ϵ^2}}{ϵ}`$, the Meissner solution in the semiinfinite interval $`[0,+\mathrm{})`$ can be obtained from (76), (77) with $`N_v=0`$ by changing the variable $`yyL`$ and proceeding to the limit $`L\mathrm{}`$. Using the formulas of section II, up to first order in $`r(T)1`$, we have

$$\varphi (y)=4\mathrm{arctan}\frac{H\mathrm{exp}\left[\frac{y}{\lambda _{J0}}\right]}{H_s+\sqrt{H_s^2H^2}},$$
(81)
$$H(y)H_1(y)=h(y)+H_1,$$
(82)
$$h(y)h_1(y)=\frac{2\lambda _{J0}H\left[H_s+\sqrt{H_s^2H^2}\right]\mathrm{exp}\left[\frac{y}{\lambda _{J0}}\right]}{\left[H_s+\sqrt{H_s^2H^2}\right]^2+H^2\mathrm{exp}\left[\frac{2y}{\lambda _{J0}}\right]},$$
(83)
$`j(y)j_{1,0}(y)=4H\left[H_s+\sqrt{H_s^2H^2}\right]`$

$$\times \frac{\left[\left[H_s+\sqrt{H_s^2H^2}\right]^2H^2\mathrm{exp}\left[\frac{2y}{\lambda _{J0}}\right]\right]\mathrm{exp}\left[\frac{y}{\lambda _{J0}}\right]}{\left[\left[H_s+\sqrt{H_s^2H^2}\right]^2+H^2\mathrm{exp}\left[\frac{2y}{\lambda _{J0}}\right]\right]^2},$$
(84)
$$J(y)J_0(y)=J_1(y)=\frac{1}{4\pi }\left[HH_1h(y)\right],$$
(85)
$$f(y)f_0(y)=f_1(y)=1\frac{r(T)}{2}\left[\lambda _{J0}^2h(y)+\frac{2}{ϵ^2}\left[h(y)+H_1H\right]^2\right].$$
(86)
#### 2 The single-vortex solution for $`(\mathrm{}<y<+\mathrm{})`$

The single-vortex solution in the infinite interval $`(\mathrm{},+\mathrm{})`$ can be obtained from (78), (79) with $`N_v=1`$ by proceeding to the limit $`L\mathrm{}`$. It has the form

$$\varphi (y)=4\mathrm{arctan}\mathrm{exp}\left[\frac{y}{\lambda _{J0}}\right],$$
(87)
$$h(y)=\lambda _{J0}\mathrm{cosh}{}_{}{}^{1}\left[\frac{y}{\lambda _{J0}}\right],$$
(88)
$$j(y)=2\mathrm{cosh}{}_{}{}^{2}\left[\frac{y}{\lambda _{J0}}\right]\mathrm{sinh}\left[\frac{y}{\lambda _{J0}}\right].$$
(89)
The quantities $`H(y)`$, $`J(y)`$ and $`f(y)`$ are given by (82), (85) and (86), respectively, with $`h(y)`$ taken from (88).

By inserting (87) into (61) with $`H=0`$, we obtain the vortex self-energy:

$$E_v=8r(T)\frac{ϵ}{\sqrt{2+ϵ^2}}.$$
(90)
The vortex flux, according to (60), is

$`\mathrm{\Phi }=\pi {\displaystyle \frac{ϵ^2}{2+ϵ^2}},`$

and the lower critical field, by (62), is

$$H_{c1}=\frac{2}{\pi }H_s=\frac{2}{\pi }\frac{\sqrt{2+ϵ^2}}{ϵ}.$$
(91)
Thus, for $`ϵ1`$, the vortex flux $`\mathrm{\Phi }\mathrm{\Phi }_0=\pi `$, and the lower critical field (91) is much larger than the corresponding field of a single junction with thick electrodes $`H_{c1}=\frac{2}{\pi }\sqrt{\frac{p}{2\lambda }}`$, in agreement with Ref. .

#### 3 The vortex solution for small Josephson screening

When the screening by Josephson currents is negligibly small, i.e., (i) for $`L\lambda _{J0}`$ and arbitrary $`H`$, or (ii) for $`H_sHH_{c2}`$ and arbitrary $`L`$, equations (76)-(80) become

$`\varphi (y)=\pi N_v+2Hy`$

$$\frac{\left(1\right)^{N_v}}{4\lambda _{J0}^2H^2}\left[\mathrm{sin}\left(2Hy\right)2Hy\mathrm{cos}\left(HW\right)\right],$$
(92)
where $`N_v=\left[\frac{HW}{\pi }\right]`$. Equations of the type (92) were first obtained for the infinite layered superconductor by means of a perturbation theory. As can be seen from (92), the overlapping of states with different $`N_v`$ now vanishes, and the period of the vortex structure is given by $`P=\frac{\pi }{H}`$, in full agreement with the general results of sections IIIA and IIIC.

### B A double-junction stack ($`N=3`$)

Owing to the symmetry (21), $`\varphi _1(y)=\varphi _2(y)\varphi (y)`$, and a double-junction stack is described by the single SG equation (72) with the single Josephson length

$$\lambda _{J0}=\frac{ϵ}{\sqrt{1+ϵ^2}}.$$
(93)
An exact solution to the boundary value problems (46) and (53) is again given by Eqs. (76)-(80), where now $`N_v`$ should be interpreted as the number of vortex planes. Thus, practically all the formulas of the previous subsection (with corresponding redefinition of the physical quantities) hold for the double-junction stack too. Note that

$`\varphi _1(y)=\varphi _2(y)\varphi (y),H_1(y)=H_2(y)H(y),h_1(y)=h_2(y)h(y),`$

$`j_{1,0}(y)=j_{2,1}(y)j(y),J_0(y)=J_2(y)J(y),f_0(y)=f_2(y)f(y),`$

$`J_1(y)=0,f_1(y)=1{\displaystyle \frac{r(T)}{\lambda _{J0}^2}}h(y),H_1=H_2={\displaystyle \frac{H}{1+ϵ^2}}.`$

According to (49), the superheating (penetration ) field is

$`H_s=\lambda _{J0}^1={\displaystyle \frac{\sqrt{1+ϵ^2}}{ϵ}},`$

which is smaller than the corresponding single-junction value (75), in agreement with Ref. .

The vortex-plane self-energy for the interval $`(\mathrm{},+\mathrm{})`$ is

$$E_v=16r(T)\lambda _{J0}=16r(T)\frac{ϵ}{\sqrt{1+ϵ^2}},$$
(94)
and the flux is

$`\mathrm{\Phi }=2\pi {\displaystyle \frac{ϵ^2}{1+ϵ^2}},`$

which immediately leads to the lower critical field:

$`H_{c1}={\displaystyle \frac{2}{\pi }}H_s={\displaystyle \frac{2}{\pi }}{\displaystyle \frac{\sqrt{1+ϵ^2}}{ϵ}}.`$

As can be seen by comparing (94) with the single-junction expression (90), the energy per vortex in the double-junction stack is higher. The vortex-plane solutions are presented schematically in Fig. 2.

## V Discussion

Within the framework of standard methods of the theory of differential equations, we have obtained a complete mathematical description of the Meissner effect and the vortex structure in periodic thin-layer $`(N1)`$-junction stacks ($`2N<\mathrm{}`$). The results of our analytical analysis of the coupled static SG equations (25), summarized in Lemmas 1,2 and the Theorem, should provide a basis for any further analytical or numerical study of these equations, not necessarily restricted to the field of superconductivity.

By proving the absence of single-vortex solutions to (25) for $`H>0`$ and specifying the actual domain of their existence ($`H=0`$, $`L=\mathrm{}`$), we have clarified a wide-spread misunderstanding: The previous estimates of $`H_{c1}`$ turn out to be irrelevant because of absolute thermodynamic instability of single-vortex configurations. In full agreement with our study of infinite ($`N=\mathrm{}`$) layered superconductors, we have shown that the only true topological (vortex) configurations that survive at $`H>0`$ are the vortex planes. Being a natural generalization of ordinary Josephson vortices in a single junction ($`N=2`$), the vortex-plane solutions for $`3N`$ inherit such properties of the former as periodicity along the layers \[Eq. (39)\] and the overlapping of states with different topological numbers $`N_v`$ \[Eq. (54)\]. A unified mathematical approach of this paper, valid for $`2N<\mathrm{}`$, allowed us to derive an exact expression \[Eq. (49)\] for the superheating (penetration) field $`H_s`$ as an explicit function of $`N`$. Within the framework of the same unified approach, we have obtained a new exact solution to the single SG equation, Eqs. (76)-(80), valid for any $`H0`$ and $`L>0`$. This solution refutes the assertions that Josephson vortices ”do not form” in single junctions with small $`L`$.

Our closed-form analytical expression for the Meissner field, Eq. (51), clearly illustrates an important feature of the Meissner effect in $`(N1)`$-junction stacks, not noticed in previous theoretical publications, namely the existence of $`\left[\frac{N}{2}\right]`$ different Josephson lengths $`\lambda _{Jk}`$ ($`k=0,1,\mathrm{},\left[\frac{N}{2}\right]1`$). This result may prove to be useful in view of the current experimental efforts to verify the interlayer tunneling model of high-$`T_c`$ superconductivity by measuring the $`c`$-axis penetration depth. \[The penetration of the parallel magnetic field with a distribution of length scales has been recently observed in the organic layered superconductor $`\kappa `$-(BEDT-TTF)<sub>2</sub>Cu(NCS)<sub>2</sub>.\]

Finally, the vortex-plane solutions of our paper provide a natural explanation for the observed coherent Josephson-vortex configurations in artificial double-junction stacks and weakly-coupled multilayers  in the presence of a static parallel external field $`H>0`$. The observation of isolated interlayer vortices in the high-$`T_c`$ superconductor Tl<sub>2</sub>Ba<sub>2</sub>CuO<sub>6+δ</sub> and in $`\kappa `$-(BEDT-TTF)<sub>2</sub>Cu(NCS)<sub>2</sub> can be explained by the fact that in these experiments the external field was set equal to zero ($`H=0`$). In this situation, isolated Josephson vortices can be pinned by structural defects, because their self-energy is lower and the spatial extension along the $`c`$-axis is smaller than those of the vortex planes at $`H=0`$. We hope that our results will stimulate further experimental investigations.

## A The solution of the finite difference equation for $`H_n(y)`$

Equations (17) can be regarded as a nonhomogeneous finite difference equation for $`H_n(y)`$ with respect to the layer index $`n`$, subject to boundary conditions (8). According to general theory of such equations, its solution can be represented in the form

$$H_n(y)=h_n(y)+H_n,$$
(A1)
where

$$h_n(y)=\frac{ϵ^2}{2}\underset{m=1}{\overset{N1}{}}G(n,m)\frac{d\varphi _m(y)}{dy}$$
(A2)
is the particular solution of (17) satisfying the boundary conditions

$$h_0(y)=h_N(y)=0,$$
(A3)
and

$$H_n=\frac{H\left(\mu ^n+\mu ^{N+n}\mu ^n\mu ^{Nn}\right)}{\mu ^N\mu ^N},$$
(A4)
$$\mu =1+\frac{ϵ^2}{2}ϵ\sqrt{1+\frac{ϵ^2}{4}},\text{ (}0<\mu <1\text{),}$$
(A5)
is the solution of the homogeneous form of (17) (with the zero right-hand side) meeting the boundary conditions

$$H_0=H_N=H.$$
(A6)
The quantities $`G(n,m)`$ in (A2) are elements of a $`\left(N1\right)\times \left(N1\right)`$ matrix $`𝐆(1n,mN1)`$. They obey the nonhomogeneous finite difference equation

$$\left(2+ϵ^2\right)G(n,m)G(n+1,m)G(n1,m)=\delta _{n,m}$$
(A7)
($`\delta _{n,m}`$ is the Kronecker index) with the boundary conditions

$$G(0,m)=G(N,m)=0.$$
(A8)
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].$$
(A9)
The following properties of $`G(n,m)`$ can be easily verified using (A7) and (A9):

$$G(n,m)=G(m,n),$$
(A10)
$$G(n,Nm)=G(Nn,m),$$
(A11)
$$G(n,m)>0\text{ for any }1n,mN1,$$
(A12)
$`{\displaystyle \underset{m=1}{\overset{N1}{}}}G(n,m)={\displaystyle \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,$$
(A13)
$$\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}.$$
(A14)
According to (A7), the inverse of the matrix $`𝐆(n,m)`$ is a Jacobian matrix $`𝐆^\mathrm{𝟏}(n,m)`$ with elements $`G^1(n,n)=2+ϵ^2`$ ($`n=1,2,\mathrm{},N1`$), $`G^1(n+1,n)=G^1(n,n+1)=1`$ ($`n=1,2,\mathrm{},N2`$), and $`G^1(n,m)=0`$ for $`\left|nm\right|>1`$. Hence all the eigenvalues $`e_j`$ of $`𝐆(n,m)`$ are positive and given by

$$e_j=\frac{\lambda _j^2}{ϵ^2},\lambda _j=\frac{ϵ}{\sqrt{2+ϵ^22\mathrm{cos}\frac{\pi j}{N}}},i=1,2,\mathrm{},N1,$$
(A15)
with the corresponding eigenvectors

$$u_{ij}=\sqrt{\frac{2}{N}}\mathrm{sin}\frac{\pi ij}{N},i,j=1,2,\mathrm{},N1.$$
(A16)
The quantities $`\lambda _j`$ determine the characteristic length scales of the stack with $`N`$ S-layers.

Note that for an infinite layered superconductor, the analog of the matrix $`𝐆(n,m)`$ is $`𝐆_{\mathrm{}}(n,m)`$, determined by the matrix elements

$$G_{\mathrm{}}(n,m)=\frac{\mu ^{\left|nm\right|}}{2ϵ\sqrt{1+\frac{ϵ^2}{4}}},\mathrm{}<n,m<+\mathrm{},$$
(A17)
that satisfy the summation rule

$`{\displaystyle \underset{m=1}{\overset{N1}{}}}G_{\mathrm{}}(n,m)={\displaystyle \frac{1}{ϵ^2}}.`$

## B Proof of Lemma 1

By introducing new functions

$`\psi _1(y)=\varphi _1(y),\psi _2(y)=\varphi _2(y),\mathrm{},\psi _{N1}(y)=\varphi _{N1}(y),`$

$$\psi _N(y)=\frac{d\varphi _1(y)}{dy},\psi _{N+1}(y)=\frac{d\varphi _2(y)}{dy},\mathrm{},\psi _{2N2}(y)=\frac{d\varphi _{N1}(y)}{dy},$$
(B1)
we convert (25) into an equivalent normal system of $`2N2`$ first-order equations

$$\frac{d\psi _i(y)}{dy}=F_i(\psi _1,\psi _2,\mathrm{},\psi _{2N2}),1i2N2,$$
(B2)
$`F_i(\psi _1,\psi _2,\mathrm{},\psi _{2N2})\psi _{i+N1},1iN1,`$

$`F_N(\psi _1,\psi _2,\mathrm{},\psi _{2N2}){\displaystyle \frac{1}{ϵ^2}}\left[\left(2+ϵ^2\right)\mathrm{sin}\psi _1\mathrm{sin}\psi _2\right],`$

$`F_i(\psi _1,\psi _2,\mathrm{},\psi _{2N2}){\displaystyle \frac{1}{ϵ^2}}\left[\left(2+ϵ^2\right)\mathrm{sin}\psi _i\mathrm{sin}\psi _{i1}\mathrm{sin}\psi _{i+1}\right],N+1iN3,`$

$`F_{2N2}(\psi _1,\psi _2,\mathrm{},\psi _{2N2}){\displaystyle \frac{1}{ϵ^2}}\left[\left(2+ϵ^2\right)\mathrm{sin}\psi _{N1}\mathrm{sin}\psi _{N2}\right],`$

subject to initial conditions

$`\psi _i(y_0)=\alpha _i,1iN1,`$

$$\psi _i(y_0)=\beta _{iN+1},Ni2N2.$$
(B3)
To prove the statement of Lemma 1, it is sufficient to observe that all $`F_i(\psi _1,\psi _2,\mathrm{},\psi _{2N2})`$ are continuous functions of their arguments for $`y(\mathrm{},+\mathrm{})`$ and $`\psi _k(\mathrm{},+\mathrm{})`$ ($`1k2N2`$). Moreover, their partial derivatives with respect to $`\psi _k`$ satisfy the relation

$$\left|\frac{F_i(\psi _1,\psi _2,\mathrm{},\psi _{2N2})}{\psi _k}\right|\frac{4+ϵ^2}{ϵ^2}$$
(B4)
for $`y(\mathrm{},+\mathrm{})`$ and $`\psi _k(\mathrm{},+\mathrm{})`$ ($`1i,k2N2`$). Thus, the Lipschitz conditions with respect to $`\psi _k`$ are met for $`y(\mathrm{},+\mathrm{})`$ and $`\psi _k(\mathrm{},+\mathrm{})`$ ($`1k2N2`$), which immediately guarantees the existence and uniqueness of a solution to (B2), satisfying arbitrary initial conditions (B3), in an arbitrary interval $`I=[L_1,L_2]`$ such that $`y_0I`$. Continuous dependence of the solution on initial data is a result of continuous dependence of $`F_i(\psi _1,\psi _2,\mathrm{},\psi _{2N2})`$ on their arguments and of the condition (B4). Infinite differentiability of the solution automatically follows from infinite differentiability of $`F_i(\psi _1,\psi _2,\mathrm{},\psi _{2N2})`$ with respect to their arguments.

## C Proof of the Theorem

Owing to the property (27), without loss of generality, we can consider the initial value problem (38) for $`y=y_00`$. To prove the Theorem, we have to show that there exists a solution to (25) in $`(\mathrm{},+\mathrm{})`$ with the properties (39)-(43) that at $`y=y_00`$ satisfies (38). By Lemma1, this solution will represent the sought unique solution to the initial value problem (38).

Consider an arbitrary finite interval $`I=[\frac{P}{2},\frac{P}{2}]`$. We start with the Meissner boundary value problem in the interval $`I`$:

$$\frac{d\varphi _n}{dy}\left(\pm \frac{P}{2}\right)=2\stackrel{~}{H}>0,\varphi _n(0)=0,1nN1.$$
(C1)
Up to a certain $`\stackrel{~}{H}=\stackrel{~}{H}_s`$, the problem (C1) has a unique solution with the properties

$$\varphi _n(y)=\varphi _n(y),yI,1nN1,$$
(C2)
$$\frac{d\varphi _n}{dy}\left(y\right)>0,yI,1nN1,$$
(C3)
$$\pi \varphi _n\left(\frac{P}{2}\right)<0,0<\varphi _n\left(\frac{P}{2}\right)\pi ,1nN1,$$
(C4)
where $`\stackrel{~}{H}_s`$ is determined by the conditions

$$\varphi _n\left(\pm \frac{P}{2}\right)=\pm \pi ,1nN1,$$
(C5)
and satisfies the inequality $`\stackrel{~}{H}_s>H_s`$, with $`H_s`$ defined via (A14). (See section IIIB.)

Now we will prove that all $`\frac{d\varphi _n}{dy}\left(0\right)`$, where $`\varphi _n`$ are the solution to (C1), for $`\stackrel{~}{H}\stackrel{~}{H}_s`$ tend to the same limiting value

$$\frac{d\varphi _n}{dy}\left(0\right)=2\sqrt{\stackrel{~}{H}_s^2H_s^2},1nN1.$$
(C6)
Indeed, relations (C6) satisfy (29) at $`y=0`$ with (30), where the substitution $`H\stackrel{~}{H}_s`$, $`L\frac{P}{2}`$ has been made. However, it is still necessary to show that the limiting value (C6) is uniquely determined by the solution to (C1). To this end, we consider the relation

$$4\left(N1\right)\frac{\stackrel{~}{H}}{H_s^2}ϵ^2\underset{n=1}{\overset{N1}{}}\underset{m=1}{\overset{N1}{}}G(n,m)\frac{d\beta _n}{d\stackrel{~}{H}}\beta _m=\underset{n=1}{\overset{N1}{}}\mathrm{sin}\varphi _n\left(\frac{P}{2}\right)\frac{d\varphi _n\left(\frac{P}{2}\right)}{d\stackrel{~}{H}}0,$$
(C7)
where $`\beta _n\frac{d\varphi _n}{dy}\left(0\right)`$. \[Relation (C6) is obtained from (29), (30) by the substitution $`H\stackrel{~}{H}`$, $`L\frac{P}{2}`$ and differentiation with respect to $`\stackrel{~}{H}`$.\] Notice that

$$\beta _n<2\stackrel{~}{H},1nN1,$$
(C8)
in the whole range $`0<\stackrel{~}{H}\stackrel{~}{H}_s`$, because all $`\frac{d\varphi _n}{dy}`$ have minima at $`y=0`$. \[See (28).\] With increasing $`\stackrel{~}{H}`$ from zero to $`\stackrel{~}{H}_s`$, all $`\varphi _n\left(\frac{P}{2}\right)`$ increase monotonously from zero to $`\pi `$. Thus, the difference on the left-hand side of (C6) first increases (for $`0<\stackrel{~}{H}H_s`$) and then decreases (for $`H_s\stackrel{~}{H}\stackrel{~}{H}_s`$). It means that the dependence $`\beta _n(\stackrel{~}{H})`$ changes from a linear one for $`0<\stackrel{~}{H}H_s`$ \[i.e., $`\beta _n(\stackrel{~}{H})=\overline{b}_n\stackrel{~}{H}`$, $`0<\overline{b}_n<2`$, $`1nN1`$\] to a nonlinear one for $`H_s\stackrel{~}{H}\stackrel{~}{H}_s`$, with

$$\frac{d\beta _n}{d\stackrel{~}{H}}>2,1nN1.$$
(C9)
At $`\stackrel{~}{H}=\stackrel{~}{H}_s`$ the right-hand side of (C6) is equal to zero, and we have

$$\frac{d\beta _n}{d\stackrel{~}{H}_s}\beta _m=b_nb_m\stackrel{~}{H}_s,$$
(C10)
where $`b_n`$ satisfy the conditions

$$ϵ^2\underset{n=1}{\overset{N1}{}}\underset{m=1}{\overset{N1}{}}G(n,m)b_nb_m=\frac{4\left(N1\right)}{H_s^2},$$
(C11)
$$0<b_{\mathrm{min}}b_nb_{\mathrm{max}},b_{\mathrm{min}}=\underset{n}{min}\left\{b_n\right\}2,b_{\mathrm{max}}=\underset{n}{max}\left\{b_n\right\}2.$$
(C12)
The integration of (C10) for $`n=m`$ under the condition $`\beta _n\left(H_s\right)=0`$ yields

$$\beta _n\left(\stackrel{~}{H}_s\right)=b_n\sqrt{\stackrel{~}{H}_s^2H_s^2}.$$
(C13)
With the help of (C8) and (C9), taken at $`\stackrel{~}{H}=\stackrel{~}{H}_s`$, we establish the following inequalities for $`b_{\mathrm{min}}`$ and $`b_{\mathrm{max}}:`$

$$2\sqrt{1\frac{H_s^2}{\stackrel{~}{H}_s^2}}<b_{\mathrm{min}}2,2b_{\mathrm{max}}<\frac{2}{\sqrt{1\frac{H_s^2}{\stackrel{~}{H}_s^2}}}.$$
(C14)
Inequalities (C14) must hold for any $`\stackrel{~}{H}_s>H_s`$. Proceeding to the limit $`\frac{H_s^2}{\stackrel{~}{H}_s^2}0`$, we get $`b_{\mathrm{min}}=b_{\mathrm{max}}=2`$, hence the relations (C6).

The Meissner solution at $`\stackrel{~}{H}=\stackrel{~}{H}_s`$ can be uniquely continued from the interval $`I`$ into the whole interval $`(\mathrm{},+\mathrm{})`$. The solution $`\varphi _n`$ ($`1nN1`$) thus obtained possesses all the required properties (38)-(43), where $`y_00`$, and $`H\sqrt{\stackrel{~}{H}_s^2H_s^2}`$. To show this, it is sufficient to prove the property (39). \[The properties (38), (41) and (42) are obvious. The property (40) results from the differentiation of (39). The property (43) is obtained from (29), (30) by integration.\]

Consider a set of functions

$$\overline{\varphi }_n(y)=\varphi _n\left(y+P\right)2\pi ,1nN1,$$
(C15)
where $`\varphi _n(y)`$ ($`1nN1`$) are the continuation of the Meissner solution at $`\stackrel{~}{H}=\stackrel{~}{H}_s`$ into the interval $`(\mathrm{},+\mathrm{}\dot{)}`$. Owing to the properties (26), (27), the functions (C15) also represent a solution to (25) in the interval $`(\mathrm{},+\mathrm{}\dot{)}`$. At $`y=\frac{P}{2}`$, the solution $`\overline{\varphi }_n`$ ($`1nN1`$) satisfies the same initial conditions as the solution $`\varphi _n`$ ($`1nN1`$). By Lemma 1, the latter means that $`\overline{\varphi }_n\varphi _n`$ ($`1nN1`$) in the interval $`(\mathrm{},+\mathrm{}\dot{)}`$, i.e.,

$$\varphi _n\left(y+P\right)2\pi =\varphi _n(y),1nN1,$$
(C16)
which accomplishes the proof of the Theorem. Note that in this proof we have not employed the symmetry relations (21) resulting from the boundary conditions (8).

FIGURE CAPTIONS

1. Fig. 1. The geometry of the problem (schematically). Here $`N=4`$, $`ap`$ ($`a`$ is the S-layer thickness), $`2L=W`$, $`H>0`$.
2. Fig. 2. The vortex-plane solutions in a double-junction stack (schematically): a) the state with a single vortex plane ($`N_v=1`$); b) the state with two vortex planes ($`N_v=2`$). The black dots conventionally denote the positions of Josephson vortices, and the arrows show the distribution of intralayer and Josephson currents.