Entanglement Properties of Disordered Quantum Spin Chains with Long-Range Antiferromagnetic Interactions
Youcef Mohdeb, Javad Vahedi, N. Moure, A. Roshani, Hyun-Yong Lee, Ravindra N. Bhatt, Stefan Kettemann, Stephan Haas
EEntanglement Properties of Disordered Quantum Spin Chains with Long-RangeAntiferromagnetic Interactions
Y. Mohdeb, ∗ J. Vahedi,
1, 2, 3, † N. Moure, ‡ A. Roshani, § Hyun-YongLee,
5, 6, 7, ¶ R. N. Bhatt,
8, 9, ∗∗ Stefan Kettemann,
1, 10, †† and Stephan Haas
4, 1, ‡‡ Department of Physics and Earth Sciences, Jacobs University Bremen, Bremen 28759, Germany Technische Univ. Braunschweig, Institut f. Mathematische Physik,Mendelssohnstr. 3, 38106 Braunschweig, Germany Department of Physics, Sari Branch, Islamic Azad University, Sari 48164-194, Iran Department of Physics and Astronomy University of Southern California, Los Angeles, CA 90089-0484 Department of Applied Physics, Graduate School, Korea University, Sejong 30019, Korea Division of Display and Semiconductor Physics, Korea University, Sejong 30019, Korea Interdisciplinary Program in E · ICT-Culture-Sports Convergence, Korea University, Sejong 30019, Korea Department of Electrical Engineering, Princeton University, Princeton, New Jersey 08544, USA School of Natural Sciences, Institute for Advanced Study ∗ ∗ ∗ , Princeton, New Jersey 08540, USA Division of Advanced Materials Science, Pohang University ofScience and Technology (POSTECH), Pohang 790-784, South Korea (Dated: September 24, 2020)Entanglement measures are useful tools in characterizing otherwise unknown quantum phasesand indicating transitions between them. Here we examine the concurrence and entanglemententropy in quantum spin chains with random long-range couplings, spatially decaying with a power-law exponent α . Using the strong disorder renormalization group (SDRG) technique, we find byanalytical solution of the master equation a strong disorder fixed point, characterized by a fixedpoint distribution of the couplings with a finite dynamical exponent, which describes the systemconsistently in the regime α > /
2. A numerical implementation of the SDRG method yields apower law spatial decay of the average concurrence, which is also confirmed by exact numericaldiagonalization. However, we find that the lowest-order SDRG approach is not sufficient to obtainthe typical value of the concurrence. We therefore implement a correction scheme which allows usto obtain the leading order corrections to the random singlet state. This approach yields a power-law spatial decay of the typical value of the concurrence, which we derive both by a numericalimplementation of the corrections and by analytics. Next, using numerical SDRG, the entanglemententropy (EE) is found to be logarithmically enhanced for all α , corresponding to a critical behaviorwith an effective central charge c = ln2, independent of α . This is confirmed by an analyticalderivation. Using numerical exact diagonalization (ED), we confirm the logarithmic enhancementof the EE and a weak dependence on α . For a wide range of distances l , the EE fits a criticalbehavior with a central charge close to c = 1, which is the same as for the clean Haldane-Shastrymodel with a power-la-decaying interaction with α = 2. Consistent with this observation, we findusing ED that the concurrence shows power law decay, albeit with smaller power exponents thanobtained by SDRG. We also present results obtained with DMRG and find agreement with EDfor sufficiently small α <
2, whereas for larger α DMRG tends to underestimate the entanglemententropy and finds a faster decaying concurrence.
I. INTRODUCTION
Long-range interactions in disordered quantum many-body systems arise in a variety of physical contexts. Forexample, randomly placed magnetic impurities in dopedsemiconductors are known to interact with each othervia exchange couplings that depend on their separationdistance . While in insulating regimes these interac-tions typically decay exponentially, J ( r ) ∝ exp( − r/ξ ),in the metallic regime they are mediated via the RKKYmechanism, and therefore decay with a power law, J ( r ) ∝ r − d , where d is the dimension of the host system. Power-law random long-range interactions have also been foundto occur in quantum glasses, where tunneling ions formlocal two-level systems which are interacting by dipole-dipole and elastic interactions . More recently, tunablelong-range Heisenberg interactions have been demon- strated in quantum simulators based on trapped ions .There have been many theoretical studies of quantumspin chains with random short-range interactions. Forexample, it is known that bond disorder drives spin-1 / , characterized by a critical ground state, i.e.a product state of singlets which are formed at randomdistances. The entanglement entropy of such critical spinchains scales logarithmically, with a central charge c ofthe corresponding conformal field theory . As the IRFPis a critical point, the entanglement entropy is logarith-mically enhanced. The corresponding central charge ˜ c is,however, smaller by a factor ln 2 than the central chargeof the corresponding clean spin chain. This has beenderived using the strong disorder renormalization group(SDRG) method . a r X i v : . [ c ond - m a t . d i s - nn ] S e p More recently, the SDRG has also been applied to dis-ordered quantum spin systems with long-range interac-tions, where it has been shown to lead to a self-consistentdescription of its thermodynamic properties in terms ofa novel strong disorder fixed point . However, quan-tum information theoretical measures, such as the en-tanglement entropy and the concurrence have not yetbeen studied in such systems. Here, we analyze thesequantities for disordered quantum spin-1 / α , using the SDRG, exact diagonaliza-tion (ED), and density matrix renormalization group(DMRG). There have been indications that such sys-tem undergo a phase transition at a critical decay ex-ponent α c from a localized regime for α > α c , to adelocalized regime for α < α c , similar to the delo-calization transition of disordered fermions with long-range hoppings . Thus, a logarithmic enhancementof the average entanglement entropy is expected at somevalue of α . Moreover, since the spin-1/2 Haldane-Shastrymodel, a Heisenberg model with power-law long-rangeinteractions, is known to be critical for α = 2 with a con-formal charge c = 1 , and as it was suggested in Ref.9 that disordered spin chains are critical when the cleanspin chain is critical with a central charge ˜ c , smaller bya factor ln 2 than the central charge of the clean criticalspin chain, one could expect a logarithmic enhancementof the entanglement entropy at α = 2 with ˜ c = ln 2 in thepresence of disorder. This motivates us to examine theentanglement properties of disordered quantum spin sys-tems with long-range power-law couplings as a functionof the power exponent α .Here, we focus on the bond disordered XX-spin chainwith long-range couplings, defined by the Hamiltonian H = (cid:88) i II. THE STRONG DISORDERRENORMALIZATION GROUP Here we describe how to apply the SDRG to thismodel, with the aim to evaluate the concurrence and theentanglement entropy .Choosing the pair with the largest coupling ( i, j ),which forms a singlet (see Fig. 1), we take the expec-tation value of the Hamiltonian in that particular singletstate within second-order perturbation in couplings withall other spins . This yields the long-range renormaliza-tion rule for the couplings between spins ( l, m ) in the XXmodel, ( J xlm ) (cid:48) = J xlm − ( J xli − J xlj )( J xim − J xjm ) J xij . (4)For the Heisenberg model, this rule differs by a numericalprefactor 1 / . Previously it was found that in long-range coupled disordered quantum spin chains there is astrong disorder fixed point . This was realized by in-specting the evolution of the width of the distributionof couplings J with the RG flow. In the short-rangecase, i.e. at the infinite randomness fixed point (IRFP),this distribution gets wider at every RG step, with width W = ( (cid:104) ln( J/ Ω ) (cid:105) − (cid:104) ln( J/ Ω ) (cid:105) ) / = ln(Ω / Ω) = Γ Ω , increasing monotonically as the RG scale Ω is lowered. Incontrast, for long-range couplings with finite α the width W saturates and converges to W = Γ, with Γ = 2 α inthe XX-limit . The convergence to a finite dynamicalexponent Γ characterizes the new strong disorder fixedpoint (SDFP). For large number of spins N, and in thelimit of a small RG scale Ω, the resulting distributionfunction of renormalized couplings J at RG scale Ω wasfound to converge to , P ( J, Ω) = 1ΩΓ Ω (cid:18) Ω J (cid:19) − / Γ Ω . (5)At the IRFP, Γ Ω increases monotonically as Γ Ω =ln Ω / Ω, when Ω, the largest energy at this renormal-ization step, is lowered. Here Ω is the initially largestenergy in the SC. FIG. 2. Definition of couplings between removed spins i, j and other spins l, m (dotted lines). At the SDFP, however, Γ Ω is found to converge to anasymptotic finite value Γ Ω → α , yielding a narrowerdistribution with finite width Γ. We note that this distri-bution is less divergent as J → P ( J, Ω),given by − ∂P ( J, Ω) ∂ Ω = P (Ω , Ω) (cid:90) (cid:32) (cid:89) i =1 dJ i P ( J i , Ω) (cid:33) × δ (cid:18) J − J + ( J − J )( J − J )Ω (cid:19) . (6)The numbering of the couplings is defined in Fig. 2.In order to proceed, let us note that the renormalizationcorrection to the bare coupling J , given by δJ = ( J − J )( J − J ) / Ω, is always smaller than J , δJ < J , as canbe checked numerically . Performing first the integralover J , we can expand the distribution function of J = J + δJ in δJ , finding to all orders in δJ − ∂P ( J, Ω) ∂ Ω = P (Ω , Ω) ∞ (cid:88) n =0 n ! (cid:104) (cid:18) ( J − J )( J − J )Ω (cid:19) n (cid:105) ∂ nJ P ( J, Ω) , where (cid:104) ... (cid:105) denotes the averaging with the distributionfunctions P ( J i , Ω) , for i = 1 , , , . Thus, we find − ∂∂ Ω P ( J, Ω) = P (Ω , Ω) (cid:32) ∞ (cid:88) n =2 A n n ! Ω n ∂ nJ (cid:33) P ( J, Ω) , (7)where A n = (cid:104) (cid:18) ( J − J )( J − J )Ω (cid:19) n (cid:105) . (8)We see that the strong disorder scaling ansatz Eq. (5) isa solution of Eq. (7) for Ω → µ = 1 / Γ = const. ,when neglecting all renormalization corrections given bythe terms n ≥ 2. To solve Eq. (7), including the FIG. 3. Probability density function p ( x = J/ Ω) on a doublelogarithmic scale for µ = 0 . µ = 0 . . p ( x = J/ Ω) on a doublelogarithmic scale for µ = 0 . µ = 0 . . 1. Inset: Magnification showing the exponentialdecay to µ as x = J/ Ω → correction terms, we make the ansatz P ( J, Ω) = 1Ω p ( x = J/ Ω) , (9)where the probability density function p ( x ) is to be de-rived with the condition that p (1) = µ , so that p ( x ) = p ( x ) = µx µ − corresponds to strong disorder scaling Eq.(5) with µ = const. . Inserting this ansatz, we find thefollowing ordinary differential equation for p ( x ),(1 − µ ) p ( x ) + x ∂∂x p ( x ) = ∞ (cid:88) n =2 A n n ! ∂ nx p ( x ) , (10)with the condition p (1) = µ .Choosing an iterative approach, we can first approx-imate the coefficients A n , by setting p ( x ) = p ( x ) = µx µ − in Eq. (9), yielding A n ≈ µ Γ ( µ )Γ ( n + 1)( n + µ ) Γ ( µ + n + 1) . (11)We note that A n ∼ /n for 0 < µ (cid:28) 1. For n = 2, thissimplifies with Γ( z +1) = z Γ( z ) to A = 4 µ / (( µ +2) ( µ +1) ), and we can solve to second order n = 2 the resulting2nd order ordinary differential equation exactly in termsof products of Hermite polynomials and Hypergeometricfunctions, as plotted in Fig. 4. For comparison, we plot p ( x ) = µx µ − in Fig. 3. We see that p ( x ) has thesame power law p ( x ) ∼ x µ − dependence for intermediate0 (cid:28) x = J/ Ω < p ( x ), but converges for x → µ approximately as p ( x → ≈ /µ / . We note that for small 0 < µ (cid:28) A n ∼ µ /n , so that higher order terms decay fast with n , andthe pdf p ( x ) obtained for n = 2 is a good approximationfor small µ < z J , where z is a real num-ber, and then integrate both sides over J . Thereby, onecan generate all moments of the distribution by partialderivatives in z and thus the entire distribution function.Since we are interested in the fixed point solution for largesystems, we set the lower integration limit to 0. The up-per integration limit is provided by Ω. We next insertthe ansatz for the distribution of renormalized couplingsEq. (5), and find a differential equation for the expo-nent µ (Ω) = 1 / Γ Ω (a detailed derivation is provided inAppendix A), ∞ (cid:88) n =1 − µ (cid:48) (Ω)(Ω t ) n ( n − µ + n ) = ∞ (cid:88) k =1 ,k (cid:48) =0 µ (Ω)Γ ( µ )(2 k )!(Ω t ) k + k (cid:48) ( k + µ ) Γ ( µ + 2 k + 1) k (cid:48) !( µ + k (cid:48) )Ω . (12)Here, Γ( x ) is the standard Gamma function. Since weare searching for the fixed point solution of the masterequation, we keep only the leading order in Ω of Eq. (12).This yields for Ω → µ (cid:48) (Ω) = 0 . (13)Thus, the solution approaches for Ω → /µ (Ω = 0) = Γ is finite, but cannot be determined by this calculationalone. To this end we need to use a scaling argument, asoutlined below. The distribution of lengths l of singlets atrenormalization scale Ω can also be derived from a masterequation, as was done in Ref. for the IRFP. However,noting that the strong disorder fixed point distribution P ( J, Ω) was obtained above to lowest order in a Taylor expansion in the difference between the renormalized andthe bare coupling J ( l ) = J ( l/a ) − α , it is clear that thereis a strong correlation between the distribution of lengthsand bare couplings. Therefore, we can make to lowestorder the ansatz, P ( J, Ω , l ) = P ( J, Ω) δ ( J − J ( l/a ) − α ) | dJ ( l/a ) − α dl | . (14)We note that in each RG step, a fraction dn/n (Ω) ofthe remaining spins n (Ω) at renormalization energy Ωare taken away. Since this is due to the formation of asinglet with coupling J = Ω, this fraction should be equalto 2 P ( J = Ω , Ω) d Ω, leading to the differential equation dnd Ω = 2 P ( J = Ω , Ω) n (Ω) . (15)Since at the SDFP P ( J = Ω , Ω) = µ/ Ω, the density ofnot yet decimated singlets at the RG scale Ω is given by n (Ω) = n (cid:18) ΩΩ (cid:19) µ , (16)where n = N/L = 1 /l is the initial density of spins.Defining l Ω as the average distance between spins at RGscale Ω, we have n (Ω) = 1 /l Ω . Thus, it follows from Eq.(16) that the RG energy scale is related to the length l Ω via Ω ∼ l − µ Ω . As we discussed above, the strongdisorder fixed point distribution is dominated by the barecouplings, which scale with distance l as J ∼ l − α . Inorder for the energy-length scaling to be consistent at allRG scales Ω, it follows necessarily that2 µ = 1 /α. (17)As we have shown above by solution of the master equa-tion that the strong disorder fixed point distribution is asolution for µ < , we can conclude that it is expected togive consistent results for α > / . We are now ready to evaluate the distribution of singletlengths by integrating over the renormalization energy Ω P s ( l ) = 2 (cid:90) Ω Ω ∗ d Ω l n (Ω) P ( J = Ω , l ; Ω) , (18)where we have used the fact that l n (Ω) is the fractionof not yet decimated spins at Ω, and the factor 2 comesfrom the normalization condition for P s ( l ), and Ω ∗ = J (( L − /a ) − α is the smallest possible energy scale.Thus, for L → ∞ we find P s ( l ) = (cid:90) Ω d ΩΩ (cid:18) ΩΩ (cid:19) µ δ (Ω − J a α l α ) J a α l α +1 (19)= al . Here we used µ = 1 / (2 α ) and the fact that Ω = J isthe largest energy scale. So far, we have assumed that l can take any value within the interval [ a, L ] . Taking intoaccount that the distance can in fact only take discretevalues l i = ia, the properly normalized probability massfunction p i for L/a → ∞ would be p i = 6 / ( π i ). How-ever, as the renormalization of the couplings will becomemore important as the RG scale Ω is lowered, correspond-ing to intermediate and long length scales l , let us checkthe actual distribution in that regime, by numerical im-plementation of the SDRG. In Fig. 5 we show resultsfor N = 200 spins, α = 1 . α = 4 . NL = 0 . , . , . 04. We find that the pref-actor scales with the inverse density, the average spacingbetween the spins l = L/N as P ( l s ) ≈ L N l − s , for a (cid:28) l s (cid:28) L, (20)for both α = 1 . α = 4 . l s < l = L/N , anda slower decay for large l s → L . Note that the numerical P ( l s ) is properly normalized, when summing over all l s ∈ [ a, L − a ]. l s P ( l s ) = 1.8 NL = 0.2 NL = 0.1 NL = 0.04 L N l s (a) l s P ( l s ) = 4.8 NL = 0.2 NL = 0.1 NL = 0.04 L N l s (b) FIG. 5. Probability distribution of singlet lengths P ( l s ),where l s is the real physical distance, on a logarithmic scale,obtained via SDRG for N = 200 and for various filling fac-tors N/L = 0 . , . , . 04. The data were generated for 10 disorder realizations with fixed α = 1 . α = 4 . L N l − s , fitting the intermediate a (cid:28) l s (cid:28) L scaling regime. III. CONCURRENCE The entanglement between any two spins of a chaincan be quantified by the concurrence between them .When the spins in the chain are in a pure state, | ψ (cid:105) , suchas the RS state, the concurrence between spins m and n is given by the correlation function C mn = |(cid:104) ψ | σ ym σ yn | ψ (cid:105)| , (21)which is the absolute value of the overlap between theoriginal state and the state obtained after spins m and n have been flipped. As reviewed above, the SDRG proce-dure yields a random singlet state (RSS) as an approxi-mation of the ground state of the system, even when long-range couplings are present. The RSS is known to be-come asymptotically exact for short-range models, char-acterized by the infinite randomness fixed point (IRFP) .This RSS is a product state that can be written in theform | ψ RS (cid:105) = (cid:79) { i,j }∈ RS | ij (cid:105) , (22)where | ij (cid:105) = ( | ↑ i ↓ j (cid:105) − | ↓ i ↑ j (cid:105) ) / √ i and j , and the directproduct extends over all singlets forming the RSS. FromEq. (22) it becomes apparent that when the system isin the RSS, the concurrence between the two spins i, j isgiven by C ij = (cid:40) i and j form a singlet in the RSS , . (23)Thus, as the RSS disregards any but the strongest cou-plings, it fails to account for corrections by residual weakcouplings. In order to include the finite amount of entan-glement which prevails between spins that do not forma singlet in the RSS, and in turn weakens the entangle-ment between the spins that do form a singlet during theSDRG procedure, we need to find a consistent strategyto include these corrections in the SDRG scheme.Before proceeding any further, let us first review theknown results for the scaling behavior of the mean andtypical concurrence at the IRFP, i.e. the fixed point forshort-range coupled random spin chains, which will serveas a baseline to compare with results obtained in theSDFP. IV. MEAN AND TYPICAL CONCURRENCEAT THE IRFP As Fisher noted in Ref. , the mean correlation func-tion between spins at long distances l is dominated byrare events. Typically, two distant spins enumerated byindices n and n = n + n will not form a singlet and willtherefore be very weakly correlated. However, in the rareevent that they do form a singlet in the RSS, they willbe strongly correlated and therefore will dominate themean correlation function, and thereby the concurrence,at large distances. As a result, the mean correlation func-tion must be proportional to the ratio of singlets formedat index distance n , P ( n ), which is related to the distri-bution function of real distances l , P ( l ), which we havereviewed above. At the IRFP, P ( l ) ∼ a/l , and we seethat the mean correlations decay faster when disorderis introduced: in the clean case they decay more slowlyas l − / . Hoyos et. al. explicitly noted that for chainswith open boundary conditions C ( n ) = 0 for even indexdistance n , which yields the more accurate result at theIRFP , , ¯ C ( n ) ∼ n × (cid:40) n is odd , n is even . (24)It was also noted in Ref. that there are additional termsin ¯ C ( n ) that decay faster for large n , the next leadingterm decaying as 1 /n / . While rare events dominate themean concurrence, we expect a different behavior for thetypical concurrence which is proportional to the typicalvalue of the coupling between two spins J at the RGscale Ω l at which they are decimated. This value can becalculated using the full IRFP distribution of J . Thereby,one finds , C typ ( l ) ∼ e − k √ l , (25)where it was used that at the IRFP the distance betweenspins is related to the RG scale Ω l by l ∼ ln (Ω / Ω l ) toobtain the scaling behavior as an extended exponential,with k being a non-universal constant of order unity .Thus, at the IRFP the typical value of the concurrencedecays exponentially fast with distance, whereas its meandecays with a power law. In the next section we introducea general approach to include corrections to the RS state,and we will, in particular, investigate how the typicalconcurrence decays with l at the SDFP. V. CORRECTIONS TO THE RANDOMSINGLET STATE To incorporate the effects of the couplings which areneglected in the RS state we define an effective Hamilto-nian as the sum of the Hamiltonion ˜ H with all effective,renormalized couplings taken into account in the SDRG,and a perturbation term ˜ H (cid:48) , which includes all couplingswhich are neglected in the SDRG,˜ H = (cid:88) { ij }∈ RS ˜ J ij S i S j (cid:124) (cid:123)(cid:122) (cid:125) ˜ H + (cid:88) { ij } / ∈ RS ˜ J ij S i S j (cid:124) (cid:123)(cid:122) (cid:125) ˜ H (cid:48) . (26)Now, we can perform perturbation theory in the term ˜ H (cid:48) to obtain the ground state of the disordered XX chain in FIG. 6. Full lines indicate strongest bare couplings. Someweaker couplings are indicated by thinner lines. Dashed linesconnect spins that form singlets in the random singlet state.Top: Spins { p, q } of Eq. (32) do not form a singlet with eachother. Bottom: { pq } of Eq. (32) do form a singlet. first order of ˜ H (cid:48) as, | ψ (cid:105) = | ψ RS (cid:105) + (cid:88) β (cid:104) ψ β | ˜ H (cid:48) | ψ RS (cid:105) E RS − E β | ψ β (cid:105) , (27)where E RS is the ground state energy of the RS state | ψ RS (cid:105) , and the sum runs over all excited states of ˜ H , | ψ β (cid:105) , as labeled by the index β , with eigen energies E β .The excited states can be obtained by combinations oftriplet states, as obtained by excitations of the singletsin the RS state. We note the useful relation˜ H (cid:48) | ψ RS (cid:105) = 14 (cid:88) { nl }(cid:54) = { mk } ( J nm + J lk − J nk − J ml ) × ( | + nl (cid:105)|− mk (cid:105) + |− nl (cid:105)| + mk (cid:105) ) (cid:79) { ij }(cid:54) = { nl }(cid:54) = { mk } | ij (cid:105) , (28)where |± nl (cid:105) ≡ | ( S = 1 , M = ± (cid:105) are two of the tripletstates formed of the spins { l, n } when the spins l and n form a singlet in the RSS. Thus, the double sum and thedirect product run over all singlet pairs in the RSS, withthe exceptions specified under the summation and directproduct signs . From this result, it becomes apparentthat the only excited states that contribute to the sumin Eq. (27) are of the form | ψ β (cid:105) = |± nl (cid:105)|∓ mk (cid:105) (cid:79) { ij }(cid:54) = { nl }(cid:54) = { mk } | ij (cid:105) , (29)whose energy difference to the RS state is given by E β − E RS = ( J nl + J mk ) / 2. With this in mind, Eq. (27)transforms into the final form of the ground state withcorrections, | ψ (cid:105) = c ψ | ψ RS (cid:105) − c ψ (cid:88) { nl }(cid:54) = { mk } J nm + J lk − J nk − J ml J nl + J mk × ( | + nl (cid:105)|− mk (cid:105) + |− nl (cid:105)| + mk (cid:105) ) (cid:79) { ij }(cid:54) = { nl }(cid:54) = { mk } | ij (cid:105) , (30)with a normalization constant c ψ given by c ψ = (cid:88) { nl }(cid:54) = { mk } (cid:18) J nm + J lk − J nk − J ml J nl + J mk (cid:19) − / . (31)Now, we are all set, and we can use the perturbedground state, given by Eq. (30), to derive the concur-rence between the spins with indices p and q . We findthe conditional expression C NSpq = c ψ (cid:12)(cid:12)(cid:12)(cid:12) J pq + J rs − J pr − J qs J ps + J qr (cid:12)(cid:12)(cid:12)(cid:12) if { pq } / ∈ RS,C Spq =1 − c ψ (cid:88) { mk }(cid:54) = { pq } (cid:18) J pm + J qk − J pk − J mq J pq + J mk (cid:19) if { pq } ∈ RS, (32)where C NSpq is the concurrence between spins p and q ifthey do not form a singlet in the RS state, and C Spq isthe one, if they do. The indexes r and s in the first linecorrespond to the spins that form a singlet in the RSstate with spins q and p , respectively, as shown in Fig.6. Eq. (32) has all the properties expected for the con-currence between two spins in a chain with long-rangecouplings. It does give a non-zero value for all pairs thatdo not form a singlet in the RS state, and it gives a con-currence smaller than 1 for spins that do. In Fig. 7(a) we show numerical results for the meanconcurrence as a function of index distance n in the disor-dered long-range XX chain with N = 800 spins randomlyplaced on L = 80000 sites, as obtained by inserting thenumerical results of the SDRG for the couplings into Eq.(32). First, we note that, in contrast with the short-range case, there is a finite concurrence for even valuesof n . This is expected by looking at the form of Eq. (32)and recalling that C ( n ) = 0 for even n was due to theimpossibility of crossing singlets in the RS state. How-ever, there is still a clear difference between even valuesof n (bottom) and odd ones (top), as indicated by theclear separation of two sets of curves. Both sets of curveshave a weak dependence on α and a regime in which itcan be fitted with a power law,¯ C ( n ) ∼ n − γ e,o , (33)where γ e,o are the decay powers for even and odd valuesof l , respectively. In fact, by using linear regression fitsin the logarithmic scale we find γ e = 1 . ± . 04 and γ o =1 . ± . 04 for both α = 0 . α = 2 . 0. Even with thecorrections to the RS state the concurrence for odd valuesof n is still dominated by rare events (singlets formed atlong distances), as indicated by the small deviation from¯ C ( n ) ∼ n − decay for all values of α . On the other hand,for even values of n , the decay is slower, the power-lawregime is smaller, and the amplitudes have a stronger α γ (c) FIG. 7. (a) Mean concurrence of disordered XX chains withpower-law, long-ranged couplings, averaged over M = 18960realizations, as function of index distance n . Power exponentsof the coupling strengths (Eq. 2): α = 0 . , ..., . N = 800spins placed randomly on one of L = 100 N sites. n = | p − q | isthe index distance. Top curves: odd n , bottom curves: even n . (b) typical concurrence as a function of index distance n .(c) Fits of linear regime in the typical concurrence to a power-law decay with power γ . The fitted exponents γ are plottedas a function of α . dependence on α , as expected from Eq. (32). For bothsets of curves, a saturation at large values of n can beseen, which is a finite size effect.The typical value of the concurrence is shown in Fig.7(b). A clear power-law behavior of the form C typ ( n ) = exp( (cid:104) ln C ( n ) (cid:105) ) ∼ n − γ ( α ) , (34)is found. Here, the power γ ( α ) has a strong dependenceon α , unlike the decay powers of the mean value. In fact,we find γ ( α ) to be linear in α , with a linear regression fitgiving γ ( α ) = 1 . α + 2 . , (35)which indicates that the typical concurrence decays fasterthan the mean concurrence for all values of α . It is alsoworth noting that since the typical concurrence decaysas a power law, it decays slower than in the IRFP case,where it has the extended exponential behavior stated inEq. (25).Now, let us see if we can use the corrections to theRS state, Eq. (32) to find the scaling behavior of thetypical concurrence analytically. The probability massfunction as a function of the index distance n decays toleading order as 1 /n , in accordance with the distributionof lengths, Eq. (20). Thus, noting that in the randomsinglet state only spins at odd index distance are pairedwe find P ( n ) = (cid:18) c n + O ( 1 n / ) (cid:19) × (cid:40) n is odd , n is even . (36)Since P ( n ) must be normalized, (cid:80) N − n =1 P ( n ) = 1, thecoefficient c depends on the weight of the faster decay-ing terms. In Ref. , c = 2 / P ( n s ) on a logarithmic scale, obtained using theSDRG for N=400 with a filling factor NL = 0 . 1. Thedata were obtained from 25000 realizations of the disor-der for each α . The red dashed line corresponds to thefunction n − s , showing good agreement for intermediateand large values of n s , in agreement with Eq. (36) with c = 2 / 3. We also observe deviations due to a faster de-caying term at small n , in agreement with Eq. (36). Thesaturation at large values of n has been checked to be afinite size effect.Now, we can calculate the mean value of the concur-rence as function of index distance n via¯ C ( n ) = (cid:104) C pq (cid:105) n = | p − q | = P ( n ) C Spq + (1 − P ( n )) C NSpq , (37)and the typical value of the concurrence via C typ ( n ) = exp (cid:0) P ( n ) ln C Spq + (1 − P ( n )) ln C NSpq (cid:1) . (38)We note that in the random singlet state without cor-rections, Eq. (37) gives with Eq. (36) C RSS ( n ) = ( c n + O ( 1 n / )) × (cid:40) n is odd , n is even , (39)when properly normalized, so that (cid:80) ∞ n =1 C RSS ( n ) = 1 , since each spin can form a singlet with only one otherspin, in which case the concurrence is exactly equal toone, giving c = 2 / 3. Now, we can find the scaling of C NSpq including the corrections to the RSS by noting thatthe distance l = | r p − r q | is always larger than the distancebetween the spins which formed singlets in the RS state, n s P ( n s ) = 0.8= 1.8= 4.8 n s FIG. 8. Probability distribution of the singlet lengths, P ( n s ),as a function of index distance, on a logarithmic scale, ob-tained via the SDRG for N=400 and a filling factor NL = 0 . α . The red dashed line corresponds to thefunction n − s , which fits the intermediate n s scaling regime. see Fig. 6 top, l = | r p − r s | , n = | r q − r r | . Thus, wecan Taylor expand Eq. (32) in l /l and l /l . Thereby wefind C NSpq = c ψ l l J l + J l ∂ l J l = c ψ α ( α + 1) l l l − α + l − α l − α − . (40)Noting that at the SDRG J l = J ( l ) − α can be relatedto the index distance n by assuming that l scales with n as dictated by the density of spins n = 1 /l , we cansubstitute l ∼ nl . As the bonding lengths l and l can take any value between 1 and l , we can average overall possible values and find as function of index distance C NSpq = k NS n − α − , with a constant k NS .Similarly, we can do an expansion C Spq in the distancebetween the spins of the singlet states in Fig. 6 bottom l = | l p − l q | and l = | l k − l m | to find C Spq = 1 + c ψ α α + 1) l − α . (41)Thus, we can insert Eqs. (40,41) into Eqs. (37,38) toget the scaling of the mean and typical concurrence withindex distance n ¯ C ( n ) = P ( n )(1 + c ψ α α + 1) ( nl ) − α )+ (1 − P ( n )) k NS n − α − , (42)and the typical value of the concurrence by C typ ( n ) = ( k NS ( nl ) − α − ) − P ( n ) (1 + c α α + 1) ( nl ) − α ) P ( n ) ∼ n − α − , (43)where we used that P ( n ) (cid:28) 1. For the typical value wethus find very good agreement of the power γ = 2+ α withthe result obtained with numerical SDRG for system size,Eq. (35), γ ( α ) = 1 . α + 2 . 02. The mean value shows amore complicated behavior with different index distanceregimes dominated by either P ( n ) or the power law decayin Eq. (42).When plotting the concurrence as function of physicaldistance l , we expect for small concentrations of spins N/L (cid:28) P ( l ) = c (cid:48) l − + O ( l − / ), Eq. (20) whereit was found that c (cid:48) = 1 / 3. Thus, we get the averageconcurrence as function of real distance l ,¯ C ( l ) = P ( l )(1 + c ψ α α + 1) l − α )+ (1 − P ( l )) k NS l − α − , (44)and the typical value of the concurrence as function ofthe physical distance l is C typ ( l ) = ( k NS l − α − ) − P ( l ) (1 + c α α + 1) l − α ) P ( l ) ∼ l − α − , (45)For comparison we show in Fig. 9 the results of numer-ical exact diagonalization for spin chains with N = 20distributed randomly of length L = 200, as plotted asfunction of the physical distance l and averaged over M = 2000 samples. We find that the average concurrencedecays with a power which increases slightly with increas-ing interaction power α but remains for all α smaller thanthe power 2, obtained in SDRG. Also, there is no evenodd effect, as expected when plotting the concurrence asfunction of the physical distance l . The typical concur-rence is found to decay with a power law, with exponent γ = 0 . α + 1 . , as obtained by a fit of all results for α = 0 . , . , . , . , . , . . , linearly increasing with α as found from SDRG, but with a smaller slope. Notethat the finite size effects seen in in SDRG, Fig. 7 areexpected to be more dominant in the smaller size used inED as presented in Fig. 9, which may explain the slowerdecay with smaller exponents observed with ED. VI. ENTANGLEMENT ENTROPY The entanglement between two segments of a spinchain, A and B , can be quantified by the von Neumannentropy of the reduced density matrix, S = − Tr( ρ A ln ρ A ) = − Tr( ρ B ln ρ B ) , (46)where ρ A = Tr B ( | ψ (cid:105)(cid:104) ψ | ) and ρ B = Tr A ( | ψ (cid:105)(cid:104) ψ | ) is ob-tained by partially tracing the complete density matrixof the system over all degrees of freedom of subchain B or A , respectively.This entanglement entropy can be used to characterizequantum phase transitions. For clean chains, it has been − − − ¯ C ( l ) (a) α = 0 . α = 1 . α = 1 . α = 1 . α = 2 . α = 2 . l − . l − l − − − − C t y p ( l ) C typ ( l ) ∼ l − γ ( α ) (b) α = 0 . α = 1 . α = 1 . α = 1 . α = 2 . α = 2 . α . . . γ γ = 0 . α + 1 . FIG. 9. Results from exact numerical exact diagonalization.Panel (a): Mean concurrence of disordered XX chains withpower-law, long-ranged couplings, averaged over M = 2000number of realizations as function of the physical distance l for various power exponents α , N = L/ 10 = 20 spins placedrandomly on one of L sites. Panel (b): typical concurrenceas a function of the physical distance l . Inset: Fits of linearregime in the typical concurrence to a power-law decay withpower γ . The fitted exponents γ are plotted as a function of α . shown that at criticality, the entropy of a subchain A oflength l scales as S ( l ) = b c + c l/a ) + k, (47)where c and c correspond to the central charges of thecorresponding 1+1 conformal field theory. In the limit ofinfinite chains with finite partitions of length l, as well asfor periodic chains with large length L (cid:29) l , b = 2 sincethere are two boundaries of the partition, while b = 1for semi infinite chains, when the partition of length l is placed on one side of the chain. k is a non-universalconstant . This scaling behavior with a logarithmic de-pendence on the segment length l is in contrast to thearea law expected for noncritical chains, where it doesnot depend on the length l of the subchain for the one-dimensional case. The simple area law is recovered awayfrom criticality, where it is found that the entropy satu-rates at large l .0In Ref. , it was shown that Eq. (47) also holds for theaverage entanglement entropy of antiferromagnetic spinchains with random short-ranged interactions. In par-ticular, using SDRG, they found that in the disorderedtransverse Ising Model, the effective central charges weregiven by ˜ c = ˜ c = ln(2) / 2, whereas in the Heisenberg andXX model ˜ c = ˜ c = ln(2). Both cases correspond to a fac-tor of ln(2) reduction of the central charge and the entan-glement entropy of their corresponding pure systems ,in accordance with a generalized c -theorem, which statesthat if two critical points are connected by a relevantRG flow, as here induced by the relevant disorder, thefinal critical point has a lower conformal charge than theinitial one .Eq. (47) applies specifically to infinite systems. InRef. , Calabrese and Cardy derived a formula valid forfinite systems of length L , S b ( l ) = b c (cid:18) Lπa sin( πl/L ) (cid:19) + k (cid:48) , (48)where k (cid:48) is a non-universal constant, and c = c has beenassumed . For periodic boundary conditions, there isan additional factor b = 2, since the subsystem is thenbounded by two boundaries, doubling the average num-ber of singlets crossing one of the boundaries, while b = 1for open boundary conditions, when there are 2 parti-tions.Refael and Moore’s method to calculate the average en-tanglement entropy in the presence of disorder is based onthe assumption that the system has been drawn to theIRFP, and the random singlet state (RSS) is a correctrepresentation of its ground state. Since the RSS corre-sponds to a product state of maximally entangled spinpairs, they note that the total entanglement entropy canbe calculated by counting the number of singlets thatcross the boundary between subsystems A and B andthen multiplying this number with the entropy of a sin-glet S = ln(2) . A schematic representation of thismethod is shown in Fig. 10 for a specific realization. Inthis example, when the boundary between A and B isdefined by line a , we obtain an entanglement entropy of3 · S , since three singlets cross over the boundary line.But if the boundary is defined by line b , the entanglemententropy is reduced to 1 · S . FIG. 10. Single realization of the random singlet state, illus-trating the entanglement entropy calculation for two differentboundaries a , b (red lines) between subsystems. Thus, for the random singlet state | ψ (cid:105) one finds S = M ln 2, where M is the number of singlets which cross the partition between subsystems A and B. Theaverage entanglement entropy of a random singlet statecan therefore be related to the probability to find a singletof a certain length n s , in the spin chain of length N , P ( n ). As Refael notes in Ref. , however, one needs totake into account the correlation with the RG history.He found that the average distance between RG steps,where a decimated singlet may cross the same partition,is exactly equal to (cid:104) m (cid:105) = 3 for the disordered nearestneighbor AFM spin chain, resulting in a correlation fac-tor 1 / (cid:104) m (cid:105) . Another way to derive this was outlined by Hoyos etal. in Ref. . The ratio of the number of singlets cross-ing the partition is the probability that a spin on oneside of the partition is entangled with another one on theother side of the partition. Thus, we can relate the en-tanglement entropy directly to the probability to have asinglet of length l , P ( l ), which we derived in the previ-ous chapter. Moreover, as we consider partial filling ofthe lattice sites with spins with density n = N/L < , we need to distinguish the EE as a function of the in-dex distances n = | i − j | between the spins from the oneplotted as function of their physical distance l = r i − r j ,where r i is the position of spin i . As a function of in-dex distance n we obtain for open boundary conditionsand 2 partitions, where one has the length l , the averageentanglement entropy, (cid:104) S n (cid:105) = ln 2 N (cid:88) n s =1 (cid:88) i P i ( n s ) | C . C . , (49)where the crossing condition (C.C.) ensures that onlysinglets are counted where spin i is on the left side ofthe partition (which contains n spins), while spin i + n s is on the right side of the partition. We first count thenumber of possible positions to place a singlet with indexdistance n s across the partition boundary, starting withthe smallest distance n s = 1, and adding successivelysinglets of larger index length. For n s < n there are, inprinciple, n s such possibilities, if the respective spins didnot yet form a singlet with another spin, each with theprobability to form a singlet of length n s , P ( n s ).In order to account for the correlation with existingsinglets, Ref. multiplied the probability P ( n s ) with afactor 1 / . This can be argued to be due to the fact thatfor every spin which may form a singlet with length n s across the boundary there is a second possibility to forma singlet, which is not crossing the boundary. Thereby,one arrives for the chain with OBC with a partition with n spins and one partition boundary at the following ex-1pression (cid:104) S n (cid:105) / ln 2 = 12 n (cid:88) n s =1 n s P ( n s )+ 12 n N − n (cid:88) n s = n +1 P ( n s )+ 12 N (cid:88) n s = N − n +1 ( N − n s ) P ( n s ) . (50)This expression is equivalent to the one given in Ref. (with the difference that they considered an embeddedpartition with two boundaries). By evaluating Eq. (50)using the result Eq. (39) P ( n ) = c /n + O (1 /n / ) forodd n , 0 for even n we find in the limit of N (cid:29) (cid:104) S n (cid:105) ≈ c ln 2 ln n + k + O (1 /n / ) . (51)We thus recover by comparison with Eq. (47) the confor-mal charge given by ˜ c = ˜ c = (6 / c ln 2. Using c = 2 / and confirmed numerically above, wethus find ˜ c = ln .For finite N we can do the summations, written interms of Polygamma functions as function of n : (cid:104) S n (cid:105) = 2 ln 2 < m > π (2 γ − − n + N + ln 16+ 2PG[0 , / n − / − , / N − / , / − ( n − / N − / n PG[1 , / n − / − N PG[1 , / N − / N − n )PG[1 , / − ( n − / N − / γ = EulerGamma = 0 . x, y ] =PolyGamma[ x, y ] is the PolyGamma function.In an attempt to take into account the correlationwith the location of other singlets in the RSS state morerigourously, we could argue that we need to multiply eachterm with the probability that the two spins did not yetform a singlet of other length with another spin, thatis (cid:81) n (cid:48) s (cid:54) = n s (1 − P ( n (cid:48) s )). Thereby, one finds for a randomsinglet state, using Eq. (39) (cid:104) S n (cid:105) / ln 2 = n (cid:88) n s =1 n s P ( n s ) (cid:89) n (cid:48) s (1 − P ( n (cid:48) s ))+ n N − n (cid:88) n s = n +1 P ( n s ) (cid:89) n (cid:48) s (1 − P ( n (cid:48) s ))+ N (cid:88) n s = N − n +1 ( N − n s ) P ( n s ) (cid:89) n (cid:48) s (1 − P ( n (cid:48) s )) . (53)By evaluating Eq. (53) we find in the limit of N (cid:29) (cid:104) S n (cid:105) ≈ c ln 2 ln n + k + O (1 /n / ) . (54)We thus recover Eq. (47) with the conformal charge givenby ˜ c = ˜ c = c ln 2 < 1. Here, for pure P ( n ) = c /n for odd n , P ( n ) = 0 for even n we have c = 8 /π = 0 . O (1 /n / ),Eq. (39), as is confirmed by our numerical results we get c = c ln 2 = 2 / / l . For small filling, N (cid:28) L , the even odd effectas function of the physical distance l is negligible and weget for intermediate range a (cid:28) l s (cid:28) L , Eq. (20), asconfirmed numerically. We thus recover Eq. (47) withthe central charge given by ˜ c = ˜ c = c (cid:48) ln 2 < . Here,using Eq. (20) we get c (cid:48) = 1.Let us next implement the SDRG numerically. Fig.11 shows results, numerically calculating the mean blockentanglement entropy using Refael and Moore’s prescrip-tion illustrated in Fig. 10, which is plotted as functionof the index distance size n of the partition, that countshow many spins are inside a partition. The XX chainwith open boundary conditions has N = 500 spins withlong-range power-law interactions Eq. (2). As mentionedabove, such a prescription implies that the RSS is a goodapproximation of the actual ground state of the chain.The result for the average entanglement entropy asfunction of index distance n is shown in Fig. 11 for vari-ous values of α with open boundary conditions. The blueline is Eq. (48) with a central charge c = ln(2), b = 1and k (cid:48) = 0 . α ,decaying by only a few percent as α is changed from α = 6 to α = 0 . l for the longranged XX-chain with N = 500 spins, open boundaryconditions for various values of α , and for a filling factor NL = 0 . 1. The average was evaluated over 20000 realiza-tions for each α . The blue dashed line corresponds to theCardy law Eq. (48) with b = 1, k (cid:48) = 0 . 05 and a centralcharge c = ln(2), in very good agreement with the an-alytical result c = ln(2). In Fig. 13 the average blockentanglement entropy as function of partition length l isshown, as obtained with numerical SDRG for the longranged XX-chain with N = 200 spins, and for variousfilling factors NL = 0 . , . , . 04, with open boundaryconditions and for two values of α . The average wasevaluated over M = 100000 realizations for each α . Thedashed lines correspond to the Cardy law Eq. (48) witha central charge ˜ c = ln(2), b = 1, in good agreement withthe analytical result.As both analytical and numerical results based onSDRG could be artefacts of the assumption that theground state remains a random singlet state, let us nextcalculate the entanglement entropy using the numeri-cal exact diagonalization method. In Figs. 14 and15 we show the results for the average block entangle-ment entropy, obtained by numerical exact diagonaliza-tion for the long ranged XX-chain for a filling factor NL = 0 . L = 140 , , , 200 and various values of α , and in Fig.15 for various sizes L = 140 , , , α = 0 . , M = 500 random samples as a functionof the real subsystem size l . Here, (cid:104)· · · (cid:105) ens denotes the en-semble average. The dashed lines correspond to Cardy’slaw, Eq. (48), with a central charge ˜ c = 1 . b = 1and k (cid:48) = 0 . 13, which confirms the weak dependence on α .The central charge is found by exact diagonalization to belarger than the one for the short ranged disordered AFMspin chain, and by 1 . α = 2 . l (physical distance), for 1 < l < L/ 4, where the criticalentanglement entropy Eq. (47) is plotted as the black linecorresponding to the approximation of the Cardy law Eq.(48) for l (cid:28) L with a central charge ˜ c = 1, b = 1. Theyellow line is Eq. (47) with ˜ c = ln 2, corresponding tothe result obtained with SDRG. We see that while for awide range of l , the central charge seems to fit to the oneof a clean critical spin chain c = 1, at small l there arestrong deviations tending rather to the SDRG result.As the averaging may diminish the dependence on α ,let us next also consider the full distribution of the en-tanglement entropy for (a) α = 0 . 1, (b) α = 0 . α = 2 . 8, where 500 random samples have been taken,and the system size N = 22 is used, and S = log 2.We note that also the distribution shows only a weak de-pendence on α . It is remarkable that the distribution ispeaked at integer multiples of ln 2, which means that evenwhen averaged over many ensembles, an integer numberof singlets crossing the partition is most likely. n S = 0.8= 1.6= 2.8= 3.6= 4.6= 6.0 FIG. 11. Average block entanglement entropy as obtainedwith the numerical SDRG for the long ranged XX-chain with N = 500 spins, open boundary conditions, and various valuesof α as function of index distance n . The average was evalu-ated over 2 × realizations for each α , The blue line is Eq.(48) with a central charge c = ln(2) b = 1 and k (cid:48) = 0 . S = 0.8= 1.6= 2.8= 3.6= 4.6= 6.0 FIG. 12. Average block entanglement entropy, obtained fromthe numerical SDRG for the long ranged XX-chain of length L = 5000, open boundary conditions for various values of α , and for a filling factor NL = 0 . l . The average was evaluated over 2 × realizationsfor each α . The blue line corresponds to the Cardy law Eq.(48) with a central charge c = ln(2), b = 1 and k (cid:48) = 0 . VII. ENTANGLEMENT MEASURESOBTAINED BY DENSITY MATRIXRENORMALIZATION GROUP We consider again the XX model in which N spins arerandomly distributed on a finite lattice with length L andinteract with each other with long-range interactions, butconsider also an external an external Zeeman magneticfield B , H xx = (cid:88) i,j
To find the ground states (GS) of Eq. (55) for differentrandom realizations, we use the density matrix renor-malization group (DMRG) method and the noise al-gorithm to avoid converging to a local minimum inDMRG calculations. As shown in Ref. 30, DMRG canbe regarded as a method for optimizing variational wavefunctions known as matrix product states (MPS) :3 S l = 1.8 NL = 0.2 NL = 0.1 NL = 0.04 (a) S l = 4.8 NL = 0.2 NL = 0.1 NL = 0.04 (b) FIG. 13. Average block entanglement entropy as a functionof the partition length l (physical distance), obtained fromnumerical SDRG for the long ranged XX-chain with openboundary conditions for N = 200 spins. Various filling factors NL = 0 . , . , . 04 were considered for both α = 1 . α = 4 . M = 100000realizations for each α . The full lines correspond to the Cardylaw Eq. (48) with a central charge ˜ c = ln(2), b = 1. | ψ (cid:105) = (cid:88) { p i } (cid:88) { α i } A p α A p α α · · · A p N α N − | p , · · · , p N (cid:105) (57)which is a representative one-dimensional tensor networkstate . In MPS representation, one can rewrite quan-tum operators in a similar tensor network, so-called thematrix product operator (MPO),ˆ O = (cid:88) { p i p (cid:48) i } (cid:88) { α i } O p p (cid:48) α O p p (cid:48) α α · · · O p N p (cid:48) N α N − | p , · · · (cid:105)(cid:104) p (cid:48) , · · · | . (58)Here, { α i } are the virtual indices which are traced out,and their dimensions are called the bond dimension ofMPO. The success of DMRG for one-dimensional sys-tems is due to the existence of an exact MPO representa-tion with finite bond dimensions for a Hamiltonian withshort-range or exponentially decaying interactions. Forexample, the Hamiltonian in Eq. (55) with α = 0 and . . . . . ¯ S v N (a) α = 0 . α = 1 . α = 2 . α = 2 . ˜c=1.4ln(2) (b) α = 0 . α = 1 . α = 2 . α = 2 . ˜c=1.4ln(2) l . . . . . ¯ S v N (c) α = 0 . α = 1 . α = 2 . α = 2 . ˜c=1.4ln(2) l (d) α = 2 . α = 2 . α = 2 . α = 2 . ˜c=1.4ln(2) FIG. 14. Average block entanglement entropy, obtainedfrom numerical exact diagonalization for the long rangedXX-chain with a filling factor NL = 0 . L =140 , , , α = 0 . , . , . , . l (physical distance). The average was evaluated over M = 2000 realizations for each α . The full lines correspondto the Cardy law Eq. (48) with a central charge ˜ c = 1 . b = 1 and k (cid:48) = 0 . η > O pp (cid:48) αβ = ( I ) pp (cid:48) S x ) pp (cid:48) η ( I ) pp (cid:48) S y ) pp (cid:48) η ( I ) pp (cid:48) S z ) pp (cid:48) η ( I ) pp (cid:48) B ( S z ) pp (cid:48) ( S x ) pp (cid:48) ( S y ) pp (cid:48) ζ ( S z ) pp (cid:48) ( I ) pp (cid:48) αβ , (59)where I is the 2 × η = 0. Unfortu-nately, an exact MPO expression of a Hamiltonian withthe power-law decaying interaction is not allowed withfinite bond dimension regardless of its exponent α . How-ever, one can decompose the power-law functions intoseveral exponential functions as follows | r j − r i | − α (cid:39) m (cid:88) l =1 λ l η | r j − r i | l , (60)where m depends on the distance | r j − r i | and α . Find-ing proper { λ l } and { η l } is not a trivial problem, butPirvu et al. in Ref. 33 have found a systematic and ele-gant way to find them. Employing the fitting procedurein Ref. 33, we have decomposed the power-law functionin Eq. (56) into 17 different exponential functions, andthey are in excellent agreement with the original power-law function in the range of α ∈ [0 . , . 8] on a lattice4 . . . . . ¯ S v N α = 0 . (a) L = 140 L = 160 L = 180 L = 200 l . . . . . ¯ S v N α = 2 . (b) L = 140 L = 160 L = 180 L = 200 FIG. 15. Average block entanglement entropy as obtained bynumerical exact diagonalization for the long ranged XX-chainwith a filling factor NL = 0 . L = 140 , , , α = 0 . α = 2 . l (physical distance). The average wasevaluated over M = 2000 realizations for each α . The yellowline corresponds to the Cardy law Eq. (48)with a centralcharge ˜ c = 1 . b = 1 and k (cid:48) = 0 . with L = 800 ( α ∈ [0 . , . 0] on L = 1200). Consideringthe random distances between neighboring spins, one canfind the following MPO tensor encoding both the long-range interaction and randomness: O p i p (cid:48) i αβ = ( I ) p i p (cid:48) i (cid:126) Γ i ⊗ S x ) p i p (cid:48) i ( (cid:126) Γ i ⊗ I ) p i p (cid:48) i (cid:126) Γ i ⊗ S y ) p i p (cid:48) i (cid:126) Γ i ⊗ I ) p i p (cid:48) i (cid:126) Γ i ⊗ S z ) p i p (cid:48) i (cid:126) Γ i ⊗ I ) p i p (cid:48) i B ( S z ) p i p (cid:48) i ( (cid:126)λ T ⊗ S x ) p i p (cid:48) i ( (cid:126)λ T ⊗ S y ) p i p (cid:48) i ζ ( (cid:126)λ T ⊗ S z ) p i p (cid:48) i ( I ) p i p (cid:48) i αβ , (61)where T denotes the transpose, (cid:126)λ = [ λ , λ , · · · λ m ] T andΓ i is a diagonal matrix with diagonal elements (cid:126) Γ i = [( η η ) r i , ( η η ) r i , · · · , ( η η m ) r i ] T , (62)and r i = | r i +1 − r i | . Performing DMRG with MPO tensorin Eq. (61), one can find the ground state of Eq. (55) fora given α , ξ and random sample. B. Results Now, let us discuss the results obtained that way withDMRG for a chain with open boundary condition andfilling factor N/L = 1 / Concurrence - First, let us begin with the concurrenceresults presented in Fig. 18. The left panel shows theaveraged concurrence as a function of the index dis-tance ( n nm = | n − m | ) between two spins at sites m and n . The distribution of concurrence at a given index dis-5 . . . . . . ¯ S v N (a) L=140 c = log 2 c = 1 (b) L=160 c = log 2 c = 1 l . . . . . . ¯ S v N (c) L=180 c = log 2 c = 1 l (d) L=200 c log 2 c = 1 FIG. 16. Same as Fig. 15 (b) plotted only for 1 < l < L/ l (cid:28) L with acentral charge ˜ c = 1, b = 1, the yellow line is Eq. (47) with˜ c = ln 2. P (a) α = 0 . P (b) α = 0 . 80 1 2 3 S vN / S P (c) α = 2 . FIG. 17. Distributions of the entanglement entropy for (a) α = 0 . 1, (b) α = 0 . α = 2 . 8, where 500 randomsamples and the system size N = 22 are used, and S =log 2. tance n nm = 10 with α = 1 and α = 2 . C ( n ) shows for some rangea power law decay, followed by an exponential decay atlarger distance l . We observe onset of the exponentialdecay at smaller lengths l for larger α . Near α = 1 wefind only a power law decay. This is consistent with theappearance of a delocalized critical state. We extract the |n-m| -4 -3 -2 -1 S m y S ny α =1.0 α =1.2 α =1.4 α =1.6 α =1.8 α =2.0 α =2.4 α =2.8 -30 -20 -10 0 log(|S S |) P α =1.0 α =2.8 FIG. 18. (Left)Concurrence averaged over 1000 random sam-ples as a function of index distance n nm = | n − m | . (Right)Distribution of the concurrence at n nm == 10. Filling factoris fixed to N/L = 0 . n for a given random sample, where S = log 2. power of the concurrence function at α = 1 and obtain γ = 1 . 56, i.e. (cid:104) C n (cid:105) ∼ n − . , which is different fromthe result obtained at the IRFP where (cid:104) C n (cid:105) ∼ n − , asreviewed above. Note that the horizontal axis of distri-bution on the right panel is on a log-scale. As expected,the concurrence is widely distributed. It follows a log-normal distribution , its center is shifted by changing α ,and its width decreases with decreasing α . Entanglement entropy - Let us next consider the entan-glement entropy Eq. (46), which one can obtain directlywith DMRG from the entanglement spectrum. In Fig. 19,an exemplary result is presented. We see a distinct struc-ture of entanglement of the ground states for a givenrandom sample at α = 0 . n stands for the length of the subsystem. At α = 2, the entanglement entropy alternates between 0and S (= log 2) throughout the entire lattice, whichindicates that the ground state is to a good approxima-tion a random singlet state, a product of local singlets,whose length is narrowly distributed and thus localized,as expected for the SDFP RS state. On the other hand,at α = 0 . 8, the S vN is non-zero everywhere. It impliesthat the ground state has a finite degree of entanglementthroughout the system. This is consistent with a delo-calization transition at α = α c . The averaged entangle-6 n S v N e n s α =0.6 α =0.8 α =1.0 α =1.2 α =1.4 α =1.6 α =1.8 α =2.0 α =2.4 α =2.8 α c datafitting ~ (1/N) + 1.0097 (a)(b) FIG. 20. (a) Entanglement entropy averaged over 1000 ran-dom samples as a function of the subsystem size n . Here, (cid:104)· · · (cid:105) ens denotes the ensemble average. The red and blacksolid lines are fitting curves to Eqs. (51) and (48), respec-tively. (b) Finite size scaling of α c at which the ground stateis most entangled for a given N . The density of spins is fixedat N/L = 0 . ment entropy [ (cid:104) S vN (cid:105) ens ] over 1000 random samples with N = 80 is presented in Fig. 20 (a) as a function of thesubsystem size n . For 1 . < α ≤ . n , i.e.area law scaling. At smaller α the entanglement entropyincreases and we find good agreement near α = 1 withCardy’s formula Eq. (48) for a finite system with openboundary conditions.Fitting curves are displayed in Fig. 20 (a) asred [Eq.(51)] and black [Eq.(48)] solid lines. Here,the extracted central charge is c fit = 0 . 93, with aconstant k = 0 . 39. After the entanglement entropyreaches a maximum at α ≈ 1, it decreases again withdecreasing α , as seen for α = 0 . , . α c with maximum entanglement entropydepends on the system size N . In analyzing its systemsize dependence, Fig. 20, we find the critical exponent tobe α c = 1 . /N → α c is estimated to be 1 . ± . M = 1000 random samples.There are peaks at integer multiples of S = log 2, S vN /S vN0 P α =0.8 B=0.0000 B=0.0001 S vN /S vN0 P α =1.0 B=0.0000 B=0.0001 S vN /S vN0 P α =2.0 B=0.0000 B=0.0001 ( a )( b )( c ) FIG. 21. Distributions of the entanglement entropy for (a) α = 0 . 8, (b) α = 1 . α = 2 . 0, where 1000 randomsamples and the system size N = L/ 10 = 80 are used, and S = log 2. both with (red lines) and without magnetic field (blacklines). Note that the effect of the magnetic field on theGS depends strongly on α . The distributions of EE arefor α < α = 1, [Fig. 21, (b)], the distribu-tion without magnetic field is almost identical to the oneat ( α, B ) = (0 . , S and 2 S are low-ered while the one at 0 is enhanced by the field. Thisindicates the emergence of free spins with no entangle-ment entropy. We see that the enhancement of probabil-ity at S vN = 0, the density of free spins increases furtherwith α see Fig. 21 (c)]. Thus, the DMRG shows thatthe entanglement entropy decreaes as α is increased be-yond the critical value α ≈ α to detect any increase with α , we can check if we cansee the increase of the entanglement entropy when theinteraction is cutoff exponentially, and the IRFP withcritical entanglement entropy is recovered. In Fig. 22the averaged entanglement entropy as a function of thecorrelation length is shown as obtained with DMRG. In-deed, while initially a decrease with decreasing correla-tion length ξ is observed at a ξ of the order of twice the7 ξ S v N e n s FIG. 22. Averaged entanglement entropy at l = 40 as a func-tion of the correlation length in units of the lattice spacing.The red solid line denotes the exact EE S vN = ln 26 log (cid:39) . . lattice spacing, we observe an increase of the EE again,in agreement with the approach to the IRFP. Susceptibility. - Let us consider next the magnetic sus-ceptibility of the ground state under a weak magneticfield B : (cid:104) χ ( α ) (cid:105) ens = (cid:104) ∆ M ( α ) / ∆ B (cid:105) ens , where M =(1 /N ) (cid:80) Ni =1 (cid:104) S zi (cid:105) . Figure 23 (a) shows the result for thatsusceptibility as function of the system size N . Notethat it converges to a single curve as increasing thesystem size, and N = 80 is large enough to see thethermodynamic limit. The susceptibility shows a verysharp increase in 0 . < α ≤ . α . We believe that the emergence of free spins leads to the significant increase of χ . We alsoinvestigate the field dependence of susceptibility withfixed N = 80, and results are shown in Fig. 23 (b) and (c).These results strongly suggest α ∗ to be close to 1 which isconsistent with the result obtained above from the con-currence and the result obtained by SDRG in Ref. 12. VIII. COMPARISON BETWEEN EXACTDIAGONALIZATION AND DMRG We have seen in the previous section that the tensornetwork extension of the DMRG yields results for en-tanglement measures which are, at least for α > N = 12spins, distributed randomly in a chain of length L = 120.We have used ED and DMRG to calculate the ground-state wave function, entanglement entropy, and spin-spincorrelations among all pairs. Results are illustrated inFig. 24 for open boundary condition. Calculations withmatrix product state (MPS) optimization have been per-formed using the ITensor C++ library . We run enoughsweeps for the entropy to converge to at least 10 − , and (()-(-)) / B = 0.00005JB = 0.00010JB = 0.00020JB = 0.00040JB = 0.00080JB = 0.00160J -6 -4 -2 () e n s B = 0.00005JB = 0.00010JB = 0.00020JB = 0.00040JB = 0.00080JB = 0.00160J -6 -4 -2 () e n s N = 30N = 40N = 50N = 60N = 70N = 80 (a)(b)(c) FIG. 23. Averaged susceptibility of the ground state as afunction of α for (a) various system sizes N with fixed density N = L/ 10 and fixed magnetic field B = 0 . J . b) forvarious magnetic fields B for N = 80. (c) The derivativeof the logarithm of the averaged magnetic susceptibility withrespect to α , where Γ( α ) = log[ (cid:104) χ ( α ) (cid:105) ens ] . a large number of states, up to 1000, was kept so thatthe truncation error is less than 10 − . Regarding theimplementation of long-range interaction, as the systemis small, we used the AutoMPO method available in ITen-sor and input all of the terms connecting sites i and j , sohere we did not consider further approximation like fit-ting interactions to a sum of exponentials as done in theprevious section or more recent SVD compression ap-proaches . The ED results were obtained with the stan-dard ARPACK diagonalization routine as implementedin SciPy .For small power exponent α = 1 both methods are inagreement for all random realizations, both in the en-tanglement entropy and the spin-spin correlations mea-surement, as seen in Fig. 24 where each row correspondsto a specific sample. When increasing α , however, onecan see that the results of DMRG gradually deviate fromED. Particularly, focusing on the entanglement entropy,it is clear that DMRG converges to a state with lowerentanglement. As seen from the amplitudes of the manybody wave function in the lower Fig. 24 (black) this8is accompanied by a breaking of the particle-hole sym-metry. In fact, it is well known that the matrix prod-uct state Ansatz of DMRG tends to prefer states oflower entanglement, when states are close in energy. In-deed, although the states obtained with ED shown inthe lower Fig. 24 (yellow) are found to have the sameenergy, the ED ground state is particle-hole symmet-ric and more strongly entangled. Thus, this is evidencethat the DMRG omits some of the singlets formed atlong distances, which therefore tends to underestimatethe entanglement while changing the energy only by anamount smaller than the numerical accuracy. It has beenreported that extensions of tree tensor networks (TTN)can capture entanglement properties in disordered sys-tems better . IX. CONCLUSIONS We find that the strong disorder fixed point, char-acterized by a fixed point distribution of the couplingswith a finite dynamical exponent, describes disorderedquantum systems of long-range coupled antiferromag-netic spin chains consistently. However, the lowest-orderSDRG, with its RS state, is found not to be sufficient toobtain the typical value of the concurrence. We thereforeproposed and implemented a correction scheme to the RSstate, allowing us to obtain the leading order corrections.These corrections yield a power law distance scaling ofthe typical value of the concurrence, which we demon-strate both by a numerical implementation of these cor-rections and by an analytical derivation. They are foundto be in agreement with each other.The entanglement entropy (EE) is calculated usingSDRG numerically and analytically and found to be log-arithmically enhanced for all α , whereas the effectivecharge is found not to depend on α and to be c = ln 2,in agreement with an analytical derivation. However,the analytical derivation uses assumptions on the corre-lation between singlets, and in a first attempt to includethese correlations, we arrived at a smaller central charge.Therefore, a more rigorous derivation is called for, whichwe leave for future research.While we confirm with numerical exact diagonalization(ED) the logarithmic enhancement of EE and a weak de-pendence on α , it fits in a wide range of distances l acritical behavior with a central charge close to c = 1, rem-iniscent of the clean Haldane-Shastry model with powerlaw decaying interaction with α = 2. Indeed, also theconcurrence, derived with numerical ED is found to de-cay with a power law, whose exponent is smaller thanthe one found by SDRG, γ = 2 and closer to the oneknown for the Haldane-Shastry model, γ = 1 . How-ever, at small distances l (cid:28) L we find strong deviations,which may indicate that the central charge converges tothe SDRG value c = ln 2 < α disorder is relevant so that the system converges to theSDRG fixed point.We also present results obtained with DMRG and findagreement with ED for sufficiently small α < 2, while forlarger α DMRG is found to tend to underestimate theentanglement entropy and finds a faster decaying concur-rence. As it is known from previous studies that DMRGunderestimates Entanglement, extensions like the treetensor network have been suggested, which also mightallow to study larger system sizes. We note that it hasbeen previously suggested that a delocalization occur ata critical value of α c . As we find a logarithmic lengthdependence for all α , as expected at a critical point, wecannot discern the delocalization transition at a specific α c in the entanglement properties within this approach. X. ACKNOWLEDGEMENTS R.N.B. acknowledges support from DOE BES GrantNo. DE-SC0002140. S.H. acknowledges support fromDOE under Grant No. DE-FG02-05ER46240. Compu-tation for the work described in this paper was supportedby the University of Southern Californias Center forHigh-Performance Computing (hpc.usc.edu). S. K. ac-knowledges support from DFG KE-807/22. H.-Y.L. wassupported by a Korea University Grant and National Re-search Foundation of Korea (NRF- 2020R1I1A3074769).9 . . . . . S v N α = 1 . α = 2 . α = 3 . . . . . . S v N α = 1 . α = 2 . α = 3 . l . . . . . S v N α = 1 . . . . . l α = 2 . l α = 3 . − − − − S y n S y m − − − − S y n S y m | n − m | − − − − S y n S y m | n − m | | n − m | − . − . − . . . . . ψ g s α = 1 . E dmrggs = − . J E edgs = − . J α = 2 . E dmrggs = − . J E edgs = − . J α = 3 . E dmrggs = − . J E edgs = − . J − . − . − . . . . . ψ g s α = 1 . E dmrggs = − . J E edgs = − . J α = 2 . E dmrggs = − . J E edgs = − . J α = 3 . E dmrggs = − . J E edgs = − . J dim H − . − . − . . . . . ψ g s α = 1 . E dmrggs = − . J E edgs = − . J dim H α = 2 . E dmrggs = − . J E edgs = − . J dim H α = 3 . E dmrggs = − . J E edgs = − . J FIG. 24. Comparison between ED and DMRG for N = 12spins distributed randomly over length L = 120 with openboundary condition: From the M = 1000 realizations, weshow the entanglement entropy, the spin-spin correlation, andthe ground-state wave function of three samples In the lowerfigure the vertical red-line is plotted at the particle-hole point,as a guide to the eyes. Appendix A: Solution of the Master Equation at theFixed Point In this appendix we show how to derive the solutionof the master equation at the fixed point. We denote µ = .First, we multiply by (cid:82) z J dJ the two sides of Eq. (7),and then plug in the Ansatz for P ( J, Ω), Eq. (9), (cid:90) Ω0 dJz J Ω µ J µ − ( µ (Ω) − Ω µ (cid:48) µ − Ω µ (cid:48) µ ln (cid:18) J Ω (cid:19) )= µ Ω − µ − (cid:90) Ω0 dJ J µ − z J (cid:90) Ω0 dJ dJ ( J J ) µ − × (cid:20)(cid:90) Ω0 dJ J µ − J J J − J (cid:90) Ω0 dJ J µ − J J J − J (cid:21) (A1)The left hand side of Eq. (A1) can be integrated and ex-pressed in terms of hypergeometric functions. One can alsointegrate over J , J and J on the L.H.S., yielding M ( µ, µ + 1 , Ω ln( z )) (cid:18) µ Ω − µ (cid:48) µ (cid:19) + µ (cid:48) µ F ( µ, µ, µ + 1 , µ + 1 , Ω ln( z ))= µ Ω − µ − ( − ln( z )) − µ γ ( µ, − Ω ln( z )) (cid:90) Ω0 dJ dJ ( J J ) µ − (cid:2) ( J − J ) − µ γ ( µ, ( J − J ) ln( z )) γ ( µ, ( J − J ) ln( z )) (cid:3) (A2)Here γ , F and M , are respectively the lower incompletegamma function, the generalized hypergeometric function,and the confluent hypergeometric function, defined by γ ( s, x ) = x (cid:82) t s − e − t dtM ( µ, µ + 1 , x ) = ∞ (cid:80) n =0 µx n n !( µ + n )2 F ( µ, µ, µ + 1 , µ + 1 , x ) = ∞ (cid:80) n =0 µ x n n !( µ + n ) (A3)Using the identity γ ( a, − z ) = − ( z a /a ) M ( a, a + 1 , z ) and Eqs.(A3) one can rewrite Eq. (A1) as M ( µ, µ + 1 , Ω ln( z )) (cid:18) µ Ω − µ (cid:48) µ (cid:19) + µ (cid:48) µ F ( µ, µ, µ + 1 , µ + 1 , Ω ln( z ))= µ Ω − µ − ( − ln( z )) − µ γ ( µ, − Ω ln( z )) ∞ (cid:88) k =0 B ( µ, k + 1)(ln( z )) k (2 k )!( k + µ ) (cid:20)(cid:90) Ω0 dJ dJ ( J J ) µ − ( J − J ) k (cid:21) (A4)with B ( µ, k + 1) being the standard Beta function.Finally we evaluate the last double integral, (cid:90) Ω0 (cid:90) Ω0 dJ dJ ( J J ) µ − ( J − J ) k = Ω µ Ω k k + µ Γ( µ )Γ(2 k + 1)Γ( µ + 2 k + 1)(A5) Allowing us to get the integrated form of the constraint equa-tion, M ( µ, µ + 1 , Ω ln( z )) (cid:18) µ Ω − µ (cid:48) µ (cid:19) + µ (cid:48) µ F ( µ, µ, µ + 1 , µ + 1 , Ω ln( z ))= µ Ω M ( µ, µ + 1 , Ω ln( z )) ∞ (cid:88) k =0 Γ ( µ )Γ(2 k + 1)(Ω ln( z )) k ( k + µ ) Γ ( µ + 2 k + 1) (A6)Which can be rewritten using Eq. (A3) in the form − µ (cid:48) (Ω) ∞ (cid:88) n =1 (Ω t ) n ( n − µ + n ) = µ (Ω)Ω × ∞ (cid:88) k =1 ,k (cid:48) =0 Γ ( µ )(2 k )!(Ω t ) k + k (cid:48) ( k + µ ) Γ ( µ + 2 k + 1) k (cid:48) !( µ + k (cid:48) ) (A7)Here we defined t = ln( z ). Appendix B: Benchmark Model Results1. Haldane-Shastry model Among the spin models in one dimension, the Haldane-Shastry chain is interesting for several reasons. It is anantiferromagnet with 1 /r exchange interactions, and it pos-sesses a Yangian symmetry which makes it integrable, there-fore, exactly solvable. This model is defined by H = J π N (cid:88) j
2. Random Heisenberg XX-model In this section, we implement the exact diagonalizationprocedure for the random short ranged AFM XX-Heisenbergmodel, defined by its Hamiltonian : H = N − (cid:88) i =1 J i (cid:0) S xi S xi + S yi S yi +1 (cid:1) (B3)Where { J i } N − i =1 are uncorrelated positive random variables,drawn from a distribution P ( J ). In Ref. 9, it was shown withSDRG method that, given an interval of length l embeddedin the infinite line, the average entanglement entropy of thisinterval, with the rest of the chain scales for large l as Eq. n . . . . . . . . S v N ¯ c = 1 , k = 0 . FIG. 25. Entanglement entropy obtained via exact diagonal-ization for the Haldane-Shastry model with N = 22 sites withperiodic boundary conditions. The solid line corresponds tothe Cardy law Eq. (48) with a central charge ˜ c = 1 . b = 2and k (cid:48) = 0 . b = 2, corresponding to the entanglement entropyof a critical system Eq (47), with an effective central charge˜ c = c × ln(2), where c = 1 is the central charge of the pure XX-Heisenberg model. Fig. (26) shows ED results for a samplewith N = 22 spins with open boundary conditions for theentanglement entropy averaged over 200 random samples asfunction of partition size n . A strong even-odd effect is seen.The yellow line is the Cardy law Eq. (48) with ˜ c = 0 . b = 1 and k (cid:48) = 0 . 68, which is in good agreement with theresult obtained for a RS state Eq. (51), when c = 0 . 8. Thepink line is Eq. (48) with ˜ c = 2 ln 2, b = 1 and k (cid:48) = 0 . N = 22 spins with periodicboundary conditions. ED reproduces the results obtained an-alytically in Ref. 9, and derived in this article, Eq. (51) for b = 2, when c = 1. Appendix C: Entanglement Entropy WithCorrections to the Random Singlet State. Having the RS state with corrections | ψ (cid:105) , Eq. (30) we canwrite the density matrix to calculate the entanglement en-tropy beyond the RSS using directly the definition of the vonNeumann entropy Eq. (46). Given that the RS state withcorrections | ψ (cid:105) is not a product state, there are no simplecombinatorical arguments, such as the counting of crossingsinglets, since the entropy of a superposition state is not thesum of the individual entropies. Moreover, a closed formulausing the definition Eq. (46) is not feasible due to the depen-dence of the sums on the specific realization of the RSS, whichmakes taking the partial trace inconceivable without consid-ering every single possible scenario, i.e., there are as manyoutcomes for the partial trace as there are possible randomsinglet states on the chain ( ∼ N ).A possible solution to this problem is to start by takinginto account in Eq. (30) the term with the largest coefficient n . . . . . . . . ¯ S v N ¯ c = 0 . , k = 0 . c = 2 . , k = 0 . FIG. 26. ED results for a sample with N = 22 spins with openboundary conditions for the entanglement entropy averagedover 200 random samples as function of partition size n . Theyellow line is the Cardy law Eq. (48) with ˜ c = 0 . b = 1and k (cid:48) = 0 . 68. The pink line is Eq. (48) with ˜ c = 2 ln 2, b = 1and k (cid:48) = 0 . n . . . . . . . ¯ S v N ¯ c = 1 . , k = 0 . FIG. 27. The ED results for the random nearest-neighbormodel with N = 22 sites on periodic chain. Data were aver-aged over 1000 random sample. The solid line corresponds tothe Cardy law Eq. (48) with a central charge ˜ c = 1 . b = 2 and k (cid:48) = 0 . | ψ (cid:105)≈ c | ψ RS (cid:105) + cδJ ( | + nl (cid:105)|− mk (cid:105) + |− nl (cid:105)| + mk (cid:105) ) (cid:79) { ij }(cid:54) = { nl }(cid:54) = { mk } | ij (cid:105) , (C1)where δJ = J nm + J lk − J nk − J ml J nl + J mk , (C2) is the maximum coefficient appearing in Eq. (30), and thecoefficient c ψ needs to be redefined as c ψ = 1 √ δJ , (C3)in order to keep the approximated state properly normalized.Now, let us consider a situation, where the RSS state is suchthat k singlets cross the partition boundary, giving the EE S ( k ) = k ln 2 . (C4)With the corrected state Eq. (C1) we find after a lengthy butstraightforward calculation, that it is possible to arrive at aconditional closed form for the entanglement entropy that de-pends where the two converted singlet pairs { nl } and { mk } are located, relative to the partition boundary. There arethree distinguishing cases that give rise to different expres-sions for the entanglement entropy as a function of the numberof crossing singlets and triplets k ( Here we set c = c ψ in Eq.(C3), which in the limit of no corrections ( δJ → , c → Case 1: Each of the two converted singlets { nl } and { mk } are at opposite sides of the cut and none of them cross theboundary: S ( k ) = − c ln (cid:18) c k (cid:19) − (1 − c ) ln (cid:18) − c k +1 (cid:19) . (C5) Case 2: Both converted singlets cross the boundary be-tween subsystems, for k ≥ S ( k ) = − c (cid:18) c k (cid:19) − c δJ + 1) ln (cid:18) c k (2 δJ + 1) (cid:19) − c δJ − ln (cid:18) c k (2 δJ − (cid:19) . (C6) Case 3: Any other relationship between the convertedpairs and the boundary, e.g. , both pairs are part of the samesubsystem or only one of them crosses the boundary. In thiscase, the approximated state brings no correction to the en-tropy, giving the same value obtained at the IRFP , Eq. (C4).Moreover, as seen in Fig. 28(a) for the specific instance k = 3, case 1 (Eq. (C5), dashed line) gives a higher entropythan the one at the IRFP (Eq. (C4), continuous line). This isexpected since the corrected state is a superposition of statesthat differ only in spin pairs { nl } and { mk } , which live on op-posite sides of the subsystem boundary, and therefore resultsin an enhancement of the quantum correlations between sub-chains. On the other hand, case 2 (Eq. (C6), dashed-dottedline) gives a lower entropy than that of Eq. (C4), also for allvalues of c (cid:54) = 1. Again, this is expected due to the fact thatthe extra correction terms are destroying the RSS, which inthis case is the maximally entangled state, given that bothpairs cross the boundary. It is worth noting, that in order toplot the entropy in Eq. (C6), Eq. (C3) was inverted in orderto obtain δJ ( c ), and the positive root was chosen. However,since Eq. (C6) is an even function of δJ , this choice becomestrivial.We observe that, even though the plot is only shown for thespecific case of k = 3 crossing singlets, the above statementsremain true for all values of k , as can be inspected via Eqs.(C5) and (C6). c S ( k = ) Case 1Case 2Case 3 (a) l › S fi α = α = α = α = α = α = (b) FIG. 28. (a) Entanglement entropy as function of c for thethree instances that occur after approximating the correctedstate to Eq. (C1). We note that for 0 < c (cid:46) . 5, pertur-bation theory is no longer valid and a different behavior inthis regime might occur. (b)Average block entanglement en-tropy calculated with the approximated state in Eq. (C1) fora chain with N = 800 spins. The power α is varied from 0 . . 00. Calabrese and Cardy’s formula (Eq. (48)) is plottedalong for reference as a black continuous line.As seen in Fig. 28(b) and (c) for a chain of length N = 800, the approximation in Eq. (C1) that gives rise toEqs. (C5,C6,C4) does not give a significant dependence on thepower α , and the entropy remains close to Cardy’s result. Bylooking at the difference of the entropy calculated with cor-rections and the one calculated solely with the RSS, we findthat they are about two orders of magnitude smaller than therespective entropy values, which is not surprising, since thecorrections to only two singlets are taken into account so far.Therefore, we can conclude that sases 1 and 2, the twocases in which corrections appear, are not frequent enoughthroughout realizations to notably affect the average entropy.In conclusion, even though the corrected state in Eq. (30)is useful to calculate the typical concurrence, it does not givea sizable correction to the average entanglement entropy gov-erned by the RSS. Next, we would have to find a way toinclude the corrections to the EE from all singlet-triplet ex-citations by taking recursively weaker and weaker correctionsinto account, which we leave for future research. ∗ [email protected] † [email protected] ‡ [email protected] § [email protected] ¶ [email protected] ∗∗ [email protected]; ***(2019-2020) †† [email protected] ‡‡ [email protected] P. W. Anderson, Phys. Rev. , 1492 (1958). N. F. Mott, Le Journal de Physique Colloques , 301(1976). C. C. Yu and A. J. Leggett, Comments Condens. MatterPhys. , 231 (1988). T. Graß and M. Lewenstein, EPJ Quantum Technology ,8 (2014). R. N. Bhatt and P. A. Lee, Phys. Rev. Lett. , 344 (1982). D. S. Fisher, Phys. Rev. B , 3799 (1994). F. Igl´oi and C. Monthus, Phys. Rep. , 277 (2005). P. Calabrese and J. Cardy, Journal of Statistical Mechan-ics: Theory and Experiment , P06002 (2004). G. Refael and J. E. Moore, Phys. Rev. Lett. , 260602(2004). J. A. Hoyos, A. P. Vieira, N. Laflorencie, and E. Miranda,Phys. Rev. B , 174425 (2007). N. Moure, S. Haas, and S. Kettemann, EPL (EurophysicsLetters) , 27003 (2015). N. Moure, H.-Y. Lee, S. Haas, R. N. Bhatt, and S. Ket-temann, Phys. Rev. B , 014206 (2018). A. D. Mirlin, Y. V. Fyodorov, F.-M. Dittes, J. Quezada,and T. H. Seligman, Phys. Rev. E , 3221 (1996). E. Cuevas, EPL (Europhysics Letters) , 84 (2004). F. D. M. Haldane, Phys. Rev. Lett. , 635 (1988). W.-J. Rao, X. Wan, and G.-M. Zhang, Phys. Rev. B ,075151 (2014). F. Igl´oi and C. Monthus, The European Physical JournalB , 290 (2018). E. Westerberg, A. Furusaki, M. Sigrist, and P. A. Lee,Phys. Rev. B , 12578 (1997). D. S. Fisher, Phys. Rev. Lett. , 534 (1992). S. Hill and W. K. Wootters, Phys. Rev. Lett. , 5022(1997). In Ref. the spin correlation function instead of the con-currence is calculated, explaining the difference in by afactor S ( S + 1) = 3 / This notation will be used throughout to simplify the longexpressions involved. Note that with the given definition of c , the concurrencefor { pq } ∈ RS is always positive. G. Vidal, J. I. Latorre, E. Rico, and A. Kitaev, Phys. Rev.Lett. , 227902 (2003). C. Holzhey, F. Larsen, and F. Wilczek, Nuclear Physics B , 443 (1994). S. R. White, Phys. Rev. Lett. , 2863 (1992). S. R. White, Phys. Rev. B , 10345 (1993). U. Schollw¨ock, Rev. Mod. Phys. , 259 (2005). S. R. White, Phys. Rev. B , 180403 (2005). D. P´erez Garc´ıa, F. Verstraete, M. M. Wolf, and J. I.Cirac, Quantum Information Computation , 401 (2007). U. Schollw¨ock, Annals of Physics , 96 (2011), january2011 Special Issue. R. Or´us, Annals of Physics , 117 (2014). B. Pirvu, V. Murg, J. I. Cirac, and F. Verstraete, NewJournal of Physics , 025012 (2010). ITensor Library(version 3.1.13) . E. M. Stoudenmire and S. R. White, Phys. Rev. Lett. ,046401 (2017). Scipy (version 1.2.1) . A. M. Goldsborough and G. Evenbly, Phys. Rev. B ,155136 (2017). J. I. Cirac and G. Sierra, Phys. Rev. B81