Model wavefunctions for an interface between lattice Laughlin and Moore-Read states

We use conformal field theory to construct model wavefunctions for an interface between lattice versions of a bosonic ν = 1/2 Laughlin state and a fermionic ν = 5/2 Moore-Read state. The properties of the resulting model state, such as particle density, correlation function and Rényi entanglement entropy are then studied using the Monte Carlo approach. Moreover, we construct the wavefunctions also for localized anyonic excitations (quasiparticles and quasiholes). We study their density profile, charge and statistics. We show that, similarly to the Laughlin-Laughlin case studied earlier, some anyons (the Laughlin Abelian ones) can cross the interface, while others (the non-Abelian ones) lose their anyonic character in such a process. Also, we argue that, under an assumption of local particle exchange, multiple interfaces give rise to a topological degeneracy, which can be interpreted as originating from Majorana zero modes.


I. INTRODUCTION
One of the characteristic features of the topological orders is the existence of nontrivial physical phenomena at the edges or interfaces with an another topological phase. While the former can be used to characterize a single topological phase [1,2], the latter can tell us how two different topological phases are related to each other (e.g. if one of them can be transformed into the other by anyon condensation [3][4][5][6][7][8][9][10]). In experiments, the interfaces can be potentially useful for example for isolating a certain edge mode to prove its existence [11,12]. From the perspective of applications, even when both sides are Abelian, the interfaces can have non-Abelian zero-energy modes similar to the Majorana zero modes [13][14][15], which can encode quantum logic gates [14][15][16][17].
The interfaces involving a non-Abelian state at least on one side are more complicated than their Abelian counterparts and our knowledge of them is less complete. An attempt at a systematic understanding of such interfaces is the topological symmetry breaking formalism (anyon condensation) [3][4][5][6][7][8][9][10]. This approach can tell us which anyons can cross the interface, as well as whether the interface is gapped or gapless and the nature of the interface modes. However, there is a number of things one cannot determine in this way, such as the nonuniversal, microscopic details of the interface and the anyons. Moreover, while being very general, it is also abstract. Thus, there is a need of concrete examples embodying these mathematical concepts and making them more accessible. There are some analytically solvable models describing interfaces, although they are so far limited only to the cases with a nonchiral topological order on both sides [4,18,19].
However, the study of interfaces between chiral topological orders is also important. In particular, the fractional quantum Hall states, which have been thoroughly studied experimentally [20], are chiral. The interfaces between Abelian quantum Hall states were already created, e.g. in Camino et al.'s interferometric experiments which aimed at demonstrating the fractional statistics [21]. There are also some proposals for creating interfaces with a non-Abelian state on at least one side [11,22,23].
There is also a considerable effort directed at creating lattice versions of FQH states. Numerical works predict the appearance of discretized FQH states in a range of simple lattice models [24][25][26][27]. Experiments have shown the creation of such topological orders in a moiré superlattice in graphene subjected to a strong magnetic field [28]. A number of alternative experimental designs were proposed, involving e.g. optical lattices [29][30][31][32][33][34][35][36][37]. In this way, one can also create discretized versions of bosonic FQH states, so far not seen in experiments. Therefore, the study of interfaces between lattice FQH states -both bosonic and fermionic -has experimental relevance.
Although, to the best of our knowledge, exactly solvable models of FQH interfaces (neither Abelian nor non-Abelian, neither lattice nor continuum) have not been found so far, such systems were studied both numerically and analytically. For Abelian states, there are systematic field theory approaches based on K matrices (related to topological symmetry breaking) [38,39] or Chern-Simons theory [40,41]. A more detailed, microscopic treatment was achieved by constructing model wavefunctions from conformal field theory, either in the continuum, by forming matrix product states from CFT vertex operators and combining matrices belonging to different FQH states [22,23], or in the lattice by combining the vertex operators directly [42]. Numerical calculations were also performed using exact diagonalization, although this method has limitations as it can be used for small system sizes [22,23,43].
In the case of non-Abelian states, apart from the topological symmetry breaking studies [5][6][7][8][9][10], there have been some works which tried to determine the properties of the interface modes at some gapless interfaces (Halperin/Pfaffian, Pfaffian/anti-Pfaffian) and/or the behaviour of anyons in the presence of these interfaces by using the "giant hole approach" [44] or by cou-pling the effective edge theories [11,[45][46][47][48][49]. Moreover, the latter approach was used for a Pfaffian/Pfaffian interface to determine the gapping conditions and the scaling of entanglement entropy for the gapped case [50]. Microscopic studies are very rare. For a Halperin/Pfaffian interface, a model wavefunction was constructed from infinite-dimensional MPS based on conformal field theory [12]. For a Pfaffian/anti-Pfaffian interface, DMRG studies was performed [51,52]. In both cases, continuum systems were investigated.
Here, we perform a microscopic study of an another example: an interface between a fermionic Pfaffian state and an Abelian bosonic Laughlin state. In contrast to the former works on non-Abelian interfaces, we focus on the lattice case. We construct model wavefunctions from conformal field theory and study their properties using Monte Carlo methods. These wavefunctions describe both the ground state and the localized anyonic excitations (both quasiholes and quasielectrons). From the ground state wavefunction (Section II), we obtain quantities such as particle density, correlation function or entanglement entropy. We note that the latter was very rarely studied for non-Abelian interfaces. Using the anyon wavefunctions (Section III) we can investigate the density profile, charge and statistics of such excitations. Such a microscopic study of localized non-Abelian anyons of a chiral topological order in the presence of an interface was, up to our knowledge, not performed before. Moreover, while the previous microscopic works on FQH interfaces concentrated on a single interface, here we also discuss the case of multiple Laughlin islands within a Moore-Read system, arguing that for local interactions across the interfaces, the parity of the particle number at each island cannot be changed locally, leading to the emergence of topological degeneracy resembling that of Majorana zero modes.

II. THE GROUND STATE
A. The construction of the wavefunction Our approach follows the CFT construction for lattice quantum Hall wavefunctions from Refs. [42,53,54], which are based on the idea devised by Moore and Read for continuum quantum Hall wavefunctions [55]. We consider a system of N sites on a plane, each one at a complex position z j = x j + iy j . The sites can be occupied by at most one particle, with n j ∈ {0, 1} denoting the occupation number of site j. That is, each site hosts a fermionic or hardcore bosonic degree of freedom (a system can contain sites of both types). We set the charge of each particle to unity. In contrast to the continuum case, we do not assume a constant magnetic field, but rather set it to zero everywhere except at the positions of the lattice sites. That is, we attach an infinitely thin solenoid, containing η i ∈ R + flux quanta, to each site. The total number of flux quanta is N φ = i η i and can be different from N , which means that we can define two kinds of filling factors. In the simplest case of a single MR or Laughlin state, the particle number M is conserved, so we can define a "lattice filling" ν lat = M/N and a "topological filling" ν = M/N φ = 1/q, q ∈ N + . The former describes the degree of discretization (the lower ν lat , the denser are the lattice points, i.e. the closer we are to the continuum), while the latter tells us which state do we discretize (i.e. it is equal to the Landau level filling factor of the corresponding continuum state).
Any wavefunction in our system can be written in the occupation number basis, where n is a vector of occupation numbers of all sites, |n is a corresponding basis state (to define it unambiguously we fix the order of fermionic creation operators in the definition of |n to be the same as the order of site indices), Ψ(n) are the unnormalized wavefunction coefficients and C is the normalization constant (which we fix to be real without loss of generality). To construct a model CFT wavefunction for a single lattice quantum Hall state, a vertex operator V i (z i , n i ) of a certain conformal field theory (depending on which kind of state we want to create) is assigned to each site. The wavefunction is then given by the correlator of the vertex operators for all sites, This method can be generalized to interfaces [42]. In general, for two given quantum Hall states, there can be many different types of interfaces, depending on the interactions between the particles on the two sides. A wavefunction for a particular type can be created by forming a correlator of the form (2), but made from the vertex operators of two different CFTs. Such a quantity is welldefined when the two CFTs can be embedded in a third one, which puts a restriction on the states for which this method can be applied. We note that in Refs. [12,22,23], model wavefunctions for interfaces (Laughin/Halperin, Pfaffian/Halperin) in continuous systems were created by patching together infinite-dimensional matrix product states, derived from conformal field theory, representing the two different fractional quantum Hall states. Our approach is similar, but it uses the CFT vertex operators directly, without the need of matrix product state representation.
So far, we employed this method only for Abelian Laughlin states [42]. Here, we use it to study an interface between a bosonic Laughlin state and a non-Abelian fermionic Moore-Read state, both at topological filling ν = 1/2. They are described by U (1) 2 and U (1) 2 × Ising CFTs, respectively. The embedding condition is satisfied, as the U (1) 2 part is the same for both states.
We assume that the system consists of two parts. The left one with, which consists of the first N L sites with η i = η L each, is described by the MR state. In the right one, consisting of next N R sites up to N = N L + N R , the particles are in the Laughlin state. The number of flux quanta per site is set to a constant within each part, η i = η R for i = N L and η i = η R for i > N L , but it can differ between the parts, i.e. we can have η L = η R . We note that in general the two parts of the system can be of any shape and can be split into disconnected regions, but we will use the L and R labels for simplicity, as this is the geometry that we will study numerically in this work.
More specifically, the planar systems considered in this work consist of sites arranged in a square lattice of size (N xL + N xR ) × N y . The interface is parallel to the y direction, as shown in Fig. 1 (a). Without the loss of generality, we set the lattice constant to unity and the position of the interface to x = 0.
The vertex operators are given by (3) Here, the V Laughlin,i (z i , n i ) and V Ising,i (z i , n i ) are the Laughlin-like and Ising-like parts of the vertex operator, respectively.
with φ(z i ) being a free chiral bosonic field and q = 2 in our case, ensuring the topological filling ν = 1/2 on both sides. We note that the same expressions can be used in the case q = 1 (an interface between a bosonic MR state and a fermionic integer quantum Hall state, both at ν = 1), but then double occupancy of the L sites has to be allowed [54]. In this work we restrict the study only to the q = 2 case. The Laughlin-like part is the same for all sites, except from the different values of η i on the two sides. The Ising part is where ψ(z i ) is a chiral Majorana field. In contrast to the Laughlin term, the Ising one is assigned only to the L sites. By evaluating the correlator (2), we obtain the unnormalized wavefunction coefficients where z i , z j denote coordinates of filled L sites, and is the charge neutrality condition (M = i n i is the total number of particles in the entire system). The charge of particles is compensated by a background charge −η i /q assigned to every site i. We note that, in contrast to the Laughlin-Laughlin interface [42], in the Pfaffian-Laughlin case the charge neutrality condition (7) enforces the conservation of the total particle number. The numbers of particles M I , I = L, R, on each side of the interface can nevertheless fluctuateprovided the same number of particles is annihilated on one side and created on the other. This means that the number of bosons and fermions in the system is not conserved. However, because the Pfaffian in (6) is nonzero only if the number of particles in the L part is even, the fermionic parity is conserved, i.e. the particles can be created and annihiliated only in pairs.
So far, we have worked on a plane. However, the numerical investigation of certain properties of our wavefunction (such as the entanglement entropy scaling) is easier for a cylinder. We start with an L x × L y rectangle on the complex plane, within which we put the sites at the positions W j =x j + iỹ j . To impose periodic boundary condition in the y (i.e. imaginary) direction, we consider the following mapping The resulting z j 's are then substituted to (6). We will consider systems on a square lattice, of size (N xL +N xR )× N y (here L y = N y ), with the interface being parallel to the periodic direction (see Fig. 1). We require that our wavefunction is scale-invariant on a plane and inversion-invariant on a cylinder. That is, when subjected to the transformations z i → cz i , c ∈ C and z i → 1/z i , it should not change (except from a multiplication by a constant). It turns out that to fulfill both conditions, one has to set η L = 3/2, η R = 1. We enforce this condition throughout this work.
We note that at these values of η, the possible system sizes can be divided into four classes, corresponding to (N L mod 8) = 0, 2, 4, 6, differing by the way the charge is distributed among L and R parts within each configuration |n . Each such a configuration has well-defined M L and thus it corresponds to the charge M L − 3 4 N L , as a background charge −η/q = − 3 4 is associated with every L site, and each particle has unit charge. The charge is not well-defined for the entire |Ψ state, as it is a linear combination of different configurations with different values of M L , because pairs of particles can be transferred across the interface. However, a charge modulo 2, ∆Q = M L − 3 4 N L mod 2 is well-defined, and equal to ∆Q = 0, 0.5, 1, 1.5 for (N L mod 8) = 0, 2, 4, 6.
In other words, only in the first class both parts of the system can be charge-neutral by themselves if the system is cut through the interface. In the third class, the charge neutrality is achieved if an integer charge (e.g. in the form of a quasiparticle) is introduced in each part. In the two remaining classes, such charge would have to be fractional.
Finally, let us remark that for some values of η, it is possible to derive a parent Hamiltonian for single Moore-Read [54,56] or Laughlin [57] lattice quantum Hall states. However, it is not straightforward to extend these calculations to the case when the system is described by two different CFTs. Another way to connect our wavefunction to a Hamiltonian is to find a short-range Hamiltonian, whose ground state is approximated by our wavefunction, following the approach from Ref. [54]: diagonalizing the Hamiltonian numerically and optimizing its coefficients to maximize the overlap between the ground state and our wavefunction. However, this would require extensive exact-diagonalization or DMRG calculations, and we leave it for future works.
Often the wavefunction itself, which may, but does not have to, be related to a Hamiltonian, can provide insights into the inner structure of given topological order. For example, the Laughlin wavefunction revealed the physical mechanism of FQHE and the nature of the fractionalized excitations [58]. The Kalmeyer-Laughlin wavefunction provided a vital example of a chiral spin liquid [59,60]. The concepts of Haldane hierarchy [61] or composite fermions [62] were also embodied in wavefunctions. Speaking of the interfaces, Regnault et al. used a model wavefuncion to study the nature of the interface modes [12,22,23], while we employed our model state to determine the properties of localized anyons in presence of the interface [42]. Therefore, we believe that the study of the wavefunction itself is important and thus, in this work we focus solely on |Ψ .

B. Numerical results -particle density
Once we have the wavefunction (6), we can study its properties numerically using Monte Carlo methods. In particular, it is straightforward to obtain the average particle density n i . On a cylinder, the density is constant in the y direction, so we define the density as a function of x We investigate this quantity for a number of systems with different sizes, which are shown in 2 (a)-(d). When the cylinder is thin, the states display large oscillations in density within the L part. In particular, at N y = 1, we have either n i = 0 or n i = 1, with no fractional values, reminiscent of the thin torus limut of the continuum FQH states. As the cylinder gets wider, the density in the L bulk becomes close to 3/4, as expected for a η = 3/2, ν = 1/2 MR state. Apart from the thin cylinders, the density inhomogenities exist mostly near the edges and the interface and get smaller as N y increases. In the R part, independently of N y , the density is close to 1/2 everywhere except from the vicinity of the cylinder. This value is expected for a η = 1, ν = 1/2 Laughlin state.
Near the interface, the density has to drop from 3/4 to 1/2. How this happens exactly depends on the size of the system. As can be seen in Fig. 2 (a)-(d), for thin cylinders, the four (N L mod 8) groups display four qualitatively different patterns of particle density. This is most striking when comparing (N L mod 8) = 4, which has an additional "step" (i.e. a local maximum of density near the interface in part R) to (N L mod 8) = 0, where such a feature is absent. For wider cylinders, the density profiles in all four groups become similar.
In the case of planar systems, the density patterns are more complicated, as the translational invariance is lost, and the density inhomogenities exist near the interface and all three edges of the L part. Two examples are shown in Fig. 2 (e),(f).
The background charge −η i /q changes abruptly from −3/4 to −1/2 at the interface. On the other hand, we have seen that the particle density changes more smoothly. This means that some excess charge is accumulated near the interface on each side. Let us investigate this more closely for the cylindrical systems we have already studied. We define the excess charge as a function of x as and a total charge accumulated in part I as Apart from thin cylinders, the excess charge is concentrated mostly near the interface x = ±0.5, with Q(−0.5) ≈ −Q(0.5). Thus, for wide enough cylinders, we have Q(−0.5) ≈ Q L . Therefore, let us study the latter quantity as a function of system size. To be precise, instead of Q L itself, we investigate the charge in part L per unit of interface length, i.e. Q L /N y . If there is a fixed density pattern in the thermodynamic limit, then this quantity should converge to a fixed value.
The results are shown in Fig. 3. The different colors correspond to different system sizes in the x direc- tion, while the marker shapes refer to the four classes (N L mod 8). It can be seen that at N y = 1 (the rightmost points from every class) the charge modulo 2 equals ∆Q. As we increase N y , Q L /N y in all four classes seems to display convergence towards a fixed, negative value of Q L /N y , lying between −0.03 and −0.07.
A similar charge accumulation at the interface was observed in Ref. [51] for the Pfaffian/anti-Pfaffian case. In this work, the accumulated charge on the two sides of the interface in the thin-torus limit was equal to ±1/4 (i.e. the charge of a non-Abelian quasiparticle). In our case, it is equal to a charge of an Abelian quasiparticle or a particle, depending on N L . Moreover, in Ref. [51] the authors observed that the dipole moment with respect to the interface is constant also outside the thin-torus limit. This dipole moment was of topological origin, as it arose from the difference of the quantized Hall viscosities on both sides of the interface. In contrast, in our work the dipole moment is not constant. This can be seen easily by noting that for large enough N y , Q(−0.5)/N y is roughly constant, therefore Q(−0.5) should grow more or less linearly with N y . Since most of the charge is concentrated at x = ±0.5, this means that the dipole moment must grow with N y . To verify this, we evaluated the dipole moment explicitly, confirming that it is indeed not conserved.

C. Numerical results -correlation function
We expect that our interface is gapless. This is because the edges of the Laughlin state are described by a chiral Luttinger liquid [1], while the MR state has also a single Majorana fermion edge mode [63,64]. In the effective interface theories considered so far for various non-Abelian interfaces [11,[45][46][47][48][49][50], the Majorana mode can be gapped only when paired with a second Majorana mode. Since there is just a single Majorana mode in the system, we expect that it cannot be gapped.
It is expected that gapless systems generated by shortrange interactions have exponentially decaying correlation functions. In our case, we do not have the parent Hamiltonian, so we cannot ensure that our wavefunction indeed can be generated by a short-range interaction. But assuming it is the case, the correlation function would give us some indication of the gaplessness of the interface.
The correlation function is given by and can be easily computed using Monte Carlo. For the ease of presentation, we choose sites with x i = x j and investigate C ij as a function of y j − y i . The results for the different values of x i = x j are shown in Fig. 4.
In the bulk, the correlation function seems to decay roughly exponentially. However, near the edges its decay seems to be slower. Near the interface, the values of the correlation function halfway across the cylinder seem to be located roughly between the results for the bulk and the edge, still showing the lack of exponential decay. This suggests that the interface is indeed gapless (provided that it is generated by a short-range interaction).

D. Entanglement entropy
The topological properties of the interface can manifest themselves in the entanglement entropy when the cut coincides with the interface. While this issue was studied using field theory for Laughlin/Laughlin [38,39] or Pfaffian/Pfaffian interfaces [50], for the Pfaffian/Laughlin interface, up to our knowledge, there were no predictions how the entropy should scale with the interface length. Thus, we are going to study the entropy numerically, using the Monte Carlo method outlined in [65,66] (see also our previous work where we used this method [42,53,54]). Within this approach, the second Rényi entropy can be obtained by sampling two independent copies of the system.
In Fig. 5, we show the entanglement entropy scaling for three series of systems, of size: (8 + 8) × N y , (10 + 4) × N y and (12 + 4) × N y . The cut is parallel to the interface, as shown in Fig. 1. We are interested especially in the four positions of the cut: in the bulks of the two sides, precisely at the interface and right next to the interface on the left (i.e. x = −1).
The scaling in the bulks in Fig. 5 corresponds to the position x = −N xL + N xL /2 + 1 for the L side and x = N xR /2 for the R side, with denoting the floor function. If the interface wavefunction indeed describes the expected topological orders, then, by applying a linear fit, we should recover the expected value of γ: γ L = ln(8)/2 and γ R = ln(2)/2, corresponding to the topological entanglement entropies of the MR and Laughlin state, respectively. These values are indicated by blue and red ticks, respectively, on the y axes of the plots. The red and blue lines denote the fits. Because for thin cylinders the linear scaling is distorted by finite-size effects, we discard these systems from the calculation. That is, in the fit we include only the data points denoted by filled symbols. The fits seem to cross the N y = 0 line relatively close to the predicted values. For γ R , the agreement is good: we obtain 0.368 ± 0.005, 0.361 ± 0.009, 0.36 ± 0.02, for (8 + 8) × N y , (10 + 4) × N y and (12 + 4) × N y , respectively, compared to ln(2)/2 ≈ 0.347. The uncertainties here are the fit uncertainties, without the inclusion of Monte Carlo uncertainties. This confirms that the R part has a ν = 1/2 Laughlin-type topological order. For γ L , we obtain −1.08 ± 0.13, −0.82 ± 0.17, −1.02 ± 0.15, respectively, compared to ln(8)/2 ≈ 1.03. That is, the agreement is worse, and the error bars are much bigger. In addition, the result for part L seems to depend strongly on the position of the cut and on which data points we take into account on the fit. Also, while the fits in Fig. 5 were performed without the inclusion of MC error bars in the weights, including them makes the result even more dependent on the number of included data points. The detailed analysis is contained in Appendix A.
Nevertheless, the fitted values oscillate around the predicted value and are clearly nonzero. Thus, we conclude that the L part is also topologically ordered, and the results are consistent with the Moore-Read topological order, although not indicating it clearly.
What happens with the entropy when the cut coincides with the interface (black markers in Fig. 5)? For almost all the investigated systems, the entropy at the interface is lower than in the bulks of both sides (excluding some thin cylinders). However, as N y increases, the interface entropy increases faster than the R bulk entropy, thus we can expect that the former will finally dominate over the latter. For large enough N y , the scaling seems to be linear. Because we do not have compelling theoretical arguments that in this case such a scaling is expected in the thermodynamic limit, we do not rule out the possibility that the percieved linear dependence is in fact nonlinear, and the nonlinearity would show up for larger N y . Nevertheless, assuming that it is linear, we perform the fit. The obtained values are close to ln(8)/2, i.e. the topological entanglement entropy of the left part. If this is indeed the case, this is similar to the case of Laughlin states at fillings 1/q L , 1/q R such that q R = a 2 q L , a ∈ N + [38,39,42]. However, the fitted parameters are subjected to the same distortions and uncertainties as γ R , thus we cannot conclude that it is indeed the case.
How far does the influence of the interface extends into the L and R parts? Next to the interface on the right (x = 1), the entanglement entropy values for large enough N y are similar to the ones in the bulk R part. However, on the left, the influence of the interface is apparent in the first column of sites next to it (x = −1). The values of the entanglement entropy (except from some low-N y systems) are lower than in the L bulk of the same system (see the violet markers in Fig. 5). The fit for large N y also yields values roughly close to the theoretical value of γ L , although again the results are uncertain due to the dependence of the fitted value on the data points included. Thus, we do not rule out the possibility that near the interface there might be some variation of γ, e.g. a similar increase as in the Laughlin-Laughlin interfaces [42].

III. THE SYSTEMS WITH ANYONS
Having determined the ground state wavefunctions, we now wonder, what are the properties of the localized anyonic excitations above the ground state.

A. The construction of the wavefunction
The wavefunctions including anyons can be obtained by inserting further vertex operators into the correlator (2). These operators depend on parameters w i , the com-plex coordinates of the anyons. The state is given by There are two differences between (14) and (1). First, now the wavefunction coefficients, as well as the normalization constant, depend on the external parameters, the anyon coordinates w. Secondly, there can be more than one degenerate state, hence we introduced the index α. Moreover, while for the ground state the fermion parity conservation was enforced by the Pfaffian factor, in the presence of anyons the correlators are nonzero both for even and odd M L . To restore the fermion parity conservation, we assume that the interaction generating our wavefunction allows to exchange particles through the interface only in pairs. Then, the Hilbert space divides into two disconnected parts, with even and odd M L . We focus on the case of even M L . While the particle coordinates are restricted to the lattice sites, the quasihole coordinates can be located anywhere on the plane/cylinder. In such a way, we will be able to move them smoothly, which will be important when evaluating their statistics.
We study two classes of anyons of the Moore-Read state. The basic non-Abelian excitations are constructed using the following vertex operator [56,67] where σ(w i ) is the holomorphic spin operator of the chiral Ising CFT, and p i = 1/2, p i = −1/2, correspond to a quasihole and a quasielectron, respectively. We note that the latter are difficult to construct in the continuum [68], whereas for the lattice their construction is simple -it requires only flipping the sign of the p i . The other group of excitations consists of Laughlin-like Abelian anyons, described by the vertex operator [69] where p i is now integer. This vertex operator describes also the excitations of the Laughlin state. Thus, these anyons are valid topological excitations of the entire system. In contrast, the ones generated by the vertex operator (15) are valid topological excitations only within the L part. Nevertheless, technically we can also attempt to put such an anyon in the R part and see what happens. We will refer to these two groups as "Abelian" and "non-Abelian" for brevity, although the reader should bear in mind that the anyonic content of the Moore-Read state is richer than the considered cases. We denote the numbers of non-Abelian and Abelian anyons as R NA and R A , and their total number as R = R NA + R A . For convenience, we will also assume that the anyons are indexed in such a way that the first R NA are non-Abelian, and the rest are Abelian.
The wavefunction coefficients for even M L are now given by the following correlator where the index α means that we take only the conformal block where the Ising fields fuse to α, and I α and J are the Ising and Jastrow parts of the wavefunction, respectively. The latter is given by where the charge neutrality is enforced by The Ising part depends on the number of non-Abelian anyons. If there are none, then there is just one conformal block, and we have If R NA = 2, then there is also one conformal block, with If R NA = 4, then there are two conformal blocks, corresponding to the situation where the pairs of σ fields fuse to α = I or α = ψ. The Ising part is given by where m α = 0, 1 for α = I, ψ, respectively,

B. Anyon charge and density distribution
Let us now verify that the anyons are well localized and that their charges agree with the theoretical prediction −p i /q. We define the excess charge at site i in the presence of anyons as the differencẽ where the index "an" means that the density is evaluated in the presence of anyons, and "GS" means the density evaluated in the ground state (i.e. without anyons). Note the difference from the definition used in Sec. II B -now we do not care for the background charge, but only for the difference of the particle density distributions with and without anyons (otherwise we would always obtain some excess charge near the interface). The charge of the localized anyon is studied by investigating the charge accumulated within some radius around the anyon position If the anyons are localized,Q k (r) should converge to a fixed value quickly as we increase r. Fig. 6 (a) shows the distribution of chargeQ i in the case of two Abelian and two non-Abelian anyons on a cylinder. Each of the anyons is located in a part where it is a valid topological excitation. It can be seen that they are indeed well localized, with most of the charge concentrated near their positions. The calculation ofQ k (r), displayed in Fig. 6 (b)-(e), shows that they indeed seem to converge to a value close to −p i /q as r increases.
As noted in Sec. III A, the definition of our wavefunction does not forbid us to put the non-Abelian anyons within the R part. The result of exchanging one Abelian and one non-Abelian anyon positions from Fig. 6 (a) is shown in Fig. 6 (f). It can be seen that still the charge is concentrated mostly in their vicinity, and approximately has the expected value −p i /q.
We note that in some cases, even with Abelian anyons only, there is some additional charge modulation at the interface. This is a finite-size effect, whose strength decreases with N y . The detailed analysis of this effect can be found in Appendix B. We note that a similar phenomenon was encountered for Laughlin-Laughlin interfaces [42].

C. Anyon statistics
To check whether the "anyons" we investigate are true anyons, we have to evaluate their statistics. We will consider the processes in which a single mobile anyon l encircles other, static anyons. The effect of anyon braiding is given by where Ψ is a vector of degenerate wavefunctions |Ψ α for all possible values of α, while γ M and γ B are the monodromy and Berry matrices. The monodromy matrix can be evaluated from the analytical continuation of the wavefunctions, while the Berry matrix can be written as γ B = exp iθ B , where the elements of θ B are given by where P is the path of the lth anyon.
To proceed further, we need to show that the conformal blocks are orthogonal if there is more than one. The overlaps can be computed using Monte Carlo, as explained e.g. in Ref. [56]. In Fig. 7 (a) and (b) we plot the overlap between the two conformal blocks for two cases of four non-Abelian anyons and two Abelian ones in four systems, depicted in Fig. 7  overlap decreasing with N is seen, with Ψ ψ |Ψ i of the order 10 −2 for the largest system. This shows that in large systems that can be studied using Monte Carlo the conformal blocks are already close to orthogonality, and we can expect that in the thermodynamic limit the orthogonality will be achieved.
It can be shown [56] that, assuming the conformal blocks are orthogonal or there is just one, we can express the Berry phase (29) solely using the normalization constant

Abelian anyons
In the case when the lth anyon is Abelian, the Berry phase can be computed analytically, under the assumption that the anyons are well-localized (which is supported by the numerical results from Sec. III B). The partial derivative ∂ ∂w l in such a case does not act on the Ising part of the wavefunction. Thus, we can easily gen-eralize the reasoning from Refs. [70]. Knowing the wavefunction coefficients (17), we can evaluate the derivative . (31) Thus, for an anticlockwise path winding at most once around each site and each anyon, the Berry phase is where S is the region of space encircled by the path P .
To deal with the first term of Eq. (32), we note that we are in fact not interested in the phase θ B αβ itself, because it contains both the statistical contribution and the Aharonov-Bohm contribution, arising from the encircled sites. To get rid of the latter, we compute the difference of phases with and without encircled anyons. For simplicity, let us assume that we compute the mutual statistics of anyons l and m, and all anyons other than m lie outside the region S. We consider two cases: θ B,in αβ , when anyon m is inside S, and θ B,out αβ , when it is outside. We have where n k in − n k out are the particle densities in the two cases. Now, the assumption of localized anyons comes into play. If the anyons are localized and far from each other, the density difference is nonzero only in the vicinity of the two locations of anyon m and is w l -independent. Thus, it can be taken out of the integral. Then, applying the residue theorem, we obtain We note that the sum of density differences within region S is equal to the charge of the anyon m, i.e. −p m /q. And thus, the Berry phase vanishes. Hence, the effect of the braiding is given by the monodromy matrix, which is equal to This recovers the Laughlin anyon statistics. The expression is valid in the entire system, i.e. in the parts L and R and for paths crossing the interface, which is consistent with the fact that the Abelian anyons are valid topological excitations of both parts. We also note that the above reasoning is valid even when p m is fractional, i.e. it yields also the mutual statistics of Abelian and non-Abelian anyons, but only if the latter is static. As far as this condition is fulfilled, there is no problem with putting a non-Abelian anyon in part R. The problems arise when it moves, as we will see in the next subsection.

Non-Abelian anyons
In the case where a non-Abelian anyon is mobile, we verify the vanishing of the Berry phase numerically. Following Refs. [56,67], we rely on the fact that the Berry phase vanishes if the normalization constant C (and hence the integrand of (29)) is lattice-periodic in w l as long as anyon l is far away from other anyons. To see this is the case, let us consider a planar system in which the anyon l moves along a rectangular path consisting of four segments P 1 , P 2 , P 3 , P 4 . We consider P 1 and P 3 being parallel to the x direction, with x increasing in the former and decreasing in the latter, and located at y = y 1 and y = y 2 . Similarly P 2 and P 4 are parallel to the y direction, with y increasing in the former and decreasing in the latter, and are located at x = x 1 and x = x 2 . Moreover, we demand that the rectangle has integer dimension in the units of lattice constants, i.e. x 2 − x 1 ∈ Z and y 2 − y 1 ∈ Z. Then, we note that P1 f (w l )dw l = x2 x1 f (x + iy 1 )dx, and P3 f (w l )dw l = x1 x2 f (x + iy 2 )dx. If f (w l ) is latticeperiodic, then f (x + iy 1 ) = f (x + iy 2 ), and the contributions of P 1 and P 3 cancel each other. Similarly, one can show that the contribution of P 4 cancels the contribution of P 2 . Thus, on this special path the statistics are determined by the monodromy. And, since the statistics are a topological property, we expect that they would not change if the path is deformed. The above reasoning can be regarded as a lattice generalization for the continuum argument that the Berry phase vanishes when C is constant.
The lattice periodicity can be demonstrated by calculating the ratio of squared normalization constants in each point of some path, C(w l ) 2 /C 2 0 , where C 0 corresponds to the starting point of the path. The method of caluclating such ratios with Monte Carlo is described e.g. in [56]. Because the system size required for the simulation of a complete braiding process is too large even for the Monte Carlo, we consider a square loop around a single lattice site (see Fig. 8 (a)), which will tell us how the normalization constant changes as we move an anyon by one lattice constant in the x and y directions (or both).
We focus on the case of two non-Abelian anyons, for which there is only one conformal block. We consider a (12+6)×12 planar system and arrange the anyons in the way shown in Fig. 8 (a). Fig. 8 (b) shows the resulting ratio of squared normalization constants while moving the quasielectron around the small square of unit length. It can be seen that the dependence is nearly periodic, with ratio being close to 1 every time the anyon is at a corner of the square. The periodicity is not perfect -there still are some discrepancies larger than the Monte Carlo error, which may be due to the insufficient separation of the quasielectron from the quasihole or from the system edge.
Therefore, we expect that the statistical phase will be determined by the monodromy. We focus on the statistical contribution to the monodromy, i.e. the monodromy after removing the Aharonov-Bohm contribution (which, in case of one anyon moving on a closed loop, can be computed by subtracting the phase with and without the second anyon within the path). This statistical contribution in the case of a single Moore-Read state is well-known. For a single anticlockwise exchange, it is equal to In the case of an interface, there are additional terms involving R sites and non-Abelian anyons, but as long as the braiding path is located in the L part, these terms do not contribute to the monodromy as no R site is encircled by the anyon, thus the result of an exchange is still given by Eq. (36). We can also ask what happens if we put the non-Abelian anyons in the R part. In such a case, the monodromy indicates that the statistics become ill-defined. To see this, we note that now non-Abelian anyons encircle the R sites. The factors (w l − z i ) p l ni for p l = ±1/2 introduce a nontrivial contribution to the monodromy every time n i = 0, i.e. when a filled site is encircled. Thus, the monodromy depends on the path, i.e. it is not statistical. Moreover, since each configuration |n corresponds 9. A schematic depiction of a system with two R islands on the L plane. The filled and empty circles denote particles/holes and anyons, respectively. The positive and negative charge is denoted by "+" and "−", respectively. The following processes are depicted here: the exchange of particles between L and R (the arrows at the bottom of each island), measurement of the charge (the arrow around the right island), and the change of parity (the arrows in the top part of the picture).
to different locations of the filled sites, each coefficient in (14) transforms in a different way and thus the effect of a braiding is no longer a phase. Therefore, we conclude that it is not possible for non-Abelian anyons to cross the interface. This conclusion does not depend on the Berry phase (for some consideration regarding the Berry phase, see Appendix C).
We note that the nontrivial monodromy of an L particle and an R anyon does not influence the boundary condition for particles going around the cylinder (see Appendix D).

IV. MULTIPLE ISLANDS AND TOPOLOGICAL DEGENERACY
So far, we discussed the properties of a single interface. Now, let us consider two disconnected islands of the R type within an L plane, as in Fig. 9. Let us also assume that the processes of exchange of particles through the interfaces are local. That is, if the islands are sufficiently far apart from each other, a pair of particles annihilated from part L should correspond to a creation of two particles in island 1 or two particles in island 2, but not one particle in each island, as shown in Fig. 9. If we fix only the total number of particles N , the Hilbert space contains configurations where the numbers of particles in the first island N R;1 is even and the ones where N R;1 id odd. There is no local process connecting the two types of configurations. Therefore, the Hilbert space fragments into two disconnected subspaces. For each of them, we can define a model ground state wavefunction using Eq. (6) or (17), but reducing the basis only to the given subspace. If k islands are introduced, then the Hilbert space fragments into 2 k−1 subspaces, which is reminiscent of the degeneracy of the Majorana modes. The appearance of topological degeneracy in a similar setting of interfaces forming several disconnected islands was already discussed in Ref. [4].
The parity of N R;1 can be measured by encircling a non-Abelian quasihole around it (see the arrow around the right island in Fig. 9). Then, the term (w l − z i ) p l ni . gives rise to a monodromy phase 0 if N R;1 is even and π if it is odd. We also propose the following procedure to switch the parities at two islands. We create a localized Abelian quasihole-quasielectron pair with p i = ±2, pinned by two pinning potentials somewhere in the L part (see the top part of Fig. 9). Then, we move the quasielectron and a quasihole into the two different islands and release the pinning potentials. While in the L part these excitations were topological, in the R part they become an ordinary particle and a hole, thus changing the parity.
In such a way, the R islands can store quantum information even though the interfaces are gapless. We note that essentially the same mechanism of creating topological degeneracy can be applied to the Laughlin-Laughlin interfaces from Ref. [42], thus connecting our model wavefunctions to earlier results, predicting the appearance of parafermion zero modes at some Laughlin-Laughlin interfaces [13].

V. CONCLUSIONS
In this work, we have constructed model wavefunctions for lattice systems at filling ν = 1/2, in which part of the system is in the fermionic Moore-Read state, and the rest is in a bosonic Laughlin state. We considered the cases in the absence and presence of localized anyonic excitations.
We have seen that the conditions of reflection and scaling invariance lead to different lattice filling factors ν lat , i.e. different particle densities on the two sides of the interface. For wide enough systems, these densities are nearly constant in the bulks of the two parts of the system, and their values are close to what is expected for the respective single quantum Hall states. Also, the constant term γ of the entanglement entropy scaling in the bulks is consistent with the values characterizing the topological order of the respective quantum Hall states.
As for the interface itself, we have found that some charge accumulates in its vicinity, due to the fact that the particle density varies more smoothly than the background charge. We observed a lack of exponential decay of the correlation function in its vicinity, consistent with the prediction that the interface is gapless. We have also shown that the scaling of the entanglement entropy at the interface is approximately linear, although the data are too noisy to determine the coefficient exactly.
We have studied the properties of the Laughlin anyons (which are valid topological excitations of the entire system) and the basic MR non-Abelian anyons. We have found that the quasiparticles of both types are welllocalized and have the expected charge irrespective of their location. However, the statistics become ill-defined if the path of a non-Abelian anyon passes through part R.
Moreover, we argued that for multiple, disconnected islands of the R part within an L system, the particle number parity at each island cannot be changed locally, i.e. it is topologically protected. However, it can be changed and measured by manipulating anyons within the L part.
The presented construction can be modified and extended in several ways. First, after allowing double occupancy, one can consider an interface between a bosonic MR state and a fermionic integer quantum Hall effect. Secondly, one can also consider different fillings ν on both sides, which would allow for all-bosonic or all-fermionic systems, at the price of enforcing different charges of the particles on the two sides. Finally, one can also use other quantum Hall states -e.g. by forming a MR/Halperin interface, studied in [11,12] for the continuous case. Since the entanglement entropy scaling on the L side and at the interface is noisy, here we provide additional results. In Fig. 10, we provide the fit parameters A, γ for all possible positions of the entanglement cut parallel to the interface. Moreover, we also study different sets of data points. Each color corresponds to a fit based on datapoints N y,min , N y,min + 1, . . . , N y,max .
It can be seen that on the L side and at the interface, the results display large fluctuations, sometimes larger than the result itself. Nevertheless, the obtained values seem to be consistent with the γ L = ln(8)/2 prediction. There is also no reason to suggest that at the interface or at its vicinity γ has a different value than on the left.
On the other hand, in the R part, the results are close to ln(2)/2 for every choice of data points. For some systems and some anyon configurations, we observe thatQ i is nonzero also far from the anyon positions. Typically, the deviation is strongest near the interface, similarly to the Laughlin/Laughlin case [42]. In show γ, i.e. minus the intercept. The results in the first (last) two rows were obtained without (with) the inclusion of error bars in the weights for the fit. The colors denote the different sets of data points included, according to the colormap in the inset of (a) (the horizontal and vertical axes correspond to Ny,min, Ny,max, respectively). of anyon configurations for a (4 + 4) × 8 cylinder. And indeed, for four of them some charge accumulates near the interface on both sides with the charge on the left being approximately equal to the charge on the right.
The presence and strength of these charge modulations depend on the anyon configuration, as well as the size of the system (especially N y and (N L mod 8)). In Fig. 12 (a) we plot the excess charge as a function of x, for a series of (6 + 6) × N y systems, with two non-Abelian quasielectrons placed at x = −4 and two non-Abelian quasiholes placed at x = 4. It can be seen that the sign of the density modulation depends on the parity of N y , which corresponds to the two possibilities (N L mod 8) = 0, 4. The magnitude of this density modulation decreases with N y . To quantify it, we define the total excess charge in part L,Q and the anyon charge in part L, 11. TheQi charge density profiles in the presence of various anyon configurations for a (4+4)×8 cylinder: (a), (b) two non-Abelian quasielectrons and one Abelian quasihole with pi = 1, (c) two non-Abelian quasiholes and one Abelian quasielectron with pi = −1, (d) an Abelian quasihole-quasielectron pair with pi = ±2 (equivalent to an R particle-hole pair), (e) an Abelian and a non-Abelian quasihole-quasielectron pair (pi = ±0.5, pi = ±1) (f) an Abelian quasihole-quasielectron pair with pi = ±4 (equivalent to two holes and two particles). The anyon positions are marked with a "×" symbol.
In Fig. 12 (b) we plot the |Q L − Q L;an |, i.e. the magnitude of the excess charge in part L not associated with anyons. It decreases exponentially with N y , suggesting that for infinitely wide cylinders the only excess charge is concentrated in the vicinity of the anyon positions. The excess charge near the interface,Q(x = −0.5), behaves very similarly to the plot in Fig. 12    We have shown that if we put the non-Abelian anyons in the R part, the monodromy in the braiding process becomes ill-defined, and thus they lose their anyonic behaviour. For completeness, here we consider the Berry phase contribution.
Let us first consider a trivial case of N L = 0 and two non-Abelian quasiholes (p 1 = p 2 = 1/2). In this case, the Pfaffian equals 1, and the only non-constant term in the Ising correlator, (w 1 − w 2 ) −1/8 , is canceled by the (w 1 − w 2 ) p1p2 term from the Jastrow part. Thus, following the same approach as in Sec. III C 1, we arrive at θ B,in − θ B,out = 2π p1p2 q = π/4. The case of N L = 2 is a little more complicated, but still tractable analytically under the assumption of localized anyons. Let us again focus on the case of two non-Abelian quasiholes (p 1 = p 2 = 1/2) and one Abelian quasielectron (p 3 = 0). The wavefunction is now given by Ψ(w, n) = 2 − n 1 +n 2 2 A n 1 +n 2 2 In Eq. (C1), A is raised to the power n1+n2 2 , because there are only two options: either two L sites are empty, and Pfaffian equals 1 ( n1+n2 2 = 0), or they are both filled and the Pfaffian is A.
The derivative of Ψ is given by Ψ, (C3) and thus, using Eq. (30), we can write the Berry phase in a process where the first quasihole encircles the second one as θ B = i 4 P ( n 1 + n 2 ) ∂A ∂w 1 dw 1 + where S is the region enclosed by the path, and δ w3 = 0, 1 for w 3 inside or outside S, respectively. Now, we need to subtract the Berry phase for the second quasihole inside and outside S. We assume that the quasielectron is outside S for both cases. The result is If the quasiholes are well-localized and far within the L part, the particle density on the two L sites is the same in the two cases. Therefore, the first term of Eq. C5 vanishes. The second term can be dealt with using the reasoning from Sec. III C 1, resulting in the phase θ B,in − θ B,out = 2π p1p2 q = π/4. In the case of N L > 2, we can attempt to prove the vanishing of the Berry phase numerically, as in Sec. III C 2. However, we observe a lack of periodicity, as seen in Fig.  13. Therefore, our current approach does not allow us to extract the value of the Berry phase. This does not rule out an appearance of periodicity for systems too large to be studied using our Monte Carlo software. The term (w i −z j ) pinj , creating nontrivial monodromy when a non-Abelian anyon encircles a filled R site, has a nontrivial effect also when an R particle encircles a non-Abelian anyon. This can occur on a cylinder. The cylindrical system mapped to a complex plane looks as in Fig. 14 (a). Thus, if an R particle goes around a cylinder, it encircles the entire L part, including the non-Abelian anyons located inside. But this can create a nonzero phase, which seems to lead to a conclusion that the boundary conditions in the R part are determined by the number of anyons in the L part. This would seem strange, as we can consider another, equivalent mapping of a cylinder to the complex plane, where the R part is inside and no L anyons are encircled, as in 14 (b), and thus the boundary conditions for R particles should be determined only by quantities related to the R part.
Let us consider a path K which encircles all N L L sites as well as k R, as well as anyons of total charge −P in /q sites. The only terms which generate nonzero monodromy of a L particle on path K are (z i − z j ) −niηj and (w l − z j ) p l nj . They give rise to a phase φ = 2π (P in − N L η L − kη R ) .
(D1) However, we can write it also as where P out is minus the charge of all anyons outside the path in the units of 1/q (i.e. their charge is −P out /q), and P = P in + P out . Due to the charge neutrality (19), the first term vanishes. The second term depends only on the sites and anyons outside the path. Thus, it is possible to express the phase on path K using only the quantities related to part R. In this way, one can see that the two mappings of the cylinder to the plane yield the same result. It is straightforward to generalize this reasoning to multiple interfaces.