Numerical investigation of gapped edge states in fractional quantum Hall-superconductor heterostructures

Abstract

Fractional quantum Hall-superconductor heterostructures may provide a platform towards non-abelian topological modes beyond Majoranas. However their quantitative theoretical study remains extremely challenging. We propose and implement a numerical setup for studying edge states of fractional quantum Hall droplets with a superconducting instability. The fully gapped edges carry a topological degree of freedom that can encode quantum information protected against local perturbations. We simulate such a system numerically using exact diagonalization by restricting the calculation to the quasihole-subspace of a (time-reversal symmetric) bilayer fractional quantum Hall system of Laughlin ν = 1/3 states. We show that the edge ground states are permuted by spin-dependent flux insertion and demonstrate their fractional 6π Josephson effect, evidencing their topological nature and the Cooper pairing of fractionalized quasiparticles. The versatility and efficiency of our setup make it a well suited method to tackle wider questions of edge phases and phase transitions in fractional quantum Hall systems.

Introduction

The fractional quantum Hall¹ (FQH) effect harbors a variety of exotic topologically ordered quantum phases, from the well understood Laughlin states all the way to states^2,3 with non-Abelian quasiparticles such as Majorana fermions and Fibonacci anyons. The experimental exploration of these systems faces two challenges: First, the desired more exotic topological orders, which could be used, e.g., for universal quantum computation,⁴ can only be accessed under extreme experimental conditions, as they are protected by very small energy gaps. Thus, despite intense efforts, definite experimental confirmation of the non-Abelian nature of a FQH phase is still lacking. Second, the topological information is encoded in degenerate ground states or the state of quasiparticles, and is therefore intrinsically hard to measure and manipulate.

To overcome both of these obstacles, several recent studies proposed focusing on more conventional FQH states, such as the Laughlin ν = 1/3 state, and to add twist defects.^5,6,7,8 Among the various physical implementations of these twist defects, one deliberately couples the edge states.^{9,10,11,12,13,14} In these proposals, the topological excitations could be localized at domain walls of differently gapped edge segments. Advantages of this approach include the comparably large gap of the Laughlin ν = 1/3 state, which protects the quantum state, and the localized nature of the topological excitations, which facilitates their measurement and manipulation. For example, refs.^15,16 use parafermion excitations at domain walls between magnetic and superconducting gapped regions of FQH edges to build topological order akin to the ${\Bbb Z}_3$ Read-Rezayi state. Barkeshli subsequently pointed out that topological information is also stored in a pair of counter-propagating ν = 1/3 edge states that are fully gapped out by a superconducting order parameter.¹⁷ The gapped edge, on which a pairing between fractionalized quasiparticles is induced, has a well defined quantized total charge that can take the values 0, 2e/3 and 4e/3 (modulo the Cooper pair charge 2e). This nonlocal observable defined along the (closed) edge distinguishes three topologically degenerate ground states of the edge. By appropriately coupling several gapped edges, one can in principle manipulate their topological ground state.¹⁷ Another approach to engineer parafermion excitations from Abelian topological order relies on lattice defects and was recently implemented numerically in refs.^18,19. The FQH edge states are also a convenient system to study the bulk-boundary correspondence in topologically ordered systems. Unlike noninteracting symmetry protected phases in two spatial dimensions, interacting integer and fractional quantum Hall states can support several distinct edge phases with different universal properties but the same symmetries.^20,21

While effective models (e.g., using a bosonized description of the edge^22,23) have permitted striking predictions at the edge of topologically ordered systems, open questions remain which can only be addressed by a microscopic approach. First, in the context of the bulk-boundary correspondence, which boundary phase is favored by certain microscopic interactions remains largely unexplored (especially when non-abelian liquids are used as the building blocks). As for localized edge modes, the braiding of non-abelian excitations relies on the possibility to tunnel quasiparticles through the bulk while keeping the edge gap open.¹⁰ This hypothesis relies on the hierarchy of energy and length scales in the system. Numerical simulations are necessary to achieve such quantitative analysis, and have the potential to identify challenges that could have been obscured by effective analytical models. They seem indispensable as experiments are undertaking the first steps to realize the ideas outlined above.^24,25

Here, for the first time, we undertake an extensive numerical calculation of a FQH system coupled to a superconductor using exact diagonalization. More explicitly, we consider a bilayer FQH system, with magnetic field perpendicular to the layers, where the orientation of the field for one layer is opposite to that for the other layer. This is equivalent to a time-reversal symmetric fractional topological insulator.²⁶ This construction permits gapping out of the edge states with singlet interlayer superconducting pairing. Our calculations are performed on a cylinder geometry in which the bilayer-FQH droplet has two edges. To make numerics feasible, we restrict our study to the subspace of zero energy bulk and edge excitations of the Laughlin ν = 1/3 state in each layer. These many-body quantum states being Jack polynomials,^27,28 we can perform an efficient evaluation of the microscopic model matrix elements in the reduced basis. Our setup can thus be straightforwardly generalized to study the edges of other FQH model states with similar properties but richer topological order (such as the Moore-Read state). It could also be cast into the matrix product state framework,^29,30,31 which should provide access to larger system sizes.

Results

Hamiltonian and effective Hilbert space

In the Landau gauge, a single-particle basis that spans the lowest Landau level on a cylinder of circumference L _y is given by

$$\phi _m(x,y) = \frac{1}{{L_y\ell _B\root {4} \of {\pi }}}e^{{\mathrm{i}}\frac{{2m\pi y}}{{L_y}}}e^{ - \frac{1}{{2\ell _B^2}}\left( {x - \frac{{2\pi m\ell _B^2}}{{L_y}}} \right)^2},$$

(1)

where $\ell _B$ is the magnetic length which we will set to unity in the rest of the manuscript. We truncate the single-particle Hilbert space of the cylinder by allowing for m to take integer values in the range −N_Φ/2 ≤ m ≤ + N_Φ/2 for some integer N_Φ. Note that m is either an integer or half integer depending on the parity of N_Φ. Here, m plays the role of both the y-momentum of the wave function and at the same time determines the location of the wave function along the x direction. This coupling of momentum and position is enforced by the lowest Landau level projection.

We now consider a bilayer system, where the two layers are distinguished by a spin index ↑, ↓ and the Hall effect in one layer is opposite in chirality to the Hall effect in the other layer. This is the case in so-called fractional topological insulators, where the labels ↑, ↓ may correspond to the physical spin and the spin-dependent magnetic field is akin to the spin-orbit coupling. An alternative scenario more relevant to traditional FQH experiments is one in which ↑, ↓ label the carriers in two adjacent quantum wells, one being electron-like and the other hole-like. A homogeneous magnetic field gives rise to edge states in both quantum wells and the direction of propagation in one quantum well is opposite to that in the other quantum well. Our numerical study applies equally well to each of these physical realizations, but we choose to describe our results using the terminology of the fractional topological insulator realization. In either case the single particle eigenstates are

$$\begin{array}{c}\psi _m^ \uparrow (x,y) = \phi _m(x,y)\cr \psi _m^ \downarrow (x,y) = \phi _m^ \ast (x,y) = \phi _{ - m}( - x,y),\end{array}$$

(2)

with −N_Φ/2 ≤ m ≤ + N_Φ/2. This ensures that the system is invariant under time-reversal symmetry ${\cal T} = K{\kern 1pt} {\mathrm{i}}\sigma _y$ for spinful fermions, where K is complex conjugation and σ _y is the second Pauli matrix acting in spin space. Note that none of the topological features we are interested in are protected by ${\cal T}$. In fact, in the electron-hole bilayer realization of the system, ${\cal T}$ is not the physical symmetry of the system, but an artificial symmetry of the model that may be broken in a microscopic realization.

In the Fock space spanned by the single-particle states $\psi _m^ \uparrow$ and $\psi _m^ \downarrow$, we consider a Hamiltonian of the form

$$\begin{array}{*{20}{l}} \hat H & = \hat H_{{\mathrm{2}} - {\mathrm{body}}, \uparrow } + \hat H_{{\mathrm{2}} - {\mathrm{body}}, \downarrow } + \mathop {\sum}\limits_{m = - s}^s {\kern 1pt} \mu _m\left( {\hat c_{m, \uparrow }^\dagger \hat c_{m, \uparrow } + \hat c_{m, \downarrow }^\dagger \hat c_{m, \downarrow }} \right) \\ & + \frac{1}{{2C}}\left( {\hat N - N_0} \right)^2 + {\mathrm{\Delta }}_0\mathop {\sum}\limits_{m = - s}^s {\kern 1pt} f_m\left( {\hat c_{m, \uparrow }^\dagger \hat c_{m, \downarrow }^\dagger + {\mathrm{h}}{\mathrm{.c}}{\mathrm{.}}} \right), \end{array}$$

(3)

where H_2−body,σ is the two-body interaction for fermions with spin σ = ↑, ↓ and $\hat c_{m,\sigma }^\dagger$ creates an electron in state $\psi _m^\sigma$. For H_2−body,σ we use the pseudo-potential Hamiltonian with V₁ being the only non-zero pseudo-potential coefficient. The coefficients μ _m describe a confining potential that is rotationally symmetric around the cylinder axis. The operator $\hat N = \hat N_ \uparrow + \hat N_ \downarrow$ measures the total number of particles. A schematic representation of this setup is sketched in Fig. 1.

Fig. 1

Schematics of the physical geometry and the one used for the numerical investigation. a Fractional topological insulator heterostructure in which carriers with spin up and down (red and blue) form a fractional quantum Hall state with opposite chirality. Proximity to superconducting reservoirs (yellow) induces a superconducting gap in their edge channels. To study the Josephson effect, relative phase φ between the left and right superconducting order parameter is included. b When imposing periodic boundary conditions along the edges, resulting in a cylinder geometry, each edge carries a topological degree of freedom. The boundary conditions can be twisted by inserting a flux ϕ into the cylinder for spin up electrons and −ϕ for spin down electrons. In the Landau gauge orbitals are localized along the cylinder, where we consider N_orb,n and N_orb,sc normal and superconducting orbitals, respectively. The typical separation between orbitals is $\frac{{2\pi \ell _B^2}}{{L_y}}$, where L _y is the cylinder perimeter. The droplet is confined by a linear potential μ _m . c With the counter-propagating edges gapped out, the bilayer FQH state on the cylinder is topologically equivalent to a single layer FQH state on a torus, where the fluxes ϕ and φ run through its two noncontractible cycles and can be used to explore its topological ground state degeneracy. It is thus topologically equivalent to the ground state degeneracy of the gapped edge modes

Due to the presence of mean-field superconductivity, the particle number is not conserved. To tune the system into a regime with finite particle number density, the charging energy of strength 1/(2C) has been added to the Hamiltonian. The two parameters C and N₀ permit tuning of the average number of particles in the system. Finally, Δ₀ is the overall strength of the superconducting coupling while f _m is the (dimensionless) variation of the superconducting order parameter along the cylinder, assuming a superconducting pairing potential that is rotationally symmetric along the cylinder axis. We will further assume a superconducting pairing potential that is nonzero only at the edge of the FQH droplet. In the electron-hole bilayer realization of the system, the term proportional to Δ₀ takes the role of a charge conserving backscattering term between the layers.

The Laughlin state edge and bulk quasihole excitations are the exact zero energy states of the model interaction $\hat H_{{\mathrm{2}} - {\mathrm{body}},\sigma }$ at filling ν = 1/3. Their corresponding wavefunctions have an analytical expression on the cylinder geometry.³² Being Jack polynomials, they can be conveniently decomposed into the occupation basis.³³ In ref. ³⁴, a careful and detailed numerical study of this state and its edge excitations was performed on the cylinder geometry using a confinement similar to the μ _m term of Eq. (3). In particular, the low energy spectrum (i.e., below the bulk energy gap) for a finite size quantum Hall droplet has the characteristic shape shown in Fig. 2a.

Fig. 2

Spectra of Hamiltonian Eq. (3) on a cylinder with circumference $L_y{\mathrm{/}}\ell _B = 15.0$, N_Φ = 21 flux quanta, 1/(2C) = 0.256271, N₀ = 30, and a linear symmetric chemical potential μ _m as shown in Fig. 1b. a Spectrum of a single layer of the cylinder with 6 fermions. Energies are shifted such that the ground state energy is zero. Shaded areas are a guide for the eye showing the ‘arcs’ in the relative momentum K _y arising from the edge states of the Hall droplet. The inset is a zoom on the low lying state within the gray box, showing the energy difference δE between the first excited state and the ground state. b and c are spectra for the full double layer with (b) Δ₀ = 0 while in (c) Δ₀ = 2.1 with N_orb,sc = 16 superconducting orbitals. Blue figures give the number of states in the circle. In (c), the superconducting gap in the edge states is apparent from the three topological “ground states” moving below the bottom of the next arc. Orange states have total spin S _z  = 0, gray states S _z  = 1, and states with higher spin only appear above this energy window. Both spectra symmetrically extend to K _y  → −K _y

To make progress in the numerical evaluation of Hamiltonian (6), we send the gap of the FQH state to infinity, i.e., we set V₁ → ∞, by projecting the Hamiltonian to the zero-energy subspace of H_2−body,↑ + H_2−body,↓. This is the space of Laughlin quasiholes in each layer. The densest state(s) in this subspace are the Laughlin FQH states with a filling fraction of 1/3 per layer and spin. With this projection to the Laughlin quasihole space in place, the Hilbert space dimension is dramatically reduced. In a second step, we diagonalize the chemical potential, superconducting, and charging energy terms in this quasihole space.

The parameters 1/2C and N₀ of the charging energy term are used to control (i) the average number of particles in the droplet and (ii) the energy difference to sectors with nearby particle numbers. We choose these parameters by diagonalizing the system in the absence of superconductivity. Note that N₀ is not equal to the particle number in the non-superconducting ground state, because the μ _m term also contributes an energy cost that depends on the particle number. We choose 1/2C and N₀ such that the ground state has a desired particle number $\tilde N$ and the lowest energy state in the sector with $\tilde N + 2$ particles is degenerate with the lowest energy state in the sector with $\tilde N - 2$ particles. The energy difference between the $\tilde N$ and the $\tilde N \pm 2$ sectors is chosen to be small enough that a moderate superconducting term can couple these sectors (as detailed in Supplementary Note 2).

The problem defined by Hamiltonian (6) has two good quantum numbers: these are the total spin measured by $\hat S_z = \frac{1}{2}\left( {\hat N^ \uparrow - \hat N^ \downarrow } \right)$ (which also encodes the fermion parity) and the relative (angular) momentum $K_y = \mathop {\sum}\nolimits_{i \in \uparrow } {\kern 1pt} m_i - \mathop {\sum}\nolimits_{j \in \downarrow } {\kern 1pt} m_j$, where m _i and m _j are the quantum numbers of the occupied ↑ -spin and ↓ -spin states, respectively, as defined in Eq. (2). In the electron-hole bilayer realization of the system, $\hat S_z$ is the particle number operator.

Spectral evidence for the topological edge states

We present the spectral features associated with the model defined by Hamiltonian (6) and sketched in Fig. 1: three low-energy states protected by an energy gap. We choose a symmetric confining potential $\mu _m = \left| m \right|$ that is linear in the orbital space index with a slope of 1 (which will serve as the unit of energy throughout this paper), as shown in Fig. 1b.

Decoupled layers

We first study the spectrum of two decoupled layers, i.e., in the absence of superconductivity when Δ₀ = 0, which is shown in Fig. 2b. It can be understood as a combination of two independent spectra of a confined Laughlin ν = 1/3 state on a cylinder with a full edge structure as displayed in Fig. 2a. In the presence of the linear confining potential $\mu _m = \left| m \right|$, the Laughlin state on the cylinder has a characteristic spectral feature as a function of K _y : focusing on a single layer with N/2 particles (assuming N even), the ground state has K _y  = 0. The lowest lying excitations appear in the momentum sector K _y  = ±N/2. δE denotes the energy difference between these states and the ground state. Further low-lying states are located in the sectors K _y  = nN/2, $n \in {\Bbb Z}$. The lowest lying states in the other momentum sectors are higher in energy, giving rise to an arc-like structure in the spectrum, as observed in Fig. 2a. These arcs are highlighted by the shaded region in Fig. 2. The lowest states at momenta K _y  = nN/2, $n \in {\Bbb Z}$ are those where the droplet has been rigidly moved by n orbitals, giving rise to an extra energy cost of about nν due to the higher chemical potential of the now occupied orbitals in comparison to the emptied ones. Since the center of mass of the wave function is moved by n orbitals, the change in center of mass momentum is K _y  = nN/2. The lowest excitations in other momentum sectors are local edge excitations or combinations thereof. Since they increase the size of the Hall droplet, they cost more energy than the rigidly moved Laughlin state.

Within these considerations, we can understand the spectrum of Fig. 2b as a finite-size representation of a collection of gapless FQH edge states. In particular we can understand the low-energy structure as superpositions of the states in the two layers ↑, ↓. We denote by $\left| {\sigma ,0} \right\rangle$, $\left| {\sigma , \pm N{\mathrm{/}}2} \right\rangle$ the three lowest states in each of the σ = ↑, ↓ sectors which occur at momenta 0 and ±N/2. The state $\left| { \uparrow ,0} \right\rangle \otimes \left| { \downarrow ,0} \right\rangle$ is then the nondegenerate ground state labelled by a blue “1” in Fig. 2b. The states $\left| { \uparrow , \pm N{\mathrm{/}}2} \right\rangle \otimes \left| { \downarrow ,0} \right\rangle$ and $\left| { \uparrow ,0} \right\rangle \otimes \left| { \downarrow , \pm N{\mathrm{/}}2} \right\rangle$ are four degenerate states at momenta K _y  = ±N/2 found at the bottom of the first arc in Fig. 2b at energy δE above the ground state. The states $\left| { \uparrow ,N{\mathrm{/}}2} \right\rangle \otimes \left| { \downarrow ,N{\mathrm{/}}2} \right\rangle$ and $\left| { \uparrow , - N{\mathrm{/}}2} \right\rangle \otimes \left| { \downarrow , - N{\mathrm{/}}2} \right\rangle$ are degenerate at K _y  = 0 and labelled by a blue “2” in Fig. 2b. They occur at exactly 2δE above the ground state.

Gapping the three-fold ground state

We now compare the low-energy structure of the system with zero (Fig. 2b) and non-zero superconducting pairing in the outer orbitals (Fig. 2c). The states that were at 2δE in the former system moved substantially below the ones formerly at δE. Thus, the spectrum cannot be decomposed into that of two independent layers anymore. Superconductivity coupled the layers. Closer inspection also reveals a tiny but nonvanishing lifting of the degeneracy between the two lowest-lying excited states at K _y  = 0. We interpret the three lowest states in the K _y  = 0 sector as the quasi-degenerate topological states of the edges and the gap above them as the superconducting gap induced in the counter-propagating FQH edge modes.

The three-fold ground state degeneracy of the gapped edge states can be understood as follows. By introducing a gap, the superconducting coupling turns the bilayer quantum Hall state with edges into a single-layer quantum Hall state on a manifold without boundary. As sketched in Fig. 1c, this manifold is topologically equivalent to a torus, where the space between the two layers becomes the interior of the torus. This is in line with the opposite sign of the Hall conductivity in the two layers, because the normal to the torus surface is also reversed. On the torus, a Laughlin state at filling ν = 1/3 has a three-fold ground state degeneracy. This degenerate ground state is thus topologically equivalent to the three ground states we observe in the superconducting bilayer system. Fractional quantum Hall ground states on the torus can be manipulated by inserting flux through the holes of the torus, see Fig. 1c. This will be demonstrated later via manipulating the flux ϕ and the Josephson phase φ, respectively.

More relevant to the physics of the bilayer heterostructure is an interpretation of the ground state degeneracy in terms of Cooper-paired Laughlin quasiparticles. Due to the mean-field superconducting order parameter, the particle number is only defined modulo 2. Assuming that the low-energy Laughlin quasiparticles (of charge e/3) that comprise the edge mode are Cooper-paired, this leaves three nonequivalent configurations for the charge of one edge: 0, 2e/3, and 4e/3—each modulo 2. Since the total number of particles of the system is quantized to integers, the two superconducting edge states cannot independently support any of these charge configurations. Rather, they either both have charge 0, or the left edge has charge 2e/3 and the right edge 4e/3, or vice versa. We thus expect a total of three nearly degenerate topological ground states from this consideration as well, in line with our numerical observation.

Beyond the two special cases shown in Fig. 2b, c, we have performed an extensive study of the spectral properties when varying the system parameters. Some results are given in Fig. 3 for the largest system size that can be reached. Another system size is discussed in Supplementary Note 1. We fix the total number of superconducting orbitals N_orb,sc, equally split between the two ends of the cylinder. The selected value is a compromise between fully covering the edge modes and a large enough non-superconducting region of N_orb,n consecutive orbitals where an incompressible liquid can develop (as depicted in Fig. 1b). For each perimeter L _y the charging energy parameter 1/2C is optimized as discussed in Supplementary Note 2 (N₀ being fixed for the full diagram to N₀ = 30). Instead of using L _y for the vertical axis, we have plotted the data as a function of the approximate width of the normal region, i.e., $\frac{{2\pi \ell _B^2N_{{\mathrm{o}}rb,n}}}{{L_y}}$. Such a quantity is more natural when comparing different system sizes.

Fig. 3

Characterization of the low-energy spectrum of Hamiltonian Eq. (3) as a function of $2\pi N_{{\mathrm{orb,n}}}\ell _B^2{\mathrm{/}}L_y$ and the strength of the pairing potential Δ₀. $2\pi N_{{\mathrm{orb,n}}}\ell _B^2{\mathrm{/}}L_y$ approximates the physical distance between the superconducting regions and can be tuned by varying L _y . The charging energy is optimized for each L _y using the procedure defined in Supplementary Note 2. Other parameters are identical to those of Fig. 2. a Gap between the three lowest states in the K _y  = 0 sector and the next excited states. Gray color indicates that the three lowest states do not have K _y  = 0. b Spread of the three lowest states in the K _y  = 0 sector as indicated in Fig. 2c divided by the gap. Gray color indicates the region in which the ratio exceeds 1 or where the three lowest states do not have K _y  = 0. c The largest eigenvalue q of the charge imbalance operator defined in Eq. (3) in the space of the tree lowest states (at K _y  = 0). Gray color is used if the three lowest states do not have K _y  = 0. d The difference r between the energy of the second lowest eigenstate at Josephson phases φ = 0 and φ = π, normalized by the spread and each time measured with respect to the ground state energy, as defined in Eq. (6). The closer r gets to 1, the better can the 6π Josephson effect be observed. The gray color has the same meaning as in (c)

In Fig. 3a, we show the energy gap above the three lowest energy states. We set the gap to zero if these three states do not have K _y  = 0. We also provide s, the ratio between the energy spread of the three lowest energy states and the gap as previously defined. We cap s to one or set it to one if the gap is zero or the three lowest energy states do not have the expected quantum numbers. To be able to claim that we have a low energy manifold made of these three states separated by a gap from the higher energy excitations, we need s < 1. The smaller s is, the closer to an exact degenerate manifold we are. As can be observed in Fig. 3b, we have a large region where s is small, beyond a critical value of Δ₀ depending on L _y .

Charge distribution

In this section we study the charge distribution between the left and right halves of the system. In a physical realization of the system, two scenarios should be distinguished. If the two edges are coupled to the same superconductor, which implies that they are phase coherent, no quantized charge can be associated with one edge alone. In contrast, if the two edges are gapped by independent superconducting reservoirs, they carry independent fractionally quantized charges. However, in this latter case the charging energy is expected to lift the ground state degeneracy and the states are not topological.

The observable that measures the charge disproportionation between the two superconducting edges is given by

$$\hat Q_{\mathrm{R}} - \hat Q_{\mathrm{L}} = \mathop {\sum}\limits_{\sigma = \uparrow , \downarrow } {\kern 1pt} \mathop {\sum}\limits_m {\kern 1pt} {\int} {\kern 1pt} {\mathrm{d}}x{\int}_0^{L_y} {\kern 1pt} {\mathrm{d}}y{\kern 1pt} {\mathrm{sgn}}(x) \times \left| {\psi _m^\sigma (x,y)} \right|^2\hat c_{m,\sigma }^\dagger \hat c_{m,\sigma },$$

(4)

where the origin of the x axis coincides with the center of the m = 0 orbital. It measures the charge difference between the left (x < 0) and right (x > 0) half of the system. We compute the expectation value of $\hat Q_{\mathrm{R}} - \hat Q_{\mathrm{L}}$ in the manifold formed by the three lowest states in the K _y  = 0 sector, yielding a 3 × 3 matrix. Since the product of left-right mirror and time-reversal symmetry leaves the system invariant, the eigenvalues of this matrix are constrained to ±q and 0, where q is an a priori unspecified real number.

Figure 4 shows the evolution of q with the strength of the superconducting pairing Δ₀. For Δ₀ = 0, we can understand the nearly quantized value q ≈ 4/3 by recalling that the three lowest K _y  = 0 states are comprised of one state for which both up- and down-spin droplets are centered around m = 0 and a pair of states in which they are both centered around m = ±1. The center of charge of the former is located exactly at x = 0, while the latter two states have an excess of charge ± q/4 right of x = 0 (i.e., x > 0) in each layer, and the opposite deficit left of x = 0 (i.e., x < 0). The operator $\hat Q_{\mathrm{R}} - \hat Q_{\mathrm{L}}$ is thus diagonal in this basis, with respective eigenvalues 0, ±q. In the thermodynamic limit, the x > 0 excess charge in the latter single layer states is equal to the quasiparticle-excitation charge ± 1/3, summing up to an expectation value q = 4/3 of $\hat Q_{\mathrm{R}} - \hat Q_{\mathrm{L}}$. In the absence of superconductivity, the relevant length scales are the perimeter L _y and the width of the droplet $2\pi N_{{\mathrm{orb}}}\ell _B^2{\mathrm{/}}L_y$ (where $N_{{\mathrm{orb}}} = 3\tilde N{\mathrm{/}}2 - 2$). The charge quantization will not be clearly observed unless both these length scales are large compared to the correlation length of the ground state (around $1.4\ell _B$ for the ν = 1/3 Laughlin state³⁵). In the inset of Fig. 4, we indeed observed a nearly quantized q≈4/3 when tuning the value of L _y to respect this criterion (obtained for $L_y{\mathrm{/}}\ell _B = 10.5$ and $2\pi N_{{\mathrm{orb}}}\ell _B{\mathrm{/}}L_y \simeq 9.6$).

Fig. 4

Largest eigenvalue q of the operator $\hat Q_{\mathrm{L}} - \hat Q_{\mathrm{R}}$, defined in Eq. (3), which measures the charge imbalance between the left half and the right half of the cylinder depicted in Fig. 1b in real space, computed in the manifold of the three lowest energy states with momentum K _y  = 0. The data was obtained with Hamiltonian Eq. (3) using the same parameters as in Fig. 2, except for $L_y{\mathrm{/}}\ell _B = 7.0$ (red, $2\pi N_{{\mathrm{orb,n}}}\ell _B{\mathrm{/}}L_y \simeq 5.4$), $L_y{\mathrm{/}}\ell _B = 10.5$ (orange, $2\pi N_{{\mathrm{orb,n}}}\ell _B{\mathrm{/}}L_y \simeq 3.6$) and $L_y{\mathrm{/}}\ell _B = 15.0$ (green, $2\pi N_{{\mathrm{orb,n}}}\ell _B{\mathrm{/}}L_y \simeq 2.5$) as indicated. Right of the respective dashed lines, the spread-to-gap ratio shown in Fig. 3b is less than one. Moreover the charge imbalance nearly vanishes in the parameter regime of optimal spread-to-spread ratio. The inset shows q as a function of L _y for Δ₀ = 0

Upon introducing the superconducting term Δ₀ ≠ 0, we observe a rapid decrease in q, reaching (and passing) zero near the value of Δ₀ that leads to an optimal spread-to-gap ratio of the nearly degenerate ground states. This can be observed by comparing Fig. 3b, c. In the latter, we present q as a function of Δ₀ and the physical distance between the superconducting edges. Two trends can be observed: (i) In the lower right corner of this parameter space, where the edges are in closest spatial proximity, q is the smallest. (ii) In contrast, the largest values of q are found in the upper left corner of this parameter space. However, in this limit of small L _y , corresponding to a thin cylinder, the charge distribution strongly varies with x even in the center of the droplet. This yields contributions to the expectation value of $\hat Q_{\mathrm{R}} - \hat Q_{\mathrm{L}}$ from the center of the droplet x = 0, such that the operator does not allow for measurement of edge properties only.

Given the limited system sizes we can study numerically, it is hard to infer the behavior of the system in the thermodynamic limit from this computation. We do, however, present data for one other system size in Supplementary Note 1 which shares the qualitative features discussed above with Fig. 3. In summary, we observed that the charge imbalance of the three lowest states in the K _y  = 0 sector evolves from being nearly quantized to 4e/3 in the limit of two decoupled layers to very small values when the states are nearly degenerate and separated by a gap from other excited states.

Spin-dependent flux insertion

To confirm the topological nature of the observed degenerate ground states, we perform a numerical charge pumping experiment. The adiabatic insertion of a magnetic flux ϕ along the cylinder axis is equivalent to changing the boundary conditions of the electronic wave functions from periodic to twisted by an angle 2πϕ/ϕ₀, where ϕ₀ is the flux quantum. In a Landau level, as ϕ is increased from 0 to ϕ₀, all single-particle orbitals are shifted by one unit of the quantum number m → m + 1. To see this, notice that changing the boundary conditions from periodic to twisted amounts to replacing m by m + ϕ/ϕ₀ in Eq. (1). In a Laughlin 1/3 state, every orbital has an average occupation 1/3, so that, in the thermodynamic limit, a fractional charge e/3 is pumped from one end of the cylinder to the opposite end in the process of the flux insertion.

We would like to utilize a flux insertion to transform the topological ground states, labeled as $\left| {{\mathrm{\Psi }}_0} \right\rangle$, $\left| {{\mathrm{\Psi }}_ + } \right\rangle$, $\left| {{\mathrm{\Psi }}_ - } \right\rangle$ by their charge imbalance 0 and ±q, into one another. As we will demonstrate, we can use charge pumping to permute these ground states. Since the two layers of our system are time-reversed partners with opposite Hall conductivities, we have to insert flux with opposite orientation for the ↑ -spin and ↓ -spin particles. Only then is a net charge pumped from one edge of the system to the other. We will refer to this as spin-dependent flux insertion (see Fig. 1b).

Suppose we start with a state $\left| {{\mathrm{\Psi }}_0} \right\rangle$ that has charge 0 on both edges. As unit ϕ₀ spin-dependent flux is adiabatically inserted, charge is transferred from the left to the right edge, so that the resulting state is $\left| {{\mathrm{\Psi }}_ + } \right\rangle$. The other ground states are expected to transform into one another analogously: $\left| {{\mathrm{\Psi }}_ + } \right\rangle \to \left| {{\mathrm{\Psi }}_ - } \right\rangle$, $\left| {{\mathrm{\Psi }}_ - } \right\rangle \to \left| {{\mathrm{\Psi }}_0} \right\rangle$. Thus, after insertion of a quantum of spin-dependent flux, we expect to obtain a permutation of the three ground states. This expectation is independent of the presence of a quantization of q, because the spectrum has to be invariant under ϕ → ϕ + ϕ₀. The observation of the state permutation under spin-dependent flux insertion could, however, be obstructed by large avoided crossings in the evolution of the energy levels. It is important to stress that the spectrum remains gapped (above the three ground states) during the entire process of spin-dependent flux insertion. This gap is provided by the superconducting order parameter that couples states of different particle number on each edge. Without the superconductivity, charge pumping would still occur between the gapless edges, but the adiabatic process would simply accumulate charge at one edge and deplete the other, mapping the eigenstates to others with ever higher energy with each quantum of spin-dependent flux inserted.

To implement the spin-dependent flux insertion, we observe that the substitution m to m + ϕ/ϕ₀ in Eq. (1) is equivalent to substituting m by m + ϕ/ϕ₀ in μ _m and f _m for an infinitely long cylinder. In a finite cylinder, such an approach is still valid for the low-energy subspace as long as the number of orbitals is larger than the number of orbitals typically covered by the incompressible liquid. In our case, this is roughly given by the number of orbitals needed by a single Laughlin state with $\tilde N{\mathrm{/}}2$, i.e., $3\tilde N{\mathrm{/}}2 - 2$.

To simulate spin-dependent flux insertion in a system with N_Φ + 1 orbitals, m = −N_Φ/2, …, N_Φ/2, we consider a system enlarged by one orbital, m = −N_Φ/2, …, N_Φ/2 + 1 and use a linear interpolation of the functions μ _m and f _m , which allows their argument to take real values and substitute in the Hamiltonian (6)

$$\mu _m \to \mu _{m + \phi /\phi _0}{\kern 1pt} {\mathrm{and}}{\kern 1pt} f_m \to f_{m + \phi /\phi _0}.$$

(5)

When tuning ϕ, the potential experiences a kink around m = 0 which would result in a kink in the energy spectrum. We have thus replaced the absolute value around 0 by a quartic polynomial interpolation that ensures the potential and its derivative are continuous. Similarly for f, we use a linear interpolation for any orbital at the boundary between a superconducting (f = 1) and a normal (f = 0) region.

The low-energy spectrum of the resulting φ-dependent Hamiltonian is plotted in Fig. 5. Up to a small avoidance, the three ground states permute as anticipated, while the spectral gap above them stays intact in the process. The overall evolution of all energy levels with a minimum at ϕ = ϕ₀/2 is a result of the specific interpolation (8). Indeed, leaving unoccupied the orbitals near the system ends (where the confining potential is more important) results in a lower total energy. This demonstrates that the superconducting coupling has indeed joined up the two layers into a topological equivalent of the torus geometry sketched in Fig. 1c. In Supplementary Note 1, we present a full phase diagram for the spin-dependent flux insertion (similar to Fig. 3) for both this system size and a slightly smaller one.

Fig. 5

Evolution of the energy levels under spin-dependent flux insertion for Hamiltonian Eq. (3) with the same parameters as in Fig. 2 except for $L_y = 7.0{\mathrm{/}}\ell _B$ and Δ₀ = 1.2. The spin-dependent flux insertion moves particles in the background of the linear onsite potential, giving rise to the overall ϕ-dependent energy shift of the eigenstates. a Evolution of the low energy spectrum. Red are the four lowest states in the K _y  = 0 sector, black are the lowest states in each of the other K _y  = 0 sectors. Thus, not all states in the gray region are shown. b Close-up of the evolution of the three lowest states corresponding to the topological edge degrees of freedom, showing how the three states are permuted (up to small anticrossings) as ϕ is changed by 2π

6π Josephson effect

As a second piece of evidence that the heterostructure realizes the topological superconducting edges, we calculate the evolution of the energy levels that corresponds to the 6π Josephson effect. In order to do so, the relative complex phase between the mean-field superconducting order parameters on the left and the right edge, φ, is varied. The Josephson effect requires quasiparticle tunneling processes between the superconducting edges. Necessarily, the spectrum of the heterostructure displays a 2π periodicity in φ. However, in the thermodynamic limit, in which the three states are degenerate, the ground state of the system does not return to itself when φ is advanced by 2π. Rather, it evolves into another degenerate ground state and only after φ is advanced by 6π does the system return to its initial state. The reason for this behavior is that the elementary excitations of the superconducting edge are Cooper paired quasiparticles of charge 2e/3, delocalized along the cylinder perimeter, which tunnel across the bulk gap.

To observe the 6π Josephson effect numerically in our finite-size setup, the energy scale associated with the tunneling must be larger than the finite size splitting between the ground states. For the system sizes accessible to exact diagonalization calculations, this is not generically the case, even when the spread-to-gap ratio shown in Fig. 3b is large. Since the tunneling amplitude is exponentially small in the distance between the edges, we expect a favorable regime for large cylinder circumference L _y (at fixed number of non-superconducting orbitals), so that the physical distance $\propto L_y^{ - 1}$ between the edges is small. Figure 6 shows the spectral evolution as a function of φ in this regime. We observe that the three low-lying states are indeed permuted as φ advances by 2π up to a residual small avoidance of the crossings between the states.

Fig. 6

Evolution of the energy levels of the 6π Josephson effect for Hamiltonian Eq. (3) with the same parameters as in Fig. 2 and Δ₀ = 2.13, varying the phase difference φ between the left and the right superconducting edge. a Evolution of the low energy spectrum. Red are the four lowest states in the K _y  = 0 sector, black are the lowest states in each of the other K _y  = 0 sectors. Thus, not all states in the gray region are shown. b Close-up of the evolution of the three lowest states corresponding to the topological edge degrees of freedom, showing a 6π periodicity with a small avoidance of the crossings between the states

To investigate in which region of phase space this type of spectral evolution can be found, we plot the ratio r of largest avoided crossing over energy spread of the ground state manifold, i.e., the quantity

$$r: = \frac{{{\mathrm{max}}\left[ {\left( {E_2(0) - E_0(0)} \right),\left( {E_1(\pi ) - E_0(\pi )} \right)} \right]}}{{E_2(\pi ) - E_0(0)}},$$

(6)

where E₀(φ), E₁(φ), E₂(φ) are the energies of the three lowest states as a function of flux φ. When the avoidance in the evolution of the energy levels vanishes, such that E₁(0) = E₂(0) and E₀(π) = E₁(π), then r → 0 and the 6π Josephson effect becomes clearly observable. In the opposite limit where the low-energy spectrum is essentially independent of φ E _i (0) = E _i (π), i = 0, 1, 2, we have r → 1 and the tunneling matrix elements are too small to overcome the finite-size induced energy splitting between the three low-lying states. Figure 3d shows r as a function of the strength of the superconducting pairing potential and the physical distance between the edges. Indeed, we find that the 6π Josephson effect is best observed when both the spread-to-gap ratio and the distance between the edges is smallest.

Discussion

We numerically studied a heterostructure of a bilayer system of FQH Laughlin states with counter-propagating edge modes that are gapped out by a mean-field superconducting order parameter. The system in the cylinder geometry with two gapped edges realizes a nonlocal topological qutrit.

Our calculations were performed using exact diagonalization and by restricting the computation to the quasihole subspace of the Laughlin state in each layer. Despite this simplification, the system size is still limited. Nevertheless, we have been able to demonstrate four key features: (i) the edges develop a spectral gap induced by the superconducting coupling, (ii) the expected number of three nearly degenerate ground states without any charge imbalance between the two halves of the system, (iii) that charge pumping can permute the ground states, and (iv) that the system exhibits a 6π-periodic Josephson effect. For each signature, we discussed the suitable parameter regime.

While the details of the phase diagram are affected by important finite-size effects, similar features can be identified in all studied systems for Δ₀ ≳ 1.5 Extrapolating from these features, we propose a physical summary and a highlight of our quantitative results on the 6π Josephson effect in Fig. 7. We focus on a given value of the pairing parameter Δ₀ = 2.0 and we explore the different regimes as a function of the size of the normal region or equivalently the distance between the two superconducting leads. This distance L_x,n is deduced from the number of orbitals without superconducting coupling as discussed previously. As seen in Fig. 7, the value of L_x,n determines the behavior of the system among three regimes. When L_x,n ≳ 3.4$\ell_B$, there is a large gap and an approximate threefold degeneracy, leading to well defined gapped edge modes. But the tunneling required for the 6π Josephson effect is exponentially suppressed. In the intermediate regime 2.4$\ell_B$ ≲ 3.4$\ell_B$ ≲ L_x,n, the 6π Josephson effect is clear and a robust gap remains. Note that the optimal value of L_x,n is roughly twice the correlation length³⁵ $\xi \simeq 1.4\ell_B$ of the Laughlin ν = 1/3 phase which allows tunneling without destroying the underlying quantum liquid. Finally, when L_x,n is small, i.e., a distance lower than the correlation length ξ, the induced gap at the edges collapses. Thus our results validate the hypothesis of the previous effective approaches and provide an estimate of the characteristic dimensions for future experimental implementations.

Fig. 7

Phase diagram along the Δ₀ = 2.0 line. The horizontal axis is the width of the normal region $L_{x,n} = \frac{{2\pi \ell _B^2N_{{\mathrm{orb,n}}}}}{{L_y}}$. The vertical axis is either the energy gap (red), the spread over gap ratio (dashed black), or the avoidance-to-spread ratio for the Josephson effect (dashed blue). We clearly discriminate three different regimes: I for small L_x,n, the backscatteting dominated regime with the breakdown of the gapped edge modes; II the intermediate region corresponding to the 6π Josephson regime; III at large L_x,n, the separated gapped edge modes

Our work provides the first quantitative study of fractional edge modes coupled to superconducting leads in a fully microscopic model. Future works, most probably relying on the density matrix renormalization group calculations, should be able to rely on our setup and results to provide new insights. In particular, it could overcome the size limitation and address the potential new emerging phases when substituting the Laughlin state with any richer topological order.

Data availability

The data have been generated using the software “DiagHam” (under the GPL license).