Model wavefunctions for interfaces between lattice Laughlin states

We study the interfaces between lattice Laughlin states at different fillings. Using conformal field theory, we derive analytical wavefunctions for the entire system and restrictions on filling factors under which they are well defined. We find a nontrivial form of charge conservation at the interface. Next, using Monte Carlo methods, we evaluate the entanglement entropy at the border, showing the linear scaling and an additional constant correction to the topological entanglement entropy. Furthermore, we construct the wavefunction for quasihole excitations and evaluate their mutual statistics with respect to quasiholes originating at the same or the other side of the interface. We show that these excitations are able to cross the border and stay localized, although their statistics may become ill-defined in such a process. Contrary to most of the previous works on interfaces between topological orders, our approach is microscopic, allowing for a direct simulation of e.g. an anyon crossing the interface. Even though we determine the properties of the wavefunction numerically, the analytical expressions allow us to study systems too large to be simulated by exact diagonalization.


I. INTRODUCTION
One of the most striking characteristics of the topological orders is the bulk-boundary correspondence -the fact that the bulk properties of the given phase can be inferred from its physics at the edge. This is, however, not a one-to-one relation, as a given bulk phase can have several different kinds of edges even if it is terminated by vacuum [1][2][3][4]. The edge is therefore richer than the bulk.
The topological systems can be characterized by their entanglement properties: the entanglement entropy and entanglement spectrum. In two-dimensional gapped systems, the entanglement entropy for a spatial bipartition has a linear scaling (area law), with a constant term (topological entanglement entropy, TEE) characterizing a given topological phase [24,25]. The entanglement spectrum, being the generalization of entanglement entropy, allows to extract the edge properties from the bulk state [26,27]. These properties were initially considered for a single topological phase. Several authors tried to extend this picture and include the entanglement between two different phases [11,12,16,20,[28][29][30]. The area law was shown to hold at the gapped borders between them, although the constant contribution depends not only on the phases involved, but also on the interaction across the interface [16,28]. The correction to the TEE can be interpreted as originating from a symmetry-protected topological phase living at the entanglement cut (which coincides with the interface in this case) [16,31]. Thus, the interfaces themselves exhibit a topological structure. As a result, the interfaces can host parafermionic modes [15,32], which can be understood as a generalization of the Majorana zero modes [33], and can in principle be used in topological quantum computation [32,34,35]. An entanglement cut can also cross through the interface, which allows to study the properties of the interface mode by comparing the entanglement entropy with the predictions for a 1D conformal field theory (CFT) [11][12][13].
The study of interfaces is also very relevant for experiments. The emblematic examples of topological phases of matter are the integer and fractional quantum Hall states (the latter will be abbreviated as FQH states throughout this work). Interfaces between different such states, exhibiting different topological orders, can be created experimentally [36,37]. For example, in an attempt to prove the existence of anyons, an interferometer was created, in which ν = 2/5 and ν = 1/3 FQH states were placed next to each other [37]. The effective theory of the interface, allowing for e/15 quasiparticle charge, was invoked in the theoretical description of this experiment [38]. Another place where the interfaces may occur is the filling ν = 5/2, where a "patched state" of Pfaffian and anti-Pfaffian was proposed to explain the thermal Hall conductivity measurements [22,39,40].
In the case of a single quantum Hall state (without interface), there are two theoretical approaches to study such a system. The first is a "top-down" one, which neglects the microscopic details of the system and focuses on the general features, described by field theories. The second is a "bottom-up" approach focusing on microscopic Hamiltonians and trial wavefunctions. The advantage of the second approach is that it is more detailed. For example, microscopically one can not only calculate the charge and statistics of the fractionalized excitations, but also their size and density profile, as well as simulate their braiding directly [41][42][43][44][45][46][47][48][49][50]. On the other hand, the bottom-up study of strongly interacting systems is challenging, as the strong correlation is not captured by mean-field techniques. As a result, the available system size is limited. In the commonly used exact-diagonalization technique one can typically study less than twenty particles. The situation is even worse in the case of interfaces, because (i) the translational symmetry of the system in the direction of the interface is lost, increasing the computational complexity of the problem, (ii) both halves of the system have to be large enough to capture the physics of respective topological phases, and (iii) the interface should be well isolated from its periodic repetition or the edges of the system. Thus, the works on the interfaces often use the top-down approach. Some of them derive the edge physics from the K matrices of the two topological phases involved [2,3,10,15,16]. Other methods involve e.g. topological symmetry breaking formalism [51,52]. In contrast, bottom-up calculations are very rare. One proposition for overcoming the limitations of exact diagonalization is presented in Refs. [11][12][13], where model wavefunctions for interfaces between continuum Halperin and Laughlin states, as well as between Abelian and non-Abelian bilayer FQH states, were formed using exact matrix product states derived from the CFT.
In this work, we focus on lattice versions of quantum Hall states [41,53,54], and propose a different microscopic approach: we derive the model wavefunctions directly from correlators of the CFT. The properties of such wavefunctions can then be investigated using Monte Carlo methods, which allow for much larger system sizes than exact diagonalization. This technique has proven successful for single FQH states [41,50,[53][54][55][56][57]. We present our method on the example of the interface between two different Laughlin states, proposing wavefunctions for the ground state as well as for localized bulk anyonic excitations.
The paper is organized as follows. In Section II we construct the ground state wavefunction. We study the properties of the wavefunctions analytically, showing a nontrivial charge neutrality condition. Also, we numerically determine the scaling of the entanglement entropy, confirming the presence of the area law at the interface and showing the existence of a correction to the topological entanglement entropy. In Section III we construct a wavefunction for the quasihole excitations. We show that they are well localized and remain so after crossing the interface. We determine the conditions under which their statistics are well-defined and evaluate the statistical phases. Section IV concludes the article.

II. THE WAVEFUNCTIONS WITHOUT ANYONS
We begin with deriving and studying the model ground state wavefunctions of a system with an interface. First (Sec. II A) we introduce the general framework of lattice CFT wavefunctions, and present an example of a single Laughlin state. Next, in Section II B, we obtain the interface wavefunction, and discuss the conditions in which it is well defined. Because these requirements enforce rather low filling, before presenting any concrete example, we study both sides of the interface separately and show that they are topological (Sec. II C). Then, in Section II D we study the entanglement between these two sides numerically.

A. Model wavefunctions from CFT -preliminaries
We start with recalling the construction of the lattice Laughlin wavefunction in planar geometry from Refs. [53,54], which builds on the framework proposed by Moore and Read for continuum FQH states [55]. The general form of the wavefunction is where n = [n 1 , n 2 , . . . , n N ] is the vector of occupation numbers of lattice sites, |n is the corresponding Fockspace basis state, N is the number of sites and C is the normalization constant. The wavefunction can describe either fermions or bosons, but we enforce the hard-core condition for the latter, i.e. n i ∈ {0, 1} in both cases. The square of the wavefunction coefficient can be expressed by a correlator of a conformal field theory, whose operator product expansion can be written as where |0 is a vacuum of this CFT, z i = x i + iy i is a coordinate of a lattice site i,z i is its complex conjugate and V (n i , z i ,z i ) is a vertex operator defined by Here, γ i is a function of the occupation of a lattice site i, and φ(z i ,z i ) is a free bosonic field. A precise form of γ i depends on the given state and will be given later. Using standard methods of the conformal field theory [58], we arrive at the following expression for the unnormalized wavefunction Here, δ γ = δ( N i=1 γ i (n i )) and χ n is an unknown phase factor, which does not affect the topological properties, so we set χ n = 1.
FIG. 1. The interface considered in this work. Blue and red colors correspond to two Laughlin fillings ν = 1/qL, ν = 1/qR, respectively, with spheres denoting the lattice sites. The green rhombus denotes the unit cell of the kagome lattice. The black plane is an example of the entanglement cut, dividing the cylinder into two cylindrical subsystems A, B.
We will work in the cylinder geometry rather than the planar one. Throughout this work, we will assume that the direction y is a periodic one. Let L be the circumference of the cylinder. Given a set of coordinates {ζ 1 , ζ 2 , . . . ζ N } on a cylinder, we relate them to the plane coordinates as z i = e 2πiζi/L . The resulting z i values can be substituted to the expression (4) and a wavefuction on a cylinder can be evaluated (see e.g. Ref. [54]).
Let us now consider an example with γ i (n i ) = qni−ηi √ q , q ∈ Z. Disregarding some factors influencing only the normalization or the gauge, we obtain a lattice analog of the Laughlin wavefunction at filling ν = 1/q, for planar geometry given by [53,54] n i is the total number of particles, and N φ = N i=1 η i is the number of magnetic flux quanta passing through the system. Thus, δ n enforces the charge neutrality (the background charge is included in the vertex operators describing sites, in contrast to the continuum case, where an additional vertex operator for background charge has to be added [55]). Because in general N φ = N , the wavefunction (5) can be characterized by two filling factors -the "Laughlin filling" ν = M/N φ = 1/q, defined as the number of particles per magnetic flux quantum, determining the topological class of the wavefunction, and the "lattice filling" ν lat = M/N , defined as the number of particles per site, controlled by η i . By tuning η i one can interpolate between continuum and lattice states [54]. One can also use CFT to derive a Hamiltonian for which (5) is a ground state [53].
B. Analytical wavefunction for the interface Now, we consider a system consisting of distinct left and right parts L and R, with the respective numbers of sites N L , N R and N = N L + N R in total. We denote the coordinates of the sites as z i = z i;L for i = 1, . . . , N L and z i = z i−N L ;R for i = N L + 1, . . . , N (analogous indices will be used for their occupation numbers). We choose the vertex operators in such a way that the first N L of them describe a Laughlin state with filling ν = 1/q L and constant η i = η L , while the next N R correspond to a similar state with filling ν = 1/q R and constant η i = η R . That is, The result is the following expression for the wavefunction coefficients where ψ L (n L ), ψ R (n R ) are the Laughlin wavefunctions (5) at the respective side of the interface (disregarding the charge neutrality), while ψ LR (n) describes the cross factors The charge neutrality is enforced by where M I = i n i;I , N φ;I = N I η I , I ∈ {L, R}. In this work, we focus only on the wavefunction, without considering the Hamiltonian generating it. However, we stress that each side of the interface can be described by model Hamiltonians with ψ L (n L ) and ψ R (n R ) as exact ground states [53].
The wavefunction (6) has several nontrivial properties. First of all, the mutual statistics of L and R particles, determined by q LR , are in general fractional. We want them to be bosonic or fermionic, so we must enforce q LR ∈ Z. Without loss of generality, we choose the Laughlin filling to be larger on the left than on the right, and set then a can be any multiple of 1/b. Interestingly, in such a case Eq. (12) coincides with a condition for a gappable Particles Quasiholes ... ...

FIG. 2.
Charge conservation in our system in the a = 2 case. A single R-type particle (red filled circle) has the same charge as a L-type particles (blue filled circles). On the other hand, an L-type quasihole (blue empty circle) has the same charge as e.g. a R-type quasiholes with the same p or one R type quasihole with a times larger p. In the latter case the two kinds of quasiholes are exactly the same object, which is signified by the "=" sign. Note that while particles are confined to their "parent" part, the quasiholes can be located anywhere. This is emphasized in the figure by placing the circles representing the particles, but not the quasiholes, on the background of respective color.
interface from Refs. [15,28]. In this article, we consider only cases with q L being prime, and a ∈ N. Gapped interfaces fulfilling (12) with natural a were considered in Refs. [10,16] using the top-down approach.
Secondly, the charge neutrality relation (11) is rather unusual. Using Eq. (12), we obtain the condition for nonzero δ n L ,n R : which means that the particles on the right have a times more charge than particles on the left. That is, when an R-type particle crosses the interface, it emerges as a Ltype particles (this rule is illustrated in Fig. 2, along with an analogous one for quasiholes, which will be derived in Sec. III A). Also, it is possible to have bosons on one side and fermions on the other choosing odd q L and even a. Due to this nontrivial behavior, it seems that the most natural realization of our wavefunctions is in spin systems, with spin flips representing the particles.
In general, the wavefunction coefficients (7) are not invariant under the scaling of coordinates z → bz, b ∈ C, in contrast to the Laughlin wavefunction (5), where the scale is arbitrary. This invariance can be restored by setting the lattice filling to be ν lat = 1/2 at both sides, which can be done by adjusting η L,R = q L,R /2. We will enforce this condition throughout this work.
The analytical results in this work are valid for all interfaces with integer a. In addition, two examples will be studied numerically: q L = 1, a = 2 and q L = 2, a = 2. The former is the interface between a fermionic ν = 1 integer quantum Hall state and a ν = 1/4 bosonic Laughlin state, the latter describes two bosonic Laughlin states with ν = 1/2 and ν = 1/8. A as a function of the cylinder circumference Ny, with the cut at the middle of the sample. The corresponding cylinder length is Nx = 2 Ny/2 , with denoting the ceiling function. Lines denote linear fits. Only the data points denoted by filled symbols are taken into account in the fitting procedure, because the data points at small Ny are more strongly affected by finite size effects. The regression is weighted based on the Monte Carlo uncertainties. The colorful ticks on the y axis denote the theoretical values ln(q)/2 of the TEE.

C. Single Laughlin wavefunctions -numerical calculations
At η = q/2 and low filling (q > 4), the lattice Laughlin states on the square lattice develop long-range antiferromagnetic correlations which destroy the topological order [54]. To prevent this from happening, we need to work on a frustrated lattice, on which a Néel ordering is impossible. We choose the kagome lattice, which, according to Ref. [54], has shorter correlation length than the triangular lattice, and thus smaller finite-size effects. Unless noted otherwise, throughout this work we consider systems on a cylinder, with (N xL + N xR ) × N y unit cells, as shown in Fig. 1. For concreteness, we set the lattice constant (i.e. the length of one of the edges of the green rhombus in Fig. 1) to 1, i.e. the nearest-neighbor distance to 0.5, although the wavefunction would not change if the coordinates are rescaled.
To show that we have a topological state on both sides, we first study single Laughlin states (q L = q R = q, N xL + N xR = N x ) before proceeding to the interfaces. The expectation value of any operatorÔ that is diagonal in the occupation number basisÔ = n O n |n n| can be written as It can be sampled using Metropolis Monte Carlo, treating |ψ(n)| 2 as the probability distribution. In such a way, we can evaluate the density-density correlation function The results for ν = 1/4 and ν = 1/8 are shown in Fig 3 (a). Although the correlation shows some antiferromagnetic-like behavior at short distances, its absolute value is decaying exponentially, showing the lack of antiferromagnetic ordering.
Next, we study the entanglement entropy. We divide the system into two subsystems A, B with a cut along the periodic direction of the cylinder (see Fig. 1 and consider a special case with just one type of Laughlin states). We choose the Rényi entropy of order 2, S where ρ A is the reduced density matrix of subsystem A. It can be calculated using the Monte Carlo method and the replica trick [59,60]. We consider two copies of the system, and write where m A , m B denote the occupation numbers within the respective subsystems for the first copy of the system, and n A , n B , analogically, for the second copy. If the total charge changes after swapping m B → n B , n B → m B , then the charge neutrality enforces ψ(m A , n B ) = ψ(n A , m B ) = 0. Eq. (16) can be evaluated numerically using the Metropolis Monte Carlo method, treating |ψ(m A , m B )| 2 |ψ(n A , n B )| 2 as the probability distribution in importance sampling. The entanglement entropy as a function of cylinder circumference is shown in Fig. 3 (b). The results, in general, show the adherence to the area law, where A is a nonuniversal coefficient and γ is the topological entanglement entropy. To obtain γ, we perform a linear fit with weights based on the Monte Carlo errors. For q = 1, 2, 4, after excluding several data points for low circumferences, which are influenced by finite-size effects, we obtain γ close to the theoretical prediction ln(q)/2. From the fits we get γ = 0.001 ± 0.007, γ = 0.334 ± 0.005, γ = 0.70 ± 0.03, close to ln(1)/2 = 0, ln(2)/2 ≈ 0.346, ln(4)/2 ≈ 0.69 for q = 1, 2, 4, respectively (the errors here are the uncertainties of the fit only). For q = 8 the situation is more complicated. The fit presented in Fig.  3 (b) yields γ = 1.10 ± 0.07, a relatively good match with the theoretical value ln(8)/2 ≈ 1.04. However, the data points from the Monte Carlo calculation exhibit some oscillations around the linear trend. This leads to a strong dependence of the fitted γ on the included data points. For example, using only N y > 3 we obtain γ = 0.81±0.05, which is further away from the expected value. Nevertheless, even though we cannot determine the value of γ accurately, the obtained values indicate that it is nonzero, and thus that the state is topological.

D. Entanglement at the interface
Next, we study the entanglement in the whole interface wavefunction. We again cut it in the y direction and investigate the scaling of the Rényi entropy as a function of the x position of the cut. Fig. 4 (a) shows the results for the q L = 1, a = 2 case for three positions of the cut: in the middle of the L or R region, and at the interface. We observe that the area law is fulfilled at all the cuts, even at the interface. In the middle of the L and R regions, the obtained value of the TEE is close to the theoretical results for single Laughlin wavefunctions. At the interface, the linear fit weighted by Monte Carlo errors yields the TEE γ LR = 0.91 ± 0.03, with the error again coming from the regression procedure. This is larger than the TEEs of the L and R parts.
Similar results are obtained for the q L = 2, a = 2 interface ( Fig. 4 (b)), although here the picture is more distorted. For the R part, the fitted γ again depends on the choice of the included data points as for the case of a single Laughlin state, probably due to finite-size effects: for N y > 2 we obtain γ R = 1.10 ± 0.06, but for N y > 3 we get γ R = 1.2 ± 0.2. Also the scaling of entropy at the interface departs further from the linear dependence than in the a = 2, q L = 1 case. Nevertheless, performing the linear fit for N y > 5, we obtain the interface TEE to be γ LR = 1.21 ± 0.03, slightly larger than the theoretical values for the L and R parts.
Let us now look closer at the dependence of the entropy on the position of the cut. Fig. 5 (a)  scaling of the entropy at all possible positions of the cut for q L = 1, a = 2. One can see that at the interface the entropy drops for thin cylinders. Fig. 5

(b) and (c)
show the fitting parameters A and γ as a function of the position of the cut. Clearly, we observe a spike in γ exactly at the interface. A similar situation is seen for q L = 2, a = 2 ( Fig. 6 (a)). However, here the spike is not located precisely at the interface, but at the nearestneighboring column of sites on the right of it. The TEE at the spike exceeds 1.5, so it is much larger than the theoretical values for the L and R parts. While the exact shape of the γ vs. x and A vs. x curves in the R part depends on the number of included data points, the spike is robust to such variations.
In summary, these results indicate that the area law still holds at the interface, as predicted in Refs. [16,28] for gapped interfaces and observed in Refs. [11,12] for a gapless one. Contrary to Refs. [11,12], we observe a spike of the TEE at the interface. Such an increase of TEE is predicted for gapped interfaces which are not fully transparent to anyons [28]. In our case, we cannot determine if there is an energy gap, as we did not construct the Hamiltonian for the entire system. However, in Section III C we will show that the statistics of some quasiholes become ill-defined when they cross the interface, suggesting that the second requirement of Ref. [28] is fulfilled. The increased TEE at the boundary between two topological orders may also be a signature of a symmetry-protected topological phase existing at the interface [16].

III. ANYONIC EXCITATIONS
Since both sides of the interface are topologically ordered, they have fractionalized excitations. In this section, we derive the wavefunction for these quasiparticles and study their properties.

A. Wavefunction with quasiholes
The model states for systems with quasiholes can be achieved by inserting further vertex operators into the correlator (2), each one corresponding to one quasihole [54,55]. Here, there will be two types of such operators, corresponding to two types of quasiholes: the left and right ones. The wavefunction coefficients are given by the following correlator (18) where Q I is the number of quasiholes of the given type (I ∈ {L, R}), V (n i , z i ,z i ) is of the same form as in Section II A, and V aI (p i;I , w i;I ,w i;I ) = : exp(i p i;I √ q I φ(w i;I ,w i;I )) : (19) are the quasihole vertex operators, with w i;I being the positions of quasiholes and p i;I being integers describing their charges (a basic quasihole has p i;I = 1). The quasihole coordinates w i;I are external parameters of the wavefunction. The quasiholes can be located anywhere on the plane, but in this work we will put them in the middle of the smallest triangles of the kagome lattice.
Since the new vertex operators have a form analogous to (3) (only with quasihole positions instead of particle ones), we can repeat the reasoning from Sec. II A, and obtain the wavefunction ψ(n) ∝ δ n L ,n R ψ L (n L , w L )ψ R (n R , w R )ψ LR (n, w) (20) where w is the collective label for all quasihole positions, and w I , I ∈ {L, R}, contains all the positions of the quasiholes of a given type. The wavefunction parts are given by (24) where all the terms not dependent on the particle positions were absorbed into the normalization.
We note that similarly to the particles, the different types of quasiholes have different charges: an L-type quasihole can be replaced for example by a R-type ones with the same p or one R-type one with a times larger p. In fact, a p i;R R-type quasihole is fully equivalent to a p j;L = p i,R /a L-type one provided that p i;R is divisible by a (one can verify that both are described by the same vertex operator). The relations between the quasiholes of different types are illustrated in Fig. 2.
The position of the quasiholes are not restricted to the L/R part of the system. However, if an R-type quasihole is not a valid topological excitation for the Laughlin filling ν = 1/q L , its statistics will become ill-defined within the L part, as we will show in Sec. III C.

B. Density profile and charge of the quasiholes
In the presence of a finite correlation length indicated by Fig. 3 (a), the quasiholes should be well localized. The The density profile and charge of a p1;R = 2 Rtype quasihole for a qL = 1, a = 2 interface in a system of size (6 + 5) × 5. The rows correspond to different positions of the quasihole: in the R part (top), at the interface (middle) or in the L part (bottom). The columns show different quantities: the deviation of the density distribution ni from half-filling (left) and the excess charge Q1;R as a function of distance r from the quasihole position (right). The horizontal dashed lines are located at Q1;R = −1. Note that the distance between nearest neighbors is r = 0.5. Monte Carlo calculations of particle density show that this is indeed the case. As an example, let us consider a p 1;R = 2 R-type quasihole in a q L = 1, a = 2 system. First, we place the quasihole in its "parent" part of the system, i.e. to the right of the interface (7 (a)). The deviation from n i = 1/2 is significant only near the quasihole position, as expected for a single FQH state. We define the excess charge within radius r from the kth I-type quasihole as where θ is the Heaviside step function. Here, we assumed that particles in part L have unit charge, and thus R-type particles have charge a. The plot of the excess charge as a function of r for the considered situation is shown in Fig. 7(b). For large r it approaches −1, i.e. its modulus is half the charge of an R particle, as expected. The quasihole is well localized also when it crosses the interface, or even when it is located precisely at the border, as seen in Figs. 7(c) and 7(e). The corresponding excess charge plots (Figs. 7(d) and 7(f)) show that while the density profile of the quasihole changes, its total charge stays the same. Note that we can also interpret the charge −1 as a lack of one L-type particle, i.e. a p 1;L = 1 L-type hole, in accordance with the charge conservation rule (24).
We can also place quasiholes in the part of the system where they are not valid topological excitaitions of a corresponding Laughlin state. An example for q L = 1, a = 2 is shown in Fig. 8. In Fig. 8 (a), two p k;R = 1 Rtype quasiholes (each having half the charge of a basic L-type hole) are placed in the R part, while in 8 (b) one of them is moved to the L part. The excess charge concentrated close to the quasihole position is −0.5, independently from this position (see Fig. 8 (c),(d)).
In some cases, fluctuations of charge density occur near the interface after a quasihole is moved across it, provided that the circumference of the cylinder is small. An example is Fig. 8 (b) where we observe that some charge has built up on the rightmost sites of the L part, and also some depletion of the charge was created on the leftmost sites of the R part. We note that, although this example above considers an R-type quasihole which is not a valid anyon at ν = 1/q L , this is not a rule. For example, similar behavior was encountered in a q L = 2, a = 2 system of size (5 + 6) × 3, with one p 1;L = 1 L-type quasihole (being a valid topological excitation of both sides) placed within the R part. Investigating several systems of different sizes and with different quasihole positions, we observed that no charge buildup appeared when the two sides of the interface fulfilled the charge neutrality rules of respective Laughlin states separately, as well as all the cases which can be achieved from this one by moving an equivalent of one L-type particle across the interface (as in Fig. 7, where the quasihole has minus the charge of one L-type particle). In all the other systems we studied, charge fluctuations occurred at the interface.
We stress that the quasihole does not "leave behind" some of its charge at the interface when crossing it, as a correct quasihole charge is observed in the vicinity of quasihole positions in Fig. 8 (b). Instead, the charge buildup probably comes from the fact that placing the quasihole to the left of the interface results in pushing some of the charge from the L part towards the R part, and some of this charge does not cross the interface.
The charge buildup does not depend on the length of the system. This can be seen in Fig. 8 (e) depicting the average particle density at given x coordinate, for systems with the same N y = 3 and different N xL = N xR + 1. Nevertheless, we observe that this effect decreases with increasing diameter of the cylinder. Fig. 9 (a) shows that the density spike near the interface gets smaller as the circumference of the cylinder grows. This is not only an effect of averaging over more sites, as can be seen in Fig. 9 (b), showing the excess charge as a function of x, (27) whose variation near the interface also decreases with increasing N y . Both the average density and excess charge near the interface seem to tend to their bulk values for wide cylinders (Fig. 9 (c), (d)). Moreover, the total excess charge accumulated on each side of the interface seems to converge to the charges of the quasiholes located in these regions ( Fig. 9 (e)), which means that the only excess charge is concentrated near the quasihole positions. Similar behavior is seen for q L = 2, a = 2, although the convergence is slower. Therefore, we expect that this effect would be negligible in the thermodynamic limit, and hence that for wide cylinders the quasiholes will be localized independently of their positions.

C. Statistics of quasiholes
Under the assumption that the quasiholes are localized, which is supported by the numerical calculations of Sec. III B, we can derive their statistics analytically, following the approach from Ref. [41]. We start from the wavefunction (20) and fix the normalization constant to be real, C = n L ,n R ψ(n L , n R , w L , w R )ψ(n L , n R , w L , w R ) (28) The total phase in the braiding process consists of two contributions: the monodromy and the Berry phase. Let us focus on the monodromy first. Since the wavefunction contains no terms depending on the positions of two quasiholes, the only term that matters is (w i;R − z j;L ) p i;R n j;L /a , for which the way the root is taken has to be defined consistently if the exponent is fractional. For a braiding process of two quasiholes in the R part, w i;R never encircles any L site (unless it goes around the cylinder -see Appendix A). Thus, in such a case (w i;R − z j;L ) p i;R n j;L /a stays in the same branch and no phase arises from this term. On the other hand, if the braiding path contains some L sites, then the contribution of (w i;R − z j;L ) p i;R n j;L /a vanishes only when p i;R is divisible by a, i.e., if the R-type quasihole is a vaild Laughlin anyon of the L side. If not, this term yields a phase 2π/a when encircling a filled L site and 0 when encircling an empty one (which suggests that the mutual statistics between L particles and basic R quasiholes is fractional). As a consequence, the phase depends on the number of encircled L particles, which is not fixed. The statistics are hence not well-defined.
Let us now proceed to the Berry phase. For concreteness, let us first assume that we move an L-type quasihole, whose position is denoted by w 1;L , around another quasihole, which can be of any type. The total Berry phase in the braiding process is given by  To get rid of the Aharonov-Bohm phase, we subtract the phase θ out , obtained when a second quasihole is outside the path of the first one, from the phase θ in obtained when the second quasihole is enclosed by the path. This cancels all the terms not depending on the density, When the quasiholes are well separated, and when there is no charge accumulation on the interface, the density difference occurs only near the two positions of the second quasihole. Moreover, it does not depend on w 1;L , so it can be taken out of the integral. Applying the residue theorem, we get where W I is the set of all I-type sites enclosed by the braiding path. Thus, the statistical Berry phase depends on the charge of the encircled quasihole, which, as we have shown in Sec. III B, is constant and quantized. For two L-type quasiholes, we obtain θ br = −2πp 1;L p 2;L /q L , as for a single Laughlin state. For L and R quasiholes, the statistical Berry phase is θ br = −2πp 1;L p 1;R /q LR . If we repeat the derivation for moving an R-type quasihole, we obtain θ br = −2πp 1;R p 2;R /q R for encircling another R-type quasihole and again θ br = −2πp 1;L p 1;R /q LR for encirlcing an L-type quasihole. Those values are the total statistical phases as long as both quasiholes are valid Laughlin anyons of the part in which they are located. Thus, we conclude that our interface wavefunction correctly reproduces the Laughlin quasihole statistics on each side, while introducing nontrivial statistics between the different types of quasiholes. Note that the statistics do not change when the anyons cross the interface, provided that they are well defined on both sides. Furthermore, the obtained values are another signature of the nontrivial mutual statistics of R-type anyons with respect to L-type particles. An L-type quasihole of charge p i;I = q I is equivalent to the absence of single L-type particle (i.e. a hole). Thus, the statistics of the R-type anyon with respect to the L-type hole is given by θ br = −2πp 1;L q R /q LR = −2πp 1;L /a, which can be fractional (provided that it is well-defined, i.e. both objects are located in the R part). This is a further example of the nontriviality of our interface.

IV. CONCLUSIONS
In this work, we have presented a class of model wavefunctions for interfaces between lattice Laughlin states, including the topological excitations. Our work is similar in spirit to Refs. [11][12][13], which derived wavefunctions for the interfaces between continuum Laughlin and Halperin states, as well as between Halperin and Pfaffian states, with the same starting point (conformal field theory) but a different method (matrix product states). Contrary to these works, we obtain a closed-form solution, which allows us to calculate some properties analytically (such as the quasiparticle statistics), and others using the Monte Carlo method (quasihole density profile, entanglement entropy).
Investigating the ground state, we have found that requiring the mutual statistics of the particles on the left and right of the interface to be bosonic or fermionic puts a constraint on the possible fillings, coinciding with a condition for gapped interfaces given in some works [10,15,16,28]. The charge neutrality condition was found to have rather nontrivial form, requiring the particles on the two sides of the interface to have different charge.
Our results indicate the existence of an entanglement area law at the interface, in alignment with the theoretical predictions for a gapless boundary [16,28] and numerical calculations for a gapless one [11,12]. In contrast to Refs. [11,12], we observe that topological entanglement entropy at the interface has a larger value than on both sides, which may suggest the presence of a symmetry-protected topological phase at the interface [16,31].
We have also introduced the wavefunctions for systems with quasihole excitations and shown that they are well localized. We found that they obey a similar charge conservation rule as for particles (one L-type quasihole being equivalent to a R-type ones). According to it, some quasiholes are valid Laughlin quasiholes of both parts of the system, while some are correct topological excitations only in the R part. We found that all quasiholes can move freely through the interface and stay localized in such a process, although if they are not valid topological excitations of the given part, their statistics become ill-defined. Otherwise, the quasiholes have nontrivial mutual statistics with respect to quasiholes of both the same and different type. Our results suggest also that the R quasiholes have nontrivial statistics with respect to L particles.
The approach taken by us has potential for further development. First, it can be extended to non-Abelian states [61], where the anyon behavior is more complex. Secondly, one can use the CFT to construct parent Hamiltonians for lattice FQH states [53], so the same can be attempted for the interface. Finally, since relatively large system sizes are available, the possibility of capturing some exotic properties, such as the presence of the parafermionic modes [15,16], can be investigated. The behavior of the (w i;R − z j;L ) p i;R n j;L /a term was covered in the main text for the paths not going around the cylinder. What happens if they do? Let us consider moving a p 1;R = 1 R-type quasihole along a closed path which winds around the cylinder once, staying in the R part throughout the process. After mapping the cylinder to the complex plane, the path looks as in Fig. 10 -it encircles the whole L region. The term in question yields a phase 2π/a for each encircled filled L site, i.e. 2πM L /a in total. Now, while M L is not well-defined as the particles can be exchanged with the R part, the exchange can only add or remove a multiple of a L-type particles. Thus, the phase is in fact well-defined and equal to