Transverse Josephson vortices and localized states in stacked Bose-Einstein condensates

Stacks of Bose-Einstein condensates coupled by long Josephson junctions present a rich phenomenology susceptible of experimental realization. We show that transverse Bloch waves excited in an array of one-dimensional condensates can carry tunneling super ows that are dynamically stable. Such states yield transverse Josephson vortices with a generic non-integer circulation in units of h=m across a stack with periodic boundary conditions forming a closed ring-shaped system. Additionally, superpositions of (linear) degenerate Bloch waves can cancel the supercurrents and give rise to families of standing-wave states with strong localization across the stack. Stable states of this type can also be found in finite size systems.


I. INTRODUCTION
Quantum tunneling is one of the most striking phenomena predicted by quantum mechanics. At a macroscopic scale it is named Josephson effect. This phenomenon is a paradigm of the phase coherence manifestation of a macroscopic quantum system. The theory of the Josephson effect was developed by B. D. Josephson for superconducting electron pairs in 1962 [1]. Since then, it has found multiple technological applications [2]. More recently, with the advent of Bose-Einstein condensates (BECs) of ultracold atoms, the Josephson effect has been demonstrated between two weakly linked condensates of neutral bosonic atoms [3,4]. As a common feature, a macroscopic quantum system is described by a complex order parameter, whose modulus and relative phase account for the density and the superfluid velocity, respectively. The particle superflow that traverses the coupled components in the Josephson effect, has a current that depends on the sine of the relative phase between the order parameters of the two weakly connected components.
The tunneling of particles between two coupled quantum systems occurs through a weak link, the Josephson junction, that can be classified into two types (short and long) [5]. When the coupling occurs at a single point, for instance a BEC in a double well potential, the link is a short Josephson junction [3,4]. Whereas when it is extended in space, the coherent transfer occurs locally along each point of the connection. Then the relative phase varies along the junction, which corresponds to a long Josephson junction. For instance, two spin components of a BEC coupled by a Raman laser [6] (internal Josephson effect) or two elongated BECs coupled along the transverse direction, as we consider in the present paper.
Most of the previous works on systems with extended coupling in ultracold gases have focused on binary BECs (see e.g. [7] and references therein). In general, the long Josephson junctions have received less attention than the short (or point-like) ones, in spite of the fact that the former presents a richer phenomenology capable of realizing superfluid equivalents of complex superconducting states [2,[8][9][10][11][12][13][14]. On the contrary, the coupled long Josephson junctions have been extensively addressed in layered superconductors [15], mainly within the model of coupled sine-Gordon equations (see, for instance [16,17]). As far as we know, only one-dimensional (1D) arrays of point-like Josephson junctions with BECs have been already experimentally realized [18]. Theoretically, the superfluid-insulator transition has been studied in twodimensional (2D) arrays of coupled 1D tubes against the absence or presence of axial periodic potentials [19]. Also, the propagation of bright solitons in arrays of BECs with negative nonlinearity has been considered [20]. Very recently, coupled atomic wires have been proposed in ultracold-gas systems for the generation of exotic phases in the presence of synthetic gauge fields [21]. In addition, the arrays of parallel one-dimensional long Josephson junctions in BECs with positive nonlinearity have been demonstrated to provide an excellent playground for the realization and stabilization of solitary waves [22].
In this work we study, analytically and numerically, a stack of linearly coupled 1D BECs with repulsive interparticle interactions that gives rise to an underlying array of coupled-parallel long Josephson junctions. In particular, we consider a stack with periodic boundary conditions forming a closed ring-shaped array of long Josephson junctions. Within the mean field description we address the dynamics of the stationary states composing a transverse, discrete Bloch band by solving the Gross-Pitaevskii equation. We will show that the periodic boundary conditions yield transverse vortex states that carry Josephson supercurrents. Interestingly, the Bogoliubov analysis reveals that these Josephson vortices can be dynamically stable, hence susceptible of experi-mental detection. From a hydrodynamical perspective, we also perform a long-wave-length linear analysis of the stationary states that allow for a somehow simpler interpretation of the system dynamics in terms of coupled wave equations for the relative phases and densities. We show that, only for particular states in the small coupling regime, the resulting hydrodynamic model resembles, although it cannot be identified with, a model of coupled sine-Gordon equations. Finally, we explore how the linear combinations of Bloch waves lead to families of nonlinear states that cancel the Josephson supercurrents and produce strongly localized density profiles that can be stable against perturbations in finite size systems.
This paper is organized as follows: in Sec. II we present the theoretical model that describes the stack of linearly coupled BECs with periodic boundary conditions in the mean-field framework. We investigate the stationary states of the system in Sec. III, in particular Bloch waves and localized states. Section IV is devoted to the analytical study of the linear stability of the stationary states in the stack, within the Bogoliubov approximation and in the hydrodynamic approach. The comparison between analytical and numerical results is presented in Sec. V, specially for transverse Josephson vortices. The nonlinear dynamics of localized states is discussed in Sec. VI. And a final discussion together with the conclusions are presented in Sec. VII.

II. MODEL OF LINEARLY COUPLED BECS
We consider a system of M coupled 1D BECs, with coherent linear coupling of frequency Ω along the junctions. Each BEC is described by the corresponding order parameter Ψ j (x, t) = n j (x, t) exp [iθ j (x, t)], where n j and θ j are the local density and phase, and j = 1, 2, ...M labels the BEC in the stack. We assume a particular arrangement with periodic boundary conditions, that is, the stack has a ring-shaped configuration with M junctions: the j-th junction lying between the j-th and (j+1)th components, and the M -th junction connecting the M -th and the 1-st BEC. The total density of the system, n T (x) = j |Ψ j | 2 , is normalized to the total number of particles, N = n T dx = j N j , which is a conserved quantity, and N j is the number of particles in the j-th BEC. A schematic representation of the system is shown in Fig. 1. Detailed proposals for the experimental realization of such systems in ultracold gases [22], including gauge dependent coupling [21], has been recently presented. At zero temperature, in the mean-field framework, the system can be described by a set of coupled Gross-Pitaevskii (GP) equations, namely for the j-th component: where the sum on the right hand side extends to the first neighbours (j − 1) and (j + 1). The 1D particle interaction strength is g = 2 ω ⊥ a, the s-wave scattering length is repulsive a > 0, ω ⊥ is the frequency of a tight transverse trap, and m is the bosonic mass. It is convenient for a later discussion, to rewrite the GP Eq. (1) in the hydrodynamic form, in terms of the densities and phases: where Q j = − 2 ∇ 2 √ n j /2m √ n j and v j = ∂ x θ j /m are the quantum-pressure energy term and the superfluid velocity, respectively, of component j. Equation (2) is the continuity equation and expresses the particle conservation in the array. The local density n j (x) varies due to either changes in the local current, J j (x) = n j (x)v j (x), within the j-th BEC or due to the Josephson current, J lj , across the junctions. We define the latter, by analogy with the local current, as J lj (x) = n lj (x) v lj (x), by means of the Josephson velocity v lj = m δy where the relative phases are θ lj = θ l − θ j , and the geometric mean densityn lj = √ n l n j . With this definition, the last term of Eq. (2) corresponds to a discrete derivative δJ lj /δy = (J j+1,j − J j,j−1 )/δy of the Josephson current. Along with the periodic boundary conditions in the y-direction, it allows for the computation of a velocity circulation around the stack On the other hand, Eq. (3) is the equation of motion of the phase, which varies locally according to the local energy content. In a stationary state, where the time variation of the phase is given by the frequency µ/ (see Eq. (6) below), where µ is the chemical potential, the local energy in the right hand side of Eq.(3) is the same (and equals µ) at every point in the system.

III. STATIONARY STATES
The BEC stack forms a discrete lattice of M sites along the y-direction, transverse to the common BEC x-axis (see Fig. 1). An effective distance of separation between two neighbour BECs can be defined provided by the linear coupling: δy = /mΩ. Thus, the discrete coordinate along the y-axis takes values y j = jδy for each j-th BEC. The characteristic length δy has to be compared with the healing length ξ = / √ m g n determined by the axial density n of the BECs. The ratio ξ/δy = Ω/gn is a measure of the tunneling-coupling strength.

Bloch states
The lattice configuration along y allows us to look for stationary states that take the form of transverse Bloch waves where ψ(x) is the axial wavefunction (with the periodicity of the discrete lattice), and K k is the transverse quasimomentum. Due to the discreteness of the system, the quasimomentum can take only M different integer values within the first Brillouin zone: where M/2 is the greatest integer less than or equal to M/2. As a result, the product space of coordinateseparated solutions Ψ j,k (x, t) = φ(x) × φ k (y j ) presents an M -fold symmetry for each wave function φ(x). All the k states but k = 0 and k = M/2 (when the latter exists), at the middle and at the end, respectively, of the Brillouin zone, are states supporting Josephson currents in the y-direction, due to the existence of nonzero and non-π relative phases between consecutive condensates.
with stationary phases θ j = k r = K x x + K k y j , where k and r are the momentum and spatial vectors, respectively. In Fig. 2 red bars represent the density and phase of Bloch waves in the simplest case of M = 3 and K x = 0 (see also Sect. V A below).
The relative phases of the Bloch waves (8), Thus, the Josephson currents are J lj = ± n sin(2πk /M )/mδy, which in the limit of large M and small k tend to J lj = n K k /m. Analogously, the transverse circulation Eq. (5) gives Γ k = ±M sin(2πk /M )/m, and tends to Γ k = 2π k /m in the same limit. The latter expression is equivalent to the circulation of a regular vortex of charge k. We will refer to these discrete vortices in the stack as Josephson vortices. Note that they are delocalized along the x-direction, and, unlike the regular vortices, do not show in general an integer-valued circulation in units of h/m.
By substituting the constant density states of the type (8) in Eq. (1), or Eq. (3), it yields the chemical potential of such Bloch waves which takes values in the range µ k ∈ gn + 2 K 2 x /2m + [− Ω, Ω], within a discrete band of energy width 2 Ω (see top panel of Fig. 2). In the limit 2πk/M 1, expanding Eq. (9) up to quadratic terms one finds µ k = gn− Ω+ 2 (K 2 x +K 2 k )/2m, which is a quadratic dispersion (as in a fully 2D continuous system) around the ground state.

Standing waves: Localized states
Along with the Bloch waves, which realize transverse current states with definite quasimomentum, the BEC stack admits standing-wave solutions without definite quasimomentum. Their existence can be tracked up to the non-interacting regime (g = 0), where the standing waves can be built from linear combination of energydegenerate Bloch waves (8) with same quasimomentum modulus |k| [23]. Due to the finite size of the stack, such states break the lattice symmetry and show differences in the particle density between neighbor BECs, which lays at the origin of the localized-density states in periodic systems known as gap solitons [23].
For each pair of energy-degenerate Bloch states sharing |k| one can also build a pair of independent linear combinations with equal weight, which we will denote by Ψ s j,|k| and Ψ a j,|k| , that have continuation in the nonlinear regime. The simplest example is the (M = 3)-stack for the families of states originated from the (real) symmetric Ψ s j,|1| = Ψ j,1 + Ψ j,−1 ∝ cos(2πj/3) and antisymmetric Ψ a j,|1| = Ψ j,1 − Ψ j,−1 ∝ sin(2πj/3) superpositions of Bloch waves with |k| = 1. Figure 2 shows the chemical potential (top panel), and the density and the phase (bottom panel) of these families in a system with zero axial momentum. Contrary to the original Bloch waves Ψ j,±1 , having equal density across the stack, the antisymmetric states contain a nodal strand Ψ a 2,|1| = 0 and two density peaks Ψ a 1,|1| = −Ψ a 3,|1| , whereas the symmetric states present a single density peak n s 2,|1| > (n s 3,|1| = n s 1,|1| ). In the nonlinear regime, for the same total density n T in the stack, the antisymmetric states have higher chemical potential than the symmetric ones.
We will focus on these types of states, having one or two density peaks, for varying M . As we will see, in the small coupling limit Ω/2gn 1 and M > 3, where n is now the maximum density in the stack, the mentioned density peaks accumulate practically all the system density. In a symmetric state with K x = 0, for instance, the nearest-neighbor BECs j ± 1 of the peak-density strand j follow GP equations which, assuming decreasing amplitudes ψ j ψ j±1 ψ j±2 , can be approximated up to first order in the neighbor amplitudes by ( Ω/2)ψ j ≈ −µ ψ j±1 , and µ ≈ g n. As a result, the l-site BECs show decreasing amplitudes in a factor for distances |l − j| away from the density peak at j. The prototypical density profile is illustrated in Fig. 3 for different values of the coupling and same chemical potential in a stack with M = 10. The bars indicate the density at each j-th BEC as calculated from the numerical solution of the Gross-Pitaevskii equations, while the open symbols mark the values given by the analytical approximation Eq. (11). The same scenario takes place for the antisymmetric states with M > 4, with decreasing amplitudes of opposite signs beyond the peaks at both sides of the nodal strand.

IV. LINEAR EXCITATIONS
In this section we focus on the analytical treatment of the linear stability of stationary states in the stack, both Bloch waves and localized states. First we make a general analysis based on the Bogoliubov equations for the linear excitations [24] in order to find stability regimes. For Bloch waves, we also consider a long-wave-length excitation approach to derive and discuss the resulting coupled sine-Gordon-like equations that are usual in superconducting Josephson junctions [17].
A. Linear stability of Bloch waves

Bogoliubov equations
Let us consider the linear excitations [u j,k (x), v j,k (x)] around the stationary states (8). After substituting (1), the excitation modes satisfy the Bogoliubov equations where where q p = 2πp/M δy is the transverse momentum of the excitation for integer p = 0, ±1, ±2, . . . M/2 , the Bogoliubov equations get decoupled for each two-dimensional wave number q = (q x , q p ): where E k,± = 2 (q 2 x ±2K x q x )/2m+gn+ Ω cos(2π k/M ). Further, we introduce the linear combinations of modes f where ζ qx = 2 q 2 x /2m is the kinetic energy of the modes along each 1D BEC, and the transverse excitations are defined through the parameters and Therefore, for each stationary state Ψ k , by solving Eqs. (16) and (17) for the frequency ω, we get the dispersion relation composed of M double branches corresponding to the different values of p. Among them, the p = 0 branch is always gapless (because α k,0 = β k,0 = 0), whereas the rest presents energy gaps given by which show up even in the non-interacting (g = 0) case. The speed of sound along each BEC (i.e. along the x coordinate) can be calculated from the gapless branch in the long wavelength limit: These quantities, ω g and c, are relevant for the energetic stability of the system. Both define energy thresholds (gaps), ω g and mc 2 , for the corresponding (energetic) instabilities in the underlying superfluid against external perturbations.
When the excitation energy (20) takes real values, we choose the usual normalization, dx (|u j | 2 − |v j | 2 ) = 1, for the corresponding modes (u j , v j ). Then, selecting real values for the Fourier amplitudes c q and d q , it follows that |c q | 2 −|d q | 2 = f The dispersion relation Eq. and k = 1, which is a discrete transverse vortex with circulation Γ 1 = 4 /m (see below).

Long wave length excitations: hydrodynamic approach
We start by introducing low energy perturbations [δn j (x, t), δθ j (x, t)] around the density and the phase of an equilibrium state Ψ j = √ n j exp(iθ j ) −→ n j + δn j exp(iθ j + iδθ j ). Then, we substitute the perturbed states in Eqs. (2) and (3), and keep terms up to first order in the perturbations. We focus on the analysis of Bloch states with n j = n. The mentioned procedure leads to 1 n −Ω α k (sin δθ j+1,j − sin δθ j,j−1 ) + β k δn j+1,j − δn j,j−1 2n , where δθ lj = δθ l − δθ j , δn lj = δn l − δn j , δn lj = δn l + δn j are the perturbations in relative phase, relative density and average density, respectively, and α k = cos(2π k/M ), β k = sin(2π k/M ). As usual in a long wavelength approximation, we have dropped the quantum-pressure term in Eq. (27). For reasons that will become apparent later, we have kept the sine functions sin δθ lj even in the linear approximation in order to track the Josephson currents, but they will be replaced by their argument (for consistency within the assumed first order approximation) at intermediate steps of the analytical derivations.
In what follows, we use the short notation ρ j = δn j+1,j /n,ρ j = δn j+1,j /n, φ j = δθ j+1,j , and alsō φ j = δθ j+1 + δθ j for the total phase. Since these quantities appear explicitly in previous expressions, we look for their equations of motion by adding and subtracting Eqs. (26) and Eq. (27) for consecutive components. For the relative quantities we get where the discrete operator δ 2 acts as δ 2 f j = f j+1 − 2f j + f j−1 . Exactly the same equations are obtained for the total quantities substituting ρ byρ and φ byφ. As can be seen, relative and total quantities are decoupled in pairs of equations (28)-(29). Within each pair, by taking the time derivative of one of the equations and making use of the others, it yields wave-like equations: whereΩ = Ω/2gn. These equations in the limit of small Josephson coupling,Ω 1, become: The above system of M pairs of equations (32)-(33) describes the dynamics of the underlying array of junctions determined by the GP Eqs. (1) in the limit of small coupling,Ω 1, and small relative phases, φ j 1. The k-dependent factor α k multiplying the transversal (discrete) derivative indicates a varying penetration length ξ J,k = δy/ |α k |. The sine functions in Eq. (33) should be formally substituted by their arguments sin φ j ≈ φ j , since we have assumed a linear regime. By keeping them, we highlight the correspondence with a sine-Gordon-like equation, where kink solutions can be found [2]. As a limiting case, the kink-type solutions have demonstrated to be useful in the search of solitonic states in the nonlinear dynamics of two-component condensates [11,14,[25][26][27].
Equations (32)-(33) admit plane wave solutions with the same phase shifts across the stack previously given in Eq. (7): with transverse momentum q p = 2πp/M δy. After substitution, one gets the amplitude relation c q = i ν d q , where ν is a real number, and the double-branched dispersion For β k = 0 (that is, when k = 0 or k = M/2, and α k = ±1) both branches coincide, in agreement with the Bogoliubov dispersion Eq. (20) in the small coupling limitΩ 1 considered here. The gaps are given by ω (1) k,p = c 2α k,p /δy. Additionally, for q p 1 the dispersion can be written as ω k,p = c q 2 x ± q 2 p , which is the usual dispersion ω = c |q| of a continuous (relativistic) 2D system when k = 0, and contains instabilities when The second dispersion branch ω (2) k,p , is not fully consistent with the Bogoliubov spectrum for β k,p = 0, unless ν = 1 and the limit Ω → 0 are considered. Such inconsistencies have been previously found in two coupled condensates (see e.g. [7]), and arise from neglecting the quantum pressure term.

B. Linear stability of localized states
We study the stability of standing-wave states having one or two density peaks in the stack of constant density BECs. The analysis of the simplest systems with M = 3, 4 components serves as an insightful starting point, with straightforward analytical solutions for the antisymmetric states. For a generic M value, we focus on the small coupling limit (Ω = Ω/2gn 1), where now n is the peak (or maximum) density in the stack.
The symmetric states involve a more cumbersome algebra. For the simplest stack with M = 3, the chemical potential is µ s |1| = 2 K 2 x /2m + (2gn s 1,|1| + gδn s 1,|1| − Ω/2 + (gδn s 1,|1| + Ω/2) 2 + 2( Ω) 2 )/2, where δn s 1,|1| = n s 2,|1| − n s 1,|1| is the density contrast between components, which increases with the chemical potential for a given coupling. An asymptotic analysis can be readily done in the small coupling limit Ω 1, where the densities fulfill (n s 2,|1| ≈ n T ) and (Ψ s 1,|1| = Ψ s 3,|1| ) ≈ −Ω Ψ s 2,|1| , and the chemical potential µ s |1| ≈ 2 K 2 x /2m + gn. In this approximation, the Fourier expansion of the excitation modes, lead to the same functional form of the spectrum (41) with the substitution of Ω by Ω/ √ 2. As a consequence, dynamical stability of the (one-peak) symmetric states is also expected in the smallΩ regime. In this case, the strong localization of the density in a single strand, with neighbor densities decreasing as powers ofΩ 2 , allows the stability prediction to be extended to arbitrary M , and even to systems with open boundary conditions. Furthermore, the same analysis also applies to the (two-peak) antisymmetric states with M > 4. As we demonstrate in Sect. VI, the numerical solutions of both the Bogoliubov equations and the nonlinear time evolution of these systems confirm these predictions.

V. TRANSVERSE JOSEPHSON VORTICES
In this section we compare the stability predictions of the linear analysis with the numerical solutions of the GP equation (1) for representative sets of parameters. Just for showing two different examples we consider one case with zero axial momentum K x = 0 and another with K x = 0. As we will see, the linear predictions coincide with the numerical results of the nonlinear dynamics, and, in particular, stable transverse Josephson vortices The Bloch waves have quasimomentum K k for the following values of k = 0, ±1. Let us analyze the solutions for each particular case: • k = 0 : All the components share the same wave function ψ 0 = √ n exp (−iµ 0 t/ ), with the chemical potential µ 0 = gn − Ω. By using the parameters α 0,p = {0 , 3/2, 3/2} and β 0,p = {0, 0, 0}, with p = 0, ±1, the dispersion curves, ω 0,p , are These expressions are plotted in the bottom panel of Fig. 4 for a system with Ω = 0.2 gn.
• k = 1 : The wave function in each BEC component is: which configure a discrete anti-vortex of circulation Γ 1 = −3 √ 3 /2m around the discrete y-direction. Analogously, the Bloch waves with k = −1 correspond to an (energetically degenerate) vortex with opposite circulation Γ −1 = −Γ 1 . The chemical potential is µ 1 = gn+ Ω/2 and the parameters α 1,p = {0 , −3/4, −3/4} and β 1,p = {0 , 3/4, −3/4}. The resulting dispersion curves, ω 1,p , are: The upper panel of Fig. 4 depicts these expressions for a system with Ω = 0.2 gn. For ω 1,±1 , the negative signs under the square root indicate the presence of instabilities on the vortex state. Even at zero axial momentum ζ qx = 0, the system is unstable at low coupling Ω/2gn < 2/3. This is the case plotted in the upper panel of Fig. 4. Its typical decay is shown in Fig. 5, which illustrates the real time evolution of the system after imprinting a small random perturbation on the stationary state. The data, axial phases (top panels) and axial densities (bottom panels) for each component, have been obtained from the numerical solution of GP Eq. (1) with periodic boundary conditions in the axial coordinate, given in units of the healing length ξ. As can be seen, the stationary configuration survives for a time lapse around t ≈ 10 ξ/c, beyond which soliton-like structures (tracing thick, dark paths on the axial density plots) appear, producing strong density and phase modulations.  Here we consider arbitrary axial momentum (K x ) states. The Bloch waves have k = 0, ±1, 2, and the excitation modes have also p = 0, ±1, 2.
• k = 1 : Now the chemical potential is µ 1 = gn + 2 K 2 x /2m, and α 1,p = {0, 0, 0, 0}, β 1,p = {0 , 1, −1, 0}. The wave functions are: which yield a discrete, transverse vortex of circulation Γ 1 = 4 /m. An interesting feature of this state is its stability, irrespective of the axial momentum, because the dispersion curves ω 1,p contain only real frequencies (see Fig. 6): We have also performed numerical simulations with the GP Eqs. (1) for the real time evolution of this state with Ω = 0.2 gn and K x ξ = 0.1 (same parameters as in Fig. 6) after being perturbed with random noise. The simulations confirm the dynamical stability of this state, since the initial configuration (53) keeps robust against the perturbations.
• k = 2 : This state lies at the edge of the Brillouin zone, having maximum chemical potential µ 2 = gn + 2 K 2 x /2m + Ω, and parameters α 2,p = {0 , −1, −1, −2} and β 2,p = {0, 0, 0, 0}. The wave function in each BEC component is: In this state the Josephson circulation vanishes, Γ 2 = 0, and the system presents a sequence of π−Josephson junctions which are unstable. This feature is captured by the linear dispersion, which shows two unstable branches Again our numerical simulations with the timedependent GP Eqs. (1), for such a state with Ω = 0.2 gn and K x ξ = 0.1, confirm the linear prediction and show the decay of the initial, constant density state.

VI. NONLINEAR DYNAMICS OF LOCALIZED STATES
We study the nonlinear dynamics of the localized states with one or two density peaks in the stack of constant Despite the small coupling, Ω = 4 × 10 −3 gn, the system is dynamically unstable. The BECs j =1, 2, 3 correspond to the panels from left to right. The initial state has been seeded with a random Gaussian perturbation (less than 1% in amplitude).
density BECs. First, we numerically solve the Bogoliubov equations in order to check the dynamical stability of the corresponding stationary state. Next, we perform the real time evolution with the GP Eq. (1) of this state after adding perturbative noise.
In the simplest stack with M = 3, the antisymmetric states are unstable. To illustrate a typical decay process, Fig. 7 shows the real time evolution of an antisymmetric state with small coupling, Ω = 4 × 10 −3 gn, and zero axial momentum, K x = 0. The data have been obtained from the numerical solution of GP Eq. (1) with periodic boundary conditions in the axial coordinate. The axial length is L x = 155 ξ. As can be seen, the initial nodal strand (middle panels in Fig. 7) remains unpopulated during the whole evolution, and its phase is essentially undefined. The decay process is qualitatively different to the Bloch wave case presented in Fig. 5. The asymmetric state shows robust features of structural stability, keeping the density pattern across the stack.  Fig. 7. Only the peak-density strand (middle panels) and the two adjacent strands (side panels) are shown, with the other strands having practically zero density.
On the contrary, the symmetric state with M = 3, for the same parameters used above, is stable against perturbations. For larger stacks (we have made calculations up to M = 11), our numerical results show that both the antisymmetric states (with two density peaks) in stacks with M ≥ 4, and the symmetric states (with one density peak) in stacks with M ≥ 3 are also stable for the mentioned small coupling. However, the stability is lost at higher coupling values (at Ω 1 × 10 −2 gn for the parameters mentioned before). As a case example of stability, the time evolution of a symmetric state in a stack with M = 10 components is shown in Fig. 8. Only the peak-density strand and its nearest neighbours are shown, since the other components have a practically null density. The initial, t = 0, state has been seeded with perturbative noise. As can be seen, the initial localized configuration is robust against the perturbations. Due to the strong density localization, the dynamics is insensitive to the change in the boundary conditions. Our results show that a one-peak state with open-boundary conditions follows a dynamics which is indistinguishable from that shown in Fig. 8.

VII. DISCUSSION AND CONCLUSIONS
The rich phenomenology presented by the stacks of parallel Josephson junctions can be readily realized in ultracold atomic gases by means of 1D or 2D optical lattices [21,22]. These systems support nonlinear states whose dynamics reflects the interplay of continuous (along the axial x-direction of the BECs) and discrete (across the stack) features, and are promising candidates for pursuing technical applications with close similarities to superconducting devices. In this work, we have contributed to this goal and have demonstrated the existence and stability of simultaneous superfluid currents flowing through both directions of a 2D stack. While the translation invariance along the x-axis allows for the excitation of axial-momentum eigenstates, the periodic arrangement of Josephson junctions induced by the linear coupling permits transverse Bloch waves carrying tunneling supercurrents. If the stack shapes a closed loop, these Josephson currents around it yield non-regular vortices whose circulation is a generic non-integer multiple of h/m.
The dispersion relations of the transverse Josephson vortices have been obtained from the analytical solution of the linear Bogoliubov equations for the condensate excitations, and compared against the nonlinear time evolution of these states as given by the numerical solution of the Gross-Pitaeskii equation. In all the cases, the subsequent nonlinear dynamics is consistent with the stability predictions of the linear anaysis.
For the sake of comparison with the usual coupled-sine-Gordon-equation model for coupled superconductors, a further linear analysis of the transverse Josephson vortices has been performed in the hydrodynamic limit. As a result, we have derived linear wave-like equations for the relative phases and densities of the BEC components that resemble the mentioned model in the limit of small coupling.
We have also shown that the Josephson supercurrents are suppressed in steady states that break the symmetry of the discrete lattice and can present a strong localization across the stack. These nonlinear states belong to continues families of solutions to the Gross-Pitaevskii equation that can be tracked up to the non-interacting regime, where they are linear superposition of degenerate Bloch waves with opposite quasimomentum. Among these families, the gap-soliton-like states showing one or two dominant density peaks find dynamical stability in finite systems within a small coupling regime.
The exploration of different topologies in the stack, or the effect of exposing the system to synthetic gauge fields [21], stand out as interesting ways of extending the present work that will be reported elsewhere. ACKNOWLEDGMENTS M. G. acknowledges financial support from Ministerio de Economía y Competitividad (Spain), Agencia Es-