Strongly interacting Majorana modes in an array of Josephson junctions

Majorana zero modes are fermions which are their own antiparticles. They are believed to exist as effective particles in the middle of the gap of a topological superconductor. The recent interest in Majorana zero modes originates in their non-Abelian exchange statistics [1], which is the basis for potential applications in quantum computation [2]. Based on the theoretical proposal [3] to realize these exotic states in semiconducting nanowires with strong spin–orbit coupling in a magnetic field and in the proximity of a conventional (nontopological) superconductor, recent experimental progress has shown signatures of Majorana zero modes in the tunnelling conductance of normal conducting–superconducting [4] and superconducting–normal conducting–superconducting systems [5]. More information and references on the fast developing field can be found in the recent reviews [6].

As the nanowire is one dimensional (1D), interaction effects become important. On the one hand, employing field theoretical methods, Gangadharaiah et al [7] have argued that strong electron–electron interactions generically destroy the topological phase by suppressing the superconducting gap. On the other hand, using a combination of analytical and numerical methods, Stoudenmire et al [8] have shown that repulsive interactions significantly decrease the required Zeeman energy and increase the parameter range for which the topological phase exists. It is believed that the origin of these effects is an interaction-driven renormalization of the Zeeman gap [8, 9]. Furthermore, it has been shown that in helical liquids the scattering processes between the constituent fermion bands open gaps, which in turn lead to a stabilization of the Majorana states against interactions [10] and that (an odd number of) Majorana zero modes are in fact stable against general interactions [11]. For further studies of the effect of electron–electron interactions on Majorana zero modes in nanowires and two-chain ladders, see [12].

Recently, a macroscopic version of the Kitaev chain [13], a toy model for a nanowire supporting noninteracting Majorana modes, was proposed [14] in a 1D array of topological superconducting islands. Its advantage over microscopic implementations is that individual parameters of the effective Hamiltonian can be potentially tuned in situ. Here, we show that additional capacitances between adjacent islands lead to an effective interaction between the low-energy Majorana degrees of freedom. We present the phase diagram of the system which demonstrates that sufficiently strong repulsive interactions will drive the system from the topologically trivial phase into the topological phase supporting Majorana zero modes before eventually leading to a Mott insulating state. We discuss how the parameters of the system can be tuned by changing the gate voltages and exploit this to propose the detection of the different phases and phase boundaries in a tunnelling experiment.

We discuss a system consisting of a 1D array of N superconducting islands (see figure 1). Because of the proximity-coupled semiconducting nanowire each of the islands has two midgap Andreev states, i.e. Majorana modes, located at the ends of the nanowire [3]. We will denote with γ_ka and γ_kb the two Majorana operators on island k associated with these zero modes. The Majorana operators are Hermitian γ_μ = γ^†_μ and fulfil the Clifford algebra {γ_μ,γ_ν} = 2δ_μν. The total Lagrangian L = T − V_J − V_M of the system consists of three terms. Coupling of the Majorana modes on nearby islands leads to the term [13, 15] $V_{\mathrm {M}} = \mathrm {i}\,E_{\mathrm {M}} \sum _{k=1}^N \gamma _{kb} \gamma _{k+1a} \cos [(\phi _{k+1} - \phi _k)/2]$ where ϕ_k is the superconducting phase of the kth island. Here and in the following, we assume for convenience periodic boundary conditions such that islands 1 and N + 1 are equivalent. Apart from the Majorana modes, the term discussed above has the additional degrees of freedom ϕ_k due to the condensate of Cooper pairs. Similar to [14], we eliminate these by connecting each superconducting island with a strong Josephson junction to a common (ground) superconductor. This fixes the superconducting phases (up to some quantum phase-slips discussed below) and is described by the Hamiltonian $V_{\mathrm {J}} = E_{\mathrm {J}} \sum _{k=1}^N (1- \cos \phi _k);$ here, E_J = ℏI_c/2e is the effective Josephson coupling of each of the Josephson junctions with critical current I_c. In the limit E_J ≫ E_M the junctions effectively pin the phases of all superconducting islands to a common value ϕ_k ≡ 0. As a result, V_M reduces to the pure Majorana coupling

$$\begin{equation} H_{\mathrm{M}} = \mathrm{i}\,E_{\mathrm{M}} \sum_{k=1}^N \gamma_{kb} \gamma_{k+1a} \quad (E_{\mathrm{J}} \gg E_{\mathrm{M}}). \end{equation} \tag{ 1 }$$

Finally, the kinetic term $T = (\hbar ^2/8e^2) \sum _{k=1}^N [ C (\dot \phi _{k+1} - \dot \phi _k)^2 + C_{\mathrm {G}} \dot \phi _k^2 ] + (\hbar /2e) \sum _{k=1}^N q \dot \phi _k$ occurs due to the capacitive couplings, where C denotes the capacitance between neighbouring islands and C_G = C_g + C_J the (total) capacitance to the ground. Here, C_J is the capacitance of the strong Josephson junction and C_g the capacitance to a common back gate at voltage V_g with respect to the ground superconductor. Apart from the charging energy, the back gate introduces a term proportional to the induced charge q = C_gV_g which is tunable with single-electron precision [16] via V_g and whose effect will become important later. The typical capacitive energy scale E_C = e²/2C_Σ of single islands depends on the total capacitance C_Σ = 2C + C_G.

Zoom In Zoom Out Reset image size

As we have seen, the strong coupling to the ground superconductor pins the superconducting phase differences to ϕ_k ≡ 0 and thus changes the energy due to the Majorana modes from V_M to H_M. The effect of the strong coupling to the ground superconductor on the charging energy T is more subtle. We first present the results for C = 0 before extending them to nonzero C: as charge and phase are conjugate variables the pinning of the phases ϕ_k greatly reduces the effect of charging [14]. An effective charging energy still arises due to quantum phase slips through the junctions. For example, changing ϕ_k from 0 to 2π leads to a charging energy (E_J ≫ E_C) [17]

$$\begin{equation} H_{Ck} = \Gamma_\Delta \cos[\pi(q/e + n_k)]= \Gamma_\Delta \cos(\pi q/e ) \mathcal{P}_k. \end{equation} \tag{ 2 }$$

Here, the tunnelling amplitude is given by Γ_Δ ≃ E^1/4_CE^3/4_Je^−S_Δ with $S_\Delta = \hbar ^{-1} \int \!\mathrm {d}\tau \,L_E =\sqrt {8E_{\mathrm {J}}/E_C}$ the dimensionless Euclidean action along the classical trajectory γ_Δ: ϕ_k∈[0,2π] in the inverted potential with L_E(τ) = −L(t = i τ).⁴ The cosine term in (2) occurs due to the Aharonov–Casher interference of the two tunnelling paths ϕ_k = 0 → 2π and ϕ_k = 0 → − 2π which lead to an indistinguishable final state and thus interfere with a phase difference depending on the total induced charge q + en_k [14, 18]. The term $n_k = \frac 12 (1- \mathcal {P}_k)\in \{0,1\}$ is the contribution to the charge due to the parity $\mathcal {P}_k = \mathrm {i}\, \gamma _{ka} \gamma _{kb} \in \{-1,1\}$ of the number of electrons on the superconducting island encoded in the state of the Majorana zero modes [14, 15]. Equation (2) is a chemical potential term: for large Γ_Δ, all the fermionic states are either filled or empty and the system is in a band-insulating state.

Going away from this special point and introducing a finite cross-capacitance parameterized by η = 2C/C_Σ with η∈[0,1], the classical path for a phase slip in ϕ_k does not only involve ϕ_k but also the other phases. To lowest nonvanishing order in η, we only need to take into account the change in ϕ_k−1 and ϕ_k+1 and obtain $S_\Delta = \sqrt {8 E_{\mathrm {J}}/E_C} [1 + (\pi ^2 -12)\eta ^2/96]$, which is accurate to more than two digits all the way up to η = 1 as numerics confirms. Additionally, the effective capacitance coupling between two islands becomes important. This term is generated by simultaneous phase slips of the two phases ϕ_k and ϕ_k+1 of neighbouring islands. Due to the symmetry between island k and k + 1, we have ϕ_k = ϕ_k+1 for the classical path γ_U: ϕ_k = ϕ_k+1∈[0,2π] with the corresponding Euclidean action $S_U= \sqrt {16 (2- \eta ) E_{\mathrm {J}}/E_C}$. The Aharonov–Casher interference in this case leads to the tunnelling amplitude Γ_U ≃ E^1/4_CE^3/4_Je^−S_U depending on the total charge 2q + e(n_k + n_k+1) of the two islands involved. Thus we obtain the interaction term

$$\begin{eqnarray} &&H_{Uk} = \Gamma_U \cos[2\pi q/e + \pi (n_k + n_{k+1})] = \Gamma_U \cos(2\pi q /e) \mathcal{P}_k \mathcal{P}_{k+1}. \end{eqnarray} \tag{ 3 }$$

We stress that (3) is an interaction term involving four Majorana operators and that the amplitude is modulated by twice the induced charge q compared to (2). We call U > 0 repulsive interaction as it prefers having an occupied site next to an empty one. For strong repulsive interactions the system is driven into a Mott-insulator or equivalently commensurate charge density wave (CDW) state.

The complete low-energy Hamiltonian in the limit E_J ≫ E_M,E_C is given by $H_{{\scriptscriptstyle \rm ANNNI}} = \sum _{k=1}^N (H_{Ck} + H_{Uk}) + H_{\mathrm {M}}$, which constitutes the transverse axial next-nearest-neighbour Ising (ANNNI) model [19] as can be seen by performing a Jordan–Wigner transformation $\mathcal {P}_k = \mathrm {i}\, \gamma _{ka} \gamma _{kb} = \sigma ^z_k$, i γ_kbγ_k+1a = σ^x_kσ^x_k+1, resulting in the spin Hamiltonian

$$\begin{equation} H_{{\scriptscriptstyle \rm ANNNI}} = \sum_{k=1}^N ( \Delta \sigma^z_k + U \sigma^z_k \sigma^z_{k+1} + E_{\mathrm{M}} \sigma^x_k \sigma^x_{k+1}). \end{equation} \tag{ 4 }$$

Here, the σ^x,y,z_k are Pauli matrices, Δ = Γ_Δcos(πq/e) and U = Γ_Ucos(2πq/e). We note that the energies Δ and U can be tuned via the charge q induced on the superconducting islands through the gate voltage V_g. As the fermionic Hamiltonian without the interaction term proportional to U has been intensively studied previously [13], we focus here on the effects of the interaction term. In our set-up, this term is most important in the case of C ≫ C_G where η ≈ 1 and $\Gamma _U \approx \Gamma _\Delta \exp [-1.25 \sqrt {E_{\mathrm {J}}/E_C}]$, which can be as large as 0.3 for E_J ≈ E_C; note that in this regime the actions S_Δ and S_U are still much larger than 1 such that the semiclassical approximation employed above is valid.

We first note that the system (4) is invariant under E_M → − E_M as well as Δ → − Δ due to the transformations σ^x,y_k → (−1)^kσ^x,y_k and σ^x,z_k → − σ^x,z_k, respectively. Thus without loss of generality, we assume in the following E_M,Δ > 0. The phase diagram of the spin model (4) contains four phases:⁵ a paramagnetic phase (PM) with a unique ground state with 〈σ^z_k〉 < 0, 〈σ^x_k〉 = 0; an antiferromagnetic phase (AFM) with doubly degenerate ground state and 〈σ^z_k〉 = 0, 〈σ^x_k〉∝(−1)^k; an 'antiphase' (AP) with a doubly degenerate ground state with 〈σ^z_k〉∝(−1)^k; and a 'floating phase' (FP) between the AFM and the AP. For Δ = 0, the duality transform of the model (4) is a sum of two quantum Ising chains, while for the noninteracting case it reduces to a single quantum Ising chain.

As we have shown above in our realization of the ANNNI model the parameters Δ and U can be easily tuned via a gate voltage. On the other hand, the coupling E_M is determined by the overlaps of the Majorana wave functions on neighbouring islands and thus fixed by the geometry of the array. Hence it is natural to consider the phase diagram as a function of Δ/E_M and U/E_M, which is shown in figure 2. For fixed Δ/E_M > 1 and sufficiently small interaction U the system is in a trivial (band-insulating) phase corresponding to the PM in the effective spin model, which is characterized by a unique ground state with $\langle \mathcal {P}_k\rangle <0$. By increasing U we cross into a topological phase (corresponding to the AFM) with $\langle \mathcal {P}_k\rangle =0$ and two degenerate ground states | ± 〉, distinguished by the (total) fermion parity $\prod _{k}\mathcal {P}_k |\pm \rangle = \pm |\pm \rangle$. The parity protection (also called topological protection) originates from the fact that any fermionic perturbation conserves the fermion parity and thus cannot mix the states | ± 〉. The phase transition between the trivial and topological phases is, for U > 0, located at

$$\begin{equation} 1-\frac{2U}{\Delta} = \frac{E_{\mathrm{M}}}{\Delta} - \frac{U E_{\mathrm{M}}^2}{2 \Delta^2 (\Delta -U)}. \end{equation} \tag{ 5 }$$

In particular, we find that Δ = E_M + 3U/2 for |U| ≪ E_M and Δ = 2U for U ≫ E_M.

Zoom In Zoom Out Reset image size

At large positive U we eventually enter incommensurate and commensurate CDW (Mott-insulator) states corresponding to the FP and the AP, respectively⁶. This region of the phase diagram cannot be reached as long as the induced charge q on the islands is homogeneous. However, replacing the common back gate by individual gates for each island yields the model (4) with site-dependent parameters Δ_i = Γ_Δ cos(πq_i/e) and U_i = Γ_U cos(π(q_i + q_i+1)/e) where q_i = G_gVⁱ_g denotes the induced charge on island i. Now using q_i = (−1)ⁱq one can enter the Mott phase for q → e/2. Note that the two degenerate ground states in the Mott insulator phase have the same fermion parity and thus are not parity protected.

The main feature of the phase diagram is its strong anisotropy under U → − U. While negative interactions suppress the topological phase, for U > 0 the ordering tendencies of the first and second terms in (4) compete with each other. In particular, starting from a noninteracting point in the trivial phase, i.e. U = 0 and Δ > E_M, competition between the band insulator and the Mott insulator will drive the system into the topological phase irrespective of the value of Δ/E_M.

As discussed above the topological phase is characterized by the existence of a doubly degenerate ground state with different fermion parities. In a finite system this degeneracy is lifted and the value of the resulting gap δ depends on the number of islands N as well as the system parameters Δ/E_M and U/E_M. On the other hand, the existence of Majorana end modes is protected by the gap between the (nearly) degenerate ground states and the first excited state, which we denote by Δ_T. This gap is given by Δ_T = 4E_M at Δ = U = 0 and decreases when going away from this point. However, as we show in the inset of figure 2, it remains of the order of E_M and significantly larger than δ for U ≠ 0 as long as one stays away from the phase boundaries.

After presenting the phase diagram, we now turn to its experimental signatures in the tunnelling conductance. The parameters Δ and U can be directly tuned through the induced charge q on the islands via a common back gate voltage V_g making it possible to choose the parameters such that the path crosses one or more phase boundaries.

Specifically, we consider the path c₁(q) = (Γ_U cos(2πq/e),Γ_Δ cos(πq/e)) with Γ_U = 0.3 Γ_Δ = 0.54 E_M such that the starting and end points at q = 0 and e/2 lie in the topological phase, while the path enters the trivial phase in between (see the red curve in figure 2). In the following, we consider an open chain of N islands. Following [20], we couple the system to an electronic lead with the tunnel Hamiltonian $H_{\mathrm { T}} = t \gamma _{1a}\sum _{k\sigma } (c_{k\sigma } - c^{\dagger} _{k\sigma })$, where c_kσ are the annihilation operators of the electrons in the lead and t is the tunnelling amplitude. We assume a constant density of states ρ₀ in the lead such that the (bare) tunnelling probability is given by Γ₀ = 2πt²ρ₀. We have calculated the differential Andreev conductance G(V ) using exact diagonalization for an open chain of length N = 24 taking only the two lowest energy states (with different fermion parities) into account⁷. In figure 3(a), the broadened conductance

$$\begin{equation} \bar{G}= \int_{-\Delta\omega/2}^{\Delta\omega/2}\frac{\mathrm{d}(eV)}{\Delta\omega}G(V) \end{equation} \tag{ 6 }$$

is plotted along c₁; the broadening $\Delta \omega \simeq \max (eV,k_{\mathrm {B}} T)$ is given by the maximum of bias voltage V and temperature T. For strong coupling Γ₀ to the lead as well as small broadening Δω, the phase boundaries are clearly visible as points where the conductance jumps from one to zero and vice versa. Weaker coupling to the lead will lead to a suppression of the conductance in the topological phase, i.e. unitary conductance cannot be observed. On the other hand, a larger broadening will eventually smear out all transitions. Thus in order to enable an experimental detection of the phase boundaries the energy broadening should be small (in units of E_M) while the coupling to the lead has to be sufficiently strong. Recent experiments on proximity coupled nanowires [4] indicate that a Majorana coupling E_M/k_B ≃ 100 mK is realistic, which is well in the range of experimental accessible temperatures. The Majorana coupling sets the topological gap Δ_T (see the inset of figure 2) as the other couplings (E_J and E_C) can be designed in a large parameter range [21].

Zoom In Zoom Out Reset image size

In order to study the regime of strong interactions we use the setup with individual back gates and induced charges q_i = (−1)ⁱq. In this setup, we consider the path c₂(q) = (Γ_U,Γ_Δ cos(πq/e)) with Γ_U = 0.3 Γ_Δ = 1.5 E_M (see the green curve in figure 2). The conductance $\bar {G}$ along c₂ is shown in figure 3(b). At q = 0 the path starts deep in the band-insulator phase and enters the topological phase at q ≈ 0.2 where the conductance becomes nonzero. When further increasing q we observe weak oscillations which are due to the finite-size lifting δ of the ground state degeneracy. For larger values of q the conductance stays zero in the incommensurate and commensurate CDW phases with the exception of the phase transition between them, where the conductance is nonzero due to finite-size effects.

Above we have shown how to realize the 1D ANNNI model using Josephson junction arrays and how to map out its phase diagram by measuring the tunnelling conductance. In this sense our proposed setup constitutes a quantum simulator for the 1D ANNNI model. In particular, the experimental control over the system parameters like the gate voltages allows us to study the effects of disorder on the phase diagram.

Interacting electrons in proximity-coupled semiconducting nanowires are described by the microscopic Hamiltonian [8]

$$\begin{eqnarray} &&\fl H_{\rm\scriptstyle NW} = - \! \int\!\mathrm{d}x\,\Psi^{\dagger}\! \left( \frac{\partial_x^2}{2m} +\mu + \mathrm{i}\, \alpha \sigma^y \partial_x + E_{\mathrm{Z}} \sigma^z \!\right) \Psi+ \! \int\!\mathrm{d}x(\Delta_{\mathrm{s}} \Psi_{\uparrow}\Psi_{\downarrow} + {\rm h.c.} + U_0 |\Psi_{\uparrow}(x)|^2 |\Psi_{\downarrow}(x)|^2)\nonumber\\ \end{eqnarray} \tag{ 7 }$$

with Ψ(x) = (Ψ_↑(x),Ψ_↓(x))^T the electron field operator, m the electron mass, μ the chemical potential, α the strength of the spin–orbit coupling, E_Z the Zeeman energy due to the applied magnetic field, Δ_s the s-wave pairing amplitude and U₀ the (short-range) Coulomb interaction. For sufficiently strong Zeeman energy (compared to the other energy scales), we only need to consider a single band similar to the ANNNI model discussed above. However, projecting the Hamiltonian (7) onto a single band strongly reduces the effect of the interaction. Specifically, we find for the effective interaction strength in the single-band model U/E_M = mU₀α²/ℏ²LE²_Z ≪ 1 with L the length of the nanowire. In this way, interacting nanowires subject to a strong Zeeman field are always in the weak coupling regime [8]. In contrast, as we showed above, the strong interaction regime for spinless fermions is readily accessible in the case of nanowires in Josephson junction arrays.

We have analysed a 1D array of Josephson junctions featuring Majorana modes, where capacitances between adjacent islands lead to interactions between the Majorana modes. We have shown that repulsive interactions generically facilitate the topological phase due to their competition with the on-site charging energies. Finally, we have proposed a tunnelling experiment to detect the phase boundaries.

We have benefited from discussions with Eran Sela and Kirill Shtengel and thank Bernd Braunecker for useful comments. This work was supported by the Alexander von Humboldt Foundation (FH) and the DFG through the Emmy–Noether program (DS).

In this appendix, we present more information about the derivation of the rates Γ_Δ and Γ_U. Starting with Γ_Δ, we are interested in the event that the phase ϕ_k of a single island changes by 2π. The relevant tunnelling matrix element t_Δ∝〈2π|e^−iHt|0〉 can be evaluated in the path-integral formalism by going to imaginary time τ = i t, cf [23–25],

$$\begin{equation} t_\Delta \propto \int \mathcal{D} [\phi_k]\, \mathrm{e}^{-\hbar^{-1} \int\!\mathrm{d}\tau\,L_E} \end{equation} \tag{ A.1 }$$

subject to the boundary condition ϕ_k(0) = 0 and $\Delta \phi _k=\phi _k(\infty ) -\phi _k(0) \in 2\pi \mathbb {Z}\setminus \{0\}$. To exponential accuracy, the path integral is dominated by the classical paths $\phi _k^{{\mathrm {\scriptstyle cl}},n}$ which minimize the action $S_E=\hbar ^{-1}\int \!\mathrm {d}\tau \,L_E$, i.e.

$$\begin{equation} t_\Delta \sim \sum_{n} \mathrm{e}^{- S_{E}[\phi^{{\rm\scriptstyle cl},n}_{k}]}, \end{equation} \tag{ A.2 }$$

where n is an index enumerating the different paths in the case when there are different minima of the action.

For η = 0, the relevant part of the action is well approximated by (keeping only terms which depend on ϕ_k)

$$\begin{equation} S_{E}[\phi_{k}]= \int_{0}^\infty\!\mathrm{d}\tau \underbrace{\left[ \overbrace{\frac{\hbar}{16 E_{C}} \phi_{k}^{\prime 2}}^{T} - \frac{\mathrm{i}\,q}{2e} \phi'_{k} + \overbrace{\frac{E_{\mathrm{J}}}{\hbar} (1-\cos\phi_{k})}^{V_{\mathrm{J}}} \right]}_{L_{E}}, \end{equation} \tag{ A.3 }$$

where ' denotes the derivative with respect to τ and we have neglected the potential proportional to E_M ≪ E_J. As the action does not depend directly on τ, the energy along the classical path minimizing the action is conserved,

$$\begin{equation} E= \frac{\partial L_E}{\partial \phi_k'} \phi_k' - L_E = T - V_{\mathrm{J}}. \end{equation} \tag{ A.4 }$$

For τ = 0 we have T = V_J = 0 such that E = 0 in our case.

We can express the kinetic energy in terms of the potential and obtain

$$\begin{equation} T= \frac{\hbar}{16 E_C} \phi_k'^2 = E+ V_{\mathrm{J}}, \end{equation}\noindent \tag{ A.5 }$$

with which we can get an alternative expression for the measure (the capacitance matrix acts as a metric)

$$\begin{equation} \mathrm{d} \tau = \sqrt{ \frac{\hbar}{16 E_C(E+V_{\mathrm{J}})} }\, \mathrm{d}\,\phi_k. \end{equation}\noindent \tag{ A.6 }$$

Due to the conservation of energy, we can go over to the Euler–Maupertuis action $S_0 = S_E+ \int \!\mathrm {d}\tau \, E$ (note that in our case E = 0 such that S₀ is in fact equal to S_E) which can be rewritten employing (A.6) as

$$\begin{equation} S_0[\phi_k] = \sqrt{\frac{\hbar}{4 E_C}} \int\!\mathrm{d}\phi_k\,\sqrt{E+V_{\mathrm{J}}} - \frac{\mathrm{i}\, q}{2e} \Delta \phi_k, \end{equation}\noindent \tag{ A.7 }$$

which is independent of imaginary time and only depends on the path chosen [26].

The action is minimized for Δϕ_k = ± 2π as each additional phase slip by 2π increases S₀. The expression corresponding to the first term in (A.7) is independent of ±. The classical path corresponds to increasing ϕ_k by 2π such that, cf [17],

$$\begin{equation} \sqrt{\frac{\hbar}{4 E_C}} \int\!\mathrm{d}\phi_k\,\sqrt{E+V_{\mathrm{J}}} = \sqrt{\frac{E_{\mathrm{J}}}{4 E_C}} \int_0^{2\pi} \!\mathrm{d}\phi_k\, (1-\cos \phi_k) = \sqrt{\frac{8 E_{\mathrm{J}}}{E_C}}. \end{equation}\noindent \tag{ A.8 }$$

We obtain the final result

$$\begin{equation} t_\Delta \sim \mathrm{e}^{-\sqrt{8 E_{\mathrm{J}}/E_C}} \sum_{\Delta\phi_k=\pm 2\pi} \mathrm{e}^{-\mathrm{i}q \Delta \phi_k/2e} \sim \mathrm{e}^{-\sqrt{8 E_{\mathrm{J}}/E_C}} \cos(q \pi/e). \end{equation}\noindent \tag{ A.9 }$$

valid up to exponential accuracy. In the main text, we use the result t_Δ = Γ_Δcos(πq/e) with $\Gamma _\Delta \simeq E_C^{1/4} E_{\mathrm {J}}^{3/4} \mathrm {e}^{-\sqrt {8 E_{\mathrm {J}}/E_C}}$. In fact, the charge q should be replaced by q + en_k. The reason is that due to the Majorana term the action is not 2π but only 4π periodic in ϕ_k. In the calculation above, we however assume the action to be 2π periodic. In fact, the action can be made 2π periodic by a gauge transformation at the expense of replacing q↦q + en_k. More information on this rather subtle point can be found in [14, 15].

The prefactor E^1/4_CE^3/4_J of Γ_Δ depends on the shape of the potential close to the turning points ϕ_k ≈ 0,±2π and cannot be obtained in our simple semiclassical analysis which only captures the physics up to exponential accuracy. However, the scaling of the prefactor can be obtained from summing up the instanton contributions [24] or by matching it to the exact solution of the Mathieu equation [17]. In our case, the potential is always given by $V_{\mathrm {J}} \simeq \frac 12 E_{\mathrm {J}} \phi _k^2$ for ϕ_k ≪ 1 such that the same prefactor E^1/4_CE^3/4_J (from the Mathieu equation) appears for all the tunnelling amplitudes.

In the case η ≠ 0, it is important to note that the phases on the different islands do not completely decouple. In lowest order in η, we need to take the phases ϕ_k±1 on the islands k ± 1 into account. Due to the symmetry of the problem, we know that ϕ_k−1(τ) = ϕ_k+1(τ). The relevant part of the action reads

$$\begin{eqnarray} &&\fl S_E = \int_0^\infty d\tau\left[\vphantom{\underbrace{3 - \cos \phi_k -2 \cos \phi_{k+1}}_{V_{\mathrm{J}}}} \frac{\hbar }{16E_C} \bigl(\phi_k'^2 + 2 (1-\eta)\phi_{k+1}'^2 \bigr) -\frac{\mathrm{i}\, q}{2e} (\phi'_k + 2 \phi'_{k+1}) + \underbrace{\frac{E_J}{\hbar}(3 - \cos \phi_k -2 \cos \phi_{k+1})}_{V_{\mathrm{J}}} \right]\!.\nonumber\\ \end{eqnarray} \tag{ A.10 }$$

Following the same line of calculation as going from (A.3) to (A.7), we obtain

$$\begin{equation} S_0 = \sqrt{\frac{\hbar}{4 E_C} } \! \int\!\mathrm{d}\phi_k \underbrace{\sqrt{[1 + 2 (1-\eta) \phi_{k+1}'(\phi_k)^2] V_{\mathrm{J}}}}_{\mathcal{L}_{\rm\scriptstyle eff}} - \frac{\mathrm{i}\, q \Delta\phi_k}{2e}, \end{equation}\noindent \tag{ A.11 }$$

where we expressed the path by giving ϕ_k+1 as a function of ϕ_k. As we are interested in processes where ϕ_k changes by Δϕ_k = ± 2π, we need to find $\phi _{k+1}(\phi _k)\colon [0,\pm 2\pi ] \mapsto \mathbb {R}$ with ϕ_k+1(0) = ϕ_k+1(±2π) such that the action is minimized. We will find the solution which corresponds to Δϕ_k = 2π below. The second solution with Δϕ_k = −2π can be obtained via the symmetry ϕ_i↦ − ϕ_i,∀i of the Lagrangian.

As the second term in (A.10) is independent of the path (it depends only on the boundary condition), we only need to minimize the first term. The extremum is attained when the function ϕ_k+1(ϕ_k) fulfils the Euler–Lagrange equation

$$\begin{equation} \frac{\mathrm{d}}{\mathrm{d}\phi_k} \frac{\partial \mathcal{L}_{\rm\scriptstyle eff}}{\partial \phi_{k+1}'} = \frac{\partial \mathcal{L}_{\rm\scriptstyle eff}}{\partial \phi_k}. \end{equation}\noindent \tag{ A.12 }$$

To first order in η (assuming that $\phi_{k+1} \in \Or(\eta)$), the equation assumes the form

$$\begin{eqnarray*} &&4 (1-\cos \phi_k)\phi_{k+1}'' +2 \sin (\phi_k) \, \phi_{k+1}' - \eta \sin \phi_k -2 \phi_{k+1} =0. \end{eqnarray*}$$

Employing the substitution ϕ_k+1(ϕ_k) = tan(ϕ_k/4)f(ϕ_k) reduces the equation to a first-order equation in f' of the form

$$\begin{eqnarray} &&\fl 2\,\cos (\phi_k/4) [8 \sin^2(\phi_k/4) f''(\phi_k) -\eta \cos (\phi_k/2) ]+ 2 [3 \sin(\phi_k/4) + \sin(3 \phi_k/4)] f'(\phi_k) =0\nonumber\\ \end{eqnarray} \tag{ A.13 }$$

which can be integrated with the solution

$$\begin{equation} f(\phi_k) = \eta \log\left[\frac12 \sin \left(\frac12\phi_k\right)\right] - \eta \frac{\log \cos (\frac14 \phi_k)}{ \sin^{2}(\frac14 \phi_k)}. \end{equation} \tag{ A.14 }$$

Plugging the solution into (A.11) and retaining the first nonvanishing term in η yields

$$\begin{equation} S_0 = \sqrt{8 E_{\mathrm{J}}/E_C} \left[ 1 + \frac{\pi^2-12}{96} \eta^2 + \Or(\eta^4) \right] - \frac{\mathrm{i}\, q \Delta \phi_k}{2e}. \end{equation} \tag{ A.15 }$$

Summing up the two contributions with Δϕ_k = ± 2π, we obtain a term proportional to cos(πq/e) as before, with the proportionality constant given by $\Gamma _\Delta = E_C^{1/4} E_{\mathrm {J}}^{3/4} \mathrm {e}^{- \sqrt {8 E_{\mathrm {J}}/E_C} [1 + (\pi ^2-12)\eta ^2/96 ]}$.

For η ≠ 0, we get additionally a next-nearest-neighbour interaction due to phase slips where both ϕ_k and ϕ_k+1 change by 2π. In fact, due to the symmetry of the problem, we can set ϕ_k+1(τ) = ϕ_k(τ). The terms of the action which change with ϕ_k and ϕ_k+1 are given by (A.3) for each of the islands and the additional contribution of the cross capacitance. Thus, we have

$$\begin{equation} S_E = \int_0^\infty\!\mathrm{d}\tau\, \left[ \frac{\hbar (2 - \eta)}{16 E_C} \phi_k'^2 + 2 E_{\mathrm{J}} (1- \cos \phi_k) \right], \end{equation} \tag{ A.16 }$$

which leads to

$$\begin{eqnarray} S_0 &=& \sqrt{ \frac{(2-\eta) E_{\mathrm{J}}}{2 E_C} } \int\!\mathrm{d}\phi_k\, \sqrt{ 1-\cos \phi_k} - \frac{\mathrm{i}\, q (\Delta \phi_k + \Delta \phi_{k+1})}{2e} \nonumber\\ &= &4 \sqrt{(2 -\eta) E_{\mathrm{J}}/E_C}- \frac{\mathrm{i}\, q (\Delta \phi_k + \Delta \phi_{k+1})}{2e} . \end{eqnarray} \tag{ A.17 }$$

As Δϕ_k+1 = Δϕ_k (because the two phases slip together), the tunnelling amplitude thus assumes the form Γ_U cos(2πq/e) with $\Gamma _U = E_C^{1/4} E_{\mathrm {J}}^{3/4} \mathrm {e}^{- \sqrt {16(2-\eta ) E_{\mathrm {J}}/E_C}}$.