Logarithmic, Fractal and Volume-Law Entanglement in a Kitaev chain with long-range hopping and pairing

Thanks to their prominent collective character, long-range interactions promote information spreading and generate forms of entanglement scaling, which cannot be observed in traditional systems with local interactions. In this work, we study the asymptotic behavior of the entanglement entropy for Kitaev chains with long-range hopping and pairing couplings decaying with a power law of the distance. We provide a fully-fledged analytical and numerical characterization of the asymptotic growth of the ground state entanglement in the large subsystem size limit, finding that the truly non-local nature of the model leads to an extremely rich phenomenology. Most significantly, in the strong long-range regime, we discovered that the system ground state may have a logarithmic, fractal, or volume-law entanglement scaling, depending on the value of the chemical potential and on the strength of the power law decay.


Introduction
In recent years, the quantum community's interest in long-range physics has steadily increased due to the emergence of promising platforms for quantum technological applications: long-range interacting quantum systems. These systems are characterized by coupling energies between pairs of microscopic constituents V i,j that decay as a power law of their distance r = |i − j|, with α > 0 [1,2]. This increased interest is largely due to the systems' stability against external perturbations, which allows for the mitigation of the detrimental effects of dynamically generated excitations [2,3]. An example of the rigidity of long-range interacting platforms against external drivings and of its utility for quantum technological applications is the possibility for such systems to host clean discrete Floquet time crystal phases [4][5][6][7]. Another example is the recently introduced advantage in the finite time performance of quantum heat-engines with a working substance hosting long-range couplings [8]. Moreover, this technological and theoretical interest is also supported from the experimental side by the possibility to implement long-range interacting systems in typical quantum simulation platforms, such as atomic molecular and optical (AMO) systems [9][10][11][12][13]. Interestingly, trapped ions setups allow tuning the power law exponent α, dictating the decay of the interaction energy with distance, from α 0 to α 3 [9].
The most important feature a system should have to be a good candidate for quantum technologies is the capability of hosting highly entangled states in its spectrum. Indeed, this crucial property is the essential ingredient to perform tasks that are classically impossible or very inefficient [14]. More precisely, entanglement is the property that makes quantum computation overtake classical one providing the computational speed-up in quantum algorithms as compared to algorithms based on the processes of classical physics [15]. Moreover, it is crucial for many quantum technological applications such as quantum teleportation [16], quantum cryptography [17] or quantum metrology [18].
A set of key quantities entering the characterization of entanglement is provided by the entanglement Rényi entropies. For their definition, one takes a partition of a given system into two subsystems A and B (the complement of A), determines the reduced density matrix of a subsystem (say, of A) ρ A by tracing out the degrees of freedom of B, and then computes its Rényi entropies: S ν = ln Tr[ρ ν A ]/(1 − ν) [19]. One of the most fundamental properties of entanglement Rényi entropies is their behavior with the size of the subsystem considered. The celebrated area law [20,21] refers to the fact that typically entanglement grows as the boundary of the subsystem considered, i.e., for a system in d dimensions and a subsystem of size L having volume ∼ L d and area ∼ L d−1 , then the entanglement entropy of the subsystem scales as ∼ L d−1 . In particular, the area law has been proven to be satisfied in the ground state of one-dimensional systems with mass gap and short-range couplings when the size of the subsystem is much larger than the correlation length [22]. At a quantum critical point, where the correlation length diverges, the area law is known to be violated by a logarithmic term proportional to the central charge of the conformal field theory (CFT) that describes the low-energy spectrum of the model [23][24][25][26][27][28]. These facts motivated initially the study of this quantity due to its similarity to the black hole entropy [20,29], and have eventually revealed the important role that entanglement plays in high-energy physics [30][31][32][33] as well as in the investigation of condensed matter systems [34][35][36].
The previous discussion changes and becomes more involved for systems with long-range couplings [2,37,38]. Indeed the prominent collective character of such non-local systems promotes entanglement spreading and leads to novel forms of equilibrium and dynamical scaling, which cannot be observed in traditional systems with local interactions [39][40][41][42][43][44]. In particular, the anomalous scaling of entanglement in the presence of long-range couplings has recently attracted great interest in the context of the so-called measurement-induced transitions [45][46][47][48][49][50][51]. In this case, the dynamical generation of entanglement is weakend by the presence of local measures applied randomly during the system evolution. More precisely, if the measurement rate is high enough, the steady state entanglement saturates to an area law value independent of the considered subsystem size, if only nearest neighbor interactions are present [3]. On the other hand, in the presence of long-range couplings, subvolume law scalings [3,[52][53][54][55], also referred to as fractal entanglement phases [56,57], appear.
These interesting dynamical phenomena have no clear equilibrium counterpart showing that their origin is directly related to the presence of long-range interactions. The entanglement properties of the ground state of a fermionic chain with long-range pairing couplings and nearest neighbors hopping amplitudes were fully characterized in Refs. [58][59][60][61][62] which reported standard logarithmic violations of the area law in the weak long-range regime.
Moreover, an anomalous logarithmic growth was found even if the mass gap is not zero, associated to the divergence of unnormalized couplings, in the strong long-range regime characterized by a power law decay exponent smaller than the system dimension. On the other hand, the authors of Refs. [63,64] considered a model of fermions with strong longrange hopping amplitudes and no pairing discovering a volume law entanglement scaling. Moreover, the entanglement properties of the Sachdev-Ye-Kitaev (SYK) model [65,66], i.e. a fully connected fermionic model with random interactions, have been extensively studied [67]. Also in this case, the eigenstates of the SYK Hamiltonian display a volume law entanglement scaling whose coefficient has been computed numerically using exact diagonalization techniques [68,69] and analytically assuming the eigenstate thermalization hypothesis [70] or using a path-integral approach which becomes exact in the large-N limit [71,72]. Finally, also in long-range bosonic systems [73,74] and in fully connected spin systems [75][76][77][78][79] only logarithmic violations of the area law were reported.
Despite the extensive amount of literature on the topic summarized above, none of the considered long-range models display a fractal entanglement scaling at equilibrium unless additional ingredients are added such as modifications of the couplings which violate time translational symmetry or the presence of a fractal Fermi surface [63]. Here, we are going to show that the subvolume law observed in measurement induced transitions [3,[52][53][54][55][56][57] is directly caused by long-range interactions and also appears at equilibrium, provided certain conditions are met.
To prove our claim, we study the ground state entanglement scaling in a prototypical model of fermions with power-law decaying hopping and pairing amplitudes, also known as the long-range Kitaev chain [2,80]. This model is sufficiently simple to allow us to perform analytic calculations but at the same time it turns out to host an extremely rich phenomenology. Using the well-known Fisher-Hartwig expansion [81,82], we were able to analytically determine the leading order dependence of the ground state entanglement on the subsystem size L in the scaling limit of an infinite chain of N → ∞ sites and infinite subsystem L → ∞ with fixed l = L/N , for different values of the available parameters. In particular, we can distinguish two main regimes: the weak long-range regime in which the coupling's power law decaying exponents are larger than the system dimension and the strong long-range regime in which they are smaller. In the former case, the system shows standard logarithmic deviations from the entanglement area law in correspondence with the quantum critical points, however, in the most interesting case of equal longrange hopping and pairing the coefficients in front of these logarithmic divergences show a nontrivial dependence on the power law decay exponent α which is not compatible with the standard scaling predicted by critical conformal field theory [23,24]. On the other hand, in the strong long-range case, the system becomes genuinely non-additive, therefore showing a logarithmic deviation from the area law even away from criticality. Most significantly, when the system chemical potential is zero, no local terms are present in the Hamiltonian (as we will see this simple fact strongly affects the nature of the ground state which becomes highly degenerate) thus resulting into a subvolume law entanglement scaling, S ∼ L 1−2α .
Summarizing, our work correctly reproduces previously known results in different limits, thus bringing several disparate results present in the literature into a coherent picture.
Moreover, we are able to detect a fractal entanglement scaling phase which is entirely due to the non-additive nature of the model and does not need the dynamical setting of measurement induced transitions to be observed.
The paper is organized as follows. In Section 2 we introduce the long-range Kitaev model and we describe its phase diagram. In Section 3 we briefly review the techniques which allow us to study the entanglement scaling of generic quadratic fermionic models (the expert reader may safely skip this part). Finally, Section 4 and 5 are devoted to the detailed characterization of the ground state entanglement scaling of the model in the weak and strong long-range regimes, respectively.

Kitaev chain with long-range couplings
We consider a generic model of spinless fermions hopping across the N sites of a onedimensional chain in the presence of pairing interactions, and with a chemical potential h. Assuming periodic boundary conditions, the system Hamiltonian reads whereĉ † j andĉ j are creation and annihilation operators for fermions at site j, while t r and ∆ r are the hopping and pairing amplitudes, respectively. We choose their dependence on the intersite distance r according to the power laws with the hopping exponent α 1 > 0, the pairing exponent α 2 > 0, and N α = N/2 r=1 r −α the Kac scaling factor [83], which guarantees extensivity of the energy in the case α i < 1, with i = 1, 2. This model, often referred to as long-range Kitaev chain [80], is emerging as a minimal model for the study of the effects of long-range couplings on a quantum system [2]. Indeed, its integrable nature makes it amenable to both analytical and numerical treatment. Moreover, as observed in Refs. [84][85][86], when the pairing and hopping power law decay exponents are equal α 1 = α 2 = α the model can be related to the quantum Ising model. In particular, in the short-range case with α → ∞, the relation becomes exact through the Jordan-Wigner mapping [87].
The quadratic nature of the Hamiltonian (2.1) allows its exact diagonalization in Fourier space via the Bogolyubov transformation where we have introduced the momentum space fermionic operatorŝ where k = 2πn/N , and n is an integer such that −N/2 + 1 ≤ n ≤ N/2 . While the Bogoliubov angles are defined by the conditions tan θ k =∆ k /(h −t k ), where Fourier transforms of the hopping and pairing amplitudes are defined as Hereafter, we set J = ∆ = 1 as the energy scale and work in units of = 1. In terms of the Bogoliubov fermions, the Hamiltonian then takes the diagonal form with the quasiparticle spectrum Since ω k (h) ≥ 0, the ground state corresponds to the Fock space vacuum for the Bogoliubov modes, defined by the conditionγ k |gs = 0, ∀k. When studying the critical properties associated with the spectrum (2.7), we must distinguish two main regimes: the weak long-range regime when α 1 , α 2 > 1, i.e., the power law decay exponents are larger than the system dimensionality, and the strong long-range regime when α 1 , α 2 < 1. In the weak long-range case, the Kac scaling is a constant in the thermodynamic limit: N α>1 → ζ(α), where ζ(α) is the Riemann zeta function. Moreover, when the system size goes to infinity, we can safely perform a continuum limit in the k variable. In particular, Eq. (2.5) may be written as where Li α (z) denotes the polylogarithm function. This leads to a continuum spectrum ω k characterized, at the critical points, by a dispersion relation that depends on α 1 and α 2 . In particular, for α 1 , α 2 > 1, the system possesses two different phases separated by two quantum critical points h c = 1, −1 + 2 1−α 1 , in correspondence of which the dispersion relation becomes gapless near to the critical mode k c = 0, π, respectively [2,88]. The critical modes of the spectrum are shown in Fig. 1a where ω 0(π) (blue(red) lines in the plot) is plotted as a function of h for different values of α 1 = α 2 . The nature of the transition is topological and the two topological phases can be distinguished by the value of the bulk topological invariant [89] Figure 1. a) Critical modes k = 0, π of the quasiparticle spectrum as a function of the chemical potential h for different values of α 1 = α 2 , two critical points emerge at h =t 0 ,t π , where in the Phase diagram of the long-range Kitaev chain in the plane (α 1 , h), for the pairing decay exponent α 2 = α 1 , α 1 is the hopping decay exponent and h is the chemical potential. The topological order parameter is q = −1 in the topological phase (blue shaded region) and q = +1 in the trivial phase (red shaded region). The phase space boundaries correspond to the solid lines h =t 0 and h =t π .
where the Bogoliubov angles are defined as θ k = arctan(∆ k /(h −t k )). Moreover, in the nontrivial phase with w = 1, the ground state is doubly degenerate, and can support Majorana edge modes [90].
In the strong long-range regime 0 < α 1 , α 2 < 1 the scenario is more complicated. Indeed, in this case, the Kac normalization factor N α diverges at large N as N α ≈ N 1−α , and the thermodynamic limit of Eq. (2.5) has to be carefully considered. In particular, as pointed out in Ref. [91], while the Fourier modes variable k = 2πn/N becomes continuous as N → ∞, the hopping and pairing amplitudest k ,∆ k , remain discrete and labeled by the integer n, reading with c α = (1 − α)2 1−α . Therefore, the presence of long-range couplings leads to a discrete spectrum ω k → ω n = 2 (h −t n ) 2 +∆ 2 n also at N → ∞. The persistence of the discrete spectrum in the thermodynamic limit does not allow us to define a continuous theory and hinders the conventional definition of quantum critical points in the Kitaev chain. In particular, the winding number in Eq. (2.9) is ill-defined as a consequence of the discontinuity in the Bogolyubov angle distribution [89]. Yet, the transition can still be characterized by the quantity This quantity has proven to be a good topological invariant in cases in which the winding number turns out to be ill-defined [89,92]. Then, also in the strong long-range regime, the behavior of the order parameter q is still consistent with a change of phase at the critical points h =t 0 ,t π [91]. However, as shown in [93], the bulk boundary correspondence turns out to be weakened by the presence of strong long-range couplings. Consequently, the change of q at the critical points is not guaranteed to be in one-to-one correspondence with the appearance of boundary topological edge states. Nevertheless, we expect bulk properties to remain consistent with a change of phase. Figure 1b shows the model phase diagram as characterized by the value of q = ±1 as a function of the chemical potential h and of the hopping power law decay exponent α 1 . Two quantum critical lines appear when varying the α 1 parameter. In particular, we notice that the location of the critical point corresponding to ω 0 = 0 is fixed to h =t 0 = 1 for any value of α 1 (blue bold line in Fig.  1b). On the contrary, the critical point corresponding to ω π = 0 (red bold line in Fig. 1b) is α 1 dependent with two different behaviors in the weak and strong long-range regimes, in particular in the thermodynamic limit we find Finally, the completely mean-field case with α 1 = α 2 = 0 needs to be treated separately. Indeed, in this case, the spectrum becomes strongly degenerate and this may alter the nature of the ground state. In particular, for completely flat couplings the sums in Eq. (2.5) can be exactly computed and, in the thermodynamic, they read t n (α 1 = 0) = δ n,0 ,∆ n (α 2 = 0) = 1 + (−1) n+1 πn . (2.14) Accordingly, the single-particle spectrum becomes where we have introduced the shortcut notation ω 0 n = ω n (α 1 = 0, α 2 = 0). It follows that an extensive number of single-particle energy levels corresponding to all the even modes become degenerate. In particular, when the chemical potential is zero h = 0 all the even modes become zero modes since at this point we have ω 0 2n (h = 0) = 0, ω 0 2n+1 (h = 0) = 2/|πn| and ω 0 0 (h = 0) = 1. This fact deeply affects the nature of the many-body ground state which is no more given by the Bogoliubov vacuum, on the contrary, it allows for a finite population of Bogoliubov fermions in an extensive number of zero modes. More precisely, the ground state for α 1,2 = 0 and h = 0 is given by a generic superposition of the form where n 0 is the number of fermions occupying the N 0 available zero modes. This ground state is highly degenerate indeed each |n 0 state can be realized in N 0 n 0 ways, leading to the exponential degeneracy As a concluding remark for this section, we stress the importance of the Kac scaling in the stabilization of the topological order in the strong long-range regime. Indeed, had we considered not properly rescaled couplings, the presence of long-range hopping α 1 < 1 would have moved the critical point to h c = O(N 1−α 1 ) → ∞, thus destroying the transition.

Entanglement scaling in free fermionic systems
We consider a bipartition of the fermionic chain described by the Hamiltonian in Eq. (2.1), in two subsystems A and B, where A is a continuous interval of chain sites of length L and B is its complementary set, see Fig. 2. Given the Hilbert spaces H A and H B associated to A and B, respectively, then the total Hilbert space of the system can be written as the tensor product H = H A ⊗ H B . If the total system is in a pure state |ψ , then the reduced density matrix, describing the state of subsystem A(B) is obtained by taking the partial trace with respect to H A(B) : ρ A(B) = Tr A(B) |ψ ψ|. The amount of entanglement between the two subsystems can be characterized by the so-called Rényi entropies of A, defined as where ν ≥ 1. These are known to provide an accurate measure for the entanglement of a bipartite system in a pure state [19]. In particular, the limit ν → 1 of the above expression corresponds to the celebrated Von Neumann or entanglement entropy The main goal of this paper is to study the Rényi entanglement entropy for the ground state of a Hamiltonian of the kind analyzed in the previous Section. In particular, we are interested in determining the dependence of S ν,L (A) on the subsystem size L in the scaling limit N → ∞, L → ∞ with fixed l = L/N and how this is affected by the presence of longrange hopping and pairing couplings in the Hamiltonian. This task may be achieved by taking advantage of the fact, that since the Hamiltonian in Eq. (2.1) is quadratic, then all its eigenstates satisfy the Wick decomposition theorem [26,94]. Accordingly, the reduced density matrix can be obtained from the two-point correlation functions. To achieve this, we introduce the 2N × 2N correlation matrix V, which is a block matrix with each 2 × 2 block defined as follows: where i and j range from 1 to N . Then, it can be shown [26,94] that this is related to the Rényi entropies through the formula It is important to notice that, from the computational point of view, this formula constitutes a dramatic simplification since the problem complexity is reduced from the diagonalization of a reduced density matrix of size 2 L × 2 L to the diagonalization of the correlation matrix (3.3) of size 2L × 2L, thus allowing to reach larger sizes L. From the analytic side, it is useful to write Eq. (3.4) as an integral on the complex plane along a contour C surrounding the eigenvalues v j ∈ [−1, 1] of V. Using Cauchy's residue theorem in order to perform the integral, one gets [95,96] S ν,L (A) = lim where we have introduced the function and the determinant Due to the translational invariance of the Hamiltonian (2.1) and given the choice of subsystem A, which is composed of contiguous sites, we can write the Fourier trasform of the correlation matrix V lj as where we have introduced the two dimensional symbol G k which, as detailed in Appendix A, can be written as where σ a , with a = x, y, z, are the Pauli sigma matrices, I is the 2 × 2 identity, and f k = γ † kγ k are the occupation numbers of the Bogoliubov fermionic modes, which for a generic state satisfy the condition 0 ≤ f k ≤ 1.
Using the techniques introduced in Refs. [26,94] the asymptotic behavior for L → ∞ of the Toeplitz determinant D L (z), entering the expression for the Rényi entropies (3.4), can be determined applying the Szegő-Widom theorem [97,98] and an extension of the Fisher-Hartwig conjecture [81,82] to non-scalar symbols [60,61]. The leading order contributions to the logarithm of D L (z) in the L → ∞ limit then read where the coefficients b p (z) of the logarithmic contribution are associated to the discontinuities of G k . More precisely, if there is a discontinuity at some k = p, this means that then the coefficient corresponding to such discontinuity can be computed as [61] b Inserting Eq. (3.10) into the integral for the Rényi entropy (3.5) one obtains where the coefficient of the logarithmic contribution can be computed as As shown in Section 2, whenever α 1,2 > 0 or α 1 = α 2 = 0 and h = 0, the many-body ground state of the system is the Bogoliubov vacuum with f k = 0 ∀k, therefore we are left with a leading order contribution given by a constant term O(1) corresponding to the standard area law in the one-dimensional case, or a logarithmic contribution which is associated to the discontinuity of the correlation matrix symbol G k . On the other hand in the specific case α 1 = α 2 = 0 and h = 0 the many-body ground state becomes highly degenerate allowing for a finite fermionic population f k = 0 for an extensive number of Bogoliubov modes, i.e., all the even modes. As a consequence, the first term in Eq. (3.13) becomes the leading contribution to the large L entanglement scaling corresponding to a volume law behavior S ν,L (α 1,2 = 0, h = 0) ≈ L. Summarizing, the machinery introduced in this section allows us to compute the leading order contribution to the scaling of Rényi entropies with the subsystem size by simply analyzing the symbol continuity properties in the different regimes.

Weak long-range regime
Let us start with the weak long-range regime corresponding to 1 < α 1 , α 2 < 2. In this case, as we have seen in Section 2, the quasiparticle spectrum is continuous in the thermodynamic limit, and the ground state is always given by the Bogoliubov vacuum with zero fermionic populations f k = 0, ∀k. Accordingly, the first term of the Fisher-Hartwig expansion (3.13) vanishes and then the leading order contribution to the entanglement scaling comes from the logarithmic term associated with the matrix symbol discontinuity.
Within the weak long-range regime, we can distinguish three different cases: α 1 > α 2 , α 1 < α 2 and α 1 = α 2 = α. Therefore, in order to proceed we must identify the location of the jumps of G k and compute the lateral limits in these three different situations. Possible sources of discontinuities for G k are the discontinuities or the zeros of the spectrum ω k (h), which appear at the two quantum critical points h = 1, −1 + 2 1−α 1 , where the spectrum becomes gapless at the soft modes k = 0, π, respectively. More precisely, G k has no discontinuities when h = 1, −1 + 2 1−α 1 , since in this case the lateral limits at the critical modes read This leads to a constant scaling of the entanglement entropy S ν,L = O(1) with the subsystem size when the system is not at quantum criticality and therefore the spectrum is gapped. This is nothing but a manifestation of the standard area law for one-dimensional gapped systems [20,21]. On the other hand, quantum criticality leads to logarithmic deviations from the area law. Let us start from the homogeneous critical point (h = 1), when the spectrum has an α 1,2 dependent dispersion relation (see Appendix C), which leads to the different lateral limits Accordingly, no discontinuity is present when the power law decay of the hopping amplitude is slower than that of the pairing, leading again to a constant entanglement entropy. In the α 1 > α 2 case instead, we have a discontinuity in the symbol, with commuting lateral limits. Inserting the expression for G ± 0 in Eq.     Fig.3, blue squares represents the numerical data, the black solid line represents our analytical prediction for the scaling in the L 1 limit, red dots are obtained from the numerics by subtracting the subleading corrections.
Then, inserting this result into the expression for the entanglement entropy, and performing the integration in Eq. (3.5) we obtain the logarithmic scaling This logarithmic scaling is analogous to the one obtained for a conformal field theory with central charge c = 1/2 [24]. This result is in agreement with previous findings [59,60] concerning the entanglement scaling in a Kitaev chain with long-range paring and nearest neighbors hopping α 1 → ∞, here we show that the same scaling holds also for finite α 1 as long as α 1 > α 2 . Figure 3b shows the numerical check of the scaling behavior of the entanglement entropy S L = S 1,L for α 1 > α 2 and h = 1. We obtain an excellent agreement once the subleading corrections are taken into account. In particular, we need to subtract from the numerical data the finite size corrections of the form where the c i = c i (α 1 , α 2 , h), i = 1, 2, 3, coefficients can be estimated from a fit with the numerical data. The most interesting case corresponds to the condition α 1 = α 2 = α which, as previously stated, is closely related to the long-range interacting quantum Ising chain. Moreover, we notice that in this regime the matrix symbol G k , hosts non-commuting lateral limits as k → 0 ± (see Eq. (4.3)). This leads to the non-trivial dependence of the logarithmic contribution coefficient on α Inserting b 0 (z) in Eq. (3.5) and performing the integration (see Appendix B), we obtain the logarithmic scaling behavior of the Rényi entropy where with z k,ν = i tan(π(2k − 1)/2ν). In particular, for ν = 2, 3, the sum in the previous expression reduces to B 2,α = 2 π 2 arctan 2 cos(απ/2) sin 2 (απ/2) + 1 , (4.10) This analytical scaling of S 2,ν at h = 1 and for α 1 = α 2 = α is compared with the numerical result in Fig. 3c. Also in this case, a good agreement is found once the subleading corrections (4.6) are taken into account.
We note that the expression for the scaling coefficients in Eq.(4.9) is valid only for integers ν > 1. Indeed, in this case ds ν /dz is a meromorphic function with poles located on the imaginary axis. This allows us to evaluate the integral in (3.5) by summing over the residues at these poles (see Appendix B for details on the calculation). On the other hand, for ν = 1, we have that which has two branch cuts from ±(1+ ) to infinity (see Appendix B). Therefore, to evaluate the integral in Eq. (3.5) for ν = 1, we perform the integration along these cuts and take into account the change in the phase of the logarithm when we go around the branch points. This reduces the integral to two real integrals, which we evaluate numerically. In the case where α 1 = α 2 = α and h = 1, the integrand still depends on α even for ν = 1, so we can still expect the coefficient for the logarithmic divergence of the von Neumann entropy S 1,L to have a nontrivial α dependence.
It is important to observe that at variance with the α 1 = α 2 cases, the scaling coefficient B ν,α cannot be written in the form where c is the central charge of some conformal field theory describing the model at the quantum critical point. This observation supports our previous claim that the case α 1 = α 2 is special and, somehow, closer to the one of a strongly interacting system such as the longrange Ising model. Indeed, while the case α 1 = α 2 continues to obey the r.h.s. of Eq. (4.13) and, so, is more likely to be described by a CFT, the case 1 < α 1 = α 2 < 2 goes beyond this description as the scaling of the ground state entanglement at the critical point cannot be related to the universal properties of a conformal field theory. A similar result is expected for the Ising model in a transverse field, where the inclusion of long-range interactions is expected to increase the effective dimension of the model and, so, disrupt any CFT description. Figure 5a shows the coefficients B ν,α for ν = 2, 3 as a function of α ∈ [1, 2], we notice that the value of the logarithmic scaling coefficients starts from zero at α = 1 and then grows with α reaching the short-range value for α = 2. Moreover, Fig. 5b shows the α dependence of the effective central charge defined as c eff (α) = 6νB ν,α /(ν + 1) as a function of α. We notice that, apart from the extrema c eff (1) = 0 and c eff (2) = 1/2, the effective charge also depends on the Rényi entropy order ν, thus confirming the fact that it cannot be considered as the proper central charge of a conformal field theory. These results are in agreement with the findings of Ref. [100], where the breakdown of conformal symmetry in a long-range fermionic chain was established.
Finally, we consider the non-homogeneous critical point h = −1 + 2 1−α 1 . In this case, the power of the dispersion relation near the soft mode k = π is not affected by the presence of long-range couplings (see Appendix C). Accordingly, also the symbol discontinuity is independent of the value of α 1,2 , in particular, we find This leads to a logarithmic contribution coefficient The corresponding scaling of the entanglement entropy is then the one obtained in Eq. (4.5), which is equivalent to the entanglement scaling in the nearest neighbor Kitaev chain, at a quantum critical point characterized by a conformal field theory with central charge c = 1/2. Figure 4 shows the entanglement scaling behavior at the non-homogeneous critical point h = −1 + 2 1−α 1 with α 1 < α 2 (Fig. 4a), α 1 > α 2 (Fig. 4b) and α 1 = α 2 (Fig. 4c). Also in this case a nice agreement with the theoretical prediction in the thermodynamic limit is found once finite size corrections are taken into account.
The results for the entanglement scaling with the subsystem size in at different critical points and for different values of the α 1 , α 2 parameters within the weak long-range regime considered in this section (1 < α 1 , α 2 < 2) are summarized in Table 1.

Strong long-range regime
The situation in the strong long-range regime is more involved. In particular previous studies on fermionic systems with strong long-range pairing interactions [59][60][61] reported logarithmic violations of the entanglement area law even away from criticality. However, in those cases, the noncritical logarithmic scaling of the ground state entanglement was associated with divergences in the long-range couplings due to the fact that no Kac scaling was introduced in the model Hamiltonian. Therefore, one may think such anomalous scalings to be trivially related to the loss of the system extensivity. On the other hand, as shown in Sec. 2, the introduction of a Kac scaling in the Hamiltonian allows us to define a model with strong long-range interaction still preserving the energy extensivity.
In particular, when a Kac scaling is introduced, the coupling divergences for α 1 , α 2 < 1 are canceled, and accordingly also the symbol discontinuity associated with them disappears. However, an infinite number of new nontrivial discontinuities arise due to the fact that the spectrum becomes discrete also in the thermodynamic limit. More precisely, as a consequence of the spectrum discontinuity, the symbol becomes discontinuous for any k = 2πn/N . Indeed, in the thermodynamic limit, G k reads Then it can be labeled by a discrete integer number n, while the k variable becomes continuous. More precisely, any real physical implementation of the model has necessarily a finite size. Therefore, the actual physical meaning of the continuum limit as N → ∞ is that the difference between two consecutive values of k is of order O(N −1 ). However, in the strong long-range case, a difference of order O(N −1 ) in the k variable results in a finite jump of the spectrum ω n which remains discrete even in the thermodynamic limits, thus resulting in a discontinuity of the matrix symbol G k for any k independently of the value of the chemical potential h. In particular, since for any α 1,2 > 0 or α 1 = α 2 = 0 and h = 0 the many-body ground state is still the Bogoliubov vacuum, then the two lateral limits corresponding to a given k ± = 2πn/N, 2π(n + 1)/N can be written as where we have introduced the angles φ n defined by the conditions cos φ n = 2(h −t n )/ω n and sin φ n = −2∆ n /ω n . Then, following the analytic procedure introduced in Section 3, for any value of h, we obtain a logarithmic scaling of the ground state Rényi entropies of the form where the B ν (h) coefficient is a function of ν, α 1 , α 2 and h. Then B ν,α 1 ,α 2 (h) is given by the sum of N contributions corresponding to the N discontinuities of the symbol, reading where, as shown in Appendix B, each contribution reads where |z l | 2 = tan 2 (π(2l − 1)/2ν), with l = 1, . . . , ν and l = (1 + ν)/2. In particular, for ν = 2 the above sum can be written explicitly as As we have already seen in the previous Sections, the most interesting situation is the one with equally long-range hopping and pairing amplitudes, i.e., with α 1 = α 2 = α, while we expect only minor differences to appear when α 1 = α 2 , as long as they are both smaller than the system dimension (here d = 1). Therefore, for the sake of simplicity, we will limit our treatment to the α 1 = α 2 = α case in the following analysis of the strong long-range regime. Figure 6a shows B 2 (h) as a function of the chemical potential h for different values of α 1 = α 2 = α. First of all, we notice that for any values of the chemical potential h = 0 and of α > 0 the scaling coefficient is of order B 2 (h = 0) = O(1), then leading to a logarithmic violation of the area law even away from the quantum critical points. Moreover, two singularities appear at the quantum critical points h =t k ,t π = 1, 0. In particular, we have a discontinuity for h = 1 and a divergence with the subsystem size for h = 0, leading to a subvolume law entanglement scaling.
These facts can be understood as follows. The spectrum is labeled by the discrete index n leading to a finite gap between the ground state and the first excited levels which are associated with discontinuities of the symbol. However, for n 1 all the modes accumulate around ω ∞ = 2|h|. This means that an extensive number of single-particle states is almost degenerate. Consequently, as long as h = 0, we may expect only the first few modes around n = 0 to provide a significant contribution to the symbol discontinuity leading to a coefficient B ν (h = 0) = O(1). Accordingly, we may expect many features of the entanglement scaling coefficients for values of the chemical potential sufficiently far from the h = 0 point, to be qualitatively reproduced by considering a single discontinuity approximation in which only the first discontinuity between the n = 0 and the first two degenerate levels n = ±1 is considered, i.e., . Then, as detailed in Appendix D within this approximation the discontinuity coefficient reads This approximation then allows us to capture the origin of the scaling coefficient discontinuity at h = 1. This originates from the fact that the zero mode gives different contributions at the two sides of the transition, indeed (see Appendix D) The single discontinuity approximation turns out to correctly reproduce the qualitative features as long as the chemical potential h is sufficiently far from h = 0 and for sufficiently large power law decay exponent α > 1/2. On the other hand, this simple approximation is no more accurate as the chemical potential approaches the h = 0 point. Indeed, in the zero chemical potential case ω ∞ = 0, and more precisely ω n ,t n and∆ n approach their asymptotic values differently if we consider the even or the odd modes (see Appendix D for more details). As a consequence, for sufficiently small α, the number of relevant symbol discontinuities grows as a power law of the subsystem size L, leading to a fractal subvolumelaw entanglement scaling. In particular, using the asymptotic expansion of ω n ,t n and∆ n in the n → ∞ limit we can extract the leading order dependence of B ν (h = 0) from L, which, as shown in Appendix D, reads Accordingly, the leading order contribution to the entanglement Rényi entropy of the system ground state at zero chemical potential takes the nontrivial form This analytic result matches the numerics in the large L limit. This is shown in Fig. 6b, where the numerical and analytical results for S 2,L are plotted as a function of ln L and for different values of α. It is important to notice that approaching the thermodynamic limit in the h = 0 case the spectrum becomes increasingly more degenerate approaching the α = 0 case. Then, for each finite N , a large number of states nearly degenerate with the ground state exists, making the estimate of the subleading corrections scaling technically challenging. Finally, as already stated in Sections 2 and 3, the mean-field case with α 1 = α 2 = 0 and h = 0 must be treated separately. Indeed, in this case the ground state degeneracy allows for a finite fermionic population of the even Bogoliubov modes, f n = 0 ∀n(even), this leads to the entanglement scaling In particular the maximal Rényi entropy is reached when f n = 1/2 ∀n(even) This tells us that the Fisher-Hartwig result, obtained as a large subsystem size expansion, actually becomes exact in this maximally entangled case.
The results for the entanglement scaling with the subsystem size for different values of the h and α = α 1 = α 2 parameters within the strong long-range regime considered in this section (0 < α < 1) are summarized in Table 2. Table 2. Summary of entanglement scaling results at different quantum critical points and for various values of α = α 1 = α 2 in the strong long-range regime.

Conclusion and outlooks
In this paper, we have further extended the understanding of the peculiar properties of entanglement in quantum systems featuring long-range interactions. At this scope, we have investigated, as a paradigmatic example, the ground state entanglement scaling of a spinless fermionic chain with long-range hopping and pairing amplitudes. The simplicity of the model and its truly non-additive nature allowed us to unveil an extremely rich and nontrivial phenomenology, which we have fully characterized both numerically and analytically in the different regions of the relevant parameters, i.e., the power law decay exponents of the hopping and pairing couplings α 1 , α 2 and the chemical potential h. In particular, two main regimes may be distinguished: the weak long-range regime with 1 < α 1 , α 2 < 2 and the strong long-range regime with 0 < α 1 , α 2 < 1.
In the weak long-range case, the system quasiparticle spectrum becomes continuous in the thermodynamic limit and the main effect of the non-local couplings is to change the dispersion relation near the gapless critical modes. Accordingly, the standard area law, typical of gapped local Hamiltonians, is satisfied in this regime apart from the logarithmic violations which appear in correspondence of the two quantum critical points located at h = 1, −1+2 1−α 1 . Such logarithmic scaling of the ground state Rényi entropies is related to discontinuities in the symbol of the correlation matrix which is a block Toeplitz matrix. The fact that the contribution to the entanglement scaling of each discontinuity only depends on the value of the symbol [60,61] at each side of the jump, allowed us to exactly compute its coefficients. Most significantly, when the hopping and pairing couplings are equally long-range, i.e., α 1 = α 2 = α, the coefficient in front of the critical logarithmic divergence at h = 1 turns out to have a non-trivial dependence on α (4.9).
Interestingly, the coefficient B ν,α is of non-universal nature, since it originates from the precise form of the spectrum in the proximity of the critical modes, and not only from the dispersion relation power law exponent. As a consequence, the critical entanglement scaling is not compatible with the result obtained from any conformal field theory and our result may be seen as a benchmark of the fact that the presence of long-range couplings explicitly breaks the critical conformal symmetry [100]. These findings demonstrate the peculiarity of the α 1 = α 2 case, whose physics is expected to be, and indeed is, closer to the one of a strongly interacting system such as the quantum Ising model, where long-range couplings are expected to increase the effective dimension and, so, disrupt integrability [101].
Moreover, for α 1 = α 2 , the critical entanglement scaling becomes α independent. In particular, when α 1 > α 2 , i.e., the pairing coupling has a slower decay with respect to the hopping, the entanglement scaling is compatible with that of conformal field theory with central charge c = 1/2. This is in agreement with the results of Ref. [60,61], where a Kitaev chain with long-range pairing and nearest neighbors hopping is considered, the validity of such results is then here extended to any long-range hopping with power law decay exponent α 1 > α 2 . The strong anisotropy between the case of dominating hopping α 1 < α 2 and the case of dominating paring α 1 > α 2 is typical of the long-range Kitaev chain [102].
In the strong long-range regime, the situation is more involved, indeed the quasiparticle spectrum can no more be considered continuous in the thermodynamic limit. Consequently, the matrix symbol of the block Toeplitz correlation matrix formally becomes discontinuous at every point of the spectrum. However, as shown in Section 5, in most situations only a few of such discontinuities truly contribute to the entanglement scaling, leading to a logarithmic dependence on the subsystem size even outside criticality. Also in this case the coefficients of such logarithmic divergence can be computed analytically for different values of the parameters α 1,2 and h.
The most interesting situation turns out to be the zero chemical potential point h = 0 in the strong long-range regime. Indeed, in this case, the coefficient in front of the critical logarithmic entanglement scaling diverges as a power law of the subsystem size, leading to a fractal subvolume-law entanglement scaling. More precisely, we were able to analytically extract the leading entanglement dependence from the subsystem size, which turns out to be of the form S ν,L ≈ L 1−2α ln L, with 0 < α = α 1 = α 2 < 1/2, where S ν,L is any ν-Renyi entropy with ν > 1. Similar sub-volume laws have already been observed in different (more complex) scenarios and, in particular, in the entanglement scaling of measurement induced phase transitions [53], where they arise due to the suppression of entanglement caused by repeated measurements in a long-range systems. Here, this phase emerges naturally in the equilibrium scaling, but it needs stronger interactions to appear with respect to the dynamical case.
Finally, in the completely mean-field case, the system presents an extensive number of degenerate modes with zero energy. These zero modes can be populated also in the many-body ground state whose degeneracy then grows exponentially with the number of zero modes. Consequently, the ground state entanglement shows a volume law behavior proportional to the size of the considered subsystem S ν,L (α = 0) ≈ L.
Our studies evidence that long-range couplings can greatly improve the scaling of entanglement at equilibrium and, therefore, that long-range interacting quantum systems represent the ideal candidate for reliable and robust quantum computation. Nevertheless, such fostered entanglement properties may not persist out-of-equilibrium, since long-range interactions have been shown to suppress the dynamical spread of entanglement in certain systems [40]. For the future, we intend to investigate these issues by performing quantum simulations of the model on actual quantum computers. This demands a careful engineering of the artificial non-local couplings on local quantum devices, a task which we are currently tackling on IBM Quantum devices [103].
The rich phenomenology hosted by the minimal long-range model we considered, already at equilibrium, suggests that many of the intriguing dynamical phenomena which are recently emerging in the quantum community, such as the non-trivial fractal entanglement scalings in the contest of measurement-induced entanglement transitions [3], can be simply ascribed to the presence of sufficiently long-range couplings among the microscopic components of the model, without any need of further complexity in the physical system under consideration. Further work is needed in order to investigate the dynamical properties of entanglement in the Kitaev chain with long-range pairing and hopping couplings subjected to a unitary or a non-unitary (measurement-like) evolution. These interesting problems are beyond the scope of this work and we leave them as an outlook for future projects.
Coarse-grained description for nonequilibrium systems and transport phenomena (CO-NEST) No. 201798CZL. AS and SS acknowledge acknowledge financial support from National Centre for HPC, Big Data and Quantum Computing (CN00000013). Access to the IBM Quantum Computers was obtained through the IBM Quantum Hub at CERN.

A Derivation of the matrix symbol
In this Appendix we provide the details for the derivation of the matrix symbol in Eq. (3.9) of the main text. We start from the definition of the correlation matrix of a stationary state |ψ , then passing to the Fourier basis we obtain Introducing the Bogoliubov transformation we can write the symbol in terms of the Bogoliubov modes as We now compute the expectation value in a stationary state associated to the fermionic populations of the Bogoliubov modes f k = γ † kγ k , so that Finally, inserting this expectation value in Eq. (A.3) and using the definition of the Bogoliubov angles tan θ k =∆ k /(h −t k ) we obtain which is the expression for the matrix symbol used in the main text.

B Coefficients of the Fisher-Hartwig expansion
The general form of the matrix symbol in Eq. (A.3) can be used to compute the different terms in the Fisher-Hartwig expansion of the Rényi entropies for large subsystem size in every situation considered in the main text. For this purpose, it is useful to rewrite G k as where we have introduced the coefficients a k = 1 − (f k + f −k ) and b k = f −k − f k and the angle φ k such that cos φ k = 2(h −t k )/ω k and sin φ k = −2∆ k /ω k . Let us start from the first term of the expansion in Eq. (3.13) this is obtained by first computing the determinant b n + a n cos δϕ n b n + a n b n − a n cos δϕ n b n − a n 1 + ϵ −1 − ϵ Figure 7. Contour of integration and cuts of the integrand in Eq. (B.8). The cuts from ±(1 + ) to the infinity correspond to ds ν (1 + , z)/dz while the cuts inside the contour, [b n − a n , b n − a n cos δφ] and [b n + a n cos δφ, b n + a n , ], are due to the other factor of the integrand.
Then, the contribution to first term in the entanglement scaling coming from each k-mode is obtained from the integral where Cauchy's residue theorem and the expression (3.6) for s ν (x, y) have been used. Finally, summing over all the modes and using the k → −k symmetry we obtain The logarithmic contribution to the entanglement scaling can be computed by considering the discontinuity coefficients. Here, we present their calculation in the general situation in which G k is discontinuous at a generic mode k = 2πn/N . We start from the definition (3.12) of the b k coefficients corresponding to each discontinuity. First of all, we consider the matrix where G ± k = lim p→k ± G p . The eigenvalues µ ± k (z) of this matrix can be written in the form From this expression we compute the coefficient B ν of the contribution of this discontinuity to the logarithmic term of the Rényi entropy. For this purpose we plug b k (z) into the contour integral for S ν,L then, performing an integration by parts, we obtain The integral over the contour C depicted in Fig. 7 can be divided into two integrals along curves enclosing respectively the cuts [b k − a k , b k − a k cos δφ k ] and [b k + a k , b k + a k cos δφ k ], which in turn can be reduced to two real integrals by performing the integration along the cuts taking into account the change in the phase of the logarithm when we go around the branch points b k ± a k and b k ± a k cos δφ k . On the other hand, we notice that for integer ν > 1, ds ν /dz is a meromorphic function with poles located at the points of the imaginary axis [60,61] z l = i tan π(2l − 1) 2ν , l = 1, . . . , ν, l = 1 + ν 2 , (B.9) and that the another factor of the integrand is analytic in the whole region outside the contour C. We can send this contour to infinity and reduce the calculation of B ν to the computation of the corresponding residues. In this way, we obtain the explicit expression This general formula can be specified in the different cases considered in the main text. In particular in weak long-range case, 1 < α 1 , α 2 < 2, the ground state corresponds to the Bogoliubov vacuum, therefore f k = 0, a k = 1 and b k = 0, ∀k. Accordingly, the first term of the expansion vanishes. Moreover the matrix symbol is continuous for generic values of the chemical potential leading to an O(1) entanglement. The only discontinuities arise at the two quantum critical points h = h c = 1, −1 + 2 1−α 1 in correspondence of the critical modes k = k c = 0, π. This leads to a logarithmic scaling with coefficient |z l | 2 + cos 2 (δφ kc /2) − i sin(δφ kc /2) where in the last step we have used the identity arctan(x) = i[ln(i + x) − ln(i − x)]/2 in order to make the expression of the coefficient explicitly real. The value of δφ kc depends on the critical point considered and the relative order of the power law decaying exponents α 1 and α 2 . In particular for h = 1 and k = 0 we find (B.12) Leading to the coefficients On the other hand, for h = −1 + 2 1−α 1 and k = π, δφ π = π independently from the values of α 1 and α 2 . This leads to the scaling coefficient In the strong-long range regime 0 < α 1 , α 2 < 1, the quasiparticle spectrum is discrete also in the thermodynamic limit, this formally leads to an infinite number of discontinuities for any mode k = 2πn/N , which are labeled by the integer n = −N/2, . . . N/2. In particular whenever α 1,2 > 0 or α 1 = α 2 = 0 and h = 0, the many-body ground state is still the Bogoliubov vacuum characterized by f k = 0, ∀k. Accordingly, the matrix symbol in the thermodynamic limit takes the form in Eq. (5.1). The coefficients of the logarithmic scaling is then given by the sum of the contributions coming from all the discontinuity, i.e., sin(δφ n /2) |z l | 2 + cos 2 (δφ n /2) with δφ n = φ n+1 − φ n . Finally, in the mean-field case α 1 = α 2 = 0 with zero chemical potential h = 0 the quasiparticle spectrum develops an extensive number of degenerate zero modes ω n = 0 corresponding to all the even modes with n = 2m. As a consequence, the ground state is characterized by a finite even mode fermionic population f 2m = 0. The leading order term in the entanglement scaling in this case is then given by the first term of the Fisher-Hartwig expansion corresponding to a volume law. In particular the maximum amount of entanglement allowed by the ground state degeneracy is obtained for f 2m = 1/2 for every even mode. In this case, the logarithmic corrections become zero since a n = b n = 0, and therefore B

D Discontinuities in the strong long-range regime
In this Appendix we provide a detailed analysis of the discontinuities of the matrix symbol G k in the strong long-range regime 0 < α 1 , α 2 < 1 for different values of the chemical potential. As discussed in the Section 5 of the main text, in this regime the matrix symbol formally develops and infinite number of discontinuities which originate from the discrete nature of the quasiparticle spetrum. However, it is important to notice that, even if the spectrum is labeled by the discrete index n leading to a finite gap between the ground state and the first excited levels, still for n 1 all the modes accumulate around ω ∞ = |h|. This means that an extensive number of single-particle states is almost degenerate. Consequently, as long as h = 0, we may expect only the first few modes around n = 0 to provide a significant contribution to the symbol discontinuity, leading to a coefficient B ν (h = 0) = O(1). Then, in order to understand the qualitative behavior of B ν (h = 0), it is useful to consider the approximation in which only the first discontinuities between the n = 0 and the first two degenerate levels n = ±1 are considered In order to compute this two contributions we have to compute the angles φ 0 and φ ±1 defined by the conditions For n = 0 we find that, independently of the value of α, the angle reads This discontinuity at the quantum critical point h = 1 is due to the fact that at this point the spectrum becomes gapless for n = 0, and it is at the origin of the discontinuity in the scaling coefficient which can be seen in Fig. 6a of the main text. The angles for n = ±1 cannot be computed exactly in close form for generic power law decaying exponent, however as a consequence of the fact thatt n =t −n , ω n = ω −n while∆ n = −∆ −n , we have that cos φ n = cos φ −n sin φ n = − sin φ −n , (D.4) and then φ n = −φ −n . Combining these properties with Eq.(5.5) we obtain number of discontinuities (see Eq. (5.5)), and the results obtained in the single discontinuity approximation. We notice that the single discontinuity approximation correctly reproduces the qualitative behavior of the scaling coefficients for sufficiently high α > 0.5 and for values of the chemical potential h which are sufficiently far from h = 0. In particular, the discontinuity of the coefficients at the quantum critical point h = 1 is captured by the approximated result.
On the other hand, when the chemical potential approaches the h → 0 limit and for sufficiently small decay exponents α < 1/2, the single discontinuity approximation turns out to be no more accurate. Indeed, in this case the number of relevant discontinuities grows with the subsystem size, leading to a subvolume law entanglement scaling. This fact can be understood by considering the h = 0 point. In this case, the spectrum accumulation point becomes ω ∞ = 0. More precisely, it is important to notice that, while at the leading order as n → ∞ the spectrum goes to zero as ω n = O(n α−1 ), independently of the parity of the mode, on the contrary next to leading order corrections differ if n is even or odd. In particular, if we perform a next to leading order expansion of the terms entering the coefficient B This result leads to the scaling of the Rényi entropy in Eq. (5.11) of the main text.