Josephson Junctions and AdS/CFT Networks

We propose a new holographic model of Josephson junctions (and networks thereof) based on designer multi-gravity, namely multi-(super)gravity theories on products of distinct asymptotically AdS spacetimes coupled by mixed boundary conditions. We present a simple model of a Josephson junction (JJ) that exhibits the well-known current-phase sine relation of JJs. In one-dimensional chains of holographic superconductors we find that the Cooper-pair condensates are described by a discretized Schrodinger-type equation. Such non-integrable equations, which have been studied extensively in the past in condensed matter and optics applications, are known to exhibit complex behavior that includes periodic and quasiperiodic solutions, chaotic dynamics, soliton and kink solutions. In our setup these solutions translate to holographic configurations of strongly-coupled superconductors in networks with weak site-to-site interactions that exhibit interesting patterns of modulated superconductivity. In a continuum limit our equations reduce to generalizations of the Gross-Pitaevskii equation. We comment on the many possible extensions and applications of this new approach.


Introduction
There has been recent interest in the potential applicability of AdS/CFT-inspired methods to traditional condensed matter problems, which are not amenable to a weak coupling quasiparticle description. This has led to the formulation and study of a large variety of models in classical gravity which aspire to capture holographically the characteristic features of some condensed matter system. A notable example is [1,2] which provides a simple gravity dual for an s-wave superconductor. More recent developments in this subject are reviewed in [3,4,5]. Ultimately one hopes that holographic techniques will provide a new efficient description of high-T c superconductors that goes well beyond the BCS theory. This simple sinusoidal relation can be substantially different in other (more general) types of junctions (see for example [7,8,9,10,11]).
Another reason to be interested in JJs is the nature of high-T c superconductivity itself.
Many high-T c superconductors enjoy a layered structure [12] that can be viewed as a natural stack of atomic scale intrinsic JJs [13] with interlayer spacing of about 15.5Å. In fact, such high quality SIS-type intrinsic JJs can be fabricated [14,15] and pose as attractive candidates of cryoelectronics.
Therefore, assuming holography can provide a new window to high-T c superconductor physics, there is an obvious interest to construct the holographic dual of Josephson junctions and more generally the dual of Josephson junction networks (JJNs).
A first step towards the construction of a holographic SNS junction has been taken recently in Ref. [16] (see also [17,18,19,20,21,22] for related setups). In this approach one is looking for solutions of the equations of motion of the standard holographic superconductor setups that are inhomogeneous across one of the field theory directions (the direction along which the SNS stack is arranged). In the present paper we will propose a distinctly different way to construct a holographic Josephson junction (one that is not necessarily restricted to SNS types).
The basic idea is to view each of the superconductors that compose the junction as a separate supergravity (or superstring) theory on its own asymptotically AdS spacetime and to model the weak link between them as a mixed boundary condition that relates the boundary conditions of the condensing symmetry-breaking field on one spacetime to the boundary conditions of the condensing symmetry-breaking field on the other spacetime.
This may look like a contrived operation on the gravitational side but it is a rather natural one from the perspective of a dual large-N quantum field theory. On the field theory side this operation amounts to a multi-trace deformation that involves products of single-trace operators from both theories. This multi-trace interaction is the only term that mediates interactions between the two theories. The single-trace operators that appear in these interactions are charged under the broken U(1) symmetries and thus mimic naturally the charge quantum tunneling effects that are present in a JJ.
Using multi-trace interactions to model the tunneling effects is natural in a regime where the mass scales of the modes that mediate the interactions between the two superconductor field theories in the full system are large compared to the typical energy scales that we consider. In that case one can integrate out these higher mass modes to obtain an effective theory at low energies. The separate gauge invariance in each of the two boundary theories implies that the effective interactions between them can only be of the multi-trace type. Abstracting from this picture the main features of the bulk description we will proceed to employ them freely in more generic situations where explicit knowledge about the boundary description is very limited or altogether absent.
The precise ingredients of our construction are presented in section 2. A simple characteristic example of a holographic JJ at zero temperature is discussed in section 3. We show that the standard sine expression for the Josephson current (1.1) is naturally reproduced in this model in a few lines and determine I max explicitly in terms of the parameters of the system. We briefly comment on the extensions and modifications that can alter this standard current-phase relation. We also discuss how this system differs from a typical Josephson junction and what kind of extensions can be used to describe more typical Josephson junctions.
Another appealing feature of the above approach is its versatility in describing very general configurations of Josephson junction networks. JJNs is a much studied subject with diverse applications. The great wealth of possibilities that they pose and the great reliability of the fabrication technologies developed for their construction makes them a prototype of complex physical systems that exhibits a variety of interesting physical behaviors. Among their many applications: they are used widely as microwave sources (see e.g. [23]), they provide controllable settings to investigate properties of granular or high-T c superconductors [24], they are frequently used as model analogs of physical systems with complex dynamics. For instance, they have been used to model biologically realistic neurons [25].
As a first simple application of our proposal in this direction we consider in section 4 a holographic Josephson junction array that can also be viewed as a special case of the honey-comb network of Ref. [26]. Using the gravitational description we find that the Cooper-pair condensates are described by a discretized Schrödinger-type equation, which has been studied extensively in the past (for a special value of one of our parameters) in radically different condensed matter and optics applications [27]. Using well known facts about this equation we show that the system in question exhibits complex behavior that includes periodic and quasiperiodic solutions, chaotic dynamics, solitons and kink solutions.
In our setup these solutions translate to one-dimensional configurations of holographic superconductor layers that exhibit interesting patterns of modulated superconductivity. In a continuum limit, the discretized Schrödinger equation becomes naturally a generalization of the Gross-Pitaevskii (GP) equation, a well-known long-wavelength description of superfluids. In this limit we recover some of the previously discrete solutions analytically.
We conclude in section 5 with an outline of further possible extensions and applications.  Figure 1: Examples of simple networks. In our context, these graphs will represent large-N CFTs coupled by double-trace interactions.
between the CFT i and the CFT j corresponds to a Lagrangian interaction of the form It will be useful to set up a convenient notation to denote graphically networks constructed in this way. For CFTs coupled by double-trace interactions, like that in eq. (2.1), we will denote the corresponding link by a single line. Such a link is undirected and joins two different CFTs, or circles back to the same CFT. The latter denotes that the corresponding CFT has itself a double-trace interaction turned on. Hence, Fig. 1(b) exhibits a simple network of two CFTs, call them CFT 1 and CFT 2 , with the Lagrangian interaction Similarly, Fig. 1(c) exhibits nine CFTs coupled pairwise by double-trace interactions in all possible ways. This example, analyzed in [28], plays a role in quenched disorder calculations using the replica trick [29]. Networks formulated in this way can have different types of links. We may use a double-line notation to denote the more general possibility of links mediated by multi-trace interactions. Figs. 2(a) and 2(b) exhibit a triple-trace and four-trace coupling respectively.
In network literature the corresponding graphs are sometimes called hypergraphs. The use of such more general couplings opens up many interesting possibilities. Most of our discussion in this paper will be focused, however, on networks with double-trace couplings only.
Being comprised by sites that correspond to interacting QFTs the above networks have in general a rich and complicated internal structure. This structure can evolve in time, vary in space, or vary from site to site. Moreover, under renormalization group flow the network graphs may change with the appearance of new links or even new sites. Indeed, it is well-known that multi-trace couplings are naturally generated under renormalization (see for instance [28]). Hence, renormalization effects in field theory can affect the quantitative features of the links. They can also change the number of sites in the following way. Many supersymmetric QFTs, like the N = 4 super-Yang-Mills theory, have a Coulomb branch.
At generic points of this branch the original gauge group is Higgsed into product gauge groups and in the far infrared one is left with a product of QFTs giving rise to a network with more sites than those of the UV theory. This is an example of an RG flow with a different number of sites in the UV and IR. In what follows, when we draw a graph representing a network we will implicitly assume that this description refers to the bare UV Lagrangian of the corresponding theory.

Networks of asymptotically AdS spacetimes
There are situations where a large-N QFT (typically at strong coupling) has a dual description in terms of a supergravity theory on an asymptotically AdS background. Accordingly, a network of such QFTs has a dual description as a multi-gravity network where each site is some supergravity theory on an asymptotically AdS background and each link is a mixed boundary condition for supergravity fields residing on different space-times [30,31,28] (see [32] for a stringy setup that involves multi-string theory networks and [33] for a brief review of the main idea and its implications for massive multi-gravity). Let us recall how such a multi-gravity network comes about in the AdS/CFT correspondence. For clarity, we will restrict to the case where the boundary QFTs are conformal.
Before the addition of multi-trace links, the CFT i s are independent field theories that do not talk to each other and the dual gravitational theory is a direct product of supergravity theories on product AdS space-times of the form k i=1 [AdS i ⊗ M i ], where M i is the internal manifold of the spacetime of the dual of CFT i .
Assume that the single-trace operators O i stated above are scalar operators 1 with scaling dimension ∆ i . The AdS/CFT correspondence maps each O i to a dual scalar field ϕ i with AdS asymptotics r i is the radial distance in the i-th AdS space that corresponds to CFT i . We use conventions where the i-th AdS boundary lies at r i → ∞. In later applications, we will assume ∆ i < d 2 , in which case the α i 's should be interpreted as the vacuum expectation values (VEVs) of the dual operators and the β i 's as the sources.
In this framework, and to leading order in the 1/N expansion, it is well known that a multi-trace coupling in the boundary CFT imposes mixed boundary conditions to the asymptotic coefficients α i , β i of the form Combined with requirements of regularity these boundary conditions fix completely the profile of the bulk solution. The inter-theory coupling induced by the relations (2.5) leads to a network of supergravity theories that can exhibit interesting collective phenomena.
We will explore these phenomena in section 4.

Designer multigravity
Let us recall how multi-trace deformations affect the profile of the bulk asymptotically AdS solution in a single theory; a subject that usually goes under the title of designer gravity [34]. In this paper we want to consider the natural extension of this framework to multi-AdS spaces, which we will suggestively call 'designer multigravity'.
As one of the simplest illustrations of the idea consider the 'network' (2.6) which contains a single vertex with a self-adjoining link. The graph (2.6) exhibits a doubletrace link, but we can equally well consider any multi-trace link. This theory has been the main focus of most previous investigations of multi-trace interactions in the AdS/CFT correspondence and designer gravity. For self-completeness and in order to set the notation, we briefly review some of the most pertinent properties of this system.
At the single site of this network resides a d-dimensional large-N CFT with a supergravity dual. Assume that the CFT has a single-trace complex scalar operator O with scaling dimension ∆. This operator is dual to a charged bulk scalar field ϕ. To further simplify the discussion we will also assume that we can consistently reduce the dynamics of the dual bulk supergravity to a (d + 1)-dimensional Einstein-Abelian Higgs model of the where θ is the phase of ϕ, namely ϕ = |ϕ|e iθ , and G, J, V are model-dependent functions of |ϕ|. A is an Abelian gauge field in the bulk and F its field strength. q is the U(1) charge of the field ϕ.
We are interested in asymptotically AdS d+1 solutions of this system. In units where the AdS radius is set to one, the potential V (ϕ) has the small-ϕ expansion When the mass m lies within the range the dual operator O can have two possible scaling dimensions (2.10) We will assume that our theory lies in the window (2.9) and will pick O to have the smaller scaling dimension ∆ ≡ ∆ − that satisfies the inequality d−2 2 < ∆ − < d 2 (the lower bound in this inequality is the standard unitarity bound in field theory).
Near the asymptotic boundary, r → ∞, the metric is that of AdS d+1 ds 2 ≃ r 2 dx µ dx µ + dr 2 r 2 (2.11) and the scalar field ϕ exhibits two independent branches The boundary condition for a generic smooth function W(α) corresponds at the boundary CFT to the multi-trace deformation [35] δL = W(O) . (2.14) When we solve the bulk equations of motion we are looking for solutions that respect the boundary conditions (2.13). In addition, we require that these solutions are everywhere regular. It turns out that regularity imposes an extra constraint on the asymptotic coefficients α, β. We will denote this additional relation as where W is a function determined by the specific dynamics of the theory, e.g. the details of the scalar potential V (|ϕ|) in the bulk action (2.7).
Combining the boundary condition (2.13) and the regularity condition (2.15) we find that α and β are completely fixed and determined as a solution to the equation which therefore can be viewed as an extremum of the function V.
For example, for boost invariant, planar configurations with vanishing gauge field one finds a solution with an acceptable naked singularity at r = 0 that has The existence of this solution and the precise value of the parameter s depends on the details of the bulk potential V (ϕ) (see [36] for additional information). For most potentials, s turns out to be a positive number. The energy density of this solution is One can also consider analogous solutions with spherical topology. It is possible to show that the energy (2.19), (2.18) provides a lower energy bound for all of these solutions [36]. Moreover, in the case of vanishing gauge field, one can argue [36] that • the theory with boundary conditions β = dW dα has a stable ground state provided the function V has a global minimum V min , and • that the minimum energy solution is the spherical soliton associated with V min .
Note that because of the presence of W, it is possible to have a stable ground state even for functions W that have no minimum. In general, this minimum involves a condensate of the charged scalar that higgses the corresponding U(1).
What we have discussed so far applies to the case of zero temperature. By studying hairy black holes in the bulk it is possible to generalize the discussion to non-zero temperature [37]. It is also possible to consider the case of non-vanishing charge density [38].
It is not hard to generalize this discussion to the case of multiple CFTs coupled together by a multi-trace interaction of the general form W(O 1 , . . . , O k ). In this case (and to leading order in the 1/N expansion), the bulk equations of motion are the same as before in each (super)gravity member of the product, but the boundary conditions change. For example, the bulk boost-invariant planar solutions in the i-th spacetime still retain the form (2.17) and each β i is still given by an equation of the form The new ingredient, responsible for the coupling between different AdS theories, lies in the mixed boundary conditions (2.5). The analog of equation (2.16) that determines the VEVs In general, these equations lead to a non-linear discrete system which can exhibit intricate behavior. This behavior includes solutions that can be periodic, quasi-periodic, chaotic or even soliton-like with energy pinned around a central site. Examples of each of these behaviors will be discussed in section 4.

A holographic superconductor with vanishing charge density
A novel type of holographic superconductor with symmetry breaking induced by doubletrace deformations was recently proposed in [38]. Since this setup will provide the basic building block of the discussion to come, it will be beneficial to recall some of its main properties.
Returning to the single-site example of (2.6) consider the case of a double-trace defor- implemented with the use of a single-trace operator O with scaling dimension ∆. The function V in (2.16) becomes in this case (for boost-invariant planar solutions) The extrema of this function obey the algebraic equation Assuming s > 0, there is one or two possible solutions to this equation depending on the sign of the constant g. For g > 0, the only solution is α = 0, which is a minimum and does not exhibit any condensate of the field ϕ. For g < 0, there are two possible solutions (2.27) The first one, α 0 , is a local maximum of V and therefore an unstable vacuum of the system.
The other solutions, α θ , are degenerate stable minima labeled by an angular variable θ.
The non-vanishing condensate of the charged scalar field ϕ in this case higgses the U(1) gauge symmetry and leads to a new type of holographic superconductor. As a concrete illustration of the idea consider a pair of two identical holographic superconductors of the type described in the previous subsection 2. 4. The parameters g, s, δ, d, ∆, which will be treated here as phenomenological parameters of the model, are chosen to be common in the two systems. Each of them has a charged scalar field O i (i = 1, 2) and a corresponding dual complex scalar field ϕ i with asymptotics (2.12). We can think of O i as the 'Cooper pair' operator in the holographic superconductor CFT i . Holographically, the VEV of these operators are given by the leading branch coefficients α i in the asymptotic expansion (2.12). We assume that both operators have the same scaling dimension ∆ < d 2 and employ the alternative quantization for the dual scalar fields ϕ i in the bulk AdS spacetime.

A holographic model of
In order to mediate an interaction that exchanges charge between the two theories we couple them via a double-trace interaction of the form h is a real number controlling the strength of the interaction and ϑ is an angular variable.
The latter is an additional tunable parameter of the interaction whose physical meaning will become clear in a moment.
In the bulk this coupling implies mixed boundary conditions. One can determine the VEVs of the dual operators O i by solving the scalar-gravity equations of motion.
Equivalently, one can extremize the function V (see eqs. (2.22), (2.23)). In the case at Extremizing with respect to α 1 and α 2 we obtain two algebraic equations which determine the α i 's uniquely in terms of the parameters g, h, s and δ up to an overall phase common in α 1 and α 2 . The relative phase between α 1 and α 2 is fixed in terms of ϑ For instance, when δ = 4 3 the solutions are In all cases Which solutions survive in cases (2) and (3) depends on the specific range of the parameters s, g, h. Assuming s > 0 we obtain • |h| < g: vacuum (1), 3 This is consistent with ∆ ∈ d 2 − 1, d 2 iff ∆ < 1, which is possible only if d = 2, 3 (the main cases of interest in this paper).
Diagonalizing the Hessian of the energy functional V (see eq. (2.19)) we find that (1) has the eigenvalues Hence, when phases with a non-zero condensate exist (g < |h|) the vacuum (1) of the normal phase is unstable (one of the eigenvalues of the Hessian is negative). In that case there is always at least one superconducting vacuum that is stable. Such a phase provides a holographic model of two weakly interacting superconducting materials at zero temperature.

Josephson current
A characteristic feature of conventional Josephson junctions is a transverse supercurrent I which is related to the condensate phase difference ϑ 12 in the following way In this section we have considered a system of two sites that represents two infinitely thin layers of superconducting material at zero charge density coupled through a weak link expressed by an interaction of the form (3.1), which can be suggestively rewritten as with the definition W E is an interaction that mediates no interlayer charge transfer. In contrast, W Jext is an interaction based on the charge-transferring operator The coupling A of this operator in W Jext , namely can be interpreted as a new interlayer background gauge potential component. 4 In this sense, for non-vanishing h, ϑ our two-site system lies in an external transverse gauge field and J ext is an externally imposed current.
The total current running across the sites of the junction can be determined in the standard way from an infinitesimal relative U(1) gauge transformation of the action. For an interaction of the form (3.9) the infinitesimal transformation gives the Lagrangian variation The second equality is the discretized version of the gradient ∂J across the interlayer direction. In the third equality we used charge conservation to set J site 1 = −J site 2 = J tot . Consequently, This is precisely what one expects. Since our system is kept at zero charge density, in the equilibrium state charge cannot flow across the junction between the two sites. As a result, irrespective of the initial configuration, once the interaction (3.9) is turned It is interesting to consider the VEV of the externally forced current operator J ext in the new vacuum. At leading order in the 1/N expansion Using eq. (3.6) we obtain J ext = I max sin ϑ = −I max sin ϑ 12 (3.18) with In the algebraically simple case of δ = 4 on the vacua (2) ± and on the vacua (3) ± (interestingly in this case I max is independent of g).
The vanishing of the total interlayer current J tot implies that the backreaction has created an equal and opposing 'Josephson current' J josephson (due to the condensate phase difference), which cancels the contribution of the externally imposed J ext . Specifically, J josephson = − J ext = I max sin ϑ 12 (3.22) in agreement with the expected sine law (3.8).
It is clear that the system we have just described is a peculiar Josephson junction unlike the typical Josephson junction commonly discussed in the literature and engineered in the laboratory. In contrast to our system in a typical junction the superconductor components have finite spatial thickness in the transverse junction direction, charge can flow in this direction, the condensate phase difference is dialed by choice (and not determined dynamically) and a Josephson current arises without having to apply a gauge field across the weak link. In section 4 we will describe how such more conventional junctions can be engineered in our framework as simple extensions of the setup we have described in this subsection.

Possible extensions and other current-phase relations
We briefly comment on a few possible extensions of this framework that may ultimately lead to different current-phase relations.
A possible alternative is to consider the same superconductor models as above, but modify the double-trace interactions. For example, one could consider other types of double-trace interactions, e.g. instead of (or in addition to) (3.1) a double-trace deformation of the form This will modify the equations (3.3) and the Josephson current (3.22) in an obvious fashion.
Another possibility is to consider higher multi-trace deformations. An example of a triple-trace interaction is Finally, it is interesting to consider coupling different types of asymptotically AdS theories. For example, one can try to couple the holographic s-wave superconductor of Ref.
[2] to the holographic p-wave superconductor of Ref. [39]. In the s-wave superconductor it is a bulk complex scalar field, dual to a complex scalar operator, that condenses. In the p-wave superconductor there is an SU(2) gauge field in the bulk and the U(1) ⊂ SU (2) is broken by the condensation of the remaining two components of the SU (2)  For simplicity, we are assuming the same constants g, h, ϑ, δ for all sites and links. Accordingly, the potential function V that we have to extremize in gravity to obtain the condensates α n is The extrema of this function are sequences of complex numbers obeying recursion relations with rich features. The solutions can be periodic, quasiperiodic (chaotic) or solitonic, and provide interesting new examples of spatially modulated, namely lattice site-dependent, superconductivity.
We will organize the discussion according to the number of boundary conditions.

No boundaries
In the case of no boundaries, the extremization of the potential function V (4.2) with respect to all the α i 's gives the recursion relations gα n + h e iϑ α n−1 + e −iϑ α n+1 + s 2 α n |α n | δ−2 = 0 , n ∈ Z . Setting we can recast (4.3) into the form gϕ n + ϕ n−1 + ϕ n+1 +s 2 ϕ n |ϕ n | δ−2 = 0 . The generic solution of these equations is parameterized by two complex numbers (these could be for instance the values of ϕ 0 , ϕ 1 at the vertices 0 and 1). There are various ways to analyze this system of equations. In fact, such equations have appeared in the past in a variety of applications and discussions of discrete dynamical systems. A notable example, that arises for the special case of δ = 4, is that of the discrete non-linear This equation is a discretized version of the Schrödinger equation with a non-linear quartic potential. Setting ψ n ≡ ϕ n e iEt (4.7) we recover our set of equations (4.5) with E =g and δ = 4.
The DNLS equation has a long history (for an extensive review and references we refer the reader to [27]). In solid state physics it appeared first in the context of the Holstein polaron model for molecular crystals [40]. In an optics context, DNLS describes wave motion in coupled nonlinear waveguides. In the basic nonlinear coupler model introduced in [41] two waveguides made of similar optical material are embedded in a different host material. DNLS, and generalizations, have also appeared in studies of nonlinear electrical lattices [42].
As a concrete algebraically simple example, in the rest of this section we will concentrate on the case of δ = 4. Qualitatively similar results are expected for generic δ. The linear stability analysis is treated for generic δ in appendix A. For most purposes we follow closely the analysis of [27] which we recommend for additional details.
In certain regimes, e.g. wheñ the solutions are regular and bounded and the Lyapunov exponent vanishes [27]. Recall that the original parameter s that appears in the AdS/CFT context (2.18) is positive, but s = s h can be both positive or negative. Moreover, one can tune s freely by adding on the field theory dual quartic-trace interactions of the form |O| 4 .
The general structure of the space of solutions is organized by a hierarchy of periodic orbits surrounded by quasi-periodic orbits [27]. The periodic orbits can be traced on the intersection of any two symmetry lines S n 0 = M n S 0 , S n 1 = M n S 1 , n = 0, 1, . . . [43,44], where the fundamental symmetry lines are defined as S 0 : z = 0 , (4.18a) For example, in the intersection of the lines S 0 and S 1 one locates the fixed point (x = x * , z = 0) with The linear stability of an orbit with period q is conveniently characterized by the value of Greene's residue [45] where DM is the linearization of the map M. The period orbit is linearly stable when 0 < ρ < 1 (elliptic periodic orbit) and unstable when ρ > 1 (hyperbolic with reflection) or ρ < 0 (hyperbolic).
In the regime of elliptic stability the fixed points (x * , 0) form the largest basins of stability among all elliptic orbits. These stable orbits, which include both periodic and quasiperiodic solutions, encircle the fixed point forming the main island on the map plane. The quasi-periodic orbits, which lie on closed curves (the Kolmogorov-Arnold-Moser (KAM) tori), densely fill the island. One can show that the map M is topologically equivalent to an area-preserving map ensuring the existence of such KAM-tori near the symmetric elliptic fixed points [46].
An illustration of the main island of elliptic orbits around the fixed point (x * , 0) for g = 1.6,s = 0.1 can be found in plot (a) of Fig. 4. Outside this island one finds regular quasiperiodic orbits and orbits that diverge.
Interesting changes in the structure of the solution space can occur as we vary the parameters of the system, here (g,s). In particular, the residue of periodic orbits can change. When the residue changes from a positive value to a negative value then a tangent bifurcation occurs where an elliptic point converts into a hyperbolic point. Whenever the residue exceeds the value of one from below a stable elliptic orbit converts into an unstable hyperbolic point with reflection accompanied by the creation of two new stable elliptic points. This kind of bifurcation is known as period-doubling bifurcation -a 1-period orbit (fixed point) converts into a 2-period orbit. The new elliptic orbits remain stable until another period-doubling bifurcation occurs. After a cascade of such bifurcations local chaos appears.
In our system, the residue ρ in eq. (4.21), is always less than one assumings > 0, In that case, only tangent bifurcations can occur and global chaos can arise through the so-called resonance overlap. Period-doubling bifurcation can instead occur wheng < 0.
Specifically, in the rangeg < −s 2 < 0 a new period-2 orbit is generated at the points The map M acts on these points by sending (x ± , 0) → (x ∓ , 0). The residue of the new solution is positive As we decreaseg further the period-2 orbit loses its stability and a new period-doubling bifurcation occurs which gives rise to a period-4 orbit. This cascade terminates at a critical parameter (see [27] and references therein) which is called the accumulation point. At that point local chaos appears.
A characteristic example of period-doubling bifurcation is depicted in plot (b) of Fig.   4. Atg = −0.62,s = 1, the originally stable fixed point of plot (a) has become unstable and two new elliptic fixed points have been created giving rise to elliptic period-2 orbits.
Besides the issue of linear stability, that was discussed above, one can also ask about the local and global thermodynamic stability of the above solutions. Local thermodynamic stability requires a positive definite Hessian of the multi-gravity energy functional (2.23) (see also (2.19)). It would be interesting to examine the extent to which local thermodynamic stability is equivalent to linear stability. On the other hand, global thermodynamic stability implies that the more stable solutions have less energy. In the following subsections 4.1.2 and 4.2 we will see that chains with a finite number of sites have a discrete set of solutions. In that case, the solutions with the minimum energy are thermodynamically favored. We hope to return to a more detailed examination of these issues in the future.

Periodic boundary conditions
From the above discussion it should be clear that, for given values of the parameters and period, one is left with at most a discrete finite set of solutions to the recursion relations

One or two boundaries
If there is a boundary, say at n = 0 with n valued only on non-negative integers, the n = 0 version of the equation (4.8) is modified tõ In that case, the whole solution is fixed by the choice of one parameter, for example ϕ 0 .
If the chain has finite size and there is also a second boundary, then the analog of (4.25) at the second boundary will fix ϕ 0 as well and the solution will be discretely unique and expressed completely in terms of the parameters of the system g, s, h. A specific example of this situation appeared in the dimer case of section 3.
A special consequence of eq. (4.25) is the fact that ϕ n are all real-valued up to a common n-independent phase. Equivalently, the phases θ n in (4.9) are all equal modulo π. This property can be deduced by using the boundary equation (4.25) to compute the conserved (i.e. n-independent) current J (4.14b). One finds J = 0 from which the above statement follows immediately.
Taking care to satisfy the conditions from the boundary equations one can proceed as before to analyze the solutions. Depending on the precise parameters and the boundary value ϕ 0 one finds again in the case of a semi-infinite chain regular or diverging solutions.
The regular solutions can be periodic or quasi-periodic. In the case of a chain with finite size only a discrete subset of the above regular solutions survives. We are not aware of any tractable analytic method that determines these regular solutions for generic sets of parameters.

Solitons and pinning of superconductivity
It is known that the DNLS equation admits also another interesting kind of solutions: soliton and kink solutions. A detailed discussion of these solutions in the real domain and related references can be found in [27] whose main points can be summarized briefly as follows. A priori one might expect that the DNLS equation does not admit such solutions. Soliton-like solutions are typically associated to integrable systems and DNLS is not integrable. It exhibits irregular chaotic behavior which in principle may prevent perfect localization. Nevertheless, it can be shown that non-integrability and discreteness appropriately combine to make such solutions possible.
The solutions of interest have the following characteristics. They are solutions where the amplitude ϕ n is exponentially localized around a single site, say at n = 0. Following the nomenclature of [27] one can distinguish between two situations: (1) Bright solitons: in this case, |ϕ n | > |ϕ n+1 | for n > 0 and |ϕ n | < |ϕ n+1 | for n < 0 with lim |n|→∞ |ϕ n | = 0.
In our context, where each site labeled by an index n, models a (1+1)-or (2+1)dimensional layer of a superconducting material, such configurations would correspond in case (1) to situations where in a chain of layers the interlayer interactions work in such a way that energy and superconductivity are strongly localized around a central site. In case (2) the opposite happens. Energy and superconductivity are modulated in such a way that they are almost uniform along the chain except around a central site where the condensate vanishes as it changes sign and superconductivity becomes very weak.
In a continuum limit (see next subsection) the above configurations appear to reconstruct a junction of three materials with one dimension higher. In case (2) we recover a configuration that is very similar to the dark soliton of [17] and reminds of an SNS junction of (2+1)-or (3+1)-dimensional superconductors. From this point of view the configuration of case (1) resembles a junction of two materials in the normal state separated by a thin superconducting layer in the middle.
Since ϕ n are now real it is convenient to view the equations (4.8) as a two-dimensional map M : R 2 → R 2 by defining a new set of R 2 coordinates (x n , y n ) = (ϕ n , ϕ n−1 ). Then, The identification of the soliton-like solutions is closely related to the structure of the fixed points of this map.
The fixed points, which by definition obey the relation x = y, are located at The x ± fixed points exist only when sgn(g + 2) = −sgn(s). Greene's residue ρ for the fixed point at the origin is [27] ρ = 1 4 (g + 2) . (4.28) Consequently, for |g| < 2, we obtain 0 < ρ < 1 which implies that the origin is a stable elliptic fixed point encircled by stable elliptic orbits. When in additions < 0 the fixed points x ± are unstable hyperbolic fixed points.
For |g| > 2 we may distinguish between the following two cases: (i)g < −2,s > 0. In that case ρ < 0 and the origin becomes an unstable hyperbolic point.
The points x ± are stable elliptic fixed points.
(ii)g > 2,s < 0. The fixed point at the origin becomes unstable and through perioddoubling bifurcation a new period-2 orbit appears located on the line x = −y.
Before proceeding to explain the main idea underlying the existence of soliton-like solutions it will be useful to introduce some language which is common in the study of dynamical systems.
A set W is called an invariant manifold of a dynamical system if for any point x ∈ W the dynamical evolution of x for any amount of time t continues to belong in W. Every fixed point p comes with its invariant manifolds. Such manifolds are called stable, and denoted as W s (p), if all points that belong on them approach asymptotically the fixed point p under dynamical evolution (namely p is an attractor on W s (p)). In contrast, an invariant manifold is called unstable, and denoted as W u (p), if all points that belong to it move asymptotically away from the fixed point p under dynamical evolution (in other words, p is a repellor on W u (p)).
For generic non-integrable maps it is known that the stable and unstable invariant found in [27] (see, for instance, Fig. 7 in [27]).
In case (ii) withg > 2,s < 0 one obtains a similar soliton-like solution, but with the new feature that adjacent amplitudes ϕ n , ϕ n+1 have alternating signs. Such solutions are known as staggered solitons [47]. The solitons of the previous paragraph are also known as unstaggered solitons.
Finally, in the case of |g| < 2,s < 0 one can consider heteroclinic orbits connecting the two unstable hyperbolic fixed points x ± . These are kink solutions of the type (2) above.

Continuum limit and the Gross-Pitaevskii equation
It is interesting to consider the continuum limit of the chain configuration described above.
In this limit with the new parameters x, G, and S kept finite. In this limit the recursion relations (4.5) turn into the second order non-linear Schrödinger differential equation where ′ ≡ d dx . For δ = 4 we recover directly a well-known equation in the context of superfluids; the Gross-Pitaevskii equation [48,49] which gives a coarse-grained description of superfluids at long wavelengths. The GP equation is typically relevant for weakly-interacting Bose-Einstein condensates or strongly bound fermionic superfluids at low temperature. It is interesting that in our formalism, which has a radically different point of departure, the same description emerges naturally (as a candidate description of strongly coupled superconductor physics) out of a framework designed specifically to deal with layered structures.
It would be worth exploring further parallels that may exist between our formalism (and the generalizations of the GP equation that it suggests) and the known applications of the GP equation to superfluidity and superconductivity.
For simplicity and concreteness let us continue to concentrate on the case of δ = 4.
Two well-known solutions of this equation are: (1) Bright solitons: for G < 0, S > 0 (2) Dark solitons: for G > 0, S < 0 For real ϕ we can find a more general class of solutions expressed in terms of the Jacobi elliptic function sn(u|m) (4.33) C and x 0 are integration constants.
For a general solution of the recursion equations (4.3) and ϑ = 0, π this is a linkdependent current. It is periodically or chaotically modulated in periodic or quasi-periodic solutions. In real solutions, where ∆θ n = 0 mod π, J n−1,n = ±2r n−1 r n sin ϑ .  In the special case where ϑ = 0, π (i.e. when no external interlayer gauge field is applied), eq. (4.41) gives the site-independent current J n−1,n = ±2J . This current vanishes for a chain with boundaries, but can be non-zero in chains without boundaries, e.g. in a circular chain with periodic boundary conditions. This is a simple example of how the topology of the network can affect the qualitative features of the configuration.

Towards a typical Josephson junction
In subsection 3.2 we discussed the similarities and differences between a two-site system and typical Josephson junctions. Here we discuss how the two-site system can be extended to look more like the typical Josephson junction.
Assume we want to describe a junction of two superconductors in three spatial dimensions linked weakly across the third direction z at a two-dimensional interface. In previous subsections we described how to deconstruct (3 + 1)-dimensional layered superconductors from an array of (2 + 1)-dimensional holographic superconductors using linear AdS/CFT arrays. To construct an (un)conventional Josephson junction of two superconductors of this type a possible strategy is described in Fig. 6.

Outline of future directions
We have proposed a novel holographic way to model the physics of Josephson junctions using networks of (super)gravity theories on asymptotically AdS spacetimes coupled via mixed boundary conditions. One of the advantages of this approach, compared to previous holographic approaches, is the versatility by which it can incorporate many different types  (e.g. s-wave, or p-wave). For example, it will be interesting to define and study interlayer transport coefficients in such models. The ultimate goal is to explore the extent to which these constructions reproduce known phenomenological features of JJNs or layered superconductor physics. For example, it would be interesting to reproduce previously observed non-sinusoidal current-phase relations in unconventional JJs.

(b) Other network architectures and complex behavior
It has been pointed by many authors (see e.g. [26]) that the architecture of a JJN can have important implications on the physical properties of the system. It is interesting to explore other configurations and examine how they affect the collective and local properties of the sites. Networks with double-trace or higher multi-trace interactions can be constructed.
An example that has been studied previously in the condensed matter literature is the Y-Josephson junction (see for instance [50]). The architecture of a Y-Josephson junction network appears in Fig. 7.
It may also be interesting to explore the existence of vortex solutions in two-or threedimensional AdS/CFT lattices. This becomes even more interesting in view of the observed connection to the GP equation in the continuum limit. A recent discussion of vortex solutions to the GP equation (and a related AdS/CFT application from a different point of view) can be found in [51].

(c) Continuous limits and deconstruction
In subsection 4.3 we considered a continuum limit of a one-dimensional holographic JJN.
In this limit the number of sites is scaled to infinity with an appropriate scaling of the other parameters of the system to zero. For a specific set of parameters we recover in this limit the Gross-Pitaevskii equation. It would be interesting to explore further relations between the generalizations of this equation suggested by our formalism and known applications of the GP methodology in superfluidity and superconductivity.
In addition, it would be interesting to explore similar continuous limits of other JJNs with more complicated topology and different ingredients. Such limits may be used to simplify some aspects of the analysis of the network or in order to attempt a novel deconstruction of one or more extra spacetime dimensions (more comments on this aspect can be found in [28]).