Strongly interacting Majorana modes in an array of Josephson junctions

An array of superconducting islands with semiconducting nanowires in the right regime provides a macroscopic implementation of Kitaev's toy model for Majorana wires. We show that a capacitive coupling between adjacent islands leads to an effective interaction between the Majorana modes. We demonstrate that even though strong repulsive interaction eventually drive the system into a Mott insulating state the competition between the (trivial) band-insulator and the (trivial) Mott insulator leads to an interjacent topological insulating state for arbitrary strong interactions.


Introduction
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 this exotic states in semiconducting nanowires with strong spin-orbit coupling in a magnetic field and in proximity to 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 conductingsuperconducting 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 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, has been proposed [14] in a one-dimensional (1D) array of topological superconducting islands. Its advantage over microscopic implementations is that the individual parameters of the effective Hamiltonian can potentially be 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. Figure 1. A 1D array of superconducting islands (light grey) coupled via strong Josephson junctions (E J ) to a common ground superconductor. Each island contains a pair of Majorana zero modes (white dots) at the end points of a semiconducting nanowire (dark grey). The tunnel coupling of individual electrons between the superconducting islands is proportional to the energy scale E M . A common gate voltage V g can be used to tune the relative strength of the different terms in the Hamiltonian. The capacitive couplings between the elements are denoted by C, C J , and C g , respectively.

Model
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] where φ k is the superconducting phase of the k-th island. Here and in the following, we assume for convenience periodic boundary condition 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 J = E J N k=1 (1 − cos φ 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 Finally, the kinetic term k ] + ( /2e) N k=1 qφ 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 g V 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 2 /2C Σ of a single islands depends on the total capacitance C Σ = 2C + C G .

Mapping on an effective spin model
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] H Ck = Γ ∆ cos[π(q/e + n k )] = Γ ∆ cos(πq/e)P k . ( Here, the tunnelling amplitude is given by (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 [18,14]. The term n k = 1 2 (1 − P k ) ∈ {0, 1} is the contribution to the charge due to the parity P k = iγ ka γ kb ∈ {−1, 1} of the number of electrons on the superconducting island encoded in the state of the Majorana zero modes [15,14]. 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 parametrized 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 ‡ Because V J ≫ V M , we take only the potential V J into account when calculating the tunnelling action S E .
in η, we only need to take into account the change in φ k−1 and φ k+1 and obtain S ∆ = 8E J /E C [1 + (π 2 − 12)η 2 /96], which is accurate to more than 2 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 : The Aharonov-Casher interference in this case leads to the tunnelling amplitude depending on the total charge 2q + e(n k + n k+1 ) of the two islands involved. Thus we obtain the interaction term 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 which constitutes the transverse axial next-nearestneighbour Ising (ANNNI) model [19] as can be seen by performing a Jordan-Wigner transformation P k = iγ ka γ kb = σ z k , iγ kb γ k+1a = σ x k σ x k+1 , resulting in the spin Hamiltonian Here, the σ x,y,z k are Pauli matrices, ∆ = Γ ∆ cos(πq/e) and U = Γ U cos(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 before [13], we focus here on the effects of the interaction term. In our set-up, this term is most , 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 one such that the semiclassical approximation employed above is valid.

Phase diagram
We first note that the system (4) The phase diagram of the resulting model including expressions for the phase transitions has been worked out in [19,22] in the parameters κ = U/∆ and E M /∆. 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 "anti phase" (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 easily be 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 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 an unique ground state with P k < 0. By increasing U we cross into a topological phase (corresponding to the AFM) with P k = 0 and two degenerate ground states |± , distinguished by the (total) fermion parity k P k |± = ±|± . 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 phase is, for U > 0, located at In particular we find ∆ = E M + 3U/2 for |U| ≪ E M and ∆ = 2U for U ≫ E M . At large positive U we eventually enter incommensurate and commensurate CDW (Mott insulator) states corresponding to the FP and the AP, respectively. This region We note that the existence of the FP, and thus the incommensurate CDW state, at small ∆/E M has 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 g V i g denotes the induced charge on island i. Now using q i = (−1) i 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 term 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-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 parity. 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 in 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.

Tunnelling conductance
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 1 (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 q = 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 T = tγ 1a kσ (c kσ − c † kσ ), 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 ρ 0 in the lead such that the (bare) tunnelling probability is given by Γ 0 = 2πt 2 ρ 0 . 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 not yet been fully established [22]. fermion parity) into account. ¶ In figure 3(a), the broadened conductancē is plotted along c 1 ; the broadening ∆ω ≃ max(eV, k B T ) is given by the maximum of bias voltage V and temperature T . For strong coupling Γ 0 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 inset of figure 2) as the other couplings (E J and E C ) can be designed in a large parameter range [21]. In order to study the regime of strong interactions we use the set-up with individual back gates and induced charges q i = (−1) i q. In this set-up we consider the path c 2 (q) = (Γ U , Γ ∆ cos(πq/e)) with Γ U = 0.3 Γ ∆ = 1.5 E M (see the green curve in figure 2). The conductanceḠ along c 2 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 ¶ The two-level approximation is appropriate in the topological phase (where we have two levels separated by ∆ T ≫ δ from the rest), as well as in the trivial phase (where there is a unique ground state) and in the Mott phase (where there are two degenerate ground state with the same fermion parity) as δ ∆ω in the latter cases such that the current vanishes.
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 set-up constitutes a quantum simulator for the 1D ANNNI model. In particular, the experimental control over the system parameters like the gate voltages allows to study the effects of disorder on the phase diagram.

Relation to nanowires
Interacting electrons in proximity-coupled semiconducting nanowires are described by the microscopic Hamiltonian [8] ) 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 0 the (shortrange) 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 0 α 2 / 2 LE 2 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.

Conclusions
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.
subject to the boundary condition φ k (0) = 0 and ∆φ k = φ k (∞) − φ k (0) ∈ 2πZ \ {0}. To exponential accuracy, the path-integral is dominated by the classical paths φ cl,n k which minimize the action S E = −1 dτ L E , i.e., where n is an index enumerating the different paths in the case that 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 ) 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, 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 with which we can get an alternative expression for the measure (the capacitance matrix acts as a metric) Due to the conservation of energy, we can go over to the Euler-Maupertuis action S 0 = S E + dτ E (note that in our case E = 0 such that S 0 is in fact equal to S E ) which can be rewritten employing (A.6) as which is independent on 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 0 . The expression corresponding to the first term in (A.7) is independent on ±. The classical path corresponds to increasing φ k by 2π such that, cf [17], We obtain the final result valid up to exponential accuracy. In the main text, we use the result t ∆ = Γ ∆ cos(πq/e) with Γ ∆ ≃ E In fact the charge q should be replace 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 on the expense of replacing q → q + en k . More information on this rather subtle point can be found in [15,14].
The prefactor E 1/4 C E 3/4 J of Γ ∆ depends on the shape of the potential close to the turning points φ k ≈ 0, ±2π and cannot be obtain 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 always is given by V J ≃ 1 2 E J φ 2 k for φ k ≪ 1 such that the same prefactor E In the case η = 0, it is important to notice 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 (A.10) Following the same line of calculation as going from (A.3) to (A.7), we obtain 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 φ k+1 (φ k ): [0, ±2π] → 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 by 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 To first order in η [assuming φ k+1 ∈ O(η)], the equation assumes the form . (A.14) Plugging the solution into (A.11) and retaining the first nonvanishing term in η yields S 0 = 8E J /E C 1 + π 2 − 12 96 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 J /E C [1+(π 2 −12)η 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 term 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 16) which leads to S 0 = (2 − η)E J 2E C dφ k 1 − cos φ k − iq(∆φ k + ∆φ k+1 ) 2e = 4 (2 − η)E J /E C − iq(∆φ k + ∆φ k+1 ) 2e . (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 Γ U = E