Parafermion braiding in fractional quantum Hall edge states with finite chemical potential

TL;DR

This paper investigates the braiding of $ ext{Z}_{2n}$ parafermions in fractional quantum Hall systems, analyzing the effects of finite chemical potential and proposing an improved protocol to ensure adiabatic braiding.

Contribution

It provides a detailed numerical study of parafermion braiding with finite chemical potential and introduces a new protocol to avoid spectral degeneracies during braiding.

Findings

Finite chemical potential causes spectral crossings that challenge adiabatic braiding.

Numerical simulations demonstrate the impact of chemical potential on energy splitting.

An improved braiding protocol is proposed to maintain adiabaticity.

Abstract

Parafermions are non-Abelian anyons which generalize Majorana fermions and hold great promise for topological quantum computation. We study the braiding of $\mathbb{Z}_{2n}$ parafermions which have been predicted to emerge as bound states in fractional quantum Hall systems at filling factor $\nu = 1/n$ ($n$ odd). Using a combination of bosonization and refermionization, we calculate the energy splitting as a function of distance and chemical potential for a pair of parafermions separated by a gapped region. Braiding of parafermions in quantum Hall edge states can be implemented by repeated fusion and nucleation of parafermion pairs. We simulate the conventional braiding protocol of parafermions numerically, taking into account the finite separation and finite chemical potential. We show that a nonzero chemical potential poses challenges for the adiabaticity of the braiding process because it leads to accidental crossings in the spectrum. To remedy this, we propose an improved braiding protocol which avoids those degeneracies.

Equations

$$\delta E\propto e^{-L/\xi }\cos (\mu L/v_F).$$

$${\chi }_k{\chi }_l$$

$${\chi }_k^{†}$$

$$H(\tau )=t_{23}(\tau ){\chi }_2{\chi }_3^{†}+t_{12}(\tau ){\chi }_1{\chi }_2^{†}+t_{24}(\tau ){\chi }_2{\chi }_4^{†}+\text{h.c.}$$

$$H=$$

$$-\int dx\frac{\Delta (x)}{n\pi a}\cos (2n\theta (x))-\int dx\frac{g(x)}{n\pi a}\cos (2n\varphi (x)),$$

$${\psi }_{\pm }(x)=\frac{e^{-in\left[\pm \varphi (x)-\theta (x)\right]}}{\sqrt{2\pi na}},$$

$$[\theta (x),\varphi (y)]=\frac{i\pi }{n}\Theta (y-x),$$

$$H$$

$$-\frac{\Delta }{n\pi a}\int dx\cos (2n\theta (x)).$$

$$\tilde{\theta }(x)$$

$$\tilde{\varphi }(x)$$

$$H$$

$$-\frac{\tilde{\Delta }}{\pi a}\int dx\cos (2\tilde{\theta }(x)).$$

$${\tilde{\Psi }}_{\pm }(x)=e^{\mp i\tilde{\mu }x/v_F}{\tilde{\psi }}_{\pm }(x),{\tilde{\psi }}_{\pm }(x)=\frac{1}{\sqrt{2\pi a}}{e}^{-i\left[\pm \varphi (x)-\tilde{\theta }(x)\right]},$$

$$H_{hl}$$

$$H_{\mu }$$

$$H_{\Delta }$$

$$\left(\begin{array}{c}{c}_{k,+}\\ c_{-k,-}^{†}\end{array}\right)=\frac{1}{\sqrt{2}}\left(\begin{array}{cc}\sqrt{1+{ϵ}_k/E(k)}& \sqrt{1-{ϵ}_k/E(k)}\\ -\sqrt{1-{ϵ}_k/E(k)}& \sqrt{1+{ϵ}_k/E(k)}\end{array}\right)\left(\begin{array}{c}{\gamma }_{k,+}\\ {\gamma }_{-k,-}^{†}\end{array}\right).$$

$$H=$$

$$E(k)=$$

$${\tilde{G}}_R(x_2,x_1,t)=$$

$${\tilde{G}}_R(x_2,x_1,\omega )=$$

$${\gamma }_{k,\sigma }(t)=e^{-iE(\sigma k)t}{\gamma }_{k,\sigma }(0).$$

$${\tilde{G}}_R(x_2,x_1,\omega )=$$

$$\frac{e^{i\tilde{\mu }L/v_F}}{2\pi v_F}\int dk$$

$$\displaystyle\tilde{G}_{R}(x_{2},x_{1},\omega)=-\frac{e^{-L/\xi_{\Delta}(\omega)}e^{i\tilde{\mu}L/v_{F}}}{2v_{F}}\Bigg{[}i+\frac{\omega\xi_{\Delta}(\omega)}{v_{F}}\Bigg{]},$$

$$-\frac{1}{\pi }\text{Im}{\tilde{G}}_R(x_2,x_1,0)\propto \exp \left(-\frac{\Delta L}{n{v}_F}\right)\cos \left(\frac{\mu L}{n{v}_F}\right).$$

$$H$$

$$-\frac{g}{K\pi a}\int dx\cos (2\tilde{\varphi }(x)).$$

$${\tilde{G}}_R(x_2,x_1,\omega )=\frac{1}{2\pi v_F}\int dk\frac{(\omega -\tilde{\mu }+k)e^{ikL/v_F}}{k^2+{\tilde{g}}^2-(\omega -\tilde{\mu }{)}^2-2i\omega \delta },$$

$$-\frac{1}{\pi }\text{Im}{\tilde{G}}_R^{(0)}(x_2,x_1,0)\propto \exp \left(-\frac{\sqrt{g^2-{\mu }^2}}{v_{F}n}L\right).$$

$$H$$

$$-\int dx\frac{\Delta (x)}{\pi a}\cos \left(2\theta (x)\right)-\int dx\frac{g(x)}{\pi a}\cos \left(4\varphi (x)\right)$$

$$-\frac{\mu }{\pi }\int dx{\partial }_x\varphi (x).$$

$$-\frac{1}{\pi }\text{Im}{G}_R^{(0)}(x_2,x_1,0)\propto e^{-\frac{\Delta L}{v_F}}\cos \left(\frac{\mu L}{v_F}\right).$$

$$-\frac{1}{\pi }\text{Im}{G}_R^{(0)}(x_2,x_1,0)\propto e^{-\sqrt{g-(\mu /2{)}^2}L/v_F},$$

Full text

I Introduction

Identical particles in two dimensions offer a more varied exchange statistics than merely the well-known bosonic and fermionic ones. The wave function of so-called anyonic particles can change by an arbitrary phase under the exchange of the particles. Moreover, if the particles span a degenerate ground state, they may even realize non-Abelian exchange statistics, where an exchange (“braiding”) implements a unitary transformation within the degenerate ground state subspace. The ground state can be used to encode quantum information, which can then be manipulated by braiding and can thus be used to implement logic gates for topologically protected quantum computation Nayak et al. (2008); Aasen et al. (2016); Alicea et al. (2011); Kitaev (2003).

The most well-known excitation believed to offer non-Abelian braiding statistics is the Majorana bound state (MBS) which has received a lot of attention in the last decade Kitaev (2001); Wilczek (2009); van Heck et al. (2012); Beenakker (2013); Franz (2013); Alicea (2012); Amorim et al. (2015); Das Sarma et al. (2015); Karzig et al. (2015, 2016, 2019). Strong experimental evidence for its existence has accumulated on several platforms Willett (2013); Nadj-Perge et al. (2014); Mourik et al. (2012); Albrecht et al. (2016); Deacon et al. (2017). MBSs are localized, but if kept at a finite distance, the exponential tails of their wavefunctions will overlap and the resulting hybridization lifts the exact ground state degeneracy Das Sarma et al. (2012); Churchill et al. (2013). In superconducting nanowires, the characteristic dependence of the energy splitting on the chemical potential and the distance has become an important tool for identifying MBS Das Sarma et al. (2012); Churchill et al. (2013); Maier et al. (2014); Zyuzin et al. (2013); Cheng et al. (2009),

[TABLE]

Here, $\xi$ is the superconducting coherence length, $L$ is the distance separating the two MBSs, $\mu$ is the chemical potential and $v_{F}$ the Fermi velocity. For the rest of this paper we will set $\hbar=1$.

Despite offering non-Abelian braiding statistics, MBSs are unable to implement all logic gates necessary for universal quantum computation in a topologically protected way Nayak et al. (2008). Therefore, non-Abelian anyons with a richer exchange statistics are needed. Recently, much attention has focused on realizing so-called $\mathbb{Z}_{p}$ parafermions in solid-state devices Lindner et al. (2012); Clarke et al. (2013); Klinovaja and Loss (2015); Alicea and Fendley (2016); Calzona et al. (2018); Chew et al. (2018); Mazza et al. (2018a); Fleckenstein et al. (2019); Mazza et al. (2018b); Snizhko et al. (2018a, b); Rossini et al. (2019); Laubscher et al. (2019). Although they do not allow for universal quantum computation, they permit to implement a richer set of quantum gates. Moreover, a network of such parafermions can be used as a building block for Fibonacci anyons, which would then indeed allow universal topological quantum computation Mong et al. (2014). An ordered set of parafermionic operators obeys the following operator algebra,

[TABLE]

In particular, MBSs can be regarded as $\mathbb{Z}_{2}$ parafermions. However, the richer exchange statistics for $p\geq 3$ makes $\mathbb{Z}_{p}$ parafermions more attractive for manipulating quantum information than MBSs Hutter and Loss (2016).

Parafermions were proposed as localized zero-modes in fractional topological insulators (FTIs). The latter can be thought of as consisting of counterpropagating quantum Hall edge states with opposite spins and fractional filling factors $\nu=1/n$ ($n$ odd) Cheng (2012); Lindner et al. (2012); Chen and Burnell (2016); Clarke et al. (2013); Snizhko et al. (2018a, b); Burrello et al. (2013a); Milsted et al. (2014). Analogously to the situation in helical edge states of 2D topological insulators (2DTIs) Fu and Kane (2009), a spectral gap in the pair of edge states can either be produced by allowing backscattering between the counterpropagating states or via the proximity effect of a superconductor (SC). These two mechanisms give rise to topologically different gaps, and the resulting interface zero-modes have $\mathbb{Z}_{2n}$ parafermions exchange statistics.

Further research has proposed their realization in different systems including 2DTIs Orth et al. (2015); Vinkler-Aviv et al. (2017); Zhang and Kane (2014); Klinovaja et al. (2014); Fleckenstein et al. (2019), fractional topological SCs Vaezi (2013); Laubscher et al. (2019), strongly interacting nanowires with spin-orbit coupling proximitized by a SC Klinovaja and Loss (2014a); Pedder et al. (2017) or in bundles of strongly interacting nanowires Klinovaja and Loss (2014b). Recently, advances in experimental platforms based on the fractional quantum Hall effect, with a potential to host parafermions, have been achieved Wu et al. (2018).

When parafermions are generated in one-dimensional systems, it is difficult to braid them by adiabatically exchanging their position without bringing them together. In nanowire-based setups this problem might be circumvented by using T-junctions Alicea et al. (2011), but for setups based on quantum Hall edge states, even this possibility does not exist. Proposals aimed at parafermion braiding are therefore based on repeated nucleation and fusion of parafermionic pairs Lindner et al. (2012); Karzig et al. (2015); Bonderson et al. (2008); Bonderson (2013); Knapp et al. (2016). One example of such a braid is shown in Fig. 1, where the two parafermionic states $\chi_{1,4}$ are effectively braided using an auxiliary pair $\chi_{2,3}$ at an intermediate step. Mathematically, the braid shown in Fig. 1 can be modeled with a simple time-dependent Hamiltonian,

[TABLE]

For $\mathbb{Z}_{n\geq 3}$ parafermions, the coupling coefficients $t_{ij}(\tau)$ can be complex and depend on the system parameters, like chemical potential, pairing gap, etc. Braiding can then be understood as a closed loop in the parameter space of the Hamiltonian. Obviously, using such a braiding protocol makes it unavoidable to bring the parafermions within a finite distance of each other, and their interactions can then be understood as due to the overlap of their wavefunctions. Alternatively, parafermions can also be coupled via Coulomb interactions rather than wavefunction overlap Snizhko et al. (2018a, b); Burrello et al. (2013a); Milsted et al. (2014).

Whenever parafermions are braided by fusion and nucleation, the notion of adiabaticity differs drastically from the case when braiding is performed by exchange of anyons in real space. In the latter case, two anyons can be physically exchanged while being kept far apart, e.g., two Majorana modes in a T-junction Alicea et al. (2011); Das Sarma et al. (2012). Adiabaticity is then related to the superconducting gap of the system which separates the ground state manifold from the quasi-particle excitations. For braiding based on nucleation and fusion, one needs to merge the parafermionic zero-modes with auxiliary parafermions, which inevitably induces finite-energy states below the superconducting gap. Hence, in this case, braiding should be performed slowly compared to these smaller gaps, rather than the larger superconducting gap of the system.

In addition to these constraints, one should note also that in both cases, the zero-modes to be braided are not exactly at zero energy because their finite separation slightly lifts the ground state degeneracy. One needs to perform the braiding fast compared to this exponentially small energy splitting, so that the ground state still appears degenerate. In principle, this sets an upper bound for the braiding time, but the latter becomes large for long systems.

Nevertheless, such a braid is still topologically protected in a weaker sense, namely in that the outcome is insensitive to small deformations of the path performed in the space of control parameters. According to Ref. Lindner et al. (2012), this path invariance holds if specific complex phases of the tunnel couplings $t_{ij}$, which would give rise to accidental degeneracies, are avoided. However, Ref. Chen and Burnell (2016) has shown using a semi-classical approach that these phases are difficult to avoid when two parafermions are separated by a SC with a finite chemical potential or a backscattering region with a magnetic field. The reasons is that both cases give rise to a distance-dependent complex phase in the coupling amplitudes and to an energy splitting similar to Eq. (1). In this work, we will assess the impact of these complex phases on braiding.

The article is organized as follows. In Sec. II, we present a short overview of our main results. In Sec. III, we use models for parafermions in fractional quantum Hall edge states or topological insulator edge states to derive the energy splitting as a function of distance and chemical potential. In Sec. IV, we simulate and compare two braiding protocols, both involving fusion and nucleation of parafermionic pairs, and compare their overlap in the presence of a nonzero chemical potential. In Sec. V, we present our conclusions.

II Summary and main results

The main subject of our study are parafermions realized in a pair of quantum Hall edge states with opposite $g$-factors at filling factor $\nu=1/n$ ($n$ odd). Such a system hosts $\mathbb{Z}_{2n}$ parafermions at the interfaces between regions gapped by either a SC or by backscattering Lindner et al. (2012); Clarke et al. (2013). We first revisit the problem of calculating the overlap between two parafermions in such a system if they are separated by a finite distance. Using a combination of bosonization and refermionization, we find the expected exponential suppression of the ground state splitting with the distance, but we also show that an oscillating term occurs if the chemical potential is nonzero. This is in accordance with a recent result based on a semiclassical analysis Chen and Burnell (2016). Moreover, we show that no such oscillating term is present in the energy splitting when two parafermions are coupled via a backscattering region with a finite chemical potential.

Similarly, we investigate $\mathbb{Z}_{4}$ parafermions emerging as zero-modes in the helical edge states of a two-dimensional topological insulator. We show that our approach captures these systems as well, and find results which are analogous to those derived for quantum Hall edge states. This suggests that the form of the energy splitting (1) is a generic property of parafermions kept at a finite distance.

From this result, we deduce the effective coupling amplitudes between spatially separated parafermions, including the complex phases, as a function of their distance. These serve as the input parameters of the parafermionic toy model (3), which allows us to numerically simulate the braiding protocol proposed by Ref. Lindner et al. (2012); Clarke et al. (2013). Since the $\mathbb{Z}_{4}$ parafermion model is the simplest example capturing all the relevant physical phenomena, we choose this as a basis. However, we have checked that the results we obtained are valid as well for $\mathbb{Z}_{2n}$ ($n$ odd) parafermions.

For a system modelled by Eq. (3) with $\mathbb{Z}_{4}$ parafermions, having at least one vanishing coupling coefficient $t_{ij}$ at any given time ensures that the ground state is always at least $4$-fold degenerate. When at least one of the other $t_{ij}$’s is non-zero at a given time, a finite gap separates the ground-state manifold from the excited states. However, the width of this gap varies as a function of time throughout the entire braiding process, and the latter can be considered as adiabatic only if it is slow compared to the smallest gap reached during the entire process. In that case, the time-dependent change of parameters induces no transitions to states outside the ground state manifold. Importantly, one should note that the gap which determines whether or not the braid is adiabatic is related to the overlap between parafermions, and is thus much smaller than the superconducting gap, which separates the parafermionic states from the quasiparticle states. A more general case of braiding of non-Abelian anyons with pairwise interactions is discussed in Ref. Burrello et al. (2013b).

In the conventional braiding protocol, pairs of parafermions are repeatedly nucleated and braided in such a way that at a given time at most two parafermions are maximally coupled to each other. For such a protocol, we show that as the chemical potential increases to values $\mu L\gtrapprox 1$, where $L$ is the length of the system, the gap between the ground states and excited states tends to become exponentially small at certain points in the braiding process. In general, the adiabatic limit can therefore be reached only for very large braiding times.

We compare this with a modified braiding protocol, which is also based on nucleation and fusion, but in which up to three parafermions can be brought together and being maximally coupled. We show that this protocol offers larger gaps throughout the entire braiding process, so the adiabatic limit can be reached more easily.

III Coupling of parafermion modes in fractionalized systems

Let us consider two fractional quantum Hall (FQH) edge states with opposite $g$-factors, such that the system contains two coupled counterpropagating fractional edge states with opposite spins. A gap can then be opened either by allowing backscattering between the two edge states or by inducing Cooper pairing via the superconducting proximity effect (Fig. 2). The bosonized theory describing the effective 1D physics is given by the following double sine-Gordon Hamiltonian Lindner et al. (2012); Clarke et al. (2013),

[TABLE]

where $v_{F}$ is the Fermi velocity, $\Delta(x)$ the superconducting pairing potential which we take to be zero outside the superconducting region and constant $\Delta(x)=\Delta$ inside, and $g(x)$ is the backscattering strength which we take to be zero outside the backscattering region and constant $g(x)=g$ inside. The chemical potential $\mu$ is assumed to be constant throughout the system, the odd integer $n=1/\nu$ is the inverse filling factor, and $a$ is a short-distance cut-off that regularizes the bosonization procedure. The bosonic fields $\theta(x)$ and $\phi(x)$ are related to the physical right-moving ($+$) and left-moving ($-$) electrons via the bosonization identity

[TABLE]

and fulfill the following fractional commutation relation

[TABLE]

where $\Theta$ is the Heaviside function. As $n$ is odd, the physical electrons satisfy the correct anti-commutation relations. The chemical potential couples to the total charge density $\rho_{C}=\sum_{\alpha=\pm}\psi^{\dagger}_{\alpha}\psi_{\alpha}=\frac{1}{\pi}\partial_{x}\phi$. The spin density is given by $\rho_{S}=\sum_{\alpha=\pm}\alpha\psi^{\dagger}_{\alpha}\psi_{\alpha}=\frac{1}{\pi}\partial_{x}\theta$. These operators are understood to be normal-ordered with respect to the ground state.

The system described by Eq. (III) has a $2n$-fold degenerate ground state and hosts $\mathbb{Z}_{2n}$ parafermions at the interfaces between the backscattering and superconducting regions Lindner et al. (2012). For strong pairing potential $\Delta$, the energy is minimized when the $\theta$ field is pinned to $\theta(x)=\pi\hat{k}_{\theta}/n$, where $\hat{k}_{\theta}$ is an integer-valued operator. Analogously, in the backscattering region the energy is minimized when the $\phi$ field is pinned to the value $\phi(x)=\pi\hat{k}_{\phi}/n$.

If neighboring parafermionic modes are kept at a finite distance, they can couple to each other either via the tunneling of quasiparticles across the superconducting region or the backscattering region. Since the braiding protocols we study involve bringing the parafermionic modes together and then separating them, our first aim is to investigate the effect of such tunneling.

The coupling between parafermions separated by a gapped region can be studied relatively easily if the coupling between a given pair of parafermions is much stronger than their coupling to all other parafermions. In the following, we will therefore study the coupling between (i) two parafermions separated by a superconducting region and (ii) two parafermions separated by a backscattering region.

III.1 Tunneling across a superconductor

We start by considering a single pair of parafermions which are separated by a short superconducting region. The Hamiltonian in the SC region reads

[TABLE]

A system described by such a Hamiltonian is depicted in Fig. 2. This fractional sine-Gordon Hamiltonian with a chemical potential can be mapped onto an integer sine-Gordon Hamiltonian without chemical potential by defining the scaling and shift transformation

[TABLE]

which ensures that $[\tilde{\theta}(x),\tilde{\phi}(y)]=i\pi\Theta(y-x)$. After defining the Luttinger parameter $K=1/n$, the Hamiltonian reads (up to a constant term)

[TABLE]

where $\tilde{\Delta}=\Delta/n$. It is known that this sine-Gordon Hamiltonian exhibits a Berezinskii-Kosterlitz-Thouless (BKT) phase transition between a gapless Luttinger liquid phase and a gapped superconducting phase Giamarchi (2013). Even tough repulsive interactions ($K<1$) will suppress the superconducting term in the renormalization group (RG) sense, for large $\tilde{\Delta}$ the superconducting term will be dominant, and the system will be in the gapped phase. Electron-electron interactions described by the Luttinger parameter $K$ affect the gapped system only weakly, so we may take $K=1$ inside the quadratic term.

Within this approximation, the effective Hamiltonian (III.1) becomes diagonal in terms of new fermionic quasiparticles defined by $\tilde{\Psi}_{\pm}\propto\exp[-i(\pm\tilde{\phi}-\tilde{\theta})]$. We then rewrite the latter in terms of the slowly oscillating fields $\tilde{\psi}_{\pm}(x)$

[TABLE]

where $\tilde{\mu}=\mu/n$, and $[\tilde{\theta}(x),\phi(y)]=i\pi\Theta(y-x)$. The procedure of first absorbing the chemical potential term in the bosonic fields before performing the approximation $K=1$ in the quadratic part of the Hamiltonian ensures that the bosonic theory can be refermionized, while keeping the correct scaling of the parameters by the Luttinger parameter $K=1/n$.

The effective Hamiltonian can now be rewritten as a fermionic Hamiltonian $H=H_{hl}+H_{\mu}+H_{\Delta}$, consisting of a helical liquid Hamiltonian $H_{hl}$, a chemical potential term $H_{\mu}$ generated by the oscillating term and a superconducting pairing Hamiltonian $H_{\Delta}$,

[TABLE]

We can diagonalize this Hamiltonian by first going to momentum space, $\tilde{\psi}_{\pm}(x)=L^{-1/2}\sum_{k}e^{ikx}c_{k,\pm}$, and then defining the Bogoliubov transformation

[TABLE]

The Bogoliubov quasiparticles obey the fermionic commutation relations $\{\gamma_{\sigma,k},\gamma^{\dagger}_{\sigma^{\prime},k^{\prime}}\}=\delta_{k,k^{\prime}}\delta_{\sigma,\sigma^{\prime}}$ and give rise to the following diagonal Hamiltonian,

[TABLE]

where $\varepsilon_{k}=v_{F}k-\tilde{\mu}$. The information about the ground state energy splitting can be extracted from the zero-frequency limit of the retarded quasiparticle Green’s function,

[TABLE]

where $\left\langle\cdots\right\rangle$ denotes the ground state expectation value with respect to the Hamiltonian (15). The time dependence of the quasiparticle operators are given by the Heisenberg equation of motion $\dot{\gamma}_{k,\sigma}(t)=i[H,\gamma_{k,\sigma}(t)]$ which is solved by

[TABLE]

By expressing the operators in terms of Bogoliubov quasiparticles we have, after performing a shift of variables $k\to k-\tilde{\mu}/v_{F}$,

[TABLE]

where $L=(x_{2}-x_{1})>0$ is the distance separating the parafermions. Since we are interested only in energies $\omega\ll\tilde{\Delta}$, the denominator is well behaved. Because $L>0$ we can integrate over momenta using the residue theorem with a contour in the upper half part of the complex plane. The Green’s function can then be written as

[TABLE]

where we defined the correlation length $\xi_{\Delta}(\omega)=v_{F}/(\tilde{\Delta}^{2}-\omega^{2})^{1/2}$. The zero-energy contribution to the nonlocal density of states is thus

[TABLE]

The oscillating term due to the finite chemical potential $\mu$ is the main result of this analysis, and is in agreement with the results for the semi-classical limit in Ref. Chen and Burnell (2016).

An investigation of a similar sine-Gordon Hamiltonian as Eq. (III.1) using the exact form factors in order to compute the tunneling density of states suggests that the exponential suppression in Eq. (20) is indeed independent of the Luttinger parameter in the quadratic part of Eq. (III.1) and depends only on $v_{F}$ and the prefactors of the cosine Pedder et al. (2016).

Our results in Eq. (20) thus generalizes the well-known exponentially decaying and oscillatory overlap of Majorana bound states in superconducting nanowires. Indeed, Majorana fermions correspond to $n=1$, in which case we recover Eq. (1) exactly.

III.2 Tunneling across a backscattering region

Next, we compare the previous results to the system depicted in Fig. 3, i.e., a pair of parafermions separated by a short backscattering region. Defining the scaling transformation $(\theta,n\phi)=(\tilde{\theta},\tilde{\phi})$ we can rewrite the Hamiltonian in the backscattering as a sine-Gordon Hamiltonian with Luttinger parameter $K=n$:

[TABLE]

As before, we can approximately use $K=1$ in the quadratic terms since we assume the system is already gapped by the backscattering term. Hence, Eq. (III.2) is quadratic in terms of the fermionic quasiparticle operators $\tilde{\psi}_{\pm}(x)$ and can be diagonalized by a similar Bogoliubov transformation.

However, the role of the chemical potential in the Green’s function turns out to be different because the chemical potential term now contains the same operator as the cosine term. Indeed, one finds for the retarded Green’s function, with the renormalized parameters $\tilde{g}=g/K=g/n$ and $\tilde{\mu}=\mu/K=\mu/n$,

[TABLE]

where $L=x_{2}-x_{1}>0$. For $\omega,\tilde{\mu}\ll\tilde{g}$, one can use $\delta\to 0$ in the denominator. The fraction can be split into a sum of two fractions with one having two poles of order one at $\pm i[\tilde{g}^{2}-(\omega-\tilde{\mu})^{2}]^{1/2}$ and the other one having a pole at $+i[\tilde{g}^{2}-(\omega-\tilde{\mu})^{2}]^{1/2}$. Since the integrals converge in the upper half part of the complex plane, we can use the residue theorem. The imaginary part of resulting retarded Green’s function reads, using the physical parameters of the theory,

[TABLE]

In contrast to Eq. (20), the chemical potential here affects only the effective gap seen by the parafermions and hence the length scale of the exponential decay. However, it does not give rise to an oscillating term. The form of the Green’s function is in agreement with Ref. Pedder et al. (2016).

To achieve an oscillatory term in the splitting of two parafermions separated by a backscattering region, one would have to add an effective Zeeman-like term $\propto B_{z}(\partial_{x}\theta)$ instead of the chemical potential term. This can be seen via the duality transformation $(\Delta,\theta,\phi,\mu)\leftrightarrow(g,\phi,\theta,B_{z})$ of the Hamiltonian (III.1).

III.3 $\mathbb{Z}_{4}$ parafermions in topological insulator edge states

Before concluding this section, we would like to apply the same approach to another important class of parafermions which do not fall into the category of $\mathbb{Z}_{2n}$ ($n$ odd) parafermions discussed so far.

As stated in the introduction, $\mathbb{Z}_{4}$ parafermionic models are among the simplest systems capturing the relevant physical phenomena. The same refermionization approach can also be used to determine the coupling strengths of $\mathbb{Z}_{4}$ parafermions appearing in helical edge states of two-dimensional topological insulator. Such time-reversal symmetric systems can be engineered in a physical setup similar to that in Figs. 2 and 3, albeit with electronic instead of fractionalized edge states. We thus also need to consider two-particle backscattering rather than one-particle backscattering: In such systems, time-reversal symmetry forbids single-particle backscattering, but two-particle backscattering, which converts two right-movers into two-left movers, is allowed and opens a gap in the helical edge state spectrum. We can describe the 1D theory of interfaces between superconducting and backscattering regions in such an edge state by the following model Hamiltonian Orth et al. (2015)

[TABLE]

Here, $K$ is the Luttinger parameter, which quantifies the strength of the electron-electron interactions inside the helical edge state. Moreover, $g(x)$ is the effective two-particle backscattering strength, $\Delta(x)$ is the superconducting pairing amplitude, and $v_{F}$ is the Fermi velocity. As before, we will assume that each region of the edge state is dominated either by $g$ or $\Delta$. The bosonic fields $\phi$ and $\theta$ obey the usual commutation relation $[\theta(x),\phi(y)]=i\pi\Theta(y-x)$. The chemical potential $\mu$ is taken to be constant everywhere and couples to the total particle density $\rho_{C}=\frac{1}{\pi}\partial_{x}\phi$.

Let us first consider two parafermions separated by a superconducting region. The superconducting pairing term is diagonal in the physical fermions. Writing the Hamiltonian in this basis inside the superconducting region, one recovers a quadratic Hamiltonian. Deep in the gapped phase, the RG-marginal electron-electron interactions can be neglected, and one finds the following zero-energy nonlocal density of states,

[TABLE]

This result is in fact identical to the corresponding result for MBS in Eq. (1).

In the opposite limit of two parafermions separated by a two-particle backscattering region, we start by defining the canonical transformation $(\tilde{\phi},\tilde{\theta})=(2\phi,\theta/2)$. As before we can then diagonalize the backscattering term by introducing fermionic quasiparticles defined by the boson fields $\tilde{\phi}$ and $\tilde{\theta}$. Deep in the gapped phase for strong backscattering strength, the theory becomes non interacting in terms of the new fermions at the Luther-Emry point $K=1/4$. The zero-energy nonlocal density of states becomes

[TABLE]

where again no oscillations are present. This corresponding result for $\mu=0$ was derived before using a different approach in Refs. Schmidt and Pedder (2016); Pedder et al. (2017).

III.4 Coupling parameter between lattice parafermions

The goal of this section is to make the connection between the coupling coefficients $t_{ij}$ in the lattice Hamiltonian (3) and the retarded nonlocal Green’s function we calculated previously. This will allow us to make the coupling constant’s dependence on the systems parameters (chemical potential, coupling strengths, length etc) explicit.

We consider the low-energy physics of parafermion modes separated by a superconductor in either a fractionalized system or a helical system. In the former case the localized modes are $\mathbb{Z}_{2n}$ parafermions, and in the latter $\mathbb{Z}_{4}$ parafermions. The physics can be described by the following effective parafermion tight binding Hamiltonian Alicea et al. (2011); Clarke et al. (2013); Lindner et al. (2012); Chen and Burnell (2016)

[TABLE]

where $t=|t|e^{i\varphi}$ denotes the effective complex coupling coefficient. As we showed in previous sections, the two parafermions are coupled via the tunneling of quasiparticles across the junction. The eigenstates of this Hamiltonian can be labeled by the integers $q=0,\cdots,p-1$ using

[TABLE]

The spectrum of the Hamiltonian can then be written as

[TABLE]

Since the nonlocal retarded Green’s function $\tilde{G}^{R}(x_{2},x_{1},\omega)$ in the limit $\omega\to 0$ contains the information about the probability amplitude for tunneling of quasiparticles between the points $x_{2}$ and $x_{1}$, as well as the phases picked up by the tunneling, we identify the Green’s function with the hopping amplitudes of the parafermion lattice model, i.e.,

[TABLE]

IV Numerical simulation of braiding

In this section we will discuss the effect of a finite chemical potential on the braiding of parafermions using the expressions for the hopping amplitudes. In particular, we will consider the braiding of two parafermions $\chi_{1}$ and $\chi_{4}$ performed by nucleating an additional parafermionic pair $\chi_{2}$ and $\chi_{3}$. The whole system is thus described by Hamiltonian (3),

[TABLE]

To fully capture the effect of a chemical potential on braiding we will only consider coupling coefficients of the form Eq. (30). A similar form can be obtained by replacing the superconducting by a backscattering region, and the chemical potential by a Zeeman-like term. The braiding scheme we are going to investigate has been proposed and studied in several papers Lindner et al. (2012); Clarke et al. (2013); Chew et al. (2018), but the detrimental effect of a finite chemical potential has so far not been assessed. Furthermore, we propose a way to mitigate this effect by a judicious choice of the time-dependence of the parameters $t_{ij}$.

During the braiding process, the coupling strengths are varied in time by controlling the distances $L_{ij}(\tau)\in[L_{min},L_{max}]$ between parafermions $\chi_{i}$ and $\chi_{j}$. The minimal model for the coupling strengths with distance-dependent phases has the form

[TABLE]

where $\varphi_{0}$ is a global phase. In the following, we will measure all length scales in units of $nv_{F}/\Delta$ and energies are in units of $\Delta$. We will focus on $\mathbb{Z}_{4}$ parafermions, as this is the simplest example which captures the relevant physics. We stress, however, that the results will be qualitatively similar for other $\mathbb{Z}_{2n}$ parafermions.

IV.1 Braiding by nucleation and fusion

Before presenting the numerical simulations, let us discuss the general features of braiding implemented by repeatedly fusing and nucleating pairs of parafermions in a 1D system Clarke et al. (2013); Lindner et al. (2012). Such a process corresponds to an adiabatic change of the parameters of the Hamiltonian (31) along a closed loop in parameter space, and acts as a non-trivial unitary transformation within the degenerate ground-state manifold. Interestingly, if one considers only specific classes of loops which obey some particular constraints, the unitary transformations which can be obtained form a $4$-dimensional representation of the braid group Orth et al. (2015); Karzig et al. (2015).

Let us be more specific by focusing on the braiding of parafermions $(\chi_{1},\chi_{4})$. At the beginning, they are far apart from each other and therefore completely decoupled, while the nucleated parafermions $\chi_{2}$ and $\chi_{3}$ are strongly coupled. During the braiding process, we want to avoid situations in which the four parafermions are all close together and we thus require that, at each step, at least one parafermion must be far away from the others. The presence of such a decoupled parafermion ensures at least a four-fold degeneracy of the ground-state manifold, since there is at least one parafermionic operator which commutes with the Hamiltonian (31). In the case where all the parafermions are coupled, both adiabatic and non adiabatic errors may occur, as studied in the case of Majorana zero-modes van Heck et al. (2012); Sekania et al. (2017); Fulga et al. (2013).

The Hilbert space of the whole system being $4\times 4=16$ dimensional, one can in principle encounter additional degeneracies, e.g., when all the couplings are switched off ($t_{12}=t_{23}=t_{24}=0$), and the ground-state manifold becomes $16$ dimensional. The appearance of such additional degeneracies must also be avoided throughout the whole braiding process, as it would make it impossible to implement an adiabatic evolution of the system.

Loops which satisfy these constraints, i.e., at least one decoupled parafermion and no additional degeneracies, are associated with unitary transformations on the ground state manifold which form a representation of the braid group. Two notable examples are provided in Fig. 4, where we plot the loops $\Gamma_{I}$ and $\Gamma_{II}$ in the $3$-dimensional parameter space spanned by $|t_{ij}|=e^{-L_{ij}}$. Pictorial representations of the main stages of loops $\Gamma_{I}$ and $\Gamma_{II}$ are given in Fig. 1 and 8 respectively. By inspecting the parafermion position in these figures, it is clear that the two oriented loops $\Gamma_{I}$ and $\Gamma_{II}$ both lead to a clockwise braid of $\chi_{1}$ and $\chi_{4}$. An anticlockwise braid can be achieved by implementing their time-reversed partners. As a clockwise braid cannot be continuously deformed into an anticlockwise one, a loop cannot be smoothly transformed in its time-reversed partner without violating the two constraints discussed above. The braiding of parafermions by nucleation and fusion can still be considered topologically protected since its outcome is insensitive to small changes of the physical parameters $L_{ij}$, provided the four-fold degeneracy is maintained throughout the braiding process.

While different loops such as $\Gamma_{I}$ and $\Gamma_{II}$ are completely equivalent in terms of the braiding outcome, in the following we will show that they can differ a lot when it comes to their actual implementation. In particular, the requirements posed by adiabaticity may vary dramatically. The reason for that traces back to the presence of complex phases in the coupling amplitudes $t_{ij}$. As already pointed out in Ref. Lindner et al. (2012), particular values of the phases $\arg t_{ij}=k\pi/4$ ($k\in\mathbb{Z})$ can indeed lead to unwanted additional degeneracies and should therefore be avoided. In presence of finite chemical potential, however, the phases of the coupling strengths are not constant but vary with the control parameter $L_{ij}$, see Eq. (32). Values of $|\mu|\gtrsim\pi(4(L_{max}-L_{min})^{-1}$ make it impossible to prevent the unwanted phases from appearing in at least one coupling strength throughout the process. As we will show in the following, braiding protocols based on different loops handle the appearance of these unwanted phases in different ways. For generic loops such as $\Gamma_{I}$, a finite chemical potential dramatically challenges the possibility to reach the adiabatic limit within reasonable braiding times. In contrast, the loop $\Gamma_{II}$ can be shown to mitigate the effects of the chemical potential, making it possible to reach adiabaticity faster and easier.

IV.2 Numerical simulation

We use the lattice Hamiltonian (31) and simulate its time evolution using the QuTiP package Johansson et al. (2012, 2013) during a total braiding time $T$. As parameters, we choose $L_{min}=0$ and $L_{max}=10$. Due to the exponential suppression of $|t_{ij}|$, this choice ensures that two parafermions separated by a length $L_{max}$ are effectively decoupled.

We implement the braiding by varying the lengths $L_{ij}$ from $L_{min}$ to $L_{max}$ (or vice versa) with a quadratic time dependence. Between two stages of the braiding protocol, we ramp the lengths up and down, respectively, between two times $\tau_{i}$ and $\tau_{f}$ as

[TABLE]

For the simulation, it is convenient to express the braiding Hamiltonian (31) as a $\mathbb{Z}_{4}$ clock model using the Fradkin-Kadanoff transformation for $j=1,2$ Fendley (2012),

[TABLE]

with $\zeta=e^{i\pi/2}$. At every lattice site, the clock operators $\hat{\sigma}$ and the flip operators $\hat{\tau}$ generalize the Pauli matrices of the $\mathbb{Z}_{2}$ Ising model, and are given by

[TABLE]

The $\mathbb{Z}_{4}$ charge $\hat{Q}=\hat{\tau}_{1}\hat{\tau}_{2}$, which shifts all the clock variables by $\zeta$, commutes with the full Hamiltonian. We can thus label the eigenstates of the Hamiltonian with its eigenvalues $\{\hat{Q}\}=\{\zeta,\zeta^{2},\zeta^{3},\zeta^{4}=1\}$. For the rest of this section, we denote by $\ket{\psi^{0}_{q}(\tau)}$ the instantaneous ground states of the Hamiltonian $H(\tau)$. They satisfy $\hat{Q}\ket{\psi^{0}_{q}(\tau)}=\zeta^{q}\ket{\psi^{0}_{q}(\tau)}$ and $H(\tau)\ket{\psi^{0}_{q}(\tau)}=E^{0}(\tau)\ket{\psi^{0}_{q}(\tau)}$. We will from here henceforth call $\ket{\psi^{0}_{q}(\tau)}$ a parity-$q$ ground state.

IV.3 Protocol A: Simultaneous exchange

Let us first focus on a braiding protocol which implements a loop like $\Gamma_{I}$. As shown in Fig. 4, it features four main steps for times $\tau=0<\tau_{1}<\tau_{2}<T$ when the Hamiltonian reads

[TABLE]

For the numerical solutions, we use a global phase $\varphi_{0}=0.4\pi$, unless specified otherwise, to minimize the effect of accidental degeneracies. The steps of the braiding protocol are depicted in Fig. 1 and they are characterized by the following values of the lengths

[TABLE]

Between any two stages, two parafermions are moved at the same time. For example, for $0<t<\tau_{1}$, parafermion $\chi_{3}$ is moved away from $\chi_{2}$ while $\chi_{1}$ is brought closer to $\chi_{2}$. We therefore refer to such a protocol as the conventional protocol with “simultaneous exchanges”.

In the upper panel of Fig. 5, we plot the coupling energies as a function of time for a total braiding time $T=50$, and in the upper panel of Fig. 6 we plot the corresponding instantaneous spectrum of the braiding Hamiltonian (31). As expected, the spectrum consists of four $4$-fold degenerate manifolds. The energy gaps between successive energy eigenstates are minimal at three times during the braiding process. These points corresponds to the times at midpoints between different stages in the protocol in Eq. (38), i.e., $\tau\in\{\tau_{1}/2,(\tau_{1}+\tau_{2})/2,(\tau_{2}+T)/2\}$. At those times, the distances between all parafermions are large, so they become approximately degenerate.

If the braiding is adiabatic, the system follows the instantaneous ground state. In the non-adiabatic case, transitions to excited states can occur. As a figure of merit for braiding protocols, we therefore consider the overlap squared between the (desired) instantaneous parity-$q$ ground state $\ket{\psi^{0}_{q}(T)}$ at the final time $T$, and the actual time-evolved state $\ket{\psi(T)}$ of the system prepared initially in the parity-$q$ ground state,

[TABLE]

The final state $\ket{\psi(T)}$ is calculated by numerically evolving the initial state in time using the time evolution operator,

[TABLE]

where $\mathcal{T}$ denotes the time-ordering operator.

The minimal energy gap $\delta$ in the upper panel of Fig. 6 sets a speed limit for adiabaticity, which can thus be reached for $T\gg 1/\delta$. The gaps are of the order of $\delta\propto e^{-L/4}$ where $L$ the maximum distance between the parafermions. In the lower panel of Fig. 6, we plot the overlaps $O_{0}(T)$ for zero chemical potential. For our choice of parameters, large oscillations of the figure of merit are visible in the short-time regime $T\lessapprox 1000$. We attribute those to Landau-Zener-Stückelberg (LZS) interferences Shevchenko et al. (2010); Knapp et al. (2016). To support this, we plot in the App. A the overlaps between the time evolved state of the system and the excited states for a given time during the braiding process. Because the energy gaps are minimal at multiple times during the process, there exist multiple paths in the energy space connecting the ground state at $\tau=0$ to the final state at time $\tau=T$. Interferences of these paths leads to oscillations in the overlap. For $T\gtrapprox 1000$ the overlap $O_{q}$ tends to unity and adiabaticity is reached.

In App. B we plot the Berry phases for this braiding process and show that they converge to the theoretically expected values Lindner et al. (2012). We also show that inverting the path in configuration space, which amounts to exchanging the zero modes $\chi_{1}$ and $\chi_{4}$ counterclockwise rather than clockwise, yields opposite Berry phases as one would expect.

For nonzero chemical potential, the scenario becomes more complex. As can be seen from the upper panel of Fig. 7, many avoided crossings occur over the entire duration of the braid. The nature of these crossings can be explained by looking at the first stage $\tau\in[0,\tau_{1}]$ of the protocol in Eq. (38). This stage involves a transition from maximal coupling between $\chi_{2}$ and $\chi_{3}$ to maximal coupling between $\chi_{1}$ and $\chi_{2}$. If we consider the corresponding Hamiltonians $H_{23}(\tau)$ and $H_{12}(\tau)$ separately, the finite chemical potential will induce oscillations in their instantaneous spectra whose amplitudes are controlled by the exponential terms $e^{-L_{23}}$ and $e^{-L_{12}}$. Crossings will then occur in both spectra, where the energy gaps between the respective ground state and first excited state go to zero. As we discuss in App. C, the full Hamiltonian $H(\tau)$ for $0\leq\tau\leq\tau_{1}$ contains both $H_{23}$ and $H_{12}$ contributions, and their combination turns these crossings into avoided crossings. Hence, the energy gap of these avoided crossings is exponentially small in the longer of the distances $L_{23}$ and $L_{12}$.

The presence of these additional avoided crossings has consequences for braiding. As can be seen from the overlap plots in the lower panel of Fig. 7, adiabaticity is not reached even at very long times. A lower bound on the braiding time can then be estimated as $\tau_{\delta(\mu)}=1/\delta(\mu)$ where $\delta(\mu)$ is the minimal gap over the total braiding time between the ground state and the first excited state as a function of the chemical potential. As explained before, as soon as those additional degeneracies occur, this minimal gap is on the order of $\delta(\mu)\propto e^{-L}$, where $L$ is the maximal distance between the parafermions. It is important to notice that not only the gaps are smaller compared to the case without chemical potential (compared to $\delta\propto e^{-L/4}$), but a chemical potential with $\mu L\gtrapprox 1$ also significantly increases the number of these gaps throughout the spectrum. Since multiple crossings cause LZS interferences, remaining in the ground state requires slower time evolution for an increased number of crossings.

IV.4 Protocol B: Sequential exchange

From the insights of the previous subsection, we propose a peculiar braiding protocol which mitigates problems related to adiabaticity and which is based on the $\Gamma_{II}$ loop. As discussed above, the latter can be obtained via a smooth deformation of $\Gamma_{I}$, meaning that they are topologically equivalent and thus implement the same braiding operation. We show this equivalence numerically in App. B. A similar protocol for Majorana zero modes has been investigated by Ref. van Heck et al. (2012), and thus our protocol can be seen as a generalization of it to parafermionic modes. The peculiarity of the new protocol, however, is that at every time at least two parafermions remain fully coupled. This is achieved by adding extra steps in between the ones listed in (38),

[TABLE]

All of these steps are shown in Fig. 8. We refer to this peculiar protocol as “sequential exchanges protocol”, as only one parafermion is moved at any time. By inspecting for example the first main stage ($0\leq\tau\leq\tau_{1})$, one can indeed see that at first $\chi_{1}$ is brought closer to $\chi_{2}$ while keeping $\chi_{2}$ coupled to $\chi_{3}$. The latter is then moved away only at later time ($\tau>\tau_{1/2}$).

Keeping one parafermionic pair fully coupled makes the gap between the ground state manifold and the excited states as large as possible, as shown in Figs. 9 and 10 (upper panels). This makes it possible to reach the adiabatic limit for shorter braiding times with respect to the conventional protocol. This can be seen by comparing the lower panels of Fig. 9 and Fig. 6 where we plot the overlap $O_{0}(T)$ for the two protocols.

Even for a non-vanishing chemical potential, these features largely persist for the sequential braiding protocol, as can be seen in Fig. 10. The adiabatic limit is still reached for much smaller time scales compared to the conventional protocol. The sequential protocol thus permits much faster braiding operations compared to the simultaneous braiding protocol, operating even for finite chemical potentials, when the coupling phases for accidental degeneracies cannot be avoided. Note that the improvements associated with our sequential protocol do not rely on the introduction of extra terms in the Hamiltonian (31), differently from what discussed in Ref. Karzig et al. (2015).

V Conclusion

In this article, we have studied the energy splitting of parafermions separated by a finite distance for parafermions realized either in fractional quantum Hall systems or interacting edge states of two-dimensional topological insulators. We found that in addition to the expected exponential decay with the distance, the energy splitting generally also features an oscillatory term if the gapped region between the parafermions is at a nonzero chemical potential. This property is analogous to a similar result for Majorana bound states and confirms a recent result for parafermions based on a semiclassical approximation.

We use the energy splitting to deduce the complex hopping amplitudes between parafermions in a 1D chain and use this model to perform a numerical simulation of the complete braiding process. Using the conventional braiding protocol, we find that a finite chemical potential induces level crossings at isolated points in the braiding process which make it very hard to reach an adiabatic limit. Therefore, we have proposed a modified braiding protocol, which offers larger gaps at all intermediate stages even in the presence of a nonzero chemical potential. We have argued that an adiabatic limit is easier to reach in this improved protocol.

Acknowledgements.

The authors acknowledge support by the National Research Fund, Luxembourg under grants ATTRACT 7556175, INTER 11223315, AFR 11224060 and PRIDE/15/10935404. AC acknowledges support by the Würzburg-Dresden Cluster of Excellence on Complexity and Topology in Quantum Matter (EXC 2147, project-id 39085490).

Appendix A Overlap between time evolved wave function and excited states

For any finite braiding time, non-adiabatic transitions outside the ground state manifold can occur. These can be quantified by computing the overlap between the time-evolved wavefunction of the system $\ket{\psi(\tau)}=U(\tau)\ket{\psi_{q}^{0}(0)}$ and the instantaneous energy eigenstates $\ket{\psi_{q}^{j}(\tau)}$ ($j=0,1,2,3$) for a given parity $q$. We define the probability to end up in the $j^{th}$ excited state at time $\tau$ by

[TABLE]

As can be seen in the Fig. 11, non-adiabatic transitions occur every time the gaps reach a minimum. As these involve multiple levels, and occur multiple times throughout the braiding process, adiabaticity is lost. These transitions can be suppressed by increasing the total braiding time as can be seen from the lower panels of Fig. 11. This scenario persists for non-zero chemical potential, where because of multiple avoided crossings in the spectrum, the adiabatic regime is only reached for much longer braiding times. As can be seen from Fig. 12, sequentially exchanging the parafermions permits much faster braiding operations.

Appendix B Berry phase

In the case of an adiabatic exchange, the final and initial state will only differ by a phase,

[TABLE]

where the dynamical phase and the Berry phase read respectively

[TABLE]

The former only depends on the instantaneous ground state energy $E^{0}(\tau)$, whereas the latter explicitly depends on the parity of the ground state. By computing the phase differences between different parity states the dynamical phase will drop out and we obtain thus information about the Berry phase. Defining $\phi_{q}=\arg\braket{\psi^{0}_{q}(0)}{\psi_{q}(T)}$, we define the Berry phase difference

[TABLE]

Adding $2\pi$ and taking the modulo $2\pi$ ensures that $\beta_{q}\in[0,2\pi]$. According to Ref. Lindner et al. (2012), the Berry phase which the parity-$q$ state picks up during a full braiding process described by the Hamiltonian (31) for $\mathbb{Z}_{p}$ parafermions is given by

[TABLE]

When the complex phases of all coupling strengths are equal, then the integer $k$ is determined by the global phase $\varphi_{0}$ of the coupling strengths of the Hamiltonian (31) by

[TABLE]

In other words, $k$ is given by the integer part of $p[\varphi_{0}+\pi]/(2\pi)$. The difference of $\pi$ with respect to Ref. Lindner et al. (2012) stems from a slightly different convention in defining the coupling coefficients. For the choice of global phase in the main text $\varphi_{0}=0.4\pi$ we have $k=2$ so that $\beta_{q}=\pi-\gamma_{q}\text{ mod }2\pi$. In Fig. 13 we plot the Berry phases as a function of total braiding time using the sequential exchange protocol corresponding to a clockwise exchange (left panel) of the zero modes ($\chi_{1},\chi_{4}$) as well as the counterclockwise exchange of the pair (right panel). This is numerically implemented by interchanging $L_{12}\leftrightarrow L_{24}$. The two scenarios correspond to opposite loops in the control parameter space and give rise to opposite Berry phase differences. Indeed, if $U_{14}$ is a representation of the braid group which exchanges the parafermions $\chi_{1}$ and $\chi_{4}$ clockwise then

[TABLE]

The operator representing the counterclockwise exchange is represented by $U_{41}=U_{14}^{-1}$. It thus follows that

[TABLE]

The Berry phases in Fig. 13 converge to the theoretical values, and the counterclockwise exchange yields opposite Berry phases.

In Fig. 14 we plot the Berry phase differences as a function of the global phase (left panel), which we compare to the expected Berry phases (45) (right panel). The Berry phase differences are computed for $T=2000$ at zero chemical potential. Note that jumps in the Berry phases occur when the global phases crosses the values $m2\pi/p$ for $m\in\mathbb{Z}$: The spectrum becomes degenerate at these points and thus adiabaticity is no longer reached. The Berry phases thus become ill-defined. Jumps between [math] and $2\pi$ are also to be expected because we are taking the modulo $2\pi$ of the phase differences.

In the lower panel of Fig. 14, we show that the sequential exchange protocol gives rise to the same Berry phases, and thus ultimately to the same braiding operators as the simultaneous exchange protocol.

Appendix C Origin of the avoided crossings

The origin of the avoided crossing in the simultaneous exchange protocol can be investigated by only considering the first stage Hamiltonian

[TABLE]

as well as the individual Hamiltonians $H_{12}$ and $H_{23}$ separately. In the left panel of Fig. 15 we plot the individual spectrum of $H_{12}(\tau)$ and $H_{23}(\tau)$. The finite chemical potential induces oscillations and degeneracies which appear throughout the spectra. The amplitude of these oscillations are controlled by the exponential prefactor $e^{-L_{ij}}$. Consider now the gap between the ground state and the first excited state, which we plot in the right panel of Fig. 15. In $H_{23}$, due to the oscillations, the gap will vanish at certain points, before converging to zero due to the exponential suppression for increasing times. When the gaps of $H_{23}$ are exactly zero, which occurs during the period where the parafermionic pair $\chi_{2,3}$ is still coupled, the energy levels of $H_{12}$ are exponentially small. The opposite argument holds for the gaps of $H_{12}$.

In the right panel of Fig. 15, we show that the gap $\delta E_{I}(\tau)$ of the full Hamiltonian has the form $\delta E_{I}(\tau)\geq\text{max}\{\Delta E_{23}(\tau),\Delta E_{12}(\tau)\}$. This means that whenever the gap of $H_{12}(\tau)$ closes due to the oscillations, $H_{23}(\tau)$ will act as a perturbation and lift the degeneracy, and vice versa. In order to minimize the non-adiabatic transitions we can increase the total gap $\delta E_{I}(\tau)$ by keeping one of the gaps of the individual Hamiltonian as large as possible. This can be performed by keeping the corresponding parafermionic pair $\chi_{2,3}$ coupled as one approaches the other parafermion $\chi_{1}$ in order to couple $\chi_{1,2}$. The exchange becomes sequential rather than simultaneous. The same argument will also hold for the remainder of the braiding process.

An analytical expression for the gap with finite chemical potential can be found by mapping the first stage Hamiltonian (49) to a $\mathbb{Z}_{4}$ clock model using the Fradkin-Kadanoff transformation (36). This gives rise to a $16\times 16$ matrix which can be block diagonalized by rotating the Hamiltonian in the parity basis. The resulting four $4\times 4$ block Hamiltonians belong then to a given parity. The spectra for are the same for all parities, so we may choose the first of these resulting $4\times 4$ Hamiltonian matrix, which corresponds to the parity 0. It reads:

[TABLE]

where $\varphi_{ij}(\tau)=\arg[t_{ij}(\tau)]=\varphi_{0}+\mu L_{ij}(\tau)$. For $\tau<\tau_{1}/2$, the parafermion pair $\chi_{2,3}$ is strongly coupled and thus $|t_{23}(\tau)|\gg|t_{12}(\tau)|$. Avoided crossings will then appear in the spectra between the ground state and first excited state at time $\tau_{c}$ when

[TABLE]

By numerically computing the eigenvalues at these times we find that the gap between the ground state and the first excited state at $\tau_{c}$ reads:

[TABLE]

where we Taylor expanded the square roots to first order in $\sqrt{2}|t_{12}(\tau_{c})|/|t_{23}(\tau_{c})|$. As expected, this gap goes to zero in the limit where $|t_{12}|\to 0$.

The smallest gap of the system $\delta(\mu)=\delta_{\mu}(\tau_{c})$ is obtained by solving Eq. (51) for $\lambda=1$. The case when $\tau>\tau_{1}/2$ is done by replacing $L_{12}\leftrightarrow L_{23}$, and the gap goes to zero when $|t_{23}|\to 0$ . As an example, the first gap in Fig. 7 appears at $\tau_{c}\approx 2.09$ with $\delta(\mu)\approx 9.5\times 10^{-4}$.