May 18, 2007 15:30 WSPC - Proceedings Trim Size: 9.75in x 6.5in joseph

1

MULTISTABILITY AND MULTI 2π-KINKS IN THE FRENKEL–KONTOROVA MODEL: AN APPLICATION TO ARRAYS OF JOSEPHSON JUNCTIONS

Karl E. Kürten^{1} and Christian Krattenthaler^{2} 1Institut für Experimentalphysik, Universität Wien, Boltzmanngasse 5, A-1090 WIEN, Austria E-mail: karl.kuerten@univie.ac.at

2Fakultät für Mathematik, Universität Wien, Nordbergstraße 15, A-1090 Vienna, Austria

                      WWW: http://www.mat.univie.ac.at/~kratt

A regular ring of Josephson junctions, connected in parallel, is studied analytically and numerically. We show that, depending on the strength of the r-well cosine potential the energy landscape of the Hamiltonian can have of the order of r^{N} /N locally stable minima separated by large barriers specified by unstable saddle points. The counting problem for the degeneracy of the total energy is equivalent to a wellknown necklace problem in combinatorial mathematics. We also demonstrate that the distribution of the phase differences as well as the energy spectrum is fractal provided that the strength of the cosine potential is sufficiently strong.

Keywords: Style file; L^{A}TEX; Proceedings; World Scientific Publishing.

1. Introduction

The Frenkel–Kontorova model,^{1,2} also known as the discrete sine-Gordon model, describes a harmonic chain mimicking a layer of atoms in a spatially periodic potential playing the role of the substrate. The particles or atoms placed on a one domensional array interact via nearest neighbours and are subjected to a periodic external multi-well potential. The model was introduced more than 60 years ago with the intention to study crystal dislocation. Meanwhile it has been applied to numerous physical systems such as charge-density waves, magnetic spirals, disordered crystals, magnetic multilayers, and adsorbed monolayers. Due to competitions between two or more length scales or frequencies, locally stable spatially modulated structures arise, where the particles are non-trivially displaced from the sites of a reference lattice. To be more precise, the Frenkel–Kontorova model has become a favourite model for disordered systems and glassy materials because of its huge number of meta-stable states. More recently, the model has also found successful applications in the investigation of fluxon dynamics in systems consisting of superconducting wires interrupted by Josephson junctions, so-called Josephson arrays.^{1,3,4} Here, a

2

discrete chain of effective particles emerges when one considers the flux flow in discrete parallel arrays of weak links between superconductors. Note that these small Josephson junctions, in contrast to long Josephson junctions associated with the continuous sine-Gordon equation, can be modelled by Kirchhoff’s circuit equations. In this case, they describe an array of N discrete Josephson elements connected by a parallel resistance/inductance combination. Moreover, Josephson devices are a candidate for the realization of models of superconducting interference devices (SQUIDS) used for experiments on quantum manipulation of qubits on which quantum computation can be realized. On the other hand, Josephson devices are also a suitable playground for fundamental investigations of soliton dynamics, chaos, and fractal behaviour.

2. The model

We consider a regular ring of small Josephson tunnel junctions connected in parallel via superconducting leads. Specifying the gauge-invariant phase difference between the macroscopic wave functions across the i-th junction by the variable ϕ_{i}, and assuming periodic boundary conditions, the total energy of the system is given by the N-particle Hamiltonian

H(ϕ_{1}, ϕ_{2}, …, ϕ_{N}) = ∑_{i=1}^{N} (\frac{1}{2}(ϕ_{i+1} − ϕ_{i})^{2} + a(1 − cos ϕ_{i})) .  (1)

The first term defines the inter-particle energy, specified by harmonic nearestneighbour interactions, while the cosine potential in the second term stems from the supercurrent through the junction given by the Josephson relation. The control parameter a can be considered as a discreteness parameter, since for a → 0 the model can be shown to be equivalent to the continuous sine-Gordon equation, where analytic solutions exist. The control parameter a > 0 also specifies the height of the barrier which separates two wells of the cosine potential. The ground state with zero energy is obtained when all units are placed in the same minimum of the cosine potential such that these configurations are ferromagnetic. In this report we mainly focus on analyzing the structure of the distribution of the phase differences as well as the structure of the energy distribution. Another macroscopic quantity, the total phase M of the chain,

M = ∑_{i=1}^{N} ϕ_{i},  (2)

is also of interest in the context of experimental investigations of fluxon dynamics,^{3} where the fluxon corresponds to a 2π-kink of a quantum phase difference ϕ_{i} between the two superconducting electrodes of the junction. Here the behaviour of travelling kinks and solitons, describing Josephson vortices, each of them carrying a magnetic flux quantum, is studied experimentally as well as theoretically.^{3,4}

3

3. Few particle systems: Two junctions

Analytic results obtained for few particle systems are often a rather useful guide for the understanding of the many-particle system, since characteristic features occurring in large systems often show up at the one- and two-particle level, respectively. In order to get some insight into the structure of the energy landscape of an Nparticle system we first consider a small system consisting of only N = 2 interacting junctions. The energy takes the form

E(ϕ_{1}, ϕ_{2}) = \frac{1}{2}(ϕ_{1} − ϕ_{2})^{2} + a ((1 − cos ϕ_{1}) + (1 − cos ϕ_{2})) ,  (3)

and the necessary condition for equilibrium is given by the set of two coupled equilibrium equations

2(ϕ_{1} − ϕ_{2}) + a sinϕ_{1} = 0,  2(ϕ_{2} − ϕ_{1}) + a sinϕ_{2} = 0.  (4)

To give an illustrative example, we first restrict the values of ϕ_{i}, i = 1, 2, to the range (−3π, 3π), where the cosine potential has only three single wells. In this case, Eqs. (4) can have up to nine distinct solutions representing locally stable minima. The homogeneous configurations −(2π, 2π), (0, 0), and (2π, 2π), specify the three degenerate ground states with zero energy. They are always stable, independent of the discreteness parameter a. Fig. 1 shows the energy landscape E(ϕ_{1}, ϕ_{2}) for the control parameter a = 25. We observe the nine valleys separated by saddle points and the corresponding maxima. The six excited states which only exist for

Fig. 1. Energy landscape E(ϕ_{1}, ϕ_{2}) for a = 25, ϕ_{i} ∈ (−3π, 3π)

4

sufficiently large values of a cannot be calculated analytically. However, for a ≫ 1, a Taylor expansion of first order in \frac{1}{a} yields the configurations

(ϕ_{1}^{∗} − \frac{2}{a}(ϕ_{1}^{∗} − ϕ_{2}^{∗}), ϕ_{2}^{∗} − \frac{2}{a}(ϕ_{2}^{∗} − ϕ_{1}^{∗})) ,  (5)

where the quantities ϕ_{1}^{∗} and ϕ_{2}^{∗} vary through all possible asymptotic values for a → ∞, that is, ϕ_{i}^{∗} ∈ {−2π, 0, 2π}. In total, these are nine configurations. The three configurations with zero energy have already been discussed above. The four configurations with ϕ_{2} = 2π − ϕ_{1} for ϕ_{1} > 0 as well as ϕ_{2} = −2π − ϕ_{1} for ϕ_{1} < 0, namely

2π (1 − \frac{2}{a}, \frac{2}{a}) , 2π (−1 + \frac{2}{a}, −\frac{2}{a}) , 2π (\frac{2}{a}, 1 − \frac{2}{a}) , 2π (−\frac{2}{a}, −1 + \frac{2}{a}) ,

correspond to 2π-kinks. The remaining two antisymmetric configurations with ϕ_{2} = −ϕ_{1}, namely

2π (1 − \frac{4}{a}, −1 + \frac{4}{a}) , 2π (−1 + \frac{4}{a}, 1 − \frac{4}{a}) ,

represent 4π-kinks. They have the highest energy. While decreasing the control parameter a, the energy landscape is changed until the excited states eventually become unstable and we only have the three degenerate ground states. Fig. 2 shows

Fig. 2. Energy landscape E(ϕ_{1}, ϕ_{2}) for a = 3, ϕ_{i} ∈ (−3π, 3π)

the energy landscape E(ϕ_{1}, ϕ_{2}) for a small value of a. The saddle points act as energy barriers. Their height is proportional to the control parameter a, while the

5

energy values of the excited states are all of the order of one. Accordingly, for large discreteness parameters a, the energy difference between these states is small compared to the energy required to overcome the barrier and switch to another state of local equilibrium. The same arguments hold for larger values of N, where we will show in the next sections that the energy landscape exhibits a highly complicated multi-valley structure with an exponentially increasing number of local minima and saddle points. All possible phase differences ϕ_{i} as a function of the control parameter a are depicted in Fig. 3. Besides the lines with ordinates 2π, −2π, and 0, arising from the three-fold degenerate ground state, we see the curves constituing the 2πkinks (second and fourth curve from above and from below) as well as the curves constituing the 4π-kinks (third curve from above and third curve from below). Inspection of the eigenvalues of the Hessian reveals that they are always stable for a > 2π and for a > 6π, respectively. Increasing the number of wells increases the number of bands.

Fig. 3. Phase differences ϕ_{i} as a function of a for N = 2 and three potential wells

For a cosine potential with five wells, that is, if we restrict the values of ϕ_{i}, i = 1, 2, to the range (−5π, 5π), all possible phase differences as a function of a are depicted in Fig. 4. In comparison with Fig. 3, here we see in addition the curves constituing the 6π-kinks, which only exist for a > 10π, and the two curves constituing the 8π-kinks, which only exist for a > 14π.

To be more general, if we restrict the values of ϕ_{i}, i = 1, 2, to the range (−rπ, rπ), a stability analysis yields that a k(2π)-kink structure is only stable for a ≥ 2(2k−1)π. These multi-2π-kinks were first discovered by Peyrard and Kruskal.^{5} Various bunched states of 2π-kinks as well as 4π-kinks were discovered later by Ustinov et al.^{4} A general analytic theory about the origin of multikinks in dispersive nonlinear systems with a piecewise parabolic potential^{6} has been presented by

6

Fig. 4. Phase differences ϕ_{i} as a function of a for N = 2 and five potential wells

Champneys and Kivshar.^{7}

4. Counting the number of locally stable energy minima

There are only few results about the concrete number of excited equilibrium configurations, though it has often been mentioned that this number can be exponentially large, growing exponentially with the length N of the chain. According to the variational principle, stationary configurations obey the set of N coupled equations

−ϕ_{i+1} + 2ϕ_{i} − ϕ_{i−1} + a sin ϕ_{i} = 0,  i = 1, …, N.  (6)

(As before, indices have to be taken mod N.) Inserting Eq. (6) into the energy function Eq. (1) the energy simplifies to

E = a ∑_{i=1}^{N} (−ϕ_{i} sin ϕ_{i} + (1 − cos ϕ_{i})) .  (7)

Let us now concentrate on large values of a, where the discreteness parameter is sufficiently strong. In the limit of large a, the stationary condition Eq. (6) reduces to

sin ϕ_{i} = 0,  i = 1, …, N.  (8)

If we restrict the ϕ_{i}’s to the range (−rπ, rπ), where r is odd, then it is obvious that the set of Eqs. (8) has exactly r^{N} distinct solutions, all describing locally stable minima in the limit a → ∞. The solution vectors take the form ϕ^{∗} = (ϕ_{1}^{∗}, ϕ_{2}^{∗}, …, ϕ_{N}^{∗}), where the components ϕ_{i}^{∗} can only take one of the r discrete values ϕ_{i}^{∗} ∈ {−(r − 1)π, …, −2π, 0, 2π, …, (r − 1)π}, which correspond to the

7

number of kinks trapped in the system. Provided that a is sufficiently large, solutions of Eq. (6) still exist near the ones in the limit a → ∞ that we just described, and they are locally stable. They can be expanded in a Taylor series in \frac{1}{a},

ϕ_{i} = ϕ_{i}^{∗} − \frac{−ϕ_{i−1}^{∗} + 2ϕ_{i}^{∗} − ϕ_{i+1}^{∗}}{a} − \frac{−ϕ_{i−2}^{∗} + 4ϕ_{i−1}^{∗} − 6ϕ_{i}^{∗} + 4ϕ_{i+1}^{∗} − ϕ_{i+2}^{∗}}{a^{2}} + · · · . (9)

Inserting this expansion into Eq. (2), we immediately find

M = ∑_{i=1}^{N} ϕ_{i}^{∗}.  (10)

Hence, for sufficiently large values of a, the total phase can only take discrete values which are multiples of 2π (M ∈ {−(r − 1)Nπ, …, −2π, 0, 2π, …, (r − 1)Nπ}).

Since the junctions reside on a ring (of size N), we consider configurations that differ by a rotation (ϕ_{i} → ϕ_{s+i})_{1≤i≤N} (for some fixed s with 1 ≤ s ≤ N) or by a mirror reflection (opposite orientations) (ϕ_{i} → ϕ_{s−i})_{1≤i≤N} (for some fixed s with 1 ≤ s ≤ N) as equivalent. Note that, because of the periodic boundary conditions, the energy is invariant under these two types of transformations. Together, rotations and reflections form the dihedral group D_{N} which acts on the particles. The dihedral group is the symmetry group of an N-sided regular polygon. Since, for large a, the stable configurations satisfying Eq. (6) can be identified with the r^{N} Ndimensional vectors whose components consist only of the value 0 and multiples of 2π in the range {−(r−1)π, …, −2π, 0, 2π, …, (r−1)π}, the number of inequivalent solutions will be equal to the number of orbits of these circular vectors under the action of the dihedral group D_{N}. In more colloquial terms, the problem of counting the inequivalent stable configurations is equivalent to the problem of determining the number D_{r}(N) of different necklaces composed of N beads, each of r possible colours. The reflection symmetry allows the necklace to be picked up out of the plane and turned over. It follows from Redfield–Pólya theory (cf. [8, Theorem 5.1] and^{9,10}) that the number of these necklaces is

D_{r}(N) = \begin{cases}\frac{1}{2}r^{(N+1)/2} + \frac{1}{2N} ∑_{d|N} φ(d) r^{N/d}, & if N is odd, \\ \frac{1}{4}r^{N/2}(r + 1) + \frac{1}{2N} ∑_{d|N} φ(d) r^{N/d}, & if N is even, \end{cases}  (11)

where φ(n) is the Euler totient function, which, by definition, is the number of positive integers ≤ n that are relatively prime to n. To give an example, for three potential wells (r = 3), corresponding to necklaces made out of three colours, the first terms of the sequence D_{3}(N) are given by

{3, 6, 10, 21, 39, 92, 198, 498, 1219, 3210, … },  (12)

N = 1, 2, 3, . . .. For N = p an odd prime number, Eq. (11) simplifies to

D_{r}(p) = \frac{r}{2} (\frac{r^{p−1} − 1}{p} + r^{\frac{p−1}{2}} + 1) .  (13)

8

Moreover, since φ(1) = 1, the number of inequivalent locally stable configurations describing all local minima of the energy landscape D_{r}(N) increases exponentially with N. More precisely, as N tends to ∞, we have

D_{r}(N) ∼ \frac{r^{N}}{2N} .  (14)

5. The energy spectrum

Let us now study the degeneracies of the energy. In the following we consider three potential wells, although the generalization to an arbitrary number r of wells is straightforward. Besides the rotations and mirror reflections described in the previous section, the Hamiltonian Eq. (1) is also invariant with respect to the transformation

ϕ_{i} → −ϕ_{i}.  (15)

If we now consider configurations associated with the same energy as equivalent, for solving the corresponding counting problem we have to regard two necklaces also as equivalent if we keep one of the three colours fixed (ϕ_{i} close to 0; colour “0”) and exchange the two remaining colours (ϕ_{i} → −ϕ_{i}; that is, we identify colours “2π” and “−2π”). Thus, the situation that we encounter is that there is not only a symmetry group acting on the configurations (in our case: the dihedral group D_{N} acting on the beads), but there is also a symmetry group acting on the set of colours, namely the symmetric group S_{2} which exchanges the two colours “2π” and “−2π”. It follows from Redfield–Pólya theory (cf. [8, Theorem 5.4, with G = D_{N} and H = S_{2}, the latter acting on three colours, that is, P_{H} = \frac{1}{2}(x_{1}^{3} + x_{1}x_{2})]) that the number of inequivalent stable configurations equals

D_{3}^{S_{2}}(N) = \begin{cases}3^{(N−1)/2} + \frac{1}{4} + \frac{1}{4N} ∑_{d|N} φ(d) 3^{N/d}, & if N is odd, \\ 2 · 3^{(N−2)/2} + \frac{1}{4·2^{v_{2}(N)}} \\ +\frac{1}{4N} ∑_{d|N} φ(d) (1 + χ(d even))3^{N/d}, & if N is even, \end{cases}  (16)

for N = 1, 2, 3, …, where v_{2}(N) is the 2-adic valuation of N (that is, the largest v such that 2^{v}|d), and where χ(A) = 1 if A is true and χ(A) = 0 otherwise. For example, the first few values of this sequence are

{2, 4, 6, 13, 22, 52, 106, 266, 630, 1652, … }.  (17)

This should be compared to Eq. (12).

Furthermore, for ϕ_{i} > 0 (0 < ϕ_{i} < 2π), the Hamiltonian is also invariant with respect to the transformation

ϕ_{i} → 2π − ϕ_{i},  (18)

as well as for ϕ_{i} < 0 (−2π < 0 < ϕ_{i}) with respect to the complementary transformation

ϕ_{i} → −2π − ϕ_{i}.  (19)

9

Fig. 5. Phase differences ϕ_{i} versus a for N = 9, and the energy spectrum versus N for a = 25

In terms of our colouring problem, these two transformations translate into the transformation exchanging the two colours “0” and “2π” (but only for those necklaces featuring only these colours), respectively into the transformation exchanging the two colours “0” and “−2π” (but only for those necklaces featuring only these colours). Necklaces which are mapped to each other under these transformations are supposed to be equivalent. It suffices to consider the transformation exchanging the two colours “0” and “2π”. If we apply this transformation to a necklace featuring only these two colours then we may obtain the same necklace (modulo rotation and mirror reflection) or not. We must subtract the half of the number of those which are not mapped to itself from D_{3}^{S_{2}}(N). In order to determine this number, we consider equivalence classes of necklaces with beads which are coloured in only two colours, where we regard two necklaces as equivalent when they arise from each other by exchanging the two colours. Redfield–Pólya theory again (cf. [8, Theorem 5.4, with G = D_{N} and H = S_{2}, the latter acting on two colours, that is, P_{H} = \frac{1}{2}(x_{1}^{2} + x_{2})]) yields that the number of these equivalence classes is

D_{2}^{S_{2}}(N) = \begin{cases}2^{(N−3)/2} + \frac{1}{4N} ∑_{d|N} φ(d) 2^{N/d}, & if N is odd, \\ 2^{(N−2)/2} + \frac{1}{4N} ∑_{d|N} φ(d) (1 + χ(d even))2^{N/d}, & if N is even. \end{cases} (20)

Then, one half of the number of necklaces with only these two colours which are not invariant under the exchange of “0” and “2π” equals D_{2}(N)−D_{2}^{S_{2}}(N). Accordingly,

10

if we write D_{3}^{E}(N) for the number of inequivalent configurations associated with distinct energies, then we arrive at

D_{3}^{E}(N) = D_{3}^{S_{2}}(N) − (D_{2}(N) − D_{2}^{S_{2}}(N)),  (21)

N = 1, 2, 3, . . .. The first few terms in this sequence are

{1, 3, 4, 11, 18, 47, 97, 254, 607, 1618, … }.  (22)

The number of inequivalent locally stable minima D_{3}^{E}(N) increases exponentially with N. More precisely, as N tends to ∞, we have

D_{3}^{E}(N) ∼ \frac{3^{N}}{4N} .  (23)

Note that this is exactly half of the asymptotic value (14) with r = 3. In fact, the quotient of (23) by (14) (with r = 3), tends to 1/2 exponentially fast (compare (12) and (22)).

The distribution of the phase differences as a function of the control parameter a is shown in Fig. 5 (left). We observe bands and “sub-bands,” as well as gap structures. The bands cluster into groups, which themselves cluster again into subgroups, and so on. In other words, each band breaks up into sub-bands, which again break up into sub-sub-bands, and so on. Similarly, between the individual bands we find gap structures as well as sub-gap structures. These structures are self-similar in the sense that each portion is quasi a reduced-scale replica of the whole.^{6,12} Moreover, exponentially many phase difference values come exponentially close to each other. Fig. 5 (right) depicts the energy spectrum for the control parameter a = 25 for N = 3, 5, 7, and 9 junctions, respectively. We clearly see how the bandand the fractal structure evolve with increasing number of junctions. According to Eq. (17) the figure depicts the 4, 18, 97, and 607 distinct energies for N = 3, 5, 7 and N = 9, respectively. Also here, a thorough analysis reveals that the inner bands exhibit features of self-similarity.

6. Discussion

We have shown that, depending on the strength of the r-well cosine potential, the energy surface of the model can have the order of \frac{r^{N}}{N} locally stable minima separated by saddle points representing large barriers. The counting problem for the degeneracy of the total energy turned out to be equivalent to a well-known problem in combinatorial mathematics, in which colourings of necklaces are enumerated, with the dihedral group acting on the beads and another symmetry group acting on the colours, the solution of the problem being provided by Redfield–Pólya theory. Another important finding is that the distribution of the phase differences as well as the energy spectrum is fractal for a sufficiently large control parameter. We observed that the distribution of the phase differences form a hierarchical band structure such that exponentially many configurations become exponentially close to each other. Preliminary investigations also show a similar hierarchical band structure for the

11

energy spectrum.^{10–12} The structure of the energy landscape is reminiscent of the appearance of glassy structures within the formation of domains, kinks, and solitons. However, in contrast to the familiar spin glass models, the model is described by a Hamiltonin without disorder.

7. Acknowledgements

The first author thanks P. Martin, S. Charpentier and M. Olivier of the Université du Littoral de Dunkerque for helpful computer assistance. We also thank F.V. Kusmartsev for useful discussions about experimental aspects. This work has been partially supported by the European Science Foundation (ESF) in the framework of the network program “Arrays of Quantum Dots and Josephson Junctions”.

References

1. O. M. Braun and Y. S. Kivshar, The Frenkel–Kontorova Model, Concepts, Methods, and Applications, Springer–Verlag Berlin Heidelberg New York

2. O. V. Zhirov, G. Casati, and D. L. Shepelyansky, Phys. Rev. E 65,026220, (2002)

3. J. Pfeiffer, M. Schuster, A. A. Abdumalikov, Jr., and A. V. Ustinov, Phys. Rev. Lett. 96,034103, (2006)

4. A. V. Ustinov, B. A. Malomed, and S. Sakai, Phys. Rev. B 57,11691, (1998)

5. M. Peyrard and M. D. Kruskal, Physica D 14,88, (1984)

6. K. E. Kürten, Fractal properties of small magnetic particles and magnetic multilayers: An exactly solvable model, Condensed Matter Theories, Vol. 20, edited by J. W. Clark, (Nova Science Publishers, New York, 2005)

7. A. Champneys and Y. S. Kivshar, Phys. Rev. E 61,3,2551, (2000)

8. N. G. de Bruijn, Pólya’s theory of counting, in: Applied Combinatorial Mathematics, E. .F. Beckenbach, ed., Wiley, New York, 1964; reprinted by Krieger, Malabar, Florida, 1981.

9. E. N. Gilbert and J. Riordan, Illinois J. Math. 5,657 (1961)

10. K. E. Kürten and C. Krattenthaler, Multistability, Phase Transitions and Fractal Structures in Magnetic Multilayers, Condensed Matter Theories, Vol. 19, edited by E. Suraud, (Nova Science Publishers, New York, 2004)

11. K. E. Kürten and F. V. Kusmartsev, Phys. Rev. B, 72, 014433, (2005)

12. M. D. Forrester, K. E. Kürten and F. V. Kusmartsev, Phys. Rev. B, 74, 1, (2006)