Entanglement entropy of random partitioning

We study the entanglement entropy of random partitions in one- and two-dimensional critical fermionic systems. In an infinite system we consider a finite, connected (hypercubic) domain of linear extent L, the points of which with probability p belong to the subsystem. The leading contribution to the average entanglement entropy is found to scale with the volume as a(p)LD, where a(p) is a non-universal function, to which there is a logarithmic correction term, b(p)LD−1 ln L. In 1D the prefactor is given by b(p)=c/3f(p), where c is the central charge of the model and f(p) is a universal function. In 2D the prefactor has a different functional form of p below and above the percolation threshold.


I. INTRODUCTION
The entanglement properties of many-body quantum systems are subjects of recent intensive theoretical studies [1][2][3][4] .The entanglement between two partitions A and B of a system being in a pure state |Ψ can be measured by the entanglement entropy 5  Most of the studies are restricted to bipartitions with a smooth, regular boundary between A and B, for example in one dimension (1D) the subsystem A contains the successive sites i = 1, 2, . . .L, and B is represented by the rest of the sites.If the system is gapped the entanglement entropy generally satisfies the so-called area law 3 : S ∼ L D−1 .In one-dimensional critical systems, with algebraically decaying correlations, the area law is supplemented by a logarithmic correction, which for conformally invariant systems is given by [6][7][8][9][10][11] S(L) c where c is the central charge of the conformal algebra.Multidimensional (D > 1) free fermion systems satisfy an area law if the spectrum is gapped 12 or the Fermi surface has high codimension 13 .In the case of a sharp D − 1 dimensional Fermi surface, there is a logarithmic correction to the are law 12,14 , and the entanglement entropy is given by the following expression [15][16][17] where the integral is over a scaled version of the spatial and the Fermi surface, in such a way that the volume of the (scaled) Fermi sea is 1, and n x and n k are unit normals to the real space boundary and the Fermi surface, respectively.
In disordered quantum spin chains the average entanglement entropy at the critical point has a logarithmic size dependence [18][19][20][21][22][23][24] , too, which can be calculated by the strong disorder RG method 25 .In higher dimensional, critical random quantum systems there is an additive logarithmic correction due to corners, the prefactor of which is universal, i.e. independent of the form of disorder 26 .
If the subsystem A is not a singly connected domain, much less (analytical) results are available.Here we mention that if A and B contain the sites of two sublattices the contact points between them scale with the volume of the system and so behaves the entanglement entropy, too 27 .The entanglement entropy of irregular subsystems with non-continuous border is also subject of research in the recent years.General upper and lower bounds have been set for fractal boundary in real space and fractal like Fermi-surface in Ref. 17 .Fractal bipartition in the topologically ordered phase of the toric code with a magnetic field was also investigated in Ref. 28 .
In the present paper we study the entanglement entropy when A is elected by a random partition, which means that points of a given domain belong to A with some probability p.This type of setting has already been used in Ref. 29 , where the low-lying part of the entanglement Hamiltonian of a random partition is calculated for the non-critical Kitaev-chain 30 .Here we consider critical fermionic models, hopping models in 1D and 2D, as well as the critical Kitaev-chain.
The rest of the paper is organised as follows.Models and the methods of calculations are presented in Sec.II.Lower and upper bounds for the entanglement entropy are calculated in Sec.III, while small p and small 1 − p expansions are performed in Sec.IV.This is followed by extensive numerical calculations in Sec.V, for different values of the occupation probability p and the linear size of the (hypercubic) domain, L. We discuss our results in Sec.VI while detailed calculations are put to the Appendices.

II. MODELS AND METHODS
We consider fermionic hopping models with half filling defined by the Hamiltonian in terms of the fermion creation, c † i , and annihilation, c j , operators at lattice sites i and j, respectively and the summation runs over nearest neighbour lattice sites.The lattice is either an infinite chain (1D) or an infinite square lattice (2D), in the latter case the components of the positions are i = (i x , i y ) and j = (j x , j y ).
Having the two-point correlation function, C(i, j) = c † i c j for i, j ∈ A, we can calculate the entanglement entropy of the system as where N A is the dimension of the correlation matrix (number of sites in the subsystem), and ζ i are the eigenvalues of the correlation matrix.As discussed in the introduction we consider finite domains of linear extent, L, (subsequent points in 1D and a square in 2D) the points of which belong to the subsystem A with probability p.We have 2 L (2 L 2 ) different subsystems in 1D (2D), for which the entanglement entropy needs to be averaged, while the average value of N A is pL (pL 2 ) in 1D (2D).
For the hopping model C(i, i) = 1/2, whereas for i = j we have in 1D and in 2D.
For the non-random partition (with p = 1) the bipartite entanglement entropy in 1D is given by Eq.(1) with the central charge c hop = 1 and the constant is c hop 1 = ln(2)/3 + (1 + γ E )/3 − 1/30 ≈ 0.723, where γ E is the Euler constant 10 .In 2D for an L × L subsystem the prefactor of L ln L in Eq.( 2) is given by 2/3, which has been verified by numerical calculations 16,34 .
Our second fermionic model is the critical Kitaev chain defined by the Hamiltonian 30 This model corresponds to the fermionic form of the critical quantum Ising chain, what is obtained after performing the standard Jordan-Wigner transformation 31 .The relevant correlation function for this model is given by 23 For the non-random partition (with p = 1) the entanglement entropy of the Kitaev chain (or the critical quantum Ising chain) is again given by Eq.( 1), with the central charge c Kit = 1/2 and the constant is c Kit 1 = c hop 1 /2 + c Kit /3 ≈ 0.528.Note, however, that if the subsystem A is not a single connected domain, as is the case for random partitions, then the entanglement entropy of the critical Kitaev chain and that of the critical quantum Ising chain is different, due to the non-local nature of the Jordan-Wigner transformation 27 .

III. LOWER AND UPPER BOUNDS FROM PARTICLE NUMBER FLUCTUATIONS
Following standard techniques 32 , we can calculate a lower bound to the entanglement entropy by using the inequality where s(x) is defined in Eq.( 4).Then, for the entropy we obtain where N = N A i=1 c † i c i is the particle number operator in the subsystem.This lower bound was used to prove the behaviour of the entanglement entropy of multidimensional free fermions 14 as well as of fractal-shaped partitions 28 .The particle number fluctuations in the case of a random partition of a uniform probability are presented in Appendix A in one and two dimensions.From these we obtain the lower bounds These bounds are plotted in Fig. 1 and Fig. 6.The leading term of the 2D lower bound in Eq. 13 is proportional to the area of 2D percolation clusters 35 , given by We note, that by shifting the parabola x(1 − x) upwards, one can also obtain upper bounds, for example with a shift of 0.08 it holds as This leads to an upper bound for the prefactor of the volume term.

IV. LIMITING BEHAVIOURS FOR p 1 AND FOR 1 − p 1
The average entanglement entropy can be calculated as a series expansion in p, performed in Appendix B for 1D, leading to the following result up to O(p 3 ) where α = 0.5335 is defined in Eq. (B7).
Similarly in two dimensions, the leading terms are The other limiting case, 1 − p 1 can be treated as follows.Here, the logarithmic corrections approach the clean system's results, but due to dilution a volume term will appear.Let us now concentrate on the infinite subsystem (denoted by B), which -in the limit 1 − p 1 -consists of two half lines (in 1D) or the whole plane without the square (in 2D), as well as some isolated points from the interval (square).By neglecting correlations between isolated points, the correlation matrix of the infinite subsystem becomes block diagonal.One block contains the correlation matrix of the infinite subsystem without the isolated sites, i.e. with p = 1.The other block corresponds to the isolated sites, and has the dimension of the number of isolated points, ÑD = (1 − p)L D , containing 1/2 in its diagonal as From this follows that the leading correction to the entanglement entropy is ) Note, that the series expansion results for the volume term agree with the lower bound for p 1 and 1−p 1.

V. NUMERICAL RESULTS
For finite subsystems of linear size, L, we have calculated the average entanglement entropy numerically.For the averaging process over different samples we have used three different methods.In the direct method a large number of samples are generated at a fixed value of p and the entanglement entropy is calculated for each sample, and this calculation is repeated for several values of p.This method generally gives an accurate average value at a given p, if the number samples is large enough (10 6 . . . 10 8 ) for each p, but comparing the results for different values of p leads to large errors since the samples are different at each p.
In the so called indirect method we use the relation where N = L (N = L 2 ) in 1D (2D), and S n denotes the average entanglement entropy of a subsystem of n sites, where the averaging is performed over all possible partitions.In practice, we have generated a large number of random samples with uniform probability for all values of n.These samples and their entanglement entropy are stored and used to calculate S (p) for different values of p.The advantage of the indirect method is that we need comparatively less samples (∼ 10 6 ) and the numerical derivation with respect to p is more smooth, compared to the direct method.
In our third, replica method we generated for each random subsystem sample one (four) replicas in 1D (2D) and fused them together.By comparing the entanglement entropy of the original and the replicated sample one can cancel the leading volume term and gain direct access to the more interesting, subleading corrections.For the best results, our calculations have generally combined the replica method with the indirect method.

A. 1D hopping model
For the 1D hopping model we used the indirect method to calculate the average entanglement entropy for domain sizes L = 16, 32, . . .1024, with 10 6 realizations in each case.In agreement with the analytical results on the lower bound in Eq.( 12) and the perturbation expansions in Eqs.(15) and (18), the leading contribution scales linearly with the number of sites in the domain, L. This is illustrated in Fig. 1, where the average numerical vale of the entanglement entropy per domain size S /L is plotted together with the analytical results.
The asymptotic behavior of the average entanglement entropy is expected to contain sub-leading terms in the form The prefactor of the volume contribution, a(p), which can be represented by extrapolating the curves in Fig. 1, is close to symmetric, a(p) ≈ a(1 − p).More interesting are the subleading terms in Eq. (20), which are conveniently analysed by the replica method using the difference 2 S (L, p) − S repl (2L, p) = b(p)(ln L − ln 2) + c 1 (p) .
(21) Here S repl (2L, p) denotes the average entanglement entropy in the replicated samples, which are obtained by joining the same sample behind another copy.By comparing results at sizes L and 2L, finite-size estimates are calculated for the prefactor, b(p), and the constant, c 1 (p), which are then extrapolated.These are plotted in Figs. 2  and 3, respectively.The prefactor, b(p), starts quadratically for small p, in agreement with the series expansion in Eq.( 15), while at p = 1 reaches the conformal result: b(1) = c hop /3 = 1/3, see in Eq.(1).The constant, c 1 (p), also appears to start quadratically at small p, while at p = 1 reaches the known result, as quoted below Eq.( 6).Interestingly, we have found an overall quadratic dependence: c 1 (p) ≈ p 2 c 1 (1), as illustrated in the inset of Fig. 3.
We have also studied the distribution of the entanglement entropy, shown in Fig. 4 for different sizes at p = 0.25.These distributions are well represented by Gaussians, as illustrated in terms of scaled distributions in the inset.Similar, Gaussian distributions are observed for other values of p as well, but close to p = 1 there is a cross-over regime, where the volume contribution, ln(2)(1 − p)L, and the logarithmic term, c/3 ln L, compete, see in Eq.( 18).

FIG. 3:
The constant of the average entanglement entropy of the 1D hopping chain (green) and that of the 1D critical Kitaev chain (red), calculated by the replica method, see text.The dashed curves on the main panel corresponds to c1(1)p 2 In the inset the difference c1(p) − c1(1)p 2 are shown for the two models.

B. Critical Kitaev chain
For the critical Kitaev chain we have calculated the entanglement entropy of random partitions, as described in the previous subsection.Here we have used finite domains of size: L = 16, 32, . . ., 512 and the averages are calculated by the indirect method over 10 6 samples.As for the 1D hopping chain, the dominant contribution to the average entanglement entropy is the volume term, as illustrated in Fig. 5 where the average entanglement entropy per domain size is shown for different values of L. The shape of the extrapolated curve is similar to that of the 1D hopping chain in Fig. 1 and it is again approxi- mately symmetric with respect to p → (1 − p).
The subleading correction terms are found to be in the same form as given in Eq. (20).Using the replica method and Eq.( 21) we have calculated estimates for the prefactor of the logarithmic term, b(p), as well as of the constant, c 1 (p), and their extrapolated values are plotted in Fig. 2 and Fig. 3, respectively.Considering the prefactor, b(p), its form is very similar to that found for the 1D hopping chain: their ratio is given by b(p) Kit /b(p) hop ≈ 1/2 = cKit/c hop .This is illustrated in the inset of Fig. 2. As seen in Fig. 3 the constant term, c Kit 1 (p), has also an approximately quadratic dependence: c Kit 1 (p) ≈ p 2 c Kit 1 (1), as illustrated in the inset of Fig. 3.

C. 2D hopping model
Here we consider the hopping model in a square lattice, in which the domain is an L × L square.We have calculated the entanglement entropy of samples having finite subsystems with linear extension L = 8, 12, 16, 24, 32, 48 and 64, while averages are obtained through the indirect method over 10 5 samples.According to the analytical results in Eqs.( 13) and ( 16) the average entanglement entropy is expected to be dominated by the surface term to which the first correction is logarithmic: This is in agreement with our numerical results in Fig. 6, showing the average entanglement entropy per domain surface.For increasing L, the curves approach the prefactor, a(p), which is approximately symmetric, a(p) ≈ a(1 − p).Comparing this figure with the one-dimensional results in Figs.1 and 5 the convergence is here slower, due to considerably smaller linear size of the domains in 2D.The prefactor of the logarithmic term is estimated through the replica method: comparing the (four times) entanglement entropy of each sample of size L, with those composed of four joint identical samples, thus having a linear size 2L Eq.( 23) defines an effective, size-dependent prefactor, b(p, L), which is plotted in Fig. 7.As seen in this figure b(p, L) starts quadratically for small p and becomes approximately linear for larger values of the probability.
To study this behaviour further we have calculated the derivative of b(p, L) with respect to p, which is shown in the first inset of Fig. 7.We note, that in the indirect method the differentiation of Eq.( 19) can be performed at each value of p, which reduces the error of the calculation.Inspecting the behaviour of ∂b(p, L) ∂p we can identify two regions.In the first regime the derivative continuously increases, while in the second regime Eq.( 23).In the first inset the derivative ∂b(p, L) ∂p is shown.
In the second inset the finite-size transition points are plotted as a function of L −3/4 , see the text.The dashed blue line is guide to the eye, the horizontal black dashed line represents the transition point for site percolation.
it becomes approximately constant.In finite subsystems there is an extended cross-over region between the two regimes, which, however, shrinks with increasing L.
We summarize these findings in the conjecture that the change in the behaviour of b(p, L) is related to the percolation transition 36 , which takes place in the random partitioning at a critical value p c = 0.592, if L → ∞.To further check this hypothesis, we have defined finite-size transition points between the two regions as the crossing point, where the linear continuation of the curve starting from p = .5for p > 0.5 reaches the value of the constant measured at p 1.These finite-size transition points are plotted in the second inset of Fig. 7 as a function of L −1/ν , with ν = 4/3 being the correlation-length critical exponent of 2D percolation, governing finite-size effects 36 .Indeed the extrapolated value of the finite-size transition point agrees with p c , within the error of the calculation.We have also checked that the volume term with a(p) shows no sign of a singularity at any value of p.

VI. DISCUSSION
We have studied the entanglement entropy of critical free-fermion models in one and two dimensions, when the sites of the subsystem were taken from a hypercubic domain of linear size L randomly, with probability p.We have investigated the average entanglement entropy by calculating lower bounds, by series expansions and performing extensive numerical calculations.When the entire system has infinite extent, the average entanglement entropy for 0 < p < 1 is found to be dominated by the volume term a(p)L D , which is supplemented by logarithmic corrections as b(p)L D−1 ln L. The volume term is non-universal, which is connected to the fact, that the distribution of the entanglement entropy is Gaussian.On the contrary, the logarithmic correction is found to contain information about the universal, critical characteristics of the system.In 1D, comparing the results of the hopping chain and that of the Kitaev chain the prefactor of the logarithm is found to scale as: b(p) = cf (p), where f (p) is a universal, model independent function and c is the central charge of the critical model.Interestingly, for both models the constant term is obtained in a pure quadratic form: c 1 (p) ≈ p 2 c 1 (1).In 2D, for the hopping model b(p) is shown to change its behaviour at the percolation transition point, p c , where the random subsystem develops an infinite cluster.According to our numerical results the derivative, ∂b(p, L) ∂p , is increasing with p for p < p c , but it saturates to a constant for p > p c .This conjectured behaviour would be interesting to justify independently by physical argumentaarguments, perhaps even with some rigorous method.
Our study can be extended to several further directions as discussed next.

A. Finite environment
We studied the case when the size of the entire system is L tot → ∞.For a finite value of L tot , the results should depend on the ratio L/L tot .In 1D, for non-random partitions with p = 1, the functional form of the entanglement entropy as a function of L/L tot is known from conformal invariance 7 .For the 1D hopping model with random partitions, we have checked that the prefactor b(p) vanishes in the case L/L tot = 1.This is illustrated in Fig. 8, where the ratio is plotted against 1/L at p = 1/2.The slope of the points, which defines b, indeed tends to zero for L → ∞, as seen in the inset of Fig. 8.

B. Position dependent selection probability
Another potential extension of our study is to consider a different type of probability distribution for selecting the points of the subsystem.For the 1D hopping model, we have also checked a position dependent probability where l i = min(i, L − i) is the distance of the point i from the nearest edge of the domain.By varying the decay exponent κ ≥ 0 one can interpolate between the The ratio in Eq.( 24) versus 1/L in the case with L/Ltot = 1 and p = 0.5.In the inset the slope of the curve calculated by two point fit is presented, which defines finitesize estimates for the prefactor b.
non-random partitioning with p i = 1 for κ → ∞ and the uniform probability partitioning with p i = 1/2 for κ = 0.The number of internal contact points between the subsystem and the environment scales as L 0 l −κ dl, which is finite for κ > 1, it scales as ln L for κ = 1 and behaves as L 1−κ for 0 ≤ κ < 1.
We have calculated the average entanglement entropy for different values of κ, ranging between 0.1 and 3.0 and the results are presented in Fig. 9.For κ > 1 the dominant contribution is S = 1/3 ln L as for the non-random partitioning case with κ → ∞.This is due to the fact, that the "volume term", which scales with the number of contact points is now O(1), thus it is subleading.In the borderline case, κ = 1, the size-dependence of S is still logarithmic, however with a different prefactor: S ≈ (1/3 + ln 2) ln L. The increase of the prefactor now is due to the "volume contribution", which scales also logarithmically.Finally, for 0 ≤ κ < 1 the average entanglement entropy scales as: S ∼ L α , with α = 1 − κ, which involves the number of contact points and now it is the dominant contribution.These results are illustrated in the inset of Fig. 9 .

C. Some open problems
It would be interesting to try to determine the function f (p) by a different method, perhaps even analytically.It could also be interesting to check, if the observed universality holds for non-fermionic models, too.In addition, one can consider models with random couplings and pose similar questions treated in this paper.No. KKP-126749 and by the Deutsche Forschungsgemeinschaft through the Cluster of Excellence on Complexity and Topology in Quantum Matter ct.qmat (EXC 2147).This publication was made possible through the support of a grant from the John Templeton Foundation.The opinions expressed in this publication are those of the authors and do not necessarily reflect the views of the John Templeton Foundation.We thank Róbert Juhász, Carsten Timm and Zoltán Zimborás for helpful discussions.
: S = −Tr B ρ B ln ρ B = −Tr A ρ A ln ρ A .Here ρ A = Tr B ρ and ρ B = Tr A ρ are the reduced density matrices of the subsystems A and B, respectively, and ρ = |Ψ Ψ|.

FIG. 1 :
FIG.1: Average entanglement entropy per domain size of the 1D hopping model calculated by the indirect method.The series expansion results for p << 1 and for 1 − p << 1 are presented by red lines, whereas the obtained lower bound is drawn by a dotted line.

FIG. 2 :
FIG.2: Prefactor of the logarithmic term of the average entanglement entropy of the 1D hopping chain (green) and twice the same the 1D critical Kitaev chain (red), calculated by the replica method, see text.The series expansion result for the 1D hopping chain at small p is shown by a dotted line.In the inset the ratio of the two prefactors are shown and the dashed line represents the ratio of the conformal charges: c Kit /c hop = 1/2.

FIG. 4 :FIG. 5 :
FIG. 4: Probability distribution of the entanglement entropy of the 1D hopping chain at p = 0.25 for different sizes: L = 64, 128, 512, 1024, from left to right.In the inset the scaled curves are shown assuming Gaussian behaviour.

FIG. 6 :
FIG.6: Average entanglement entropy per domain volume of the 2D hopping model calculated by the indirect method.The perturbative result for p 1 and 1 − p 1 are presented by red dotted line, whereas the lower bound is drawn by a black dotted line.

4 FIG. 7 :
FIG.7: Effective, size-dependent prefactor of the logarithmic correction term in the average entanglement entropy of the 2D hopping model calculated through the replica method in FIG. 8:The ratio in Eq.(24) versus 1/L in the case with L/Ltot = 1 and p = 0.5.In the inset the slope of the curve calculated by two point fit is presented, which defines finitesize estimates for the prefactor b.

FIG. 9 :
FIG.9: Average entanglement entropy with site dependent probability in one dimension.The logarithm of the average entanglement entropy ln S as a function of ln L, for 1 ≤ κ ≤ 3. The black lines are guides to the eye, the slope of the lower and upper lines are 1/3 and 1/3 + ln(2), respectively.In the inset the α exponent of the power-law increase of the average entanglement entropy is shown as a function of 0.1 ≤ κ ≤ 0.5, see text.