Probability distribution of the entanglement across a cut at an infinite-randomness fixed point
PProbability distribution of the entanglement across a cut at an infinite-randomnessfixed point
Trithep Devakul , Satya N. Majumdar and David A. Huse Department of Physics, Princeton University, NJ 08544, USA Laboratoire de Physique Theorique et Mod´eles Statistiques, Universit´e Paris-Sud, Orsay, France.
We calculate the probability distribution of entanglement entropy S across a cut of a finite one-dimensional spin chain of length L at an infinite-randomness fixed point using Fisher’s strong ran-domness renormalization group (RG). Using the random transverse-field Ising model as an example,the distribution is shown to take the form p ( S | L ) ∼ L − ψ ( k ) , where k ≡ S/ log [ L/L ], the largedeviation function ψ ( k ) is found explicitly, and L is a nonuniversal microscopic length. We discussthe implications of such a distribution on numerical techniques that rely on entanglement, such asmatrix product state (MPS) based techniques. Our results are verified with numerical RG simula-tions, as well as the actual entanglement entropy distribution for the random transverse-field Isingmodel which we calculate for large L via a mapping to Majorana fermions. I. INTRODUCTION
Entanglement has emerged as a key ingredient in thestudy of quantum many-body systems. In particular, thematrix product state (MPS) or density matrix renormal-ization group (DMRG) methods for numerically studyingquantum states work directly with the bipartite entan-glement of a pure many-body quantum state across each“cut”. In a disordered system, this entanglement mayvary widely between different locations of the cut. Thusit is of interest to understand the probability distributionof the entanglement for such states. We are particularlyinterested here in highly-excited eigenstates that are ofinterest in studies of many-body localization (MBL).At quantum criticality, the ground-state entanglementof certain translationally-invariant systems is known todisplay universal behavior characteristic of the associatedconformal field theory. In one spatial dimension, the vonNeumann entropy (in bits) for a subsystem of length L at such a conformally-invariant quantum critical pointscales for large L as S = − Tr { ρ log ρ } = c L + const , (1)where c is the central charge and ρ is the reduced densityoperator for the subsystem.A similar scaling of entanglement as the logarithm of L applies for a number of disordered one-dimensionalsystems whose ground states are governed by infinite-randomness fixed points of strong-randomness renormal-ization groups . This class of systems include the disor-dered Heisenberg and XXZ spin chains, and the trans-verse field Ising model (TFIM) . For the disorderedTFIM, as long as it can be transformed to noninteract-ing fermions, this scaling of the entanglement also appliesto all excited eigenstates .With the surge of interest in the many body localiza-tion (MBL) phase transition, which might be describedby an infinite randomness fixed point , a deeper under-standing of this entanglement scaling becomes even morerelevant. Also, critical points within the MBL phase, should they exist, are governed in one-dimensional sys-tems by infinite-randomness fixed points and exhibit thislogarithmic scaling of the entanglement (Eq 1) with aneffective (irrational) central charge .The practical importance of quantum entanglement ishighlighted by the success of the density matrix renor-malization group (DMRG) algorithm, which relies onthe low entanglement nature of certain quantum statesto accurately represent them as matrix product states(MPS) . There has been progress in obtaining highly ex-cited eigenstates of MBL systems using this MPS frame-work . A variant of DMRG called DMRG-X has beendeveloped to treat highly-excited states, with the specificgoal of treating MBL systems . The accuracy of thesealgorithms depends on the bond dimension χ allowed oneach of the internal bonds of the MPS. Representing astate to a certain accuracy requires the underlying MPSto have a bond dimension χ on a bond which will gen-erally scale exponentially with the entanglement entropy χ ∼ e aS , where a depends on the structure of the en-tanglement spectrum (eigenvalues of the reduced densitymatrix). The logarithmic growth of S with L thus leadsto a polynomial growth in the necessary χ and there-fore of the computational complexity of the DMRG. Ina disordered system, S and thus the necessary bond di-mension χ varies along the spin chain, so to evaluate thescaling of the DMRG algorithm’s computation time weneed to understand the distribution of the entanglement.In this paper, we derive an explicit expression for theprobability distribution of the entanglement entropy atthe infinite randomness fixed point of the critical ran-dom TFIM for a finite open system, like those studiedby DMRG. From the perspective of Fisher’s strong dis-order renormalization group (RG) for the random crit-ical TFIM , there is no difference between a ground orexcited state , and so our results may also be applied toDMRG-X calculations. This expression is used to derivethe scaling behavior and probability of finding entangle-ments from various points in the distribution. A find-ing is that the computational time is not dominated atlarge sample length L by the local maximum of the en- a r X i v : . [ c ond - m a t . s t a t - m ec h ] D ec tanglement. We check our findings with numerical sim-ulations of the RG itself, as well as calculations on thedisordered transverse-field Ising model, and find excellentagreement. II. FISHER’S RENORMALIZATION GROUP
We first begin with a rough overview of Fisher’s RGanalysis for the random transverse-field Ising chainHamiltonian of the form H = (cid:80) i J i σ xi σ xi +1 − (cid:80) i h i σ zi .At criticality, J i and h i are independent random cou-plings drawn from the same distribution (which, for con-venience, we assume has equal weights for positive andnegative couplings). By performing a Jordan-Wignertransformation in the σ z basis, this Hamiltonian canbe expressed in terms of the conventional Majoranafermions γ j , with only nearest neighbor couplings: H = i (cid:80) L − j =1 g j γ j γ j +1 , with g i − = h i and g i = J i .The RG proceeds by always treating the strongest cou-pling, in exactly the same way as is done for a randomsinglet phase. The energy scale Ω is defined to be thestrongest of all the couplings Ω = max j | g j | , beginningat Ω = Ω for the unrenormalized bare model. An RG“step” begins by finding the coupling with | g j | = Ω.Since the distribution of the g ’s is very broad, we willalmost always have | g j − | (cid:28) | g j | (cid:29) | g j +1 | . As a result,the two Majoranas γ j and γ j +1 form a two-level systemwith eigenenergies ± g j that is only weakly coupled toneighboring Majoranas. We put this two-level systemin one of its local eigenstates (Fisher always chose theground state, since that is what he was focused on) andthen treat the coupling to its neighbors perturbatively.This results, from leading order in degenerate perturba-tion theory, in an effective coupling g (cid:48) = ± g i − g i +1 /g i between the neighboring Majoranas γ j − and γ j +2 , withthe sign of the new coupling depending on whether theground or excited state was chosen. After each such deci-mation, the energy scale Ω is decreased accordingly. Notethat the entanglement structure and the magnitudes ofthe renormalized couplings are independent of whetherthe ground or excited state was chosen. At this quantumcritical point, this RG flows to “infinite randomness”,meaning the probability distribution of log | g | becomesarbitrarily broad, and the approximation of keeping onlythe leading-order perturbative coupling becomes asymp-totically exact .It is convenient to use the scaled log couplings β i =log(Ω /g i ) ≥ / Ω),which is the RG flow parameter. After a decimation,the new β is then simply given by β (cid:48) = β i − + β i +1 ,greatly simplifying the flow equations. Solving the flowequations results in the fixed point distribution p ( β/ Γ) = e − β/ Γ . (2)In the original spin language, a field h i being “inte-grated out” corresponds to the local state of the renor-malized spin ˜ σ xi being fixed to |→(cid:105) or |←(cid:105) , which →
2) resulting in the two spins on eitherside being fully correlated, but not yet entangled. The cutnow lies “within” that combined spin. When the field on thissite is decimated (2 → S across the cut is half the number N of decimations across it, S = N/ more microscopically represents an equal-amplitude lin-ear combination of some particular pattern of the bare σ zi ’s within the cluster and its opposite under flipping allspins. This step thus introduces one bit of entanglementacross any cut within this cluster. A bond J i being “in-tegrated out” corresponds to the two renormalized spins˜ σ zi and ˜ σ zi +1 being fixed to be either parallel or antipar-allel and thus combined to be one renormalized spin, soit introduces new microscopic correlations but no newentanglement.In this picture, an entanglement cut obtains one newbit of entanglement whenever a field h i is decimatedacross the cut, which is precisely half of all decimations.The decimations across the cut alternate between h i ’sand J i ’s: When the cut is between two renormalized spinsit is on a bond J i . When this J i is decimated these twospins are fully correlated and the cut is then within thenew renormalized spin and thus “on” a field h i . Whenthat h i is decimated the entanglement increases by onebit and the cut returns to being on a bond between tworenormalized spins. This process is shown in Fig 1. Thuswhen the number of decimations across a cut is N , theentanglement across that cut is S = N/ III. DERIVATION OF THE ENTANGLEMENTDISTRIBUTION
To obtain the probability distribution for the the num-ber of decimations N across a single cut for a system offixed total length L , we first instead consider the problemof having a system at a fixed log energy cutoff Γ. Thisis significantly simplified due to the fact that the “steps”in log Γ between two decimations are, at the fixed point,independent and identically distributed (i.i.d.) events .Suppose on running this RG the most “recent” deci-mation across our cut was at log energy cutoff Γ . Oncethe cutoff has increased to a higher Γ, what is the prob-ability R (Γ , Γ ) that no more decimations across our cuthave happened? At the critical fixed point, this proba-bility depends only on (cid:96) ≡ log Γ / Γ , so R (Γ , Γ ) = R ( (cid:96) ).Solving the flow equations , one finds R ( (cid:96) ) = (cid:32) √ √ e − −√ (cid:96) − − √ √ e − √ (cid:96) (cid:33) . (3)And therefore, the probability r ( (cid:96) ) that the next decima-tion occurs at (cid:96) is given by r ( (cid:96) ) = − dR ( (cid:96) ) d(cid:96) = 1 √ (cid:16) e − −√ (cid:96) − e − √ (cid:96) (cid:17) . (4)Now we may proceed towards the distribution of totaldecimations N between Γ and a final cutoff Γ. We beginwith a decimation that occurred initially at Γ , followedby N decimations happening at { Γ , Γ , . . . , Γ N } , withall Γ i − < Γ i < Γ. Taking advantage of the i.i.d. natureof the steps in log Γ, we define (cid:96) n = log Γ n / Γ , and (cid:96) =log Γ / Γ , so that the probability of having N decimationsis obtained by integrating over all possible { (cid:96) , . . . , (cid:96) N } with their respective probabilities, p ( N | (cid:96) ) = (cid:90) (cid:34) N (cid:89) n =1 r ( (cid:96) n − (cid:96) n − ) d(cid:96) n (cid:35) R ( (cid:96) − (cid:96) N ) . (5)This integral can be expressed as a convolution in thevariables ∆ (cid:96) n = (cid:96) n − (cid:96) n − , which we can then apply themethod of Laplace transforms to, resulting in˜ p ( N | a ) = (cid:90) ∞ p ( N | (cid:96) ) e − (cid:96)a d(cid:96) = [˜ r ( a )] N ˜ R ( a ) (6)where ˜ · denotes the Laplace transformation, and a is theLaplace conjugate variable. ˜ r ( a ) can be calculated ex-plicity from Eq. 4 to be ˜ r ( a ) = (1 + 3 a + a ) − .To invert the Laplace transform, we employ theBromwitch integral p ( N | (cid:96) ) = 12 πi (cid:90) γ + i ∞ γ − i ∞ da e a(cid:96) ˜ p ( N | a ) (7)= 12 πi (cid:90) γ + i ∞ γ − i ∞ da ˜ R ( a ) exp [ (cid:96)H ( a ; N/(cid:96) )] (8)where H ( a ; x ) = a + x log ˜ r ( a ) (9)and γ is chosen such that the contour of integration isto the right of all singularities of the integrand in thecomplex plane. In the limit of large (cid:96) , with N/(cid:96) kept fi-nite, this integral can be well approximated by the saddle x φ ( x ) FIG. 2: Plot of the function φ ( x ) (Eq. 11). point method. By analyticity of H , the saddle point for x > a axis at the minimum of H ( a ; x ) tothe right of both of its singularities. Then, up to poly-nomial prefactors, we have that p ( N | (cid:96) ) ∼ exp [ − (cid:96)φ ( N/(cid:96) )] (10)where φ ( x ) = − min a H ( a ; x ).Minimizing H ( a ; x ) we get φ ( x ) = 3 − √ x − x x log (cid:104) x + x (cid:112) x (cid:105) , (11)which is shown in Fig 2. This function has a finite valueat x = 0 of φ (0) = (3 − √ /
2, reaches a minimum at φ (1 /
3) = 0, and goes as φ ( x ) ≈ x (log(2 x ) −
1) + 3 / O (1 /x ) for large x .The connection to a system of fixed size can now bedone by relating L = Γ , up to some proportionalityconstant (which we set to 1 for now as it does not affectany of our immediate conclusions). This works becausein the large (cid:96) limit, the fluctuations of log lengths areof order one: log L = 2 log Γ + O (1). Thus, after thesubstitution (cid:96) = log L , for the critical random TFIMthe distribution of the bipartite entanglement S = N/ p ( S | L ) ∼ L − ψ ( k ) (12)for a finite sample of length L with open ends, with ψ ( k ) = φ (4 k ) and k ≡ S/ log L .There are a few interesting regions in this function thatdeserve mentioning. The fraction of cuts with zero en-tanglement is non-zero but vanishingly small, scaling as L − θ , with θ = ψ (0) = ∼ = 0 . k = 1 / ∼ = 0 . S ≈ log L for a single cut . Also ofinterest is the typical largest entanglement cut to appearin a sample. This is the S which appears with probabil-ity scaling p ( S | L ) ∼ L − , of which there are O (1) of ineach sample. This happens when ψ ( k ) = 1, which is at k ∼ = 0 . k .The cumulant generating function for this distribution g ( t | (cid:96) ) = log (cid:104) e Nt (cid:105) has been obtained analytically previ-ously . To see the relation between these two results,one can express in the limit of large (cid:96) (using the samesaddle point approximation), g ( t | (cid:96) ) = log (cid:90) ∞ p ( N | (cid:96) ) e Nt dN (13) ≈ − (cid:96) min x { φ ( x ) − tx } + . . . (14)= − − √ e t (cid:96) + . . . (15)in agreement with Ref. 17 at large (cid:96) up to additive loga-rithmic corrections. Thus, the connection between φ ( x )and g ( t | (cid:96) ) is via a Legendre transform. IV. IMPLICATIONS FOR DMRG
What does this distribution of entanglement mean fornumerical techniques such as DMRG that rely on en-tanglement? To accurately represent a state as a ma-trix product state (MPS), the number of states kept (orbond dimension) χ is related to the entanglement acrossa cut on that bond. Having a high bond dimension al-lows DMRG to capture states more accurately, but atthe cost of increased computation time. At the infiniterandomness fixed point, the entanglement comes only inthe form of S maximally entangled pairs, so the entan-glement spectrum is therefore simply 2 S equal nonzeroeigenvalues, followed by zeros, for which a bond dimen-sion of χ = 2 S is needed. Using an algorithm that allotsbond dimension independently for each bond as neededto attain a certain accuracy, we can ask the question ofhow the distribution of entanglement affects the scalingof the computation time for such a numerical technique.The DMRG-X algorithm for finding highly excitedMBL eigenstates relies on the diagonalization of a d χ × d χ effective Hamiltonian, where d = 2 is thelocal physical number of degrees of freedom. This diago-nalization can be done in full, or using alternatives suchas shift-invert Lanczos. The computational complexityof these algorithms scale as some power c of the bonddimension. For full exact diagonalization, c = 6, while c would be smaller for methods that scale more favorablyat large bond dimension. Letting T denote the compu-tation time, we have T ∼ χ c ∼ L ck log 2 , where we havedefined k ≡ S/ log L .We can now ask what the scaling of the mean compu-tational time (cid:104) T ( L ) (cid:105) is with system size L . The averagecomputation time per bond scales as (cid:104) T ( L ) (cid:105) ∼ (cid:90) ∞ p ( S | L ) L ck log 2 dS . (16) Up to logarithmic corrections, for large L (cid:104) T ( L ) (cid:105) ∼ L max k [ ck log 2 − ψ ( k )] . (17)For full exact diagonalization, c = 6 and one finds thetime to be dominated by bonds with k ∗ = 2 / √ ∼ =0 . (cid:104) T ( L ) (cid:105) ∼ L α , where α = ck ∗ log 2 − ψ ( k ∗ ) ∼ = 0 . L α . These powers willonly be smaller if a method that scales more favorablythan full exact diagonalization can be used. This ignoresthe strong correlation of entanglement entropy betweenbonds within a sample, which are important but don’taffect these exponents for the mean time.The maximum entanglement bonds were found ear-lier to occur near k ∼ = 0 . V. NUMERICAL RESULTS
To check our analytic findings, we numerically run theinfinite randomness RG on many samples. The systemis treated as an array of blocks which are the bonds be-tween the Majorana modes in the renormalized model.Each block i has a coupling strength g i , a length l i anda normalized distribution p i ( N ) for the number of deci-mations N across all internal cuts within the block. Theblocks are initialized with a coupling g i from the fixedpoint distribution (Eq 2) with Γ = 1, an initial length l i = 1, and a trivial initial internal distribution of N of p i ( N ) = δ N, . Upon decimation of a block i , the blocks i − i , and i + 1 are merged to a single block, and theprobability distribution is updated as p (cid:48) ( N ) = l i − p i − ( N ) + l i p i ( N −
1) + l i +1 p i +1 ( N ) l i − + l i + l i +1 . (18)The coupling strength is updated as prescribed by theRG rules, and the new length is simply the summation.To obtain the number of decimations given a fixed flowinterval (cid:96) = log Γ / Γ , we simply run the RG until thelog energy cutoff has reached the desired value. Then,the probability distribution can be sampled from all theremaining blocks in the system weighted by their length.The probability distribution for a system of fixed size L ,on the other hand, can be obtained by initializing thesystem with an odd number of blocks L (each of initiallength l = 1) and simply running the RG until only oneblock remains, which is guaranteed to have length L . -2 -4 -6 p ( N | l og Γ ) Γ = 20Γ = 40Γ = 100 N -2 -4 -6 p ( N | l og L ) L = 10 + 1 L = 10 + 1 L = 10 + 1 FIG. 3: Probability distribution of N from numerical RGsimulations obtained by (top) running the RG until a fixedlog energy cutoff Γ, and (bottom) by running the RG to com-pletion for a fixed system size L . Dashed lines show the nor-malized analytical results (Eq 10) with (cid:96) = log Γ in the fixedΓ case, and (cid:96) = log L in the case with fixed L . There isexcellent agreement between the analytical and numerical re-sults. Figure 3 shows the numerical results for fixed Γ, andfixed L , along with the (normalized) analytical predic-tions from Eq 10. In the case of fixed L , the substitution (cid:96) = log Γ is used, taking advantage of the relationshipbetween the mean L and Γ. In both cases, the agreementbetween numerical results and the approximate analyti-cal expression is excellent. VI. TRANSVERSE FIELD ISING MODEL
Finally, we compare with results for the random trans-verse field Ising model, which can be obtained by a map-ping to noninteracting majorana fermions and the entan-glement entropy obtained from the single particle corre-lation matrix . Due to self duality, this system is criticalwhen the distribution of J and h are the same. We pickthe disorder distribution to be the fixed point distribu-tion Eq 2. Since we can always rescale the energy, weset Ω = 1 so the cumulative distribution for J is givenby P ( J < j ) = j / Γ for 0 ≤ j ≤
1, and similarly for theon-site fields h . Γ is the parameter that flows towardsinfinite randomness, so we expect the Fisher RG to holdfor large Γ. At small Γ, the distribution becomes moreand more narrow arround J = h = 1, and so realisticallyour results cannot apply for very small Γ. On the otherhand, very high Γ runs into problems with machine preci-sion at large system sizes. We will focus on Γ = 1, which S/‘ l og [ p ( S | ‘ ) ] / ‘ L = 20 L = 50 L = 100 L = 200 L = 400 L = 800 FIG. 4: Distribution of entanglement entropy S sampled fromcuts of random excited eigenstates of the random transversefield Ising model of varying L . The random fields and cou-plings are taken from the fixed point distribution (Eq 2)with Γ = 1, corresponding to a flat distribution. With (cid:96) = log L/L for L ≈ . × − , all these curves col-lapse well towards a single curve. The dashed line shows theanalytic prediction (Eq 19) with a fixed constant shift to alignwith the numerical data. corresponds to a flat disorder distribution as is commonlyused in studies of localization (cf. the Anderson model).In comparing with actual microscopic models, there isa proportionality constant L , Γ / Γ = (cid:112) L/L , whichdepends on the microscopic details. There is therefore afree parameter L in defining (cid:96) = log L/L . However,log [ p ( S | (cid:96) )] /(cid:96) = − φ (2 S/(cid:96) ) + const( (cid:96) ) (19)is a universal function of
S/(cid:96) for all L (with the constantproviding the normalization having some dependence on L ).Figure 4 shows the entanglement distribution for theTFIM with Γ = 1 in this manner. For all L , the entan-glement distribution shows a peak at an entanglementof S = 1 bit, and collapses very nicely to a single curvein the tail, in very good agreement with the analyticalcurve. Notice that the fitted L is typically much smallerthan 1, this is a result of the fact that for L (cid:28) ξ the lo-calization length, entanglement will grow much quickerthan log L , and so at large L , S ∼ log L/L will appearto have a small L . VII. CONCLUSION
We have successfully obtained, using Fisher’s strongrandomness RG, an analytic expression for the prob-ability distribution of entanglement entropy for a ran-dom TFIM chain at the infinite randomness fixed point.The distribution is found to have a large deviation form p ( S | L ) ∼ L − ψ ( S/ log[ L/L ]) with ψ ( x ) given explicitly. Al-though the results were obtained for a system at a finitelog energy cutoff Γ, they can equally be applied to themore applicable case of a finite length chain by the sub-stitution L ∼ Γ . These results can also be applied toexcited states as well, and are verified by numerical cal-culations.The distribution of entanglement entropy is particu-larly relevant for DMRG studies. The typical entangle-ment to appear in a sample grows logarithmically with L and is in agreement with previous calculations for meanentanglement. We find the typical maximum entangle-ment entropy to appear in a sample of length L , anddiscover that it is higher than the entanglement whichdominates the computational time for DMRG using ex-act diagonalization of an effective Hamiltonian.We thank Vedika Khemani and Shivaji Sondhi for stim-ulating discussions. G. Vidal, J. I. Latorre, E. Rico and A. Kitaev, Phys. Rev.Lett. , 227902 (2003). G. Refael and J. E. Moore, Phys. Rev. Lett. , 260602(2004). D. S. Fisher, Phys. Rev. B , 6411 (1995). D. S. Fisher, Phys. Rev. B , 3799 (1994). D. Pekker, G. Refael, E. Altman, E. A. Demler and V.Oganesyan, Phys. Rev. X , 011052 (2014). R. Vosk, D. A. Huse and E. Altman, Phys. Rev. X ,031032 (2015). A. C. Potter, R. Vasseur and S. A. Parameswaran, Phys.Rev. X , 031033 (2015). L. Zhang, B. Zhao, T. Devakul and D. A. Huse, Phys. Rev.B , 224201 (2016). R. Vosk and E. Altman, Phys. Rev. Lett. , 067204(2013). R. Vasseur, A. C. Potter and S. A. Parameswaran, Phys. Rev. Lett. , 217201 (2015). U. Schollwock, Annals of Physics , 96 (2011). D. Pekker and B. K. Clark, arXiv:1410.2224. X. Yu, D. Pekker and B. K. Clark, arXiv:1509.01244. M. Serbyn, A. A. Michailidis, D. A. Abanin and Z. Papi´c,Phys. Rev. Lett. , 160601 (2016). V. Khemani, F. Pollmann and S. L. Sondhi, Phys. Rev.Lett. , 247204 (2016). We only consider a finite system with one cut, while Re-fael and Moore consider a finite subsystem of an infinitesystem and thus two cuts, so our results differ by a factorof 2. M. Fagotti, P. Calabrese and J. E. Moore, Phys. Rev. B , 045110 (2011). I. Peshel, J. Phys. A: Math. Gen.36