A tunable macroscopic quantum system based on two fractional vortices

According to Leggett [11], the first experiments produced a quantum superposition of macroscopically distinct states. In order to validate this situation experimentally, Leggett proposed [11] a 'real-time' experiment concerning the two classical ground states |L〉 and |R〉 of a double-well potential: 'one starts the system off in, say, |L〉, switches one's measuring apparatus off, then switches it on again at time t and detects whether it is in |L〉 or |R〉, and by making repeated runs of this type with varying values of t, plots a histogram of the probability P_L(t) that the system is in |L〉 at time t'. A quantum-mechanical calculation would yield

$$\begin{equation} P_L(t) = [1 + \cos ( \Delta_{01} t)]/2 , \end{equation} \tag{ 1 }$$

where Δ₀₁ = (E₁ − E₀)/ℏ and E₁ − E₀ is the energy splitting separating the two lowest quantum-mechanical energy levels of the double-well potential.

This experiment was realized with superconducting charge [12], phase [13] and flux [14] devices, where quantum oscillations were based on Cooper pairs, the Josephson phase and the flux generated by a superconducting ring, respectively. But until now, no experiment has reported coherent quantum oscillations for Josephson vortices. These are vortices of electric current appearing in superconducting Josephson junctions (JJs). Their macroscopic quantum nature was only demonstrated in escape experiments [15], where the macroscopic ground state tunnels out of a potential well.

In this paper, we propose a 'real-time' experiment for Josephson vortices based on the experimental setup described in [16]. Furthermore, we show how our system can be tuned from a classical regime, where no quantum oscillations are observable, into the quantum regime.

There is also a different concept for observing such vortex oscillations [17]. It does not involve injector currents, which are sources of noise in our device. On the other hand, this alternative concept requires a bias current and an external magnetic field which also brings fluctuations into the system. These two components are absent in our setup.

Both approaches are based on long Josephson junctions (LJJs). These devices have attracted the attention of researchers for several decades because an LJJ is a relatively clean and well-controlled nonlinear system, which allows one to study fundamental physics. For example, Cherenkov radiation from a fast moving fluxon [18] was observed, the ratchet effect [19–21] was demonstrated and the Zurek–Kibble scenario of phase transitions [22] was tested. Furthermore, LJJs have applications as a flux flow oscillator in submillimetre receivers [23]. New possibilities have opened up with the advent of π JJs and technologies allowing us to combine conventional 0 and π junctions into one device [24, 25].

In a one-dimensional 0–π LJJ, vortices carrying half of the magnetic flux quantum Φ₀ ≈ 2.07 × 10⁻¹⁵ Wb appear spontaneously at the boundaries between 0 and π regions [26, 27]. These semifluxons [28, 29] have a degenerate ground state of either positive or negative polarity. The orientation ↑ or ↓ depends on the direction of the supercurrent circulating around the 0–π boundary. The classical behaviour of semifluxons has been studied extensively theoretically [26, 28–30] and experimentally [31–34].

One can also construct 'molecules' consisting of two or more semifluxons, for example, a molecule with antiferromagnetically (AFM) arranged semifluxons and antisemifluxons 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 the ↑↓ and ↓↑ states by applying a small uniform bias current to the LJJ [35–37].

Such an AFM semifluxon molecule was suggested to observe quantum oscillations of semifluxons [38]. 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 double-well potential the height of the barrier 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 the 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), where μ(x) is a smooth function without discontinuities and H(x) is a Heaviside step function. This κ discontinuity of the Josephson phase ϕ(x) is produced by a pair of tiny current injectors situated next to each other and contacting one of the superconducting electrodes of the junction [39]. The current I_inj is injected into one injector and collected from the other one. It causes a strong gradient of the phase of the superconducting wave function between the two injectors. This in turn leads to a strong gradient of the Josephson phase which can be considered as a 'discontinuity' if the distance between the current injectors is small compared to the Josephson penetration depth λ_J. This phase discontinuity κ∝I_inj provokes the appearance of a pinned Josephson vortex, whose magnetic flux can have arbitrary values controlled by the injected current and is not limited to multiples of the flux quantum Φ₀.

Such a device was proposed and successfully tested [39–41]. For example, at κ = π there appear semifluxons and the junction is equivalent to the 0–π LJJ described above. By creating discontinuities with κ ≠ π (|κ| ⩽ 2π) one can spontaneously form and study vortices carrying an arbitrary topological charge of $\wp =-\kappa$ or $\wp =-\kappa +2\pi \,{\mathrm {sgn}}\,{\kappa }$. The fractional flux associated with a vortex carrying a topological charge $\wp$ is $\Phi =\Phi _0\,\wp /(2\pi )$. Fractional fluxons can form a variety of ground states [42], have characteristic eigenfrequencies [43] and get depinned by overcritical bias currents [44, 45].

In this paper we investigate a 'molecule' [16] 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 macroscopic quantum system that is suitable for the experimental observation of coherent quantum oscillations.

This paper is organized as follows. In section 2 we describe the concept of our tunable system and introduce the model of an LJJ with two discontinuities κ₁ and κ₂ separated by a distance a. Based on this model we derive in section 3 analytical expressions for stationary phases and their corresponding energies. In section 4 we address the question of the preparation of the initial state and show that the energy barrier can be reduced to reach the quantum regime discussed in more detail in section 5. The distinguishability of the two states of the molecule is the topic of section 6. In section 7 we qualitatively discuss the sensitivity of the system to fluctuations, and section 8 summarizes our main results. Finally, to keep the paper self-contained, an appendix provides a collection of the relations for elliptic functions.

Previously [38], we considered molecules consisting of two semifluxons in a 0–π–0 LJJ situated at the beginning and the end of the π segment of length a. The two states of the molecule were composed of an AFM arranged pair of a semifluxon and an antisemifluxon. Such a molecule exhibits two degenerate ground states. We found that macroscopic quantum effects can be observed if a is restricted to values [38] a_c < a ≲ a_c + 0.02, where lengths are measured in units of the Josephson penetration depth λ_J and a_c = π/2 is a bifurcation point. For a < a_c there exists only one static solution with a constant phase (flat phase state), while for a > a_c there are two degenerate solutions (AFM states).

It is rather difficult to fabricate the junction with the value of a that lies 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. For a 0–π–0 junction, the value of a_c is fixed. For a 0–κ–2κ junction, however, the value of a_c can be tuned electronically over a large range as described in the following.

Consider the 0–κ–2κ LJJ of figure 1(a) with two discontinuities κ, see figure 1(b). For κ = π this reduces to the previous case of a 0–π–2π LJJ or, as the Josephson phase is 2π periodic, to a 0–π–0 LJJ. For a > a_c the two classically degenerate ground states are molecules consisting of vortices with topological charges (+π,−π) or (−π,+π). Now, let us vary κ around κ = π. The two corresponding states shown in figures 1(c) and (d) then are (κ,κ − 2π) and (κ − 2π,κ) [42]. Note that the critical distance [42]

$$\begin{equation} a_{\rm c}(\kappa)=2\,{F}\left(\frac{\pi}{4},\sqrt{1-\sin\frac{\kappa}{2}}\right), \end{equation} \tag{ 2 }$$

where F(ϕ,k) is the elliptic integral defined in (A.1), is not constant but depends (weakly) on κ as shown in figure 2. We find that the value of a_c may change from a_c(π) = π/2 ≈ 1.57 to $a_{\mathrm { c}}(0)=a_{\mathrm { c}}(2\pi )=2\ln (1+\sqrt {2})\approx 1.76$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Thus, we can fabricate an 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 figure 2. Note that the required precision for a is less restrictive by one order of magnitude than in the previous proposal [38]. In addition, we obtain a system that can be tuned into the quantum regime. In the following sections we analyse this system in detail.

The dynamics of the Josephson phase μ(x,t) in infinitely long LJJs without bias current and dissipation is described by the sine-Gordon equation [29]

$$\begin{equation} \mu_{xx}(x,t)-\mu_{tt}(x,t)-\sin[\mu(x,t)+\theta(x)]=0. \end{equation} \tag{ 3 }$$

We use dimensionless quantities, where lengths are written in units of the Josephson penetration depth λ_J and time is written in units of ω⁻¹_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

$$\begin{equation} \theta(x) \equiv -\kappa_1\,{H}\left(x+\frac{a}{2}\right) - \kappa_2\,{H}\left(x-\frac{a}{2}\right) \end{equation} \tag{ 4 }$$

describes two discontinuities −κ₁ and −κ₂ located at x = −a/2 and +a/2, where H(x) is the Heaviside step function. The values of κ₁ and κ₂ can be controlled [40] by injector currents I_inj1∝κ₁ and I_inj2∝κ₂. The two discontinuities divide the x-axis into three parts, which we will refer to as the left, middle and right parts of the junction, see figure 1(a).

The boundary conditions for μ(x,t) are given by

$$\begin{equation} \mu_x(-\infty,t)=\mu_x(+\infty,t)=0 . \end{equation} \tag{ 5 }$$

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 (3) can be derived from the Lagrangian density

$$\begin{equation} \mathcal{L}\equiv {\textstyle\frac{1}{2}}{\mu_t^2}(x,t)-\mathcal{U} , \end{equation} \tag{ 6 }$$

where

$$\begin{equation} \mathcal{U}\equiv {\textstyle\frac{1}{2}}{\mu_x^2}(x,t)+1-\cos\left[\mu(x,t)+\theta(x)\right] \end{equation} \noindent \tag{ 7 }$$

is the potential energy density.

To avoid infinite energies, the restriction

$$\begin{equation} \cos\left[\mu(\pm\infty,t)+\theta(\pm\infty,t)\right]=1 \end{equation} \tag{ 8 }$$

ensures that the potential energy density (7) vanishes at x = ± ∞. This condition restricts the phases at x = ± ∞ to the values

$$\begin{equation} \mu_{\mathrm{L}}\equiv\mu(-\infty,t)=2\pi n_{\mathrm{L}}, \end{equation} \tag{ 9 }$$

$$\begin{equation} \mu_{\mathrm{R}}\equiv\mu(+\infty,t)= 2\pi n_{\mathrm{R}} + \kappa_1 + \kappa_2, \end{equation} \noindent \tag{ 10 }$$

where n_L and n_R are integer numbers.

In order to find stationary solutions of (3) we integrate its static version

$$\begin{equation} \mu_{xx}-\sin(\mu+\theta)=0 \end{equation} \tag{ 13 }$$

in each region of the junction. As a result we obtain the conservation relation

$$\begin{equation} {\textstyle\frac{1}{2}}\mu_x^2+\cos(\mu+\theta)=2k^2-1 \end{equation} \tag{ 14 }$$

or

$$\begin{equation} \left(\frac{1}{2k}\mu_x\right)^2=1-\frac{1}{k^2}\sin^2\left(\frac{\mu+\theta}{2}+\frac{\pi}{2}\right), \end{equation} \tag{ 15 }$$

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 (A.6) we see that the phase

$$\begin{equation} \mu(x)=(2n-1)\pi-\theta(x)\pm2\,{\mathrm{am}}\left[k(x+c),k^{-1}\right] \end{equation} \tag{ 16 }$$

solves (13) in every region. Here n is an integer number and c is an arbitrary constant. Furthermore, we have introduced the Jacobi amplitude function am(u,k) defined by (A.4) and (A.5).

From the boundary conditions (5) and (8) in (14), we find that

$$\begin{equation} k=1 \end{equation} \tag{ 17 }$$

in the left and right parts of the junction.

By matching (14) for the outer parts (k = 1) and the inner part (k ≠ 1) of the junction at x = ± a/2, we obtain

$$\begin{eqnarray} k^2&=&1+\sin\frac{\kappa_1}{2}\sin\left(\mu_{\mathrm{l}}-\frac{\kappa_1}{2}\right) \nonumber\\ &=&1{-\sin\frac{\kappa_2}{2}\sin\left(\mu_{\mathrm{r}}-\kappa_1-\frac{\kappa_2}{2}\right)} . \end{eqnarray} \tag{ 18 }$$

Therefore, k can vary between the values 0 and $\sqrt {2}$ in the middle of the junction. Since the behaviour 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) the inner region with k > 1 and (c) the inner region with k < 1.

Outer regions with k = 1.

With the help of (A.7) we obtain from (16) the expression

$$\begin{equation} \mu(x)=\mu_{\mathrm{L,R}}+4\arctan\Big(\mathrm{e}^{a/2-|x|} \tan\frac{\mu_{\mathrm{l,r}}-\mu_{\mathrm{L,R}}}{4}\Big) \end{equation} \tag{ 19 }$$

for k → 1, where μ_l and μ_r are restricted to

$$\begin{equation} |\mu_{\mathrm{l,r}}-\mu_{\mathrm{L,R}}|<2\pi. \end{equation} \tag{ 20 }$$

Since we have to match these two solutions and their derivatives to their counterparts in the inner region, we define the signs

$$\begin{equation} \sigma_{\mathrm{l}}\equiv\mathrm{sign}[\mu_x(-a/2)]=\mathrm{sign}(\mu_{\mathrm{l}}-\mu_{\mathrm{L}}), \end{equation} \tag{ 21 }$$

$$\begin{equation} \sigma_{\mathrm{r}}\equiv\mathrm{sign}[\mu_x(+a/2)]=\mathrm{sign}(\mu_{\mathrm{R}}-\mu_{\mathrm{r}}) \end{equation} \tag{ 22 }$$

of μ_x(x).

Inner region with k > 1.

With the help of (A.4) it can easily be verified that

$$\begin{equation} \mu(x)=\kappa_1-\pi+2\,\sigma\,{\mathrm{am}}\left[k\left(x+{a}/{2}\right)+\sigma \,{F}(\psi_{\mathrm{l}},k^{-1}),k^{-1}\right] \end{equation} \tag{ 23 }$$

matches the solution for the left part of the junction, that is, μ(−a/2) = μ_l and σ = σ_l. Here we have used the abbreviation

$$\begin{equation} \psi_{\mathrm{l,r}}\equiv{\textstyle\frac{1}{2}}\left({\mu_{\mathrm{l,r}}+\pi-\kappa_1}{}\right). \end{equation} \tag{ 24 }$$

To match the solution for the right part of the junction, that is, μ(a/2) = μ_r and σ = σ_r, we have the additional constraint

$$\begin{equation} a={\sigma}k^{-1}\left[\,{F}(\psi_{\mathrm{r}},k^{-1})-\,{F}(\psi_{\mathrm{l}},k^{-1})\right]. \end{equation} \tag{ 25 }$$

Inner region with k < 1.

Similar to the case k > 1, we find that

$$\begin{equation} \mu(x)={(2n-1)\pi}+\kappa_1+2\sigma_{\mathrm{l}}\,{\mathrm{am}} [k(x+{a}/{2})+\sigma_{\mathrm{l}}\,k\,{F}(\varphi_{\mathrm{l}},k),k^{-1}], \end{equation} \tag{ 26 }$$

where we have used

$$\begin{equation} \varphi_{\mathrm{l,r}}\equiv \arcsin\left(k^{-1}\sin\bar\psi_{\mathrm{l,r}}\right) \end{equation} \tag{ 27 }$$

and

$$\begin{equation} \bar\psi_{\mathrm{l,r}}+n\pi \equiv\psi_{\mathrm{l,r}},\quad|\bar \psi_{\mathrm{l,r}}|\leqslant \pi/2. \end{equation} \tag{ 28 }$$

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

$$\begin{equation} a=\left[\sigma_{\mathrm{r}}\,{F}(\varphi_{\mathrm{r}},k)-\sigma_{\mathrm{l}}\,{F}(\varphi_{\mathrm{l}},k)\right] +(1-\sigma_{\mathrm{l}}\sigma_{\mathrm{r}})\,{K}(k)+4\nu\,{K}(k). \end{equation} \tag{ 29 }$$

Due to the properties of elliptic integrals, for a given value of a the integer number ν is limited to the values

$$\begin{equation} 0\leqslant \nu < \frac{a}{2\pi}+\frac{1}{2}. \end{equation} \tag{ 30 }$$

We now have expressed the stationary solutions of the sine-Gordon equation (3) in terms of μ_L, μ_R, μ_l and μ_r. The possible values of μ_L and μ_R are restricted by (9) and (10) and depend on the topological charge of the system. The values of μ_l and μ_r have to be determined from (18), (25) and (29). Examples will be presented in section 4.2.

In this section we use the results of the previous section to calculate the energy of the stationary solutions.

Integrating (7) and using (14) we obtain the energy

$$\begin{equation} U_{\mathrm{l}}=\int_{-\infty}^{-a/2} \mu_x^2\, \mathrm{d} x= \int_{\mu_{\mathrm{L}}}^{\mu_{\mathrm{l}}} \mu_x\, \mathrm{d} \mu \end{equation} \tag{ 31 }$$

of the solution in the left part of the junction. From (19) we find that

$$\begin{equation} \mu_x=2 \sin \left[ {\textstyle \frac{1}{2}}\left(\mu(x)-\mu_{\mathrm{L}}{}\right) \right] \end{equation} \tag{ 32 }$$

and arrive at

$$\begin{equation} U_{\mathrm{l}}=8\sin^2 \left[\textstyle \frac{1}{4} \left( {\mu_{\mathrm{l}}-\mu_{\mathrm{L}}}{} \right) \right]. \end{equation} \tag{ 33 }$$

Similarly, we obtain the expression

$$\begin{equation} U_{\mathrm{r}}=8\sin^2 \left[\textstyle \frac{1}{4} \left( {\mu_{\mathrm{r}}-\mu_{\mathrm{R}}}{} \right) \right] \end{equation} \tag{ 34 }$$

for the right part of the junction.

For the middle part of the junction we again combine (7) and (14) and obtain the energy

$$\begin{equation} U_{\mathrm{m}}=\int_{-a/2}^{a/2} \left[ \mu_x^2+2(1-k^{2}) \right]\,\mathrm{d} x \end{equation} \tag{ 35 }$$

of the phase in this region of the junction. From (16) we find that

$$\begin{equation} \mu_x^2=4k^2 \,{\mathrm{dn}}^2 \left[k(x+c),k^{-1}\right] \end{equation} \tag{ 36 }$$

and

$$\begin{equation} \pm \,{\mathrm{sn}} \left[k(x+c),k^{-1}\right] = (-1)^n \cos \left[\textstyle \frac{1}{2} \left( {\mu(x)-\kappa_1}{} \right) \right]. \end{equation} \tag{ 37 }$$

Equations (36) and (A.11) allow us to evaluate the integral. Finally, we obtain with the help of (A.15) and (37) the energy

$$\begin{equation} \fl U_{\mathrm{m}}=2a(1-k^{2})+4k\mathcal{E}\left(ka,k^{-1}\right) - {4} k^{-1} \cos\frac{\mu_{\mathrm{l}}-\kappa_1}{2} \cos\frac{\mu_{\mathrm{r}}-\kappa_1}{2}\,{\mathrm{sn}}(ka,k^{-1}) , \end{equation} \tag{ 38 }$$

where $\mathcal {E}$ is the Jacobi epsilon function defined by (A.11). The total energy

$$\begin{equation} U=U_{\mathrm{l}}+U_{\mathrm{m}}+U_{\mathrm{r}} \end{equation} \tag{ 39 }$$

of a stationary solution μ(x) is the sum of the energies of these three domains.

To study the stability of a stationary solution μ(x) we insert the ansatz

$$\begin{equation} \mu(x,t)=\mu_{}(x)+\psi(x)\mathrm{e}^{-\mathrm{i}\omega t} \end{equation} \tag{ 40 }$$

into the sine-Gordon equation (3) and linearize it around the stationary solution μ(x). Since μ(x) solves the stationary sine-Gordon equation (13), we obtain the differential equation

$$\begin{equation} -\psi''(x)+\cos\left[\mu_{}(x)+\theta(x)\right]\psi(x)=\omega^2\psi(x) \end{equation} \tag{ 41 }$$

for the eigenmodes ψ(x) with eigenvalues ω², which, after substituting (16), takes the form

$$\begin{equation} \psi''(x)+[\,1+\omega^2-2\,{\mathrm{sn}}^2(k(x+c),k^{-1})]\psi(x)=0. \end{equation} \tag{ 42 }$$

Note that (42) is the Lamé equation [46]. The boundary conditions and the matching conditions at x = ± a/2 for ψ(x) follow from the boundary and matching conditions for μ(x) discussed in section 2.2.

As long as the lowest eigenvalue ω²₀ of the Lamé equation (42) is positive, the stationary solution μ(x) is stable. When it becomes negative, the stationary solution μ(x) is unstable. In section 4.2 we solve (42) numerically.

In a previous publication [38] it was shown that a configuration with the discontinuities κ₁ = π and κ₂ = −π (0–π–0 LJJ) provides an effective double-well potential with degenerate ground states where quantum tunnelling can be observed.

An intuitive way to arrive at this state from the κ₁ = κ₂ = 0 state is to use

$$\begin{equation} \kappa_1=\kappa_{\rm i}, \quad \kappa_2=-\kappa_{\rm i} \end{equation} \tag{ 43 }$$

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.

$$\begin{equation} \kappa_1=\pi+\kappa_{\mathrm{t}}, \quad \kappa_2=-\pi+\kappa_{\mathrm{t}}, \end{equation} \tag{ 44 }$$

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 (43) and another one creating κ_t and wired to produce two discontinuities of the same sign like in (44). The advantages of this approach are: instead of using the first current source, one could replace the junction by a conventional 0–π–0 LJJ and then only apply injector currents to tune κ_t according to (44); (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. A 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

$$\begin{equation} \kappa_1=\kappa_2=\kappa . \end{equation} \tag{ 45 }$$

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 $\kappa _{\mathrm { c}}^{\upharpoonleft \upharpoonleft }(a)$ [42], where $\pi <\kappa _{\mathrm { c}}^{\upharpoonleft \upharpoonleft }(a)<2\pi$, 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; and (ii) one can operate the quantum system 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,⁴ and (ii) we do not know whether 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.

We consider the particular injector-tuning technique described by (45) for characterizing the states, its energies and stabilities during this process.

4.2.1. Determination of phase parameters

After the preparation 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 (9) and (10).

When we tune κ no (anti)fluxon is emitted or absorbed. Therefore, n_L and n_R do not change and the topological charge

$$\begin{equation} \mu_{\mathrm{R}}-\mu_{\mathrm{L}}=2\kappa-2\pi \end{equation} \tag{ 46 }$$

is induced in the system when the discontinuities have the value κ.

In the stationary solutions presented in section 3 we still have to find μ_l and μ_r. In our examples we restrict ourselves to the values of a covering the grey region in figure 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, (18) restricts the values of μ_l and μ_r to the intervals

$$\begin{equation} \frac{\kappa}{2}-\pi<\mu_{\mathrm{l}}+2\pi n_1<\frac{\kappa}{2} \end{equation} \tag{ 47 }$$

and

$$\begin{equation} \frac{3\kappa}{2}<\mu_{\mathrm{r}}+2\pi n_2<\pi+\frac{3\kappa}{2}. \end{equation} \tag{ 48 }$$

Here we use instead of (20) the more restrictive condition |μ_l,r − μ_L,R| < π which requires n₁ = 0 and n₂ = 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 centre.

By using the constraints (47) and (48) together with (18) we express μ_r in terms of μ_l and obtain

$$\begin{equation} \mu_{\mathrm{r}}=-\mu_{\mathrm{l}}+2\kappa-2\pi \end{equation} \tag{ 49 }$$

and

$$\begin{equation} \mu_{\mathrm{r}}=\mu_{\mathrm{l}}+\kappa-\pi \end{equation} \tag{ 50 }$$

as the only two possible solutions.

Next we substitute these two expressions into (29) and solve the resulting equations numerically for μ_l for a given distance a. We denote the solution corresponding to (49) by an index m, and the two solutions corresponding to (50) by indices + and −. Figure 3(a) illustrates this process for a = 2.

Zoom In Zoom Out Reset image size

When we insert the three values of μ_m and μ_± into (19) and (26) we find the phases μ_m(x) and μ_±(x) depicted in figure 3(b). These solutions satisfy the symmetry relations

$$\begin{equation} \mu_{\rm m}(-x) = -\mu_{\rm m}(x) +2\kappa -2\pi \end{equation} \noindent \tag{ 51 }$$

and

$$\begin{equation} \mu_{\pm}(-x) = -\mu_{\mp}(x) +2\kappa -2\pi . \end{equation} \noindent \tag{ 52 }$$

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 section 4.2.2.

4.2.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:

$$\begin{equation} \omega_{\mathrm{p}}=\sqrt{\frac{2\pi j_{\rm c}}{\Phi_0 C}},\quad E_{\mathrm{J}}=\frac{j_{\rm c} w \Phi_0}{2\pi} ,\quad \lambda_{\mathrm{J}}=\sqrt{\frac{\Phi_0}{2\pi \mu_0 2 \lambda_{\mathrm{L}} j_{\rm c}}}, \end{equation} \tag{ 53 }$$

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 Φ₀ ≈ 2.07 × 10⁻¹⁵ Wb. Typical parameters are λ_L = 90 nm, w = 1 μm and C = 4.1 μF cm⁻² with j_c = 4 kA cm⁻² [48]. This corresponds to λ_J = 6.0 μm, ω_p = 2π × 274 GHz and E_Jλ_J = 496 meV.

Without loss of generality, we restrict ourselves to 0 < κ < π. The results for the range π < κ < 2π are obtained by substituting κ → 2π − κ.

In order to calculate the energies of the solutions we insert the values for μ_l obtained in the previous section into (39) and obtain the energies U_m(κ) of the state μ_m(x) and U₊(κ) = U₋(κ) of the states μ_±(x). In figure 4(a) we depict the corresponding energy difference

$$\begin{equation} \Delta U(\kappa)\equiv U_{\mathrm{m}}(\kappa)-U_{\pm}(\kappa) \end{equation} \tag{ 54 }$$

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.

Zoom In Zoom Out Reset image size

The stability of the solutions μ_m(x) and μ_±(x) is determined by the eigenvalues of (42). 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 figure 4(b). Positive eigenvalues characterize stable solutions, whereas negative eigenvalues characterize unstable solutions.

In the lower part of figure 2, that is, for a < a_c(π) ≈ 1.57, there is only one stationary solution μ_m(x), so that ΔU(κ) ≡ 0 in figure 4(a). Figure 4(b) shows that the corresponding eigenvalue ω²₀ is positive; therefore this solution is stable.

For a_c(π) < a < a_c(0) ≈ 1.76 we have to distinguish between κ > κ_c(a) and κ < κ_c(a), where the function κ_c(a) is the inverse of a_c(κ) defined by (2). For κ > κ_c(a) there are three stationary solutions μ_m(x) and μ_±(x). As shown in figure 4(a) the energy barrier ΔU(κ) becomes lower when we decrease κ and vanishes at κ = κ_c(a). Our numerics shows that ΔU(κ)∝(κ − κ_c)² for κ → κ_c. From the eigenvalues ω²₀(κ) displayed in figure 4(b) corresponding to these stationary solutions, we conclude that μ_m(x) is unstable whereas μ_±(x) are stable.

At κ = κ_c(a) all three solutions have the same energy and ΔU(κ) vanishes. The eigenvalues ω²₀(κ) join at the bifurcation point κ = κ_c(a) and vanish, see figure 4(b). For κ < κ_c(a) the two stable solutions μ_±(x) disappear, while the unstable solution μ_m(x) becomes stable, see figure 4(b).

Finally, for a > a_c(0) there are always three stationary solutions μ_m(x) and μ_±(x). According to figure 4(b), μ_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 figure 4(a). The eigenvalues ω²₀(κ) corresponding to the three solutions vanish too at κ → 0.

In this limiting case, 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 [49]. 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.

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

$$\begin{equation} \mu(x,t)=\mu_{\rm m}(x)+\sum_{n=0}^\infty q_n(t)\psi_n(x) \end{equation} \tag{ 55 }$$

in terms of the stationary solution μ_m(x) and the corresponding eigenmodes ψ_n(x) defined by (42). By inserting this expansion into the Lagrangian density (6) 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 ω²₀ of (42) is sufficiently separated from the next higher eigenvalue we can expect that for low energies the particle only moves along the 'q₀ direction'. Motivated by this simple picture, we omit the higher modes in (55) and use the approximation

$$\begin{equation} \mu(x,t)\approx\mu_{\rm m}(x)+q_0(t)\psi_0(x), \end{equation} \tag{ 56 }$$

where ψ₀(x) is the eigenfunction corresponding to the lowest eigenvalue ω²₀. We insert this equation into our Lagrangian density (6) and take into account only terms up to the fourth order in q₀ and obtain the Lagrangian

$$\begin{equation} L={\textstyle \frac{1}{2}} \dot{q_0}^2-U(q_0), \end{equation} \tag{ 57 }$$

where

$$\begin{equation} U(q_0)\equiv {\textstyle \frac{1}{2}}{\omega_0^2}{}q_0^2+{\textstyle\frac{1}{24}}Kq_0^4 \end{equation} \tag{ 58 }$$

is the potential energy. Here the positive parameter K is defined by

$$\begin{equation} K\equiv-\int_{-\infty}^{+\infty}\mathrm{d} x\,\cos[\mu_0(x)+\theta(x)]\,\psi_0^4(x) \end{equation} \tag{ 59 }$$

and ψ₀ is normalized according to $\int \psi _0^2\,\mathrm {d} x=1$. Due to the symmetry relation (51) for μ_m(x) there is no third-order term in (58).

Note that for ω²₀ < 0, where μ_m(x) is unstable, U(q₀) describes a double-well potential with a maximum at q₀ = 0 and two minima at $q_0=\pm \sqrt {6/K}|\omega _0|$. For ω²₀ > 0 it only has one minimum at q₀ = 0. In the first case the oscillation frequencies around the minima are $\omega _\pm =\sqrt {2}\,|\omega _0|$. In the second case ω₀ is the oscillation frequency around the minimum.

For the Lagrangian (57) we can derive the stationary Schrödinger equation

$$\begin{equation} \left(\frac{-\hbar_{\mathrm{eff}}^2}{2}\frac{\partial^2}{\partial q_0^2} +U(q_0)\right)u_j(q_0)=\varepsilon_j u_j(q_0), \end{equation} \tag{ 60 }$$

where

$$\begin{equation} \hbar_{\mathrm{eff}}\equiv \frac{\hbar \omega_{\mathrm{p}}}{E_{\mathrm{J}}\lambda_{\mathrm{J}}} \end{equation} \tag{ 61 }$$

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 (60) numerically for the potential U(q₀) given by (58). Additionally, we compare our numerical results for δε₀₁ ≡ ε₁ − ε₀ with the scaled energy splitting

$$\begin{equation} \Delta=8\hbar_{\mathrm{eff}}\sqrt{\frac{2\Delta U}{\pi\hbar_{\mathrm{eff}}\omega_\pm}}\exp\left(-\frac{16\Delta U}{3\hbar_{\mathrm{eff}}\omega_\pm}\right) \end{equation} \tag{ 62 }$$

in the semiclassical limit [50].

Figure 5 shows our numerical results for δε₀₁ and the semi-classical expression Δ for different values of a. We note that the semi-classical expression Δ is a good approximation for δε₀₁ as long as δε₀₁ is small. To have a good quantum-mechanical two-level system the two lowest eigenvalues ε₀ and ε₁ have to be well separated from the higher eigenvalues. Therefore, we additionally compare δε₁₂ ≡ ε₂ − ε₁ to δε₀₁.

Zoom In Zoom Out Reset image size

To establish a two-level system at a temperature T three conditions have to be fulfilled: (i) ΔU ≫ k_BT to suppress thermal hopping between the two classical ground states; (ii) δε₀₁ ≫ k_BT to observe coherent oscillations; and (iii) δε₁₂ ≫ δε₀₁ to have approximately a two-level system. For large values of ΔU, the energy splitting δε₀₁ becomes small. Therefore, we have to find parameters where ΔU and δε₀₁ have reasonable values.

From figure 5, 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 δε₀₁ corresponding to approximately 100 mK look promising. To check condition (i) we calculate the energy difference ΔU from (54) and find values corresponding to temperatures between 0.9 and 1.2 K. Therefore, we may expect to observe coherent quantum oscillations as defined by (1) with a frequency Δ₀₁/(2π) ≈ 2 GHz for temperatures below 100 mK which is experimentally accessible.

Two typical examples are shown in figure 6; one for a = 1.8 and κ = 0.08π (a), and one for a = 1.6 and κ = 0.68π (b), indicated by two arrows in figure 5. 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 figure 6(a) we find that ΔU = 1.66 × 10⁻⁴ (957 mK), δε₀₁ = 0.17 × 10⁻⁴ (95 mK) and δε₁₂ = 1.28 × 10⁻⁴ (736 mK), while for the parameters of figure 6(b) we find that ΔU = 1.94 × 10⁻⁴ (1118 mK), δε₀₁ = 0.17 × 10⁻⁴ (95 mK) and δε₁₂ = 1.44 × 10⁻⁴ (829 mK). From figure 5 we conclude that for larger values of a the energy splitting δε₀₁ becomes more sensitive to κ. The advantages and disadvantages of both regimes in terms of read-out are discussed in the next section.

Zoom In Zoom Out Reset image size

The elliptic integral of the first kind is defined by

$$\begin{equation} F(\phi,k)\equiv\int\nolimits_0^\phi \frac{\mathrm{d}\theta}{\sqrt{1-k^2\sin^2\theta}}, \end{equation} \tag{ A.1 }$$

where the argument k is called modulus.

For ϕ = π/2 it reduces to the complete elliptic integral of the first kind

$$\begin{equation} {K}(k)\equiv\,{F}(\pi/2,k). \end{equation} \tag{ A.2 }$$

The elliptic integral F(ϕ,k) obeys the symmetry relation

$$\begin{equation} {F}(n\pi\pm\phi,k)=2n\,{K}(k)\pm\,{F}(\phi,k), \end{equation} \tag{ A.3 }$$

where n is an integer.

For k < 1 the elliptic integral (A.1) is a monotonically 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

$$\begin{equation} {\mathrm{am}}(u,k)\equiv\phi. \end{equation} \tag{ A.4 }$$

The relation

$$\begin{equation} {\mathrm{am}}(u,k)\equiv \arcsin\{k^{-1}\sin[\,{\mathrm{am}}(ku,k^{-1})]\} \end{equation} \noindent \tag{ A.5 }$$

extends this definition for k > 1.

The Jacobi amplitude function is monotonically increasing for k < 1 and periodic for k > 1. Furthermore, it solves the differential equation

$$\begin{equation} \left[\frac{\mathrm{d}}{\mathrm{d} u}\,{\mathrm{am}}(u,k)\right]^2=1-k^2\sin^2[\,{\mathrm{am}}(u,k)] \end{equation} \tag{ A.6 }$$

and for k → 1 it reduces to

$$\begin{equation} \,{\mathrm{am}}(u,1)=2\arctan(\mathrm{e}^{u})-\pi/2. \end{equation} \tag{ A.7 }$$

The three Jacobi elliptic functions

$$\begin{equation} {\mathrm{sn}}(u,k)\equiv\sin[\,{\mathrm{am}}(u,k)] \end{equation} \tag{ A.8 }$$

and

$$\begin{equation} {\mathrm{cn}}(u,k)\equiv\cos[\,{\mathrm{am}}(u,k)], \end{equation} \tag{ A.9 }$$

as well as

$$\begin{equation} {\mathrm{dn}}(u,k)\equiv\frac{\mathrm{d}}{\mathrm{d} u}[\,{\mathrm{am}}(u,k)], \end{equation} \tag{ A.10 }$$

together with the Jacobi epsilon function

$$\begin{equation} \mathcal{E}(u,k)\equiv\int\nolimits_0^u\mathrm{d} t\ \,{\mathrm{dn}}^2(t,k), \end{equation} \tag{ A.11 }$$

are defined in terms of the Jacobi amplitude function am(u,k).

These functions satisfy the relations

$$\begin{equation} {\mathrm{sn}}(u,k)= k^{-1}\,{\mathrm{sn}}(ku,k^{-1}) \end{equation} \tag{ A.12 }$$

and

$$\begin{equation} {\mathrm{cn}}(u,k)={\mathrm{dn}}(ku,k^{-1}) \end{equation} \tag{ A.13 }$$

as well as

$$\begin{equation} {\mathrm{dn}}(u,k)={\mathrm{cn}}(ku,k^{-1}), \end{equation} \tag{ A.14 }$$

and the addition theorem

$$\begin{equation} \fl \mathcal{E}(u_1+u_2,k)=\mathcal{E}(u_1,k)+\mathcal{E}(u_2,k)-k^2\,\,{\mathrm{sn}}(u_1,k)\,{\mathrm{sn}}(u_2,k)\,{\mathrm{sn}}(u_1+u_2,k). \end{equation} \tag{ A.15 }$$