Overlap of parafermionic zero modes at a finite distance

Parafermion bound states (PBSs) are generalizations of Majorana bound states (MBSs) and have been predicted to exist as zero-energy eigenstates in proximitized fractional quantum Hall edge states. Similarly to MBSs, a finite distance between the PBS can split the ground state degeneracy. However, parafermionic modes have a richer exchange statistics than MBSs, so several interaction terms are allowed by the underlying $\mathbb{Z}_{2n}$ symmetry, rendering the effective Hamiltonian governing a pair of PBSs at a finite distance nontrivial. Here, we use a combination of analytical techniques (semiclassical instanton approximation) and numerical techniques (quantum Monte Carlo simulations) to determine the effective coupling Hamiltonian. For this purpose, we go beyond the dilute one-instanton gas approximation and show how finite-size effects can give rise to higher-order parafermion interactions. We find excellent agreement between the analytical results and Monte Carlo simulations. We estimate that these finite-size corrections should be observable in some of the recently proposed experiments to observe PBSs in strongly correlated systems.


I. INTRODUCTION
Recent developments in the field of fractional quantum Hall (FQH) systems have allowed the creation of devices containing regions which can be subject to induced superconductivity or interedge scattering [1][2][3]. Such systems are theoretically expected to be suitable for hosting bound states with non-Abelian statistics at the interfaces between superconducting and scattering regions [4][5][6][7][8][9]. Such states, referred to as Z n parafermion bound states (PBSs), are generalizations of Majorana bounds states (MBSs), the latter corresponding to PBSs with n = 2 [10].
Being one-dimensional and subject to interactions, systems hosting parafermions are usually modeled using bosonization methods [28]. Parafermionic modes emerge at the interfaces between regions where FQH states are coupled by backscattering and superconducting pairing, thus generating nontrivial band gaps. Both effects give rise to cosine terms in the bosonization language, so parafermion systems can be modeled as inhomogeneous sine-Gordon Hamiltonians. While a homogeneous sine-Gordon model is among the rare examples of an exactly solvable interacting model, no such exact solution is known for the inhomogeneous model. Therefore, the theoretical modeling typically rests on approximations.
One such method is the instanton gas approximation, which is based on a semiclassical treatment of the sine-Gordon Hamiltonian [29][30][31][32]. Single instantons, time-like kinks of the phase field, are the solutions of the corresponding classical Euler-Lagrange equation and can be used to compute the ground-state energy [33,34]. Multiple instantons can easily be accounted for in the dilute one-instanton gas approximation, where different instantons are assumed to be far apart, so that their interaction can be neglected. This method can be extended towards multiple interacting instantons, also called molecular instantons [35][36][37]. Used in the context of resurgence theory [38,39], multi-instantons can offer a path to go beyond the first-order approximation. As we will show, such interacting multi-instanton configurations can become important for PBSs located at a finite distance from each other.
In this work, we go beyond the dilute one-instanton gas approximation [32] and include bi-instantons [38], i.e., configurations including pairs of correlated instantons, in a system that hosts Z 2n parafermion bound states. This second-order process adds an extra ground-state splitting term that exponentially decays twice as fast as the one-instanton term [29] and contains a distinct oscillatory length dependence in the case of a nonzero chemical potential. These approximate results are then compared to a path-integral quantum Monte Carlo simulation and we find good quantitative agreement.
We proceed to analyze this ground-state splitting correction in terms of higher-order parafermion interactions, i.e., interaction terms containing more than two parafermion operators. We show that the bi-instanton is associated with a four-parafermion interaction, a result which can be generalized to arbitrary higher orders and associated with 2k-parafermion interactions. We briefly discuss the prospects of an experimental observation of bi-instanton corrections.
The paper is organized as follows: in Sec. II we introduce the sine-Gordon model and the basic assumptions utilized to simplify the action. Moreover, we introduce the effective low-energy model governing the parafermionic modes. In Sec. III, we obtain the finitesize energy splitting of the ground state, first by explain-ing the dilute one-instanton gas approximation, and then by extending it to include bi-instantons. We continue by comparing the analytical results with Monte Carlo simulations in Sec. IV. Finally, we present our parameter estimates and conclusions in Sec. V.

II. MODEL
The model system for studying parafermion bound states consists of a pair of FQH edges containing two counter-propagating modes. To generate parafermion bound states, it is necessary to engineer two topologically distinct spectral gaps in different regions such that a subgap parafermion bound state appears at the interfaces between these gapped regions [4,6]. Many proposals consider an interface between ferromagnetic (FM) and superconducting (SC) regions [7,9,32]. Denoting by ψ L,R (x) the annihilation operators for left-and right-moving electrons in the FQH state, the ferromagnet induces a backscattering gap corresponding to a term ∆ FM ψ † L ψ R +h.c.. On the other hand, the proximity effect from a nearby superconductor creates and annihilates Cooper pairs and gives rise to a term ∆ SC ψ † L ψ † R + h.c.. Hence, a pair of FQH edge states with such an FM-SC-FM junction (see Fig. 1) should host parafermion bound states at the interfaces.

A. Sine-Gordon Hamiltonian
We consider a FQH state with filling factor 1/n and group velocity v. We first bosonize the left-and rightmoving electrons in terms of chiral bosonic fields ϕ L,R (x) such that where ξ is the correlation length, which is the inverse of the high-energy cutoff, ξ = v/E cutoff (using = 1). We continue by defining the fields φ = (ϕ R + ϕ L )/2 and θ = (ϕ R − ϕ L )/2, which satisfy the commutation , where Θ is the Heaviside function. Hence, φ(x) and (n/π)∂ x θ(x) are canonically conjugate variables. The backscattering and pairing terms can then be simplified to ∆ SC [sin(2nφ) + 1] and ∆ FM [sin(2nθ) + 1], respectively.
The Euclidean action for the whole system thus becomes an inhomogeneous sine-Gordon model, where ∆ FM (x) and ∆ SC (x) vanish, respectively, outside the FM and SC regions and are constant inside * *

FM SC FM
FIG. 1. The system comprises a pair of FQH edges with two counter-propagating edge modes at filling factor 1/n. The FQH edge is subject to induced superconductivity and ferromagnetic coupling leading to an FM-SC-FM junction. The parafermion zero modes appear at the interfaces between the FM and SC regions, which is illustrated by the asterisks.
those regions. We assume the chemical potential µ to be constant along the system. We focus on the limit where L FM is the length of the FM region. In this limit, the field θ(x) is pinned to a minimum of sin(2nθ) inside the FM region, while the field φ is allowed to fluctuate. By completing the square and considering θ(x) as constant inside the FM region, the effective action of the system becomes that of a simple sine-Gordon model with a Berry-phase term due to the chemical potential, Here and in the following, we use ∆ ≡ ∆ SC and L ≡ L SC to denote the pairing strength and the length of the superconducting region, respectively. Due to the sine potential, the field φ is allowed to fluctuate about the minima of the sin(2nφ) potential. However, static solutions do not capture the whole picture and instanton solutions that interpolate between different minima are necessary for the ground-state energy. In the next section, we review the dilute one-instanton gas approximation and extend it to include bi-instanton processes.

B. Effective parafermion Hamiltonian
At low energies, the predictions resulting from the sine-Gordon model can be translated into an effective parafermion Hamiltonian, which is valid at energies below the band gaps, |E| ∆, ∆ FM . To describe the effective parafermion interaction, we consider two parafermion modes described by Z 2n parafermion operators α L,R . These operators satisfy the parafermionic commutation relations α 2n L,R = 1, α † L,R = α 2n−1 L,R and α L α R = e iπ/n α R α L . A generic interaction Hamiltonian that preserves the Z 2n charge can be written as The phase diagrams of parafermion Hamiltonians, with similar building blocks, have been studied using different methods [10,[40][41][42][43]. This Hamiltonian can be written in the eigenstate basis of the first-order interaction term, α † L α R |q = −e iπ(q−1/2)/n |q , where the eigenstates |q are labeled by an integer q ∈ {0, 1, 2, ..., 2n − 1} which corresponds to the total charge inside the superconducting region in units of the fractional charge e/n [32]. As a consequence, powers of the first-order interaction term satisfy such that the energy eigenvalues E(q) of the Hamiltonian H Z2n can be expressed as a sum, In the following, we will use a semiclassical instanton calculation as well as a Monte Carlo simulation to determine the effective coupling strengths t 1,2 and the phases λ 1,2 .

III. INSTANTON CALCULATION
A. Review of the dilute one-instanton gas In this section, we briefly review the dilute oneinstanton gas approximation [29][30][31][32] to find the energy splitting between different ground states. In order to obtain the energies, it is necessary to compute the transition rates between different configurations of the field φ(x, τ ). In imaginary time, such transition amplitudes correspond to the matrix elements, between two stationary states |j ± , in which φ(x, τ ) is pinned at a minimum of the sine potential and which are eigenstates, in the limit ∆L/v 1, of the Hamiltonian corresponding to the action given by Eq. (3). Moreover, T ξ/v is a large time. This transition amplitude can be conveniently calculated with the action S[φ].
We focus on classical solutions of the action S[φ] with constant spatial profile, ∂ x φ = 0, [32] because solutions with a nonzero ∂ x φ have a larger action and only contribute subleading corrections to the transition amplitude for a given number of instantons. Using this simplification, the equation of motion corresponding to the action (3) becomes and allows us to define the states |j which correspond to the stationary solutions φ j = −π/4n + jπ/n (j ∈ {0, 1, .., 2n − 1}) which minimize the classical action. Solitons correspond to the non-stationary classical solutions and can be found by straightforward integration of Eq. (8). They take the form of a classical field that interpolates between two stationary solutions and is centered at an imaginary time τ 0 (for = ±), Here, φ + cl refers to an instanton with final state |j + 1 , while φ − cl denotes an anti-instanton with final state |j − 1 . Moreover, we defined ω = 2 ∆v/ξ. The action of the soliton field has two different contributions, i.e., a kinetic term (S 0 ) and a Berry-phase term (S Bp ). While the former is responsible for the system's overall energy scale, the latter will induce an oscillatory behavior in the energy splitting, where we integrated over the position and used Eq. (8) to simplify the action. Calculating the integral for the soliton, one finds As the classical instanton (anti-instanton) interpolates between |j and |j + 1 (|j − 1 ), this allows us to identify the transition rate between these two stationary states. By calculating the path integral with fluctuations η around the classical solution, the quantum amplitude of a transition starting in state |j at imaginary time τ = 0 and ending in state |j + 1 (|j − 1 ) at τ = T is given by where, in the first line, we expanded the action around classical solitons between states |j and |j ± 1 , such that we have a Gaussian integral over all fluctuations η(x, t) which vanish at the boundaries. This remaining path integral contains a translation-invariant zero mode, which can be integrated out and creates a factor T πvS 0 /n [29,30]. The remaining integral without zero modes (denoted by the prime in the path-integral measure) is Gaussian, andF = −∂ 2 τ −v 2 ∂ 2 x +ω 2 sin(2nφ cl ) is the differential operator associated with the fluctuations. By computing the Gaussian integral in the second line, we obtain the determinant without zero modes, det [F ], multiplied by N , the normalization constant from the measure D [η] [29]. Following Refs. [29,32], the determinant can be calculated by multiplying and dividing by det x + ω 2 is the differential operator of a harmonic oscillator, and computing the ratio of determinants using the zeta regularization method for a Neumann boundary condition [44], As a result, the transition rate becomes where K = ω/(π √ n). We identify in the results for where we summed over all combinations of n I and nĪ with n t = n I − nĪ. In the second line, we used the summation form of the Dirac delta function, expanding n t and distributing the exponential e ±iπq together with the instanton amplitude.
To obtain an expression for the energy splitting, we expand j + | e −HT |j − in terms of the eigenstates of a general two-parafermion Hamiltonian and we compare both expressions. The eigenstates |q of the effective lowenergy parafermion Hamiltonian (4) are similar to Bloch waves and can be written as linear combinations in the |j basis [44]. By noticing that |j = |j + 2n due to the periodicity of the sine function, we can write By considering the completeness of the basis |q , we can expand |j ± to obtain Neglecting the constant energy ω/2, we can set Eq. (17) and Eq. (15) equal, to arrive at We note that E(q) decays exponentially as ∝ e −Lω , whereas the chemical potential gives rise to a typical oscillation [32]. A final equality between Eq. (18) and Eq. (6) relates the parafermion coupling parameters to the microscopic model constants, In the dilute one-instanton gas approximation, all higherorder coupling processes are absent, i.e., t k = 0 for k > 1.

B. Beyond the dilute one-instanton gas
Moving beyond the dilute one-instanton gas approximation, we now consider a system composed of instantons and bi-instantons, the latter corresponding to a correlated two-instanton event [37][38][39]. In such a biinstanton, the change of φ(τ ) during the transition is still fast compared to the distance between two instantons, but in contrast to the instanton gas limit, this distance is not infinite [38].
In contrast to a single instanton, the bi-instanton is not an exact solution of the classical Euler-Lagrange equation and has a different winding number. Nonetheless, it is a solution up to a correction which is exponentially small in the distance between the instantons, so it can have a significant contribution to the quantum mechanical path integral. Bi-instantons will have an energy scale of the order of e −2S0 [29], but interactions between the two instantons produce corrections to this energy [36]. We consider the action for a configuration of two instantons at positions ±τ 0 given by (for 1 , 2 = ±) n arctan e ω(τ +τ0) + 2 2 n arctan e ω(τ −τ0) .
For 1 = 2 such a configuration describes a pair of instantons ( 1,2 = 1) or anti-instantons ( 1,2 = −1), whereas for 1 = − 2 , the field describes an instanton anti-instanton pair. The action of such a bi-instanton is twice the classical action of a single instanton plus a positive energy due to the repulsive interaction between them (see the Appendix. A for a derivation), In the following, we focus on the bi-instantons with 1 = 2 ≡ because these contribute to the energy splitting. In contrast, an instanton-anti-instanton pair ( 1 = − 2 ) does not allow additional transitions between different |j , so it just adds a constant value to the energies.
By direct application of the concepts previously introduced for the dilute one-instanton gas, we compute the transition rate due to a single bi-instanton, where we followed the same procedure as in Eq. (12), with the translational invariance of the bi-instanton center of mass now being responsible for the zero mode. Here, z is the distance between the centers of the two instantons and the integral over z comes from the quasi zero modes of D[η] [36].M = −∂ 2 τ − v 2 ∂ 2 x + (4∆v/ξ) sin(2nφ ) is the differential operator associated with the bi-instanton. The eigenvalues ofM are two fold degenerate so det [M ] = (det [F ]) 2 [37].
The correction to the bi-instanton amplitude due to the instanton interactions, [II] = [I] 2 dze −4S0e −ωz , can be simplified by using the semiclassical approximation S 0 1 and by introducing a regularization parameter c [39] in the instanton-instanton interaction given by Eq. (21), where γ E ≈ 0.5772 is the Euler number and we need to subtract the divergent term O(1/c) which corresponds to non interacting instantons. We proceed by generalizing Eq. (15) to a gas consisting of both instantons and biinstantons, which we call a dilute two-instanton gas approximation. Between two normalized stationary states |j + and |j − , we consider all possible combinations of (anti-)instantons and (anti-)bi-instantons where n t = n I − nĪ + 2n II − 2nĪĪ is the net number of instantons. By comparing Eq. (24) with Eq. (17), it is possible to obtain the energy splitting including the subleading correction due to bi-instantons, In Fig. 2, we show the energy splitting for a Z 4 parafermion (n = 2), comparing the result for the dilute oneinstanton gas (18) (dashed line) with the corrected result including bi-instantons (25) (solid line) for a value √ 2e −S0 = π and [γ E + ln(4S 0 )]/2 = 0.1. While energy crossings between two consecutive q mod 4 occur at the same energy, the crossings between q and q + 2 mod 4 are shifted. By setting Eq. (25) and Eq. (6) equal, we recover the leading order result (18) and, in addition, the second-order coupling amplitude and phase,

IV. MONTE CARLO
We start with a discretization of the action (3), using a lattice constant a x and time step a τ . In the following, we set v = 1, which leads to where τ i = ia τ and x j = ja x . Moreover, N x and N τ are chosen such that N x a x = L and N τ a τ = T . We choose a rectangular grid, N τ = 100, and set a x = a τ = a. To simulate instanton configurations which connect states |j − and |j + , we assume a twisted boundary condition in the τ direction, such that φ τ Nτ +1 ,xi = φ τ1,xi + πδQ/n [45], where δQ = j + − j − is an integer. This condition is enforced in the simulation via an energy penalty nax 2πaτ (φ τ1,xj −φ τ Nτ ,xj +πδQ/n) 2 , and causes the simulated configuration to have a net number δQ of instantons.
The spatial boundary conditions are open, allowing the field to fluctuate freely in the x direction. However, fields which fluctuate in the x direction lead only to subleading contributions compared to configurations which are spatially constant.
The Monte Carlo algorithm is initialized with a vanishing configuration φ τi,xj = 0 for all i, j. In each Monte Carlo step, we propose a field update φ → φ , which is accepted with probability P = min(| exp(−α∆S)|, 1) that is determined by the change of the action ∆S = . We introduced an additional parameter α ≥ 1 which will be explained shortly. Hence, if the updated action S[φ ] is smaller than the original action S[φ], then ∆S < 0 and the update will be accepted (P = 1). On the other hand, if ∆S > 0, we draw a uniformly distributed random number X ∈ [0, 1] and accept the proposed update if X < P . Otherwise, the update is rejected and the original field configuration is kept.
A "bead update" which randomly displaces the field along all imaginary-time coordinates has a small acceptance rate, which results in an inefficient simulation. To achieve an improved convergence towards fluctuations around the classical saddle-point solutions of the path integral, we implemented a local field update scheme at the level of individual nodes, i.e., φ τi,xj = φ τi,xj + δφ, at a random position (i, j), where δφ ∈ [−δ, δ] is a uni-form random displacement and δ a dynamical interval width. After 100 Monte Carlo steps, we update δ depending on the acceptance rate. If the acceptance rate is below 0.6, we reduce δ to 0.8δ. Typically, this means δ converges to 0.1. Since the action difference after a single local update is independent of the imaginary time T (interpreted as an inverse temperature), an artificial parameter α with αa τ ∝ T is introduced to obtain a local importance sampling scheme corresponding to an average inverse temperature. We found optimal performance and fast convergence to fluctuations around the classical saddle-point solutions for the choice α = N τ .
Every 200 steps, we randomly attempt center-of-mass moves of the field along τ : For a fixed x j , we propose a uniform shift of the field for all τ i , i.e., φ xj → φ xj + δ, where φ xj = (φ τ1,xj , φ τ2,xj , . . .). The move is accepted according to the same criterion as for the individual node updates. By swapping τ and x, a center-of-mass shift is attempted along the spatial direction. This routine allows slight changes in the center of mass of instantons and improves convergence.
After the system equilibrates (we typically consider 10000 local updates, although the system usually converges in less than 5000), we store the final configuration. The simulation is restarted by creating either an instanton or an anti-instanton, with the same likelihood, at a random imaginary time τ after which we let the system equilibrate again. The creation or annihilation of instantons followed by equilibration is repeated 60 times before the next configuration is stored. This procedure allows us to sample configurations in the neighborhood of different saddle points which otherwise cannot be reached in reasonable computation time. For each δQ, we collect N c = 3000 configurations based on which we compute the observables explained in the following. To simulate the action, given by Eq. (27), we use n = 2, ω = √ 2, a = 0.6, and N τ = 100. For each length N x and chemical potential µ, we simulated configurations with 0 ≤ δQ ≤ 2.
The twisted boundary conditions effectively decouple different instanton sectors and allow us to simulate quantities such as Eqs. (12) and (22). Since all δQ sectors contribute to physical observables, we define the following average: where φ δQ are configurations containing δQ instantons or anti-instantons (negative δQ). We denote the transition rates of δQ instantons as The Monte Carlo result depends on a normalization constant N MC which is independent of δQ. To check the results we plot G 2 MC and G 1 MC [see Fig. 3(a)]. Fitting the expressions for G for different choices of parameters (ω, T and n = 2). Note that even for small lengths, G 2 MC agrees with the analytical result even though it is near the region where the theory is no longer valid. Indeed we later show that the energy splitting deviates from the expected value in this region.
From the transition rates, we follow the derivation of Eq. (18) to compute the energy splitting. In particular, we note that Eq. (17) can be written as a Fourier transform of transition rates, i.e. δQ|e −HT |0 = 1 2n q e −iπqδQ/n q|e −E(q)T |q , (30) from which follows The energy splitting is plotted in Fig. 3(b) for L = 12 and shows a good agreement between simulation and theory. We isolate the subleading contributions in Fig. 3(c), which demonstrates that the corrections are two orders of magnitude lower than E(q) for the chosen parameters, but still in very good agreement with the theory. The bi-instanton corrections are more visible when we consider E(q) as a function of length for a fixed chemical potential; see Fig. 4 for µ = 0.5. According to the theory, the level crossing between q = 0 and q = 2 at L ≈ 6 is shifted, which we confirm by the simulations. For small lengths, the semiclassical approximation does not fully capture the behavior of the system  4. |E(q)| as a function of length for µ = 0.5 with the same color scheme as Fig. 3. Note that for L < 6, it is possible to clearly observe an asymmetry between the q = 0 and 2 sectors.
and we see differences between the analytical results and Monte Carlo simulations. Even so, it is still possible to find regions in which the semiclassical regime describes bi-instantons and the effects are not negligible.

V. CONCLUSION
In this work, we have analyzed the relevance of subleading finite-size effects to the effective description of parafermion bound states existing in a pair of FQH channels with FM-SC-FM interfaces.
We derived a dilute bi-instanton gas approximation and demonstrated that the outcomes are compatible with high-order parafermion interactions in the effective model. We found a relation between the coupling parameters and the microscopic constants, which reveals that many-body parafermion interaction terms can be observed in small systems and manifest themselves in characteristic oscillations of the ground-state energy level splitting as a function of chemical potential or system size.
We estimate that the subleading correction due to fourparafermion interaction can be of the order of 10% of the leading contribution. Indeed, using typical experimental values for the FQH edges [1,19,46], we can estimate v = 10 5 m/s, E cutoff ∼ 2 meV the bulk gap of FQH, L ∼ µm and ∆ ∼ 1 meV, from which we find that the correction is of the order of 0.1 ∼ 10%, as shown in Fig. 5.
The bi-instanton approximation is also compatible with Monte Carlo simulations. We highlight that the agreement between analytical and numerical methods indicates the existence of high-order interactions that cannot be explained by previous results. Far from being a mere hindrance, these subleading terms, which do not exist in the case of Majorana bound states, can give rise to novel phases and mechanisms that lie beyond the usual description of parafermion chains.
The first and third terms are part of the classical action of the instantons centered at τ 0 and −τ 0 . This can be seen by rewriting 1 2 V ( π 4n + f + )f 2 − ≈ V ( π 4n + f − ). The linear term in f − of S + can be integrated by parts using the equation of motion, We can do the same procedure for S − , such the classical action of two instantons is S(φ + 2 ) = 2S 0 + 4 2 S 0 e −2ωτ0 . (A6)