Entanglement negativity between separated regions in quantum critical systems

We study the entanglement between disjoint subregions in quantum critical systems through the lens of the logarithmic negativity. We work with systems in arbitrary dimensions, including conformal field theories and their corresponding lattice Hamiltonians, as well as resonating valence-bond states. At small separations, the logarithmic negativity is big and displays universal behavior, but we show non-perturbatively that it decays faster than any power at large separations. This can already be seen in the minimal setting of single-spin subregions. The corresponding absence of distillable entanglement at large separations generalizes the 1d result, and indicates that quantum critical groundstates do not possess long-range bipartite entanglement, at least for bosons. For systems with fermions, a more suitable definition of the logarithmic negativity exists that takes into account fermion parity, and we show that it decays algebraically. Along the way we obtain general results for the moments of the partially transposed density matrix.


I. INTRODUCTION
Quantum critical groundstates, such as the ddimensional transverse-field Ising model at its transition, possess algebraically decaying correlation functions at large separations.It is thus natural to ask whether this translates to the amount of entanglement between separated subregions.Remarkably, for general bosonic conformal field theories (CFTs) in one spatial dimension, it was shown that this is not the case: the logarithmic negativity (LN) [1], which is a proper entanglement monotone [2] contrary to the entanglement entropy, decays faster than any power [3,4].In higher dimensions, this was shown to hold for the non-interacting scalar field [5,6] as well as for resonating valence-bond (RVB) states [7] on arbitrary lattices (in agreement with continuum results [8]).But can it be true for all highly correlated quantum critical systems in any dimension?
We show that the answer is yes by using (i) exact properties of the 2-spin reduced density matrix in the d-dimensionsal quantum Ising model, and (ii) a nonperturbative expansion of twist operators in general CFTs.However, it is not the end of the story for systems that contain local fermion operators, such as massless Majorana or Dirac fermions.For such fermionic theories another definition of the LN was introduced [9] that takes into account the fact that physical density matrices must respect fermion parity.In fermionic quantum critical systems, we find that this alternative LN, a proper entanglement monotone for fermions [10], decays algebraically.

II. ENTANGLEMENT NEGATIVITY
Consider the groundstate of a quantum critical system (such as a CFT), and A = A 1 A 2 is a subregion made of two disjoint parts, see Fig. 1.Tracing out the complement of A gives the reduced density matrix ρ A .The bosonic (b) [1,11,12] and fermionic (f ) [9] logarithmic negativities read where ∥X∥ 1 = Tr √ XX † is the trace norm.The bosonic partial transpose (PT) T b 1 swaps the entries of the density matrix pertaining to subregion A 1 only, where the subscript denotes the subregion supporting the state.One can use (2) for fermionic states, but this leads to many shortcomings.For one, if the initial state is Gaussian, the partially transposed one is generally not Gaussian [13], making the calculation of E b arduous even in the simplest cases.Second, the LN defined via (2) fails to detect topological phases [9].One reason behind these issues is that T b 1 does not account for the anticommutation relations between fermions.One can define an alternative fermionic partial transpose T f 1 that removes the above shortcomings [9].In the numberoccupation basis, the fermionic partial transpose acts as (2) but it incorporates an additional phase factor that depends on the occupation numbers, see example below and the Supplemental Material (SM).

III. SCALING WITH SEPARATION
Let us begin our discussion with the simplest case: A 1 and A 2 are adjacent.For simplicity, we can take the subregions to be hyper-cubes (such as squares in d = 2, see Fig. 1) of linear size ℓ.In that case, both LNs obey the boundary law where ℓ d−1 is the size of the shared boundary, δ is a UV cutoff such as the lattice spacing, and c 1/2 is a positive coefficient.Equation (3) follows from the fact that when A is the entire space, both LNs reduce to the 1/2-Rényi entropy [1,9], which obeys a boundary law Being a UV property, the scaling (3) continues to hold as the subregions shrink so that they are no longer complementary, but |∂A 1 | is replaced by the shared boundary.For d = 1, one gets a logarithm instead [4,14].
Next, we increase the separation.For finite separations, the LNs are scaling functions of r/ℓ, and we first study what happens when this ratio approaches zero.In the limit r → δ we should recover the boundary law (3), so that for small but finite distances, we replace the cutoff by the separation r to obtain the universal scaling where κ b/f is a positive coefficient determined by the data of the CFT.For instance, for the 2d free relativistic scalar field, we numerically verified the scaling (4) by working with a lattice of coupled quantum harmonic oscillators, and obtained κ b = 0.0226 ± 0.0004.In 1d, the scaling becomes κ b/f log(ℓ/r), where ℓ is the length of each of the two intervals.In a fermionic CFT, we now argue that κ f = κ b .It is simplest to make the argument on the infinite cylinder, where A 1,2 are infinite strips separated by a strip of width r.If κ b ̸ = κ f , as we take the distance r → 0, E b −E f would diverge as 1/r.This contradicts the fact that for A 1 and A 2 adjacent, the LNs are equal.The difference between the LNs comes from matrix elements that hop into the separating strip, and so these contributions decrease as r → 0. For the 1d free Dirac fermion, we can show equality.Indeed, for 1d CFTs E b = c 4 log(ℓ/r) at small r/ℓ [4].The free-fermion central charge is c = 1.
To evaluate E f , we use the moments of ρ T f 1 obtained for even integers [14], which we analytically continue to get E f = 1 4 log(ℓ/r), matching E b .Taking r → δ, we recover the adjacent case [3].
In the opposite limit of large separations, one possibility is power-law decay since the correlation length is infinite, and correlation functions of all local operators decay algebraically.This is what happens for the mutual information I(A 1 , A 2 ) [15], which is a measure of both entanglement and classical correlations [16].However, contrary to the naive expectation, for 1d CFTs, E b decays faster than any power [3].It suggests that sufficiently separated subregions share little (but non-zero [17,18]) distillable entanglement, at least as quantified by E b .Beyond 1d, the LN was observed to decay exponentially for the free scalar CFT in 2d and 3d [5,6].We wish to answer the questions: How general is this rapid suppression of entanglement?Does it hold in fermionic systems when one uses appropriate entanglement measures?

IV. SKELETAL APPROACH TO NEGATIVITY
We begin by considering the simplest case where A 1 and A 2 are two spins (qubits) that belong to the groundstate of a d-dimensional Hamiltonian.This extreme skeletal regime [19,20] faithfully captures the absence of power-law scaling in the thermodynamic limit.Indeed, for large separations, the regions become point-like relative to their separation, which justifies fusing operators: this is the celebrated operator-product expansion (OPE).Working with the minimal Hilbert space is tantamount to keeping the dominant operators in the OPE, and we shall show that it gives non-perturbative access to their contributions.An advantage of the skeletal approach is that we obtain the entire density matrix, from which one can extract entanglement measures other than the LN, which becomes essential as the number of spins in A increases.
By tracing out the complement of A, one obtains the reduced density matrix in the S z -basis {| ↓↓⟩, | ↓↑⟩, | ↑↓⟩, | ↑↑⟩}, where m z = ⟨S z i1 ⟩, and For simplicity, we have assumed that (i) sites i 1 and i 2 are related by symmetry, and (ii) the groundstate is real in the S z -basis.These conditions are not essential, but are respected for instance in the transverse-field Ising model on a general d-dimensional lattice, such as hypercubic, H = − ⟨i,j⟩ S x i S x j − h i S z i .Regarding the first condition, one can work on the infinite lattice or with periodic boundary conditions; for open boundary conditions, the two sites should be related by reflection.The second condition is respected due to the anti-unitary symmetry given by complex conjugation in the S z -basis, T = i K i , which we call time-reversal.One sees that ρ T b 1 A is given by ρ A | ϱ y r →−ϱ y r , coherent with the fact that time-reversal on A 1 yields T 1 S y i1 T 1 = −S y i1 , whereas the x, z-components remain invariant and the transformation does not act on i 2 .We now examine the fate of entanglement at the quantum critical point h c .For r ≫ 1, ϱ z r → m 2 z approaches a finite constant, while ϱ x,y r → 0 since the 1-point functions ⟨S x,y i ⟩ vanish at criticality due to the Ising Z 2 symmetry.Further, we have the strict inequalities 0 < m z < 1/2 since the 1-point function does not vanish because S z i is even under the Z 2 symmetry, and the groundstate is not fully polarized, which only occurs when h → ∞.We can therefore accurately approximate ρ A at sufficiently large separations as a diagonal matrix with strictly positive eigenvalues, implying that the LN vanishes.Thus, there exists a sudden-death distance r sd at which E b (r ⩾ r sd ) = 0.In particular, the decay is not algebraic.Interestingly, this distance is small: r sd = 3 for d = 1 as was shown using (5) [19,21,22].We thus have that A 1 and A 2 rapidly become un-entangled with separation for all d since the positivity of the PT is a necessary [23] and sufficient [24] condition for the separability of two qubits.Furthermore, we show that these properties hold across the entire phase diagram, (SM).
As a second example, we consider nearest-neighbor SU (2)-symmetric RVB states of spins on the square lattice [25].These states describe critical phases of matter but do not correspond to CFTs in the continuum limit.Nonetheless, we can also argue that there exists a suddendeath distance in this case.In the skeletal regime, the two-spin density matrix is also given by Eq. ( 5), with ϱ x r = ϱ y r = ϱ z r because of the SU (2) symmetry.Moreover, the spin-spin correlation functions decay exponentially fast [26][27][28][29], |ϱ z r | ∼ exp (−r/ξ) with ξ = 1.35(1) [28].The eigenvalues of the partially transposed density matrix become all positive for a finite values of ϱ z r sd , corresponding to a finite sudden-death distance r sd .For larger subregions of two or more spins, the density matrices also depend on dimer-dimer correlations, which decay algebraically [28,29].However, the argument still applies and the eigenvalues of the partially transposed density matrix are guaranteed to be positive for small but finite correlations, corresponding to a sudden-death distance.In fact, the argument readily generalizes for RVB states defined on arbitrary lattices, including triangular lattices where the system describes a gapped phase.It is known that the negativity in RVB states on arbitrary lattices is at most exponentially decaying with the separation [7].Our present analysis is consistent but stronger than those results, and further shows that distillable long-range entanglement exactly dies at finite distance in these states.

V. NEGATIVITY MOMENTS
In our analysis below, we establish the rapid decrease of E b in the thermodynamic limit in all CFTs for all d.The CFT analysis employs the replica trick [30], where one needs to evaluate the normalized moments [3,4] with n = 3 being the first integer that yields a nonzero answer.The n → 1 analytic continuation of the even sequence yields the LN E b .For our 2-qubit case, we can exactly evaluate the moments.We find at large separations , where the ellipsis denotes subleading terms.We have the asymptotic relation , where ∆ σ is the scaling dimension of the leading Z 2 -odd operator σ, since S x i is the microscopic ferromagnetic order parameter.S y i is also Z 2 -odd, but it maps to a CFT operator with a higher scaling dimension, consistent with the strong spin anisotropy of the Ising model.More precisely, ϱ y r = A y /r 2(∆σ+1) , so that the corresponding low-energy operator is the temporal component of the descendant of the order parameter, ∂ 0 σ.This can be analytically verified with the exact solution obtained by fermionizing in 1d [31].We thus have that the moments scale as where we have written the expression in a form that holds true in the thermodynamic limit for a large family of CFTs including the d = 1, 2 Ising CFTs, with ℓ being the linear size of A 1 and A 2 .The power law that applies to all CFTs, irrespective of their spectra, is given in the next section.For the d = 1 Ising CFT, we have ∆ σ = 1/8, so that (7) scales as (ℓ/r) 5/2 , in agreement with the continuum result [32].Moreover, we verified that the scaling (7) holds for the free scalar and free Dirac fermion in 1d and 2d.
For the RVB states on the square lattice, the moments scale as E b n (r ≫ ℓ) ∼ exp (−2r/ξ) for the 2-spin case.For larger subsystems, the moments also depend on the dimer-dimer, or 4-spin correlation functions, which decay algebraically.These correlations thus dominate the spinspin ones, and we expect the negativity moments to decay algebraically with the separation.In either cases there is no sudden death for the negativity moments.This argument generalizes to arbitrary lattices.

VI. FIELD-THEORETICAL ANALYSIS
Consider generic (d+1)-dimensional CFTs, where the subsystems A 1 , A 2 are regions of R d with linear size ℓ and separated by a distance r, see Fig. A1A2 having support on A. For r ≫ ℓ, the twist operators can be factorized into a product of a twist operator acting on A 1 and another on A 2 : Tr As in Ref. [33], the large-r limit allows us to further expand each twist operator in terms of local operators in the direct product of the CFT, which is the OPE fusion of extended operators, where k j labels a complete set of operators on the j th copy, and r Ai j is an arbitrary point in A i on that copy.
We have the analogous equation for Σ (n)T1 A1 , but the coefficients C {kj } are replaced by C {kj } .These coefficients depend on the region they pertain to, but we leave this dependence implicit.For scalar operators we have [33] and similarly for C {kj } with a time-reversed manifold C (n)T1 A1 . The argument generalizes to fields with spin, such as fermions and vector fields, see [34].Crucially, since the coefficients depend on manifolds related by time-reversal symmetry, we have C {kj } = ±C {kj } depending on the parity of {k j }.We thus find When n → 1, the C, C coefficients reduce to groundstate one-point functions which vanish on R d+1 .Therefore, E b decays faster than any power, generalizing the d = 1 results [3,4] to CFTs in arbitrary dimensions, including fermionic ones.In contrast, the moments E b n decay algebraically.The leading contribution corresponds to the combination of operators with the lowest sum of conformal dimensions which is odd under time reversal.Also, the non-vanishing terms in the expansion are simply 2(1−n) times the corresponding ones in the expansion of the Rényi mutual information [33].Finally, for regions of different shapes, ℓ 2 /r 2 in (7) becomes ℓ 1 ℓ 2 /r 2 , and the overall prefactor depends on the geometry (ℓ i is the linear size of A i ).Such shape dependence also holds for the fermionic LN discussed below, and its analysis represents an avenue for future research.
As an application of ( 10), let us consider the Ising CFTs in 1d and 2d.The first primary operators are I, σ, ϵ with respective dimensions 0 < ∆ σ < ∆ ϵ .These fields are time-reversal symmetric; the first odd operator is the time component of the descendant ∂ µ σ, with scaling dimension ∆ σ +1.The first contribution in the expansion (10) thus comes from the configuration {σ, ∂ µ σ} where the Z 2 symmetry requires an even number of σ fields, and the moments scale precisely as in (7).

VII. FERMION ENTANGLEMENT
As we did for spin Hamiltonians, we first work in the skeletal regime where each subregion contains a site hosting a single fermion mode.As we explained above, this skeletal regime faithfully yields the long-distance scaling by virtue of the OPE.Our approach will capture the non-perturbative contributions of the leading fermionic and bosonic operators in the spectrum.To obtain subleading contributions, one would increase the number of modes kept in each subregion.The creation operators are c † i1 , c † i2 , and the reduced density matrix ρ A and its fermionic PT read in the occupation-number where n = ⟨n i1 ⟩ is the average occupation, ϱ n r = ⟨n i1 n i2 ⟩ is the occupation 2-point function, g r = ⟨c † i1 c i2 ⟩ is the fermion Green's function, and r = |i 2 − i 1 | the separation (see SM).For ρ A , one sets the elements in parentheses (14 and 41) to zero, while to get ρ T f 1 A we set the g-terms not in parentheses (23 and 32) to zero.The fermionic PT is seen to correspond to the regular PT T b 1 accompanied by an additional multiplication by a phase [9], see also SM.For simplicity, we have assumed that (i) sites i 1 and i 2 are related by symmetry and (ii) the state conserves the fermion number, such that ⟨c i1 c i2 ⟩ vanishes.As before, we work with Hamiltonians on the infinite lattice or with periodic boundary conditions; for open boundary conditions, the two sites should be related by reflection.The above two conditions hold in a great variety of cases, including Hamiltonians that lead to CFTs such as the Su-Schrieffer-Heeger (SSH) model at its topological phase transition [35], or more interestingly the interacting Hubbard model of spinless fermions on the square lattice, Here, V ⩾ 0 is a repulsion between nearest neighbors, while the phase (±1) in the hoppings generates a π-flux on each square

FIG. 2.
Fermionic LN E f for Dirac fermions in 1d (red) and 2d (blue) versus r/ℓ in log-log scale.Data obtained by numerical diagonalization of the correlation matrix for ℓ = 120 (1d) and ℓ = 11 (2d).The solid lines correspond to the prediction (12) with ∆ f = d/2.For d = 1 we find B = 1/4, in agreement with Ref. [14], whereas for d = 2 it is obtained from a fit.plaquette so that the V = 0 band structure has two Dirac cones.At half-filling, V c corresponds to a Gross-Neveu-Yukawa quantum critical point separating a Dirac semimetal from a gapped charge density wave (CDW) phase where the occupations for the two sites in the unit cell differ [36,37].Entanglement was previously studied for V ⩽ V c but for adjacent subregions [38].
Using (11) we can analytically evaluate the fermionic LN.Let us work at half-filling n = 1/2.From the definition (1), we get , where ρ n r = ϱ n r − n2 is the connected correlation function.If the lattice model describes a quantum critical system, we have the following power-law decay at large separations: g r = A f /r 2∆ f and ρ n r = A n /r 2∆n .In turn, this implies the power-law decay of the fermionic LN, where we have written the expression in a form that we conjecture holds true in the thermodynamic limit for a large family of CFTs, with ℓ being the linear size of A 1 and A 2 .In contrast, the bosonic LN of ( 11) is readily seen to vanish for r ⩾ r sd , in agreement with our findings of the previous section.The leading term in (12) is always given by the fermionic scaling dimension ∆ f .Interestingly, the scaling dimension ∆ n can be smaller than ∆ f since n i will generally have overlap with a scalar in the CFT, and the unitary bound for scalars, (d − 1)/2, is lesser than for fermions, d/2.So E f (r ≫ ℓ) is not determined by the operator with the lowest dimension, but rather by the lowest fermionic operator.Indeed, the first terms in the expansion are We observe that no term coming solely from bosonic operators appears.In contrast, for the mutual information between the two sites in the same limit, we get I = 4|g| 2 + 8(ρ n r ) 2 + • • • .The second term is purely bosonic, and would dominate in theories with ∆ n < ∆ f , where it would characterize non-entangling, i.e., classical, correlations.A case where this occurs is the model H f at V c where the number operators n i1 and n i2 are not even under the Z 2 symmetry exchanging the two sites within a unit cell.As such ρ n r will have overlap with the CFT correlator of the order parameter field σ, which has a lower scaling dimension than the fermion [36,37,39].We also find that the moments E f n , defined analogously to (7), scale with the same power as E f in (12), see SM.
In Fig. 2, we test (12) for free Dirac fermion CFTs in d = 1, 2. The Gaussian groundstates allow us to reach large system sizes ℓ d ≫ 1 on the lattice.The scaling agrees perfectly with (12) for the dimensions ∆ f = d/2.

VIII. OUTLOOK
We showed that at large separations the bosonic LN E b decays faster than any power, whereas the fermionic LN E f decays algebraically for fermionic quantum critical systems.There is thus no distillable long-range entanglement in spin (bosonic) systems, in stark contrast with fermionic ones whose entanglement is infinitely more robust.It would be interesting to understand the scaling function for general CFTs, but the answer escapes the present expansion of twist operators.Further, it would be of interest to adapt the twist expansion to show the algebraic decay of the fermionic LN (12) in the continuum.For boson/spin systems, a fundamental question arises: since quantum critical groundstates do not possess such long-range entanglement, what physically relevant quantum states do?A natural setting would be to study nonequilibrium states.A result in this direction was obtained for quantum circuits at the 1d measurement-driven phase transition, where algebraic decay was observed [40,41].Our study motivates charting the space of long-range entangled states.
with the deviation being parametrically small in the separation, Here, λ ∞ i denotes the r = ∞ value of the ρ A eigenvalues λ i in Eq. (S4).It can be seen that the eigenvalues λ i at r = ∞ are strictly positive for 0 , where m x = ⟨S x i ⟩ is the 1-point function of the order parameter.The Heisenberg uncertainty principle constrains the 1-point functions to obey m 2 x + m 2 z ⩽ 1/4.Thus, this eigenvalue only vanishes at the two extreme points in the phase diagram: h = 0 where m x = ±1/2, and h = ∞ where m z = 1/2.For instance, at the critical point the order parameter vanishes m x = 0, while m z < 1/2 since the state is not fully polarised along the transverse field direction.We conclude that there exists a finite critical separation r sd at which λi become strictly positive, and so E b (r ⩾ r sd ) = 0, (S8) implying that the 2-spin state ρ A becomes separable beyond a critical separation.We note this separability also holds at the extremes h = 0, ∞ since the state ρ A is then manifestly separable.Indeed, one can check that while some eigenvalues λi vanish, none become negative.

Bosonic negativity moments
Although the density matrix becomes separable it is not invariant under PT, and we can quantify this non-invariance via the moments (S11)

B. TWO-FERMION DENSITY MATRIX
We compute the matrix element [ρ A ] ij of the two-fermions density matrix (11) in the basis given in the main text.We begin with the off-diagonal elements, and similarly we have [ρ A ] 32 = g * r .On the second line, each of the three last terms vanishes because fermionic operators satisfy the exclusion principle c 2 = (c † ) 2 = 0.Moreover, we used the cyclic property of the trace from the second to the third line.We compute the diagonal elements using similar arguments.We start with [ρ A ] 44 ,

FIG. 1 .
FIG. 1. a) SubregionA is the union of two disjoint parts A1 and A2, here two squares in d = 2. b) Replicated spacetime used to evaluate the moments of the partially transposed density matrix E b n for n = 3.We show the flow of time for a trajectory that traverses the three copies of A1 (top), and A2 (bottom).The partial transposition on A1 time-reverses the trajectory relative to A2.

1 . 1 A1A2n
The trace Tr ρ n A1A2 corresponds to a normalized partition function in theory where n copies of space are sewed along the crosssection corresponding to imaginary time τ = 0 and position x ∈ A 1 A 2 , as illustrated in the bottom-right panel of Fig. 1.For the trace Tr ρ T b , the sewing for x ∈ A 1 connects replicas in the opposite order because of the partial transposition, see the top-right panel of Fig. 1.The resulting n-sheeted spacetimes are denoted C (n) A1A2 and C (n)T1 A1A2 , respectively.The sewing can be implemented by non-local defect/twist operators Σ

3 = 768 m 2 z( 1
the above results for the eigenvalues, we find for large separations at criticalityE b n = b n xy + • • • = b n A x A y r 2(2∆σ+1) + • • • (S10)in agreement with the general CFT scaling obtained in the main text.For n = 3, 4, 5, the b n coefficients read b