A tunable macroscopic quantum system based on two fractional vortices

We propose a tunable macroscopic quantum system based on two fractional vortices. Our analysis shows that two coupled fractional vortices pinned at two artificially created \kappa\ discontinuities of the Josephson phase in a long Josephson junction can reach the quantum regime where coherent quantum oscillations arise. For this purpose we map the dynamics of this system to that of a single particle in a double-well potential. By tuning the \kappa\ discontinuities with injector currents we are able to control the parameters of the effective double-well potential as well as to prepare a desired state of the fractional vortex molecule. The values of the parameters derived from this model suggest that an experimental realisation of this tunable macroscopic quantum system is possible with today's technology.


I. INTRODUCTION
In a one-dimensional 0-π long Josephson junction (0-π LJJ) vortices carrying half of the magnetic flux quantum Φ 0 ≈ 2.07 × 10 −15 Wb appear spontaneously at the boundaries between 0 and π regions. 1 These semifluxons 2,3 have a degenerate ground state of either positive or negative polarity. The orientation ↑ or ↓ depends on the direction of supercurrent circulating around the 0-π boundary. The classical behavior of semifluxons has been studied extensively theoretically 1-5 and experimentally. [6][7][8] One can also construct "molecules" consisting of two or more semifluxons, for example, a molecule with antiferromagnetically (AFM) arranged semifluxon and antisemifluxon situated at a distance a from each other in a 0-π-0 LJJ. Such a molecule exhibits two degenerate ground states ↑↓ or ↓↑. One can switch the molecule between ↑↓ and ↓↑ states by applying a small uniform bias current to the LJJ. 9,10 Such an AFM semifluxon molecule was suggested as a semifluxon based qubit. 11 The low energy dynamics of the system was reduced to the dynamics of a point-like particle in a double-well potential and the parameters of this potential and the mass of the particle were related to the parameters of the 0-π-0 LJJ. The parameters for which quantum effects emerge were estimated. However, in such a qubit the height of the barrier in a doublewell potential depends only on geometrical parameters of the 0-π-0 LJJ, in particular, on the length a of the π part, which makes it impossible to tune the barrier during experiment.
In a conventional LJJ with the Josephson phase φ and the first Josephson relation j s = j c sin φ one can create artificially an arbitrary, electronically tunable κ discontinuity of the phase φ, so that φ(x) = µ(x) + κ H(x) and µ(x) is a smooth function without discontinuities and H(x) is a Heaviside step function. For κ = π this junction is equivalent to the 0-π LJJ described above. Such a device was proposed and successfully tested. [12][13][14] By creating discontinuities with κ = π (|κ| ≤ 2π) one can spontaneously form and study vortices carrying an arbitrary topological charge of ℘ = −κ or ℘ = −κ + 2π sgnκ that automatically appear pinned at the κ-discontinuity. The fractional flux associated with a vortex carrying a topological charge ℘ is Φ = Φ 0 ℘/(2π). Fractional fluxons can form a variety of ground states, 15 have characteristic eigenfrequencies 16 and get depinned by overcritical bias currents. [17][18][19] In the present paper we investigate an AFM molecule 20 consisting of two fractional vortices in which the topological charge of the vortices can be tuned electronically by changing κ. For appropriate parameters of the LJJ we obtain an electronically tunable qubit.
The paper is organized as follows. In Sec. II we describe the concept which underlies our idea to observe macroscopic quantum effects. In Sec. III we introduce the model of a LJJ with two discontinuities κ 1 and κ 2 separated by a distance a. Based on this model we derive in Sec. IV analytical expressions for stationary phases and their corresponding energies. In Sec. V we show how the state of the system can be manipulated. In particular, we discuss how the energy barrier can be reduced to reach the quantum regime, which is discussed in Sec. VI. In Sec. VII we analyze how the states of the molecule can be distinguished. Finally, to keep the paper self-contained an appendix provides a collection of relations for elliptic functions.

II. CONCEPT
Previously 11 we considered an AFM molecule of two semifluxons in a 0-π-0 LJJ and found that macroscopic quantum effects can be observed if the length a of the π segment is restricted to values 11 a c < a a c + 0.02. Here the scaled length a is measured in units of the Josephson penetration depth λ J . The critical distance a c is a bifurcation point: for a < a c there exists only one static solution with constant phase (flat phase state), while for a > a c there are two degenerate solutions (AFM states) that can be used as two classical states of a qubit. In the present paper we only consider infinitely long LJJs, where a c = π/2.
It is rather difficult to fabricate the junction with the value of a that lays within the narrow range mentioned above. Therefore, it is highly desirable to make a or a c tunable. The scaled length a is fixed by design and can vary only slightly due to variations of λ J , e.g., with temperature. The value of a c , however, can be varied in a larger range as described in the following.
Thus, we can fabricate a LJJ with the scaled distance 1.57 < a < 1.76 between the discontinuities. Then, we can decrease κ starting from κ = π to bring a c (κ) as close as needed towards a, as shown in Fig. 2. Note that the required precision for a is less restrictive by one order of magnitude than in the previous proposal 11 . In addition we obtain a system which can be tuned into the quantum regime. In the following sections we analyze this system in detail. The critical distance ac for asymmetric AFM molecules as a function of κ. It separates the flat phase state (white region) from the AFM state (gray region). As an example we consider a LJJ with a = 1.7 (straight line). By decreasing κ we can tune ac towards a, e.g., κ = 0.3π corresponding to ac = 1.66 as indicated by the point on the straight line.

III. MODEL
The dynamics of the Josephson phase µ(x, t) in infinitely long LJJs without bias current and dissipation is described by the sine-Gordon equation 3 We use dimensionless quantities, where lengths are written in units of the Josephson penetration depth λ J and time is written in units of ω −1 p , where ω p is the plasma frequency. Energies are measured in units of E J λ J , where E J is the Josephson energy per length. The subscripts x and t denote partial derivatives with respect to coordinate and time, accordingly.
The function describes two discontinuities −κ 1 and −κ 2 located at x = −a/2 and x = +a/2, where H(x) is the Heaviside step function. The values of κ 1 and κ 2 can be controlled 13 by injector currents I inj1 ∝ κ 1 and I inj2 ∝ κ 2 . The two discontinuities divide the x-axis into three parts to which we will refer to as the left, middle and right part of the junction, see Fig. 1(a). The boundary conditions for µ(x, t) are given by Although there are discontinuities at x = ±a/2 the phase µ(x, t) and its derivative µ x (x, t) with respect to x have to be continuous at x = ±a/2 at any time t.
The sine-Gordon equation (2) can be derived from the Lagrangian density where is the potential energy density.
To avoid infinite energies the restriction cos [µ(±∞, t) + θ(±∞, t)] = 1 (7) assures that the potential energy density (6) vanishes at x = ±∞. This condition restricts the phases at x = ±∞ to the values where n L and n R are integer numbers.

IV. STATIONARY SOLUTIONS
In this section we first derive analytical expressions for the stationary Josephson phase µ(x) and the corresponding energies. These expressions depend on the values and µ r ≡ µ (+a/2) (9b) of µ at the boundaries (the left and the right) between the regions which have to be determined selfconsistently. Furthermore, we analyze the stability of the stationary solutions.

A. Phases
In order to find stationary solutions of Eq. (2) we integrate its static version µ xx − sin(µ + θ) = 0 (10) in each region of the junction. This procedure leads to a conservation relation where k is a positive real number and has different values in the inner and outer parts of the junction. By comparing this equation with Eq. (A.6) we see that the phase solves Eq. (10)  From the boundary conditions (4) and (7) in Eq. (11) we find k = 1 (14) in the left and right part of the junction. By matching Eq. (11) for the outer part (k = 1) and the inner part (k = 1) of the junction at x = ±a/2 we obtain Therefore, k can vary between the values 0 and √ 2 in the middle of the junction. Since the behavior of the Jacobi amplitude function depends on the value of k as discussed in the appendix, we consider three different cases: (a) outer regions with k = 1, (b) inner region with k > 1 and (c) inner region with k < 1.
a. Outer regions with k = 1. With the help of Eq. (A.7) we obtain from Eq. (13) the expression for k → 1, where µ l and µ r are restricted to Since we have to match these two solutions and their derivatives to their counterparts in the inner region we define the signs b. Inner region with k > 1. With the help of Eq. (A.4) it can easily be verified that (19) matches the solution for the left part of the junction, that is µ(−a/2) = µ l and σ = σ l . Here we have used the abbreviation To match the solution for the right part of the junction, that is µ(a/2) = µ r and σ = σ r , we have the additional constraint c. Inner region with k < 1. Similar to the case k > 1 we find where we have used Note that there is only one integer number n. Therefore, the values of ψ l,r and µ l,r are restricted to |ψ r − ψ l | ≤ π and |µ r − µ l | ≤ 2π. The additional constraint to match the solution for the right part reads Due to the properties of elliptic integrals, for a given value of a the integer number ν is limited to the values We now have expressed the stationary solutions of the sine-Gordon equation (2) in terms of µ L , µ R , µ l and µ r . The possible values of µ L and µ R are restricted by Eqs. (8) and depend on the topological charge of the system. The values of µ l and µ r have to be determined from Eqs. (15), (21) and (25). Examples will be presented in Sec. V B.

B. Energies
In this section we use the results of the previous section to calculate the energy of the stationary solutions.
Integrating Eq. (6) and using (11) we obtain the energy of the solution in the left part of the junction. From Eq. (16) we find and arrive at Similarly, we obtain the expression for the right part of the junction.
For the middle part of the junction we again combine Eqs. (6) and (11) and obtain the energy of the phase in this region of the junction. From Eq. (13) we find and (32) where E is the Jacobi epsilon function defined by Eq. (A.11). The total energy of a stationary solution µ(x) is the sum of the energies of these three domains.

C. Stability analysis
To study the stability of a stationary solution µ(x) we insert the ansatz into the sine-Gordon equation (2) and linearize it around the stationary solution µ(x). Since µ(x) solves the stationary sine-Gordon equation (10), we obtain the differential equation for the eigenmodes ψ(x) with eigenvalues ω 2 , which, after substituting Eq. (13), takes the form Note that Eq. (38) is the Lamé equation 21 . The boundary conditions and the matching conditions at x = ±a/2 for ψ(x) follow from the boundary and matching conditions for µ(x) discussed in Sec. III. As long as the lowest eigenvalue ω 2 0 of the Lamé equation (38) is positive, the stationary solution µ(x) is stable. When it becomes negative, the stationary solution µ(x) is unstable. In Sec. V B we solve Eq. (38) numerically.

V. TAILORING MOLECULE STATES
Our ultimate goal is to have a tunable qubit based on an AFM molecule. For this purpose we first bring the system into one of its degenerate ground states and choose κ such that the barrier between the two states is high. By tuning κ we then reduce this energy barrier to reach the quantum regime.

A. Initializing the AFM state
In a previous publication 11 it was shown that a configuration with the discontinuities κ 1 = π and κ 2 = −π (0-π-0 LJJ) provides an effective double-well potential with degenerate ground states where quantum tunneling can be observed.
An intuitive way to arrive at this state from the κ 1 = κ 2 = 0 state is to use and sweeping κ i from 0 to π. Then, we want to deal with a molecule consisting of a direct and a complementary vortex. Therefore, we have to change the discontinuities appropriately, e.g., where the (de)tuning κ t starts at 0 and changes in the range between −π and +π. In this approach one can employ two injector current sources: one creating κ i wired to produce two discontinuities of opposite signs like in Eq. (39) and another one creating κ t and wired to produce two discontinuities of the same sign like in Eq. (40). The advantages of this approach are: (i) instead of using the first current source, one could replace the junction by a conventional 0-π-0 LJJ and then only apply injector currents for the tuning κ t according to Eq. (40); (ii) this technique works also in an annular JJ; (iii) after κ i has reached the value π we know that we have prepared the (+π, −π) state rather than the (−π, +π) state. The disadvantage of this approach is that one always injects a current ∝ κ = π or more into the system, which brings additional noise.
An alternative way is to use only one injector current source, which is coupled to injectors as Starting from the κ = 0 state and no vortices in the JJ, one slowly increases the value of κ and finds the symmetric ferromagnetic vortex configurations (κ, κ). When κ is increased above κ ↿↿ c (a), 15 where π < κ ↿↿ c (a) < 2π, the molecule emits one fluxon and reconfigures into one of the degenerate asymmetric antiferromagnetic states (κ, κ − 2π) or (κ − 2π, κ). To reach the (+π, −π) or (−π, +π) state we decrease κ to π. Then we further decrease κ to reach the quantum limit. This technique has two advantages: (i) one needs only one injector current source; (ii) one can operate the qubit at small κ values, which brings less noise. This technique, however, has two disadvantages: (i) it does not work in annular LJJs, where the emitted fluxon cannot leave easily; 22 (ii) we do not know if we prepared the (κ, κ − 2π) or (κ − 2π, κ) state.
For our further calculations we choose the second approach where κ, and therefore, the additional noise in the system is small in the quantum limit.

B. Tuning the barrier height
We consider the particular injector-tuning technique described by Eq. (41) for characterizing the states, its energies and stabilities during this process.

Determination of phase parameters
After the initialization procedure discussed in the previous section the two discontinuities have the value κ = π and the AFM state has the topological charge µ R − µ L = 0. For simplicity we choose n L = 0 and n R = −1 in Eqs. (8).
When we tune κ no (anti)fluxon is emitted or absorbed. Therefore, n L and n R do not change and the topological charge is induced in the system when the discontinuities have the value κ.
In the stationary solutions presented in Sec. IV we still have to find µ l and µ r . In our examples we restrict ourselves to the values of a covering the gray region in Fig. 2 near the boundary a c (κ), that is, 1.5 ≤ a ≤ 2 and 0 ≤ κ ≤ 2π. It turns out that for these parameters there are only solutions for k < 1. For this limitation, Eq. (15) restricts the values of µ l and µ r to the intervals κ 2 − π <µ l +2πn 1 < κ 2 (43a) and 3κ 2 <µ r +2πn 2 < π + 3κ 2 . (43b) Here we use instead of Eq. (17) the more restrictive condition |µ l,r −µ L,R | < π which requires n 1 = 0 and n 2 = 1.
This condition is reasonable because even in the case when we have a whole integer fluxon at one of the discontinuities the phase changes at most by π when one goes from infinity up to the fluxon center. By using the constraints Eqs. (43) together with Eq. (15) we express µ r in terms of µ l and obtain and as the two only possible solutions. Next we substitute these two expressions into Eq. (25) and solve the resulting equations numerically for µ l for a given distance a. We denote the solution corresponding to Eq. (44a) by an index m, and the two solutions corresponding to Eq. (44b) by an index + and −. Figure 3 illustrates this process for a = 2.
When we insert the three values of µ m and µ ± into Eqs. (16) and (22)   and The solution µ + (x) describes a molecule consisting of a κ vortex at x = −a/2 and the complementary κ−2π vortex at x = +a/2. The solution µ − (x) is a complementary molecule (κ − 2π, κ), which has the same energy. These two solutions are stable whereas the solution µ m (x) is unstable, see Sec. V B 2.

Energies and stability analysis
Depending on the parameters of the LJJ there is either a single stable stationary solution µ m (x), or two stable solutions µ ± (x) and one unstable solution µ m (x). The two stable solutions are separated by an energy barrier which is governed by the unstable solution. To reach the quantum limit we want to reduce this energy barrier by tuning κ. Therefore, we investigate how the energies and stabilities of the stationary solutions depend on κ.
For our examples we use the following JJ parameters: (46) where λ L is the London penetration depth, j c is the critical current density of the JJ, C is the capacitance of the JJ per unit of area, w is the JJ's width, and Φ 0 ≈ 2.07 × 10 −15 Wb. Typical parameters are λ L = 90 nm, w = 1 µm, and C = 4.1 µF/cm 2 with j c = 100 A/cm 2 . This corresponds to λ J = 38 µm, ω p = 2π × 42.8 GHz, and E J λ J = 78.4 meV.
In order to calculate the energies of the solutions we insert the values for µ l obtained in the previous subsection into Eq. (35) and obtain the energies U m (κ) of the state µ m (x) and U + (κ) = U − (κ) of the states µ ± (x). In Fig. 5 we depict the corresponding energy difference for different values of a. This figure shows that the energy barrier can be tuned down to zero by decreasing the discontinuity κ from π to 0. The stability of the solutions µ m (x) and µ ± (x) is determined by the eigenvalues of Eq. (38). These eigenvalues are calculated numerically. We have used a junction of length l = 20 to emulate an infinitely long JJ. The results are depicted in Fig. 6. Positive eigenvalues characterize stable solutions, whereas negative eigenvalues characterize unstable solutions.
In the lower part of Fig. 2, that is for a < a c (π) ≈ 1.57, there is only one stationary solution µ m (x), so that ∆U (κ) ≡ 0 in Fig. 5. Figure 6 shows that the corre-  sponding eigenvalue ω 2 0 is positive, therefore this solution is stable.
Finally, for a > a c (0) there are always three stationary solutions µ m (x) and µ ± (x). According to Fig. 6, µ m (x) is unstable and µ ± (x) are stable. All three solutions reach the same energy at κ → 0 and the energy barrier vanishes as ∆U (κ) ∝ κ, see Fig. 5. The eigenvalues ω 2 0 (κ) corresponding to the three solutions vanish too at κ → 0.
In fact, all three solutions are just single fluxons. The two stable solutions µ + (x) and µ − (x) are fluxons weakly pinned at −a/2 and +a/2, respectively. The unstable solution µ m (x) is a fluxon at x = 0. It is rather interesting that ∆U (κ) vanishes for κ → 0 even for a > a c (0). Therefore, one can always make the barrier as small as needed for arbitrary a > a c (π) = π/2 and is not limited by the interval of 1.57 a 1.76! This limit κ → 0 is similar to the heart-shaped qubit, in which a fluxon is trapped in a double-well potential created by a non-uniform magnetic field 23 . In this situation, residual pinning in the system can play a major role and make the parasitic pinning potential larger than the one we construct here.

VI. QUANTUM REGIME
We expect to observe coherent quantum oscillations between the two degenerate states corresponding to the classical states µ + (x) and µ − (x) in the region where the energy barrier in Fig. 5 becomes sufficiently small, i.e., the coupling between these two states is sufficiently large. In order to quantify this effect we map the dynamics of the system to the dynamics of a single particle in a double-well potential and calculate the energy splitting δε separating the two lowest energy levels. Because the oscillation frequency between the two states µ + (x) and µ − (x) is given by δε/ we are then able to determine the values of the discontinuity κ ∝ I inj for a given value of a where the quantum oscillations become observable.  Fig. 5 the right vertical axis shows the corresponding temperature for typical parameters mentioned after Eq. (46). For the parameters indicated by two arrows more details are shown in Fig. 8.

A. Single-mode approximation
In order to map the dynamics of the complete system to the dynamics of a single particle in a double-well potential we express the phase in terms of the stationary solution µ m (x) and the corresponding eigenmodes ψ n (x) defined by Eq. (38). By inserting this expansion into the Lagrangian density (5) and integrating over x we obtain a Lagrangian for the mode amplitudes q n (t) which describes the motion of a fictitious particle in many dimensions. If the lowest eigenvalue ω 2 0 of Eq. (38) is sufficiently separated from the next higher eigenvalue we can expect that for low energies the particle only moves along the "q 0 direction". Motivated by this simple picture, we omit the higher modes in Eq. (48) and use the approximation where ψ 0 (x) is the eigenfunction corresponding to the lowest eigenvalue ω 2 0 . We insert this equation into our Lagrangian density (5) and take into account only terms up to the fourth order in q 0 and obtain the Lagrangian where is the potential energy. Here the positive parameter K is defined by and ψ 0 is normalized according to ψ 2 0 dx = 1. Due to the symmetry relation (45a) for µ m (x) there is no thirdorder term in Eq. (51).
Note that for ω 2 0 < 0, where µ m (x) is unstable, U (q 0 ) describes a double-well potential with a maximum at q 0 = 0 and two minima at q 0 = ± 6/K|ω 0 |. For ω 2 0 > 0 it only has one minimum at q 0 = 0. In the first case the oscillation frequencies around the minima are ω ± = √ 2 |ω 0 |. In the second case ω 0 is the oscillation frequency around the mimimum.

B. Energy splitting
For the Lagrangian (50) we can derive the stationary Schrödinger equation is the dimensionless Planck constant and the energy eigenvalues ε j are given in units of E J λ J . In order to calculate the energy splitting we solve Eq. (53) numerically for the potential U (q 0 ) given by Eq. (51). Additionally, we compare our numerical results for δε 01 ≡ ε 1 − ε 0 with the scaled energy splitting in the semiclassical limit. 24 Figure 7 shows our numerical results for δε 01 and the semi-classical expression ∆ for different values of a. We note that the semi-classical expression ∆ is a good approximation for δε 01 as long as δε 01 is small. To have a good quantum-mechanical two-level system (qubit) the two lowest eigenvalues ε 0 and ε 1 have to be well separated from the higher eigenvalues. Therefore, we additionally compare δε 12 ≡ ε 2 − ε 1 to δε 01 .
To establish a two-level system at a temperature T three conditions have to be fulfilled: (i) ∆U ≫ k B T to suppress thermal hopping between the two classical ground states; (ii) δε 01 ≫ k B T to observe coherent oscillations; (iii) δε 12 ≫ δε 01 to have approximately a twolevel system. For large values of ∆U , the energy splitting δε 01 becomes small. Therefore, we have to find parameters where ∆U and δε 01 have reasonable values.
From Fig. 7 we find that condition (iii) is violated if κ is tuned too close to κ c (a) whereas condition (ii) is violated if κ becomes too large. Energy splittings δε 01 corresponding to approximately 25 mK look promising.
To check condition (i) we calculate the energy difference ∆U from Eq. (47) and find values corresponding to approximately between 100 mK and 135 mK. Therefore, we may expect to observe coherent quantum oscillations for temperatures below 25 mK which is experimentally accessible.
Two typical examples are shown in Fig. 8; one for a = 1.8 and κ = 0.08π (a), and one for a = 1.6 and κ = 0.68π (b), indicated by two arrows in Fig. 7. The first example corresponds to the case a > a c (0) while the second example corresponds to the case a c (π) < a < a c (0). The results for the two regimes look similar: For the parameters of Fig. 8 (a) we find ∆U = 1.23 × 10 −4 (111 mK), δε 01 = 0.27 × 10 −4 (25 mK) and δε 12 = 1.19 × 10 −4 (109 mK) while for the parameters of Fig. 8 (b) we find ∆U = 1.47 × 10 −4 (135 mK), δε 01 = 0.27 × 10 −4 (25 mK) and δε 12 = 1.34 × 10 −4 (122 mK). From Fig. 7 we conclude that for larger values of a the energy splitting δε 01 becomes more sensitive to κ. The advantages and disadvantages of both regimes in terms of read-out are discussed in the next section. In both cases the energy splitting δε01 corresponds to 25 mK.

VII. READ-OUT
To observe coherent quantum oscillations one has to distinguish between the (κ, κ − 2π) state described by µ + (x) and the (κ − 2π, κ) state described by µ − (x) of the molecule. Therefore, it is necessary to readout its state. One possibility is to readout the flux associated with each fractional vortex in a molecule, similar to the earlier experiments. 10 The magnetic fluxes measured in these states should be distinguishable.
In order to see if this is possible we calculate the flux measured on the left half of the LJJ for the two different states µ − (x) and µ + (x) and the results are shown in Fig. 9.
For a < a c (0) the flux difference ∆Φ ≡ Φ + − Φ − vanishes at the bifurcation points κ c (a). For a > a c (0) the value of ∆Φ remains finite for κ → 0, as in this limit the two states correspond to two integer fluxons weakly pinned at x = ±a/2. For larger values of a the flux difference ∆Φ increases. Therefore, in junctions with larger a the two states are easier to distinguish. In this case, however, the qubit is harder to control because the relevant range of κ becomes smaller, see Fig. 7.

VIII. CONCLUSIONS
We have presented a concept of a qubit based on a molecule consisting of two fractional Josephson vortices.
Two degenerate ground states are separated by an energy barrier, which can be tuned during the experiment by changing simultaneously the values of the discontinuities. The concept may work in both linear and annular geometries.
In particular, we have obtained analytical solutions for the stationary phases and their energies in unbiased LJJs with two discontinuities. Furthermore, we have analyzed the stability of the stationary solutions by calculating the corresponding eigenmodes numerically. We have used these eigenmodes to map the low-energy dynamics of the system to the dynamics of a particle in a one-dimensional double-well potential and have solved the corresponding Schrödinger equation.
Finally we have shown how the energy barrier can be tuned with the help of injector currents to reach the quantum limit. Our results indicate that for typical parameters a quantum-mechanical two-level system (qubit) can be established for temperatures below 25 mK which is at the limit of modern dilution refrigerators. Experiments with such fractional vortex molecules are in progress in the Tübingen group.

Appendix: Elliptic integrals and elliptic functions
In this appendix we introduce the elliptic integral of first kind and the Jacobi amplitude function and provide a collection of formulas which are used in the present paper. These relations are in accordance with Ref. 21.

Elliptic integrals
The elliptic integral of first kind is defined by where the argument k is called modulus. For φ = π/2 it reduces to the complete elliptic integral of first kind K(k) ≡ F(π/2, k). (A. 2) The elliptic integral F(φ, k) obeys the symmetry relation F(nπ ± φ, k) = 2n K(k) ± F(φ, k), where n is an integer.

Jacobi amplitude function
For k < 1 the elliptic integral (A.1) is a monotonously increasing function of φ. The Jacobi amplitude function is the inverse of the elliptic integral. For a given value u ≡ F(φ, k) we have to find the corresponding integration limit φ. Then the Jacobi amplitude function is defined by am(u, k) ≡ φ. (A.4) The relation am(u, k) ≡ arcsin{k −1 sin[ am(ku, k −1 )]} (A.5) extends this definition for k > 1. The Jacobi amplitude function is monotonously increasing for k < 1 and periodic for k > 1. Furthermore, it solves the differential equation and the addition theorem E(u 1 + u 2 , k) = E(u 1 , k) + E(u 2 , k) − k 2 sn(u 1 , k) sn(u 2 , k) sn(u 1 + u 2 , k). (A. 15)