Introduction

As a measure of the quantum correlations in many-body systems, the entanglement entropy (EE) has become an indispensable tool in modern physics. The EE quantifies the amount of non-local correlation between a subsystem and its compliment. The ability of the EE to expose the non-local features of a system offers a way to characterize topological orders1,2, solve the black hole information paradox3, and quantify information scrambling in a thermalization process under unitary time evolution4,5,6.

Recently, the EE was measured in quantum many-body systems for the first time7,8. Specifically, the second Rényi EE (2REE), one of the variants of the EE, of a pure quantum state was measured in quantum quench experiments using ultra-cold atoms. For such pure quantum states, it is believed that the EE of any small subsystem increases in proportion to the size of the subsystem just like the thermodynamic entropy9. This is called the volume law of the EE. However, when the size of a subsystem becomes comparable to that of its complement, it is observed experimentally that the EE starts deviating from the volume law, and eventually decreases (Fig. 1). The curved structure of the size dependence of the EE, which first increases linearly and then decreases, universally appears in various excited pure states, for example, energy eigenstates10 and states after quantum quenches11,12. We call this curved structure a Page curve, after D. Page13, and it is of fundamental importance in explaining these examples to reveal a universal behavior of the Page curve.

Fig. 1
figure 1

A schematic picture of our setup. The second Rényi Page curve for pure states, \(S_2(\ell )\), follows the volume law when \(\ell\) is small, but gradually deviates from it as \(\ell\) grows. At the middle, \(\ell {\mathrm{ = }}L{\mathrm{/}}2\), the maximal value is obtained, where the deviation from the volume law is ln 2 (see the Results section). Past the middle \(\ell {\mathrm{ = }}L{\mathrm{/}}2\), it decreases toward \(\ell {\mathrm{ = }}L\) and becomes symmetric under \(\ell \leftrightarrow L - \ell\)

Full size image

Despite their ubiquitous appearance, the theoretical understanding of the Page curves is limited to the case of a random pure state13,14, which is a state at infinite temperature and can be defined only in a finite-dimensional Hilbert space. By contrast, the cold-atom experiments address finite temperature systems and an infinite-dimensional Hilbert space; therefore, it is important to develop a theory of the Page curve applicable to these situations. Additionally, there are practical needs for the estimation of the slope of the volume law. The slope is often employed, for example, to calculate the corresponding thermodynamic entropy8,10 and to detect a transition between the energy eigenstate thermalization hypothesis (ETH) phase and the many-body localized (MBL) phase15,16. However, since the experimentally or numerically accessible sizes of the systems are small, the estimation of the volume-law slope is deteriorated by the curved structure of the Page curve.

In this work, we show that Page curves in broad classes of excited pure states exhibit universal behaviors. We first derive the function of the Page curve for canonical thermal pure quantum (cTPQ) states, which are pure states representing thermal equilibrium at a temperature β−117,18. In particular, the Page curve of the 2REE is controlled only by two parameters, an effective dimension ln a and an offset ln K, for any Hamiltonian and at any temperature. We then conjecture and numerically verify that this feature of the Page curve universally appears in any sufficiently scrambled pure states representing equilibrium states; that is, our function fits the 2REE of the energy eigenstates of a non-integrable system and the states after quantum quenches including the state realized in the above-mentioned experiment8. By contrast, in the case of the energy eigenstates of an integrable system, which are not scrambled at all, we find that their Page curves deviate from our function. Since our function enables us to estimate the slope of the volume law from small systems with high accuracy and precision, our result is also numerically effective in detecting the ETH-MBL transition15,16,19,20 and improves the estimation of the critical exponent.

Results

Derivation of the Page curve in cTPQ states

Let us consider a lattice Σ containing L × M sites (Fig. 1), equipped with a translation-invariant and local Hamiltonian H. We divide Σ into two parts, A and B, each containing \(\ell \times M\) and \((L - \ell ) \times M\) sites. The n-th Rényi EE of a pure quantum state |ψ〉 is defined as

$$S_n(\ell ){\mathrm{ = }}\frac{1}{{1 - n}}{\mathrm{ln}}\left( {{\mathrm{tr}}\rho _{\mathrm{A}}^n} \right),$$
(1)

where ρA ≡ trB|ψ〉〈ψ|. We call \(S_n(\ell )\) as a function of \(\ell\) the n-th Rényi Page curve (nRPC). We note that we use the term differently from how it is used in the context of quantum gravity, where it denotes the temporal dependence of the entanglement during the formation of a black hole. To simplify the calculation, we assume \(\ell ,L \gg 1\).

To derive the behaviors of nRPCs for any Hamiltonian, we utilize cTPQ states, which are proposed along with studies of typicality in quantum statistical mechanics21,22,23,24,25,26,27.The cTPQ state at the inverse temperature β is defined as

$$\left| \psi \right\rangle {\mathrm{ = }}\frac{1}{{\sqrt {Z_\psi } }}\mathop {\sum}\limits_j z_je^{ - \beta H/2}\left| j \right\rangle ,$$
(2)

where \(Z_\psi \equiv \mathop {\sum}\nolimits_{i,j} z_i^ \ast z_j{\langle}i|e^{ - \beta H}|j\rangle\) is a normalization constant, {|j〉} j is an arbitrary complete orthonormal basis of the Hilbert space \({\cal H}_\Sigma\), and the coefficients {z j } are random complex numbers \(z_j \equiv (x_j + iy_j){\mathrm{/}}\sqrt 2\), with x j and y j obeying the standard normal distribution \({\cal N}(0,1)\). For any local observable, the cTPQ states at β reproduce their averages in thermal equilibrium at the same inverse temperature18. As a starting point, we here derive the following exact formula of the 2RPC at a temperature β−1 for any Hamiltonian (see the Methods section for the calculations and results for the nRPCs),

$$S_2(\ell ){\mathrm{ = }} - {\mathrm{ln}}\left[ {\frac{{{\mathrm{tr}}_{\mathrm{A}}({\mathrm{tr}}_{\mathrm{B}}e^{ - \beta H})^2 + {\mathrm{tr}}_{\mathrm{B}}({\mathrm{tr}}_{\mathrm{A}}e^{ - \beta H})^2}}{{({\mathrm{tr}}e^{ - \beta H})^2}}} \right].$$
(3)

We also give several simplifications of Eq. (3). The first step is to decompose the Hamiltonian H as H = HA + HB + Hint, where HA,B are the Hamiltonians of the corresponding subregion and Hint describes the interactions between them. Since the range of interaction Hint is much smaller than \(\ell\) and \(L - \ell\) due to the locality of H, we obtain the simplified expression

$$S_2(\ell ){\mathrm{ = }} - {\mathrm{ln}}\left( {\frac{{Z_{\mathrm{A}}(2\beta )}}{{Z_{\mathrm{A}}(\beta )^2}} + \frac{{Z_{\mathrm{B}}(2\beta )}}{{Z_{\mathrm{B}}(\beta )^2}}} \right){\mathrm{ + }}{\mathrm{ln}}R(\beta ),$$
(4)

where R(β), coming from Hint, is an O(1) constant dependent only on β and \(Z_{{\mathrm{A,B}}}(\beta ) \equiv {\mathrm{tr}}_{{\mathrm{A,B}}}(e^{ - \beta H_{{\mathrm{A,B}}}})\).

Further simplification occurs through the extensiveness of the free energy, −ln ZA,B/β, which is approximately valid when \(\ell ,L - \ell \gg 1\). In the region in question, ln ZA,B is proportional to the volume of the corresponding subregion, and thus, we replace ZA(2β)/ZA(β)2 and ZB(2β)/ZB(β)2 with \(Q(\beta )a(\beta )^{ - \ell }\) and \(Q(\beta )a(\beta )^{ - (L - \ell )}\), respectively. Here, a(β) and Q(β) are O(1) constants dependent only on β, and 1 < a(β) holds because of the concavity and monotonicity of the free energy. Finally, we reach a simple and universal expression for the 2RPC:

$$S_2(\ell ){\mathrm{ = }}\ell {\mathrm{ln}}a(\beta ) - {\mathrm{ln}}\left( {1 + a(\beta )^{ - L + 2\ell }} \right){\mathrm{ + }}{\mathrm{ln}}K(\beta ),$$
(5)

where K ≡ R/Q. This is our first main result. The same simplifications can be applied to a general nRPC (see the Methods section). For example, concerning the 3RPC, we have

$$S_3(\ell ){\mathrm{ = }}\ell \frac{{{{\mathrm {ln}}}b}}{2} - \frac{1}{2}{\mathrm{ln}}\left( {1{\mathrm{ + }}K_1^\prime \frac{{b^\ell }}{{a^L}}{\mathrm{ + }}b^{ - L + 2\ell }} \right){\mathrm{ + }}{\mathrm{ln}}K_2^\prime ,$$
(6)

where b, \(K_1^\prime\) and \(K_2^\prime\) are O(1) constants that depend only on β.

The significance of Eq. (5) is apparent: it tells us that the 2RPC is determined by only two parameters: a(β) and K(β). The first term represents the volume law of entanglement for \(\ell \le L/2\), and thus, its slope, ln a(β), is a density of the 2REE in the thermodynamic limit L → ∞. The third term, ln K, represents an offset of the volume law. The second term gives the deviation from the volume law, which stems from the highly non-local quantum correlation between subsystems A and B. Here, we see that the way that quantum correlations appear in cTPQ states is completely characterized by the volume-law slope, a(β). As \(\ell\) approaches L/2, the quantum correction to the volume law becomes stronger and eventually becomes exactly ln 2 at \(\ell {\mathrm{ = }}L{\mathrm{/}}2\), independent of the inverse temperature β and the Hamiltonian. This is a unique feature of the 2RPC, as we do not observe such universal behaviors in the nRPCs for n ≥ 3. With regard to the third term, a similar offset term appears at zero temperature. It comes from the degeneracy of a quantum state at zero temperature, and is referred as the topological EE for topological states1. Similarly, ln K contains the degeneracy term, but it also contains other terms, e.g., ln Q. It is an interesting future direction to decode the topological EE from ln K.

In addition, by using Eq. (5), the mutual information is straightforwardly obtained. Suppose that the state is the cTPQ state (Eq. (2)) and the system is divided into three parts, A, B, and C. The (second Rényi) mutual information between A and B is defined as \(I_2 \equiv S_2^{\mathrm{A}} + S_2^{\mathrm{B}} - S_2^{{\mathrm{A}} \cup {\mathrm{B}}}\). It becomes and calculated as

$$I_2{\mathrm{ = }}{\mathrm{ln}}\left( {\frac{{a^{ - \ell } + a^{ - L + \ell }}}{{(a^{ - \ell /2} + a^{ - L + \ell /2})^2}}} \right) + q{\mathrm{ln}}K,$$
(7)

where \(S_2^X\) is the 2REE in X, q = 1 when \({\mathrm{A}} \cup {\mathrm{B}}\) is connected and q = 0 when A∪B is disconnected, \(\ell\) is the sum of the length of A and B, and, for simplicity, we take the both lengths to be \(\ell /2\). Eq. (7) explains the observed size dependence of the mutual information in ref. 8. I2 grows exponentially with \(\ell\) for \(\ell < L/2\), and shows a linear growth for \(L/2 < \ell\). See Supplementary Note 1 for the detailed explanations.

Finally, to confirm the validity of the approximations and clarify the advantages of our formula (5), we present numerical simulations of the 2RPC of cTPQ states for the S = 1/2 XY chain under a periodic boundary condition,

$$H{\mathrm{ = }}\mathop {\sum}\limits_{i = 1}^L \left( {S_i^xS_{i + 1}^x + S_i^yS_{i + 1}^y} \right).$$
(8)

This system is mapped to the free fermion system by the Jordan-Wigner transformation28, and the quantities \({\mathrm{tr}}_{\mathrm{A}}({\mathrm{tr}}_{\mathrm{B}}e^{ - \beta H})^2\) and \({\mathrm{tr}}_{\mathrm{B}}({\mathrm{tr}}_{\mathrm{A}}e^{ - \beta H})^2\) can be efficiently calculated in a large system (L ~ 100) by the correlation functions of the system29. We numerically calculate the 2RPC of the cTPQ states at the inverse temperature β = 4 by evaluating Eq. (3). As Fig. 2 shows, the numerical data of the 2RPC are well-fitted by our formula (5) for all system sizes L and subsystem sizes \(\ell\) (details of the fitting is described in the Methods section). In addition, we compare several estimates of the density of the 2REE from numerical data in the inset: ln a from the fits by our formula, the density of the 2REE for half of the system, S2(L/2)/(L/2), and the average slope of the curve between \(\ell {\mathrm{ = }}1\) and \(\ell = 5,\,(S_2(5) - S_2(1)){\mathrm{/}}4\). It is clear that ln a does not contain any systematic error compared with the other two estimates, which represents one of the advantages of our formula (5). We also numerically check the validity of the approximations in deriving the formula (5) in Supplementary Figure 1.

Fig. 2
figure 2

Second Rényi Page curve in cTPQ states. The dots represent the second Rényi Page curves in the cTPQ states of the spin system (9) at an inverse temperature β = 4 calculated by Eq. (3) for various system sizes L. The lines are the fits by Eq. (5) for the numerical data. The inset shows the fitted values of ln a, S2(L/2)/(L/2), and the average slope of the curve between \(\ell {\mathrm{ = }}1\) and \(\ell {\mathrm{ = }}5\). The dotted lines are the extrapolations to L → ∞ by 1/L scaling for ln a and S2(L/2)/(L/2) and by 1/L2 scaling for the average slope

Full size image

General conjecture for scrambled states

So far, we have focused on the cTPQ states, but they are merely a canonical example of pure states locally reproducing the Gibbs ensemble. Here, we pose a conjecture for other scrambled pure states: Eq. (5) works as a fitting function for generic scrambled pure states.

In the following two subsections, we numerically check this conjecture. We calculate the 2RPCs of various pure states, namely, the excited energy eigenstates of (non-)integrable models, and the states after a quantum quench. We show that the conjecture holds for the eigenstates of the non-integrable model but not for those of the integrable model. We also numerically reveal that Eq. (5) well fits the 2RPC averaged over the time evolution after a quantum quench.

Numerical results for energy eigenstates

As the ETH claims, in a wide class of models, the energy eigenstates look thermal—the expectation values of the local observables reproduce those of the Gibbs ensemble30,31,32. From the viewpoint of ETH, its extension to non-local quantities is interesting10,33. We test whether the formula (5) applies to the 2RPC, which is highly non-local at \(\ell {\mathrm{ = }}O(L)\), in particular \(\ell \simeq L/2\).

As an example, we take the S = 1/2 XXZ spin chain with/without next-nearest neighbor interactions under the periodic boundary condition,

$$H{\mathrm{ = }}\mathop {\sum}\limits_{i = 1}^L \left( {S_i^xS_{i + 1}^x + S_i^yS_{i + 1}^y + \Delta S_i^zS_{i + 1}^z + J_2{\mathbf{S}}_i \cdot {\mathbf{S}}_{i + 2}} \right),$$
(9)

where we set Δ = 2 and J2 = 4 for non-integrable cases and J2 = 0 for integrable cases34.

Figure 3a, b show the 2RPC of the eigenstates of this model with various energies, which are obtained by exact diagonalization. We see that the fit by the formula (5) works quite well for the non-integrable cases, although not for the integrable cases. Moreover, as Fig. 3c, d clearly indicate, the residuals of the fits per site for all eigenstates decrease with respect to L for the non-integrable cases but increase for the integrable cases. We therefore numerically conclude that our formula (5) is applicable to non-integrable models but not to integrable models. We provide a discussion of the physics behind this result in Supplementary Note 2.

Fig. 3
figure 3

Second Rényi Page curve for general energy eigenstates. a 2RPCs of several energy eigenstates of the non-integrable Hamiltonian, Eq. (9) with Δ = 2 and J2 = 4 (dots), and the fits by our formula (5) (lines). The inset shows the energy spectrum of the Hamiltonian, and the arrows indicate the eigenstates presented in the figure. b Same as figure a for the integrable Hamiltonian (Δ = 2, J2 = 0). c Residuals of fits per site \(r_i \equiv L^{ - 1}\mathop {\sum}\nolimits_{\ell = 0}^L (S_2(\ell )_{i,{\mathrm{data}}} - S_2(\ell )_{i,{\mathrm{fit}}})^2\), where \(S_2(\ell )_{i,{\mathrm{data}}}\) is the 2REE of the i-th eigenstate and \(S_2(\ell )_{i,{\mathrm{fit}}}\) is a fitted value of it, for all eigenstates of the non-integrable Hamiltonian (10) with Δ = 2 and J2 = 4 (we consider only the sector of a vanishing total momentum and magnetization). The eigenstates are sorted in descending order in terms of the residuals, and the horizontal axis represents their percentiles. The fits become better as the size of the system increases. d Same as figure c for the integrable Hamiltonian (Δ = 2, J2 = 0). The fits become worse as the size of the system increases

Full size image

We comment on our results of the energy eigenstates from the viewpoint of ETH. First, the success of our formula in the non-integrable case is important for the following reasons: the corresponding thermal “ensemble” is not the usual microcanonical ensemble (mixed state) but is rather a thermal pure state. Although the Page curve necessarily deviates from the volume-law slope in a pure state, the way how it deviates always exhibits a universal behavior. Furthermore, \(S_2(\ell )\) is a highly non-local and complicated observable. We thus expect that our results will bring the studies of ETH to the next step, i.e., its non-local extension. Second, with regard to the extension of ETH to non-local quantities, there are new proposals in which the effect of the energy fluctuation is incorporated10,35. This is called subsystem ETH35. In the subsystem ETH, the authors suggest that the volume-law slopes of the higher-order REE for the energy eigenstates may be different from those of the Gibbs ensembles or the cTPQ state. We provide a discussion on this deviation in Supplementary Note 2. In our numerical calculations on the energy eigenstates, however, we do not see the deviation of the 2RPC from (5), which is derived for the cTPQ states. This might be because of the limitation of the system sizes, and it would be interesting to see the deviation indeed occurs in larger systems. Indeed, this extension was recently analyzed in ref. 36. The EE of energy eigenstates in a non-integrable model was studied there by substituting ρdiag in Supplementary Eq. (11) by the microcanonical density matrix (ensemble), and the result supports our generalized formula (Supplementary Eq. (11)). Third, the failure of our formula in the integrable cases is surprising because ETH for local observables was proved to hold for almost all eigenstates, even in integrable systems37. By contrast, our results in Fig. 3c, d clearly show that almost all eigenstates indeed violate formula (5), or non-local ETH. A similar observation was made in ref. 38, where the EE of eigenstates in a different integrable model from ours was explicitly calculated. They found that the subsystem size dependence of the EE is completely different from the result by Page, where the EE is close to its maximum value13. Their result is consistent with our finding that the 2REE of an integrable model is qualitatively different from that of a non-integrable model.

Numerical results for quenched states

Next, as a second example of thermal quantum states where our formula (5) applies, we consider pure states after quantum quenches in closed quantum systems, where some parameters of their Hamiltonians are abruptly changed11,39,40. When such states become stationary after a long time, they are considered to represent a thermal equilibrium corresponding to the Gibbs ensemble in non-integrable systems, while they do not in integrable systems because infinitely many conserved quantities block thermal behaviors41,42.

Here, we numerically simulate the dynamics of the 2RPC after a quantum quench from pure states and find that our formula (5) explains well the 2RPC of the stationary states. We note that experimental measurement of the dynamics of the 2RPC was already realized in ref. 8. Indeed, our results explains the experimental data well: Fig. 4a in ref. 8 is well-fitted by our formula (5), and the parameters are estimated to be ln(a) = 0.974 and ln K = 0.162 (Supplementary Figure 2).

Fig. 4
figure 4

Dynamics of the second Rényi Page curve after quantum quench. a Time evolution of the 2RPC in a non-integrable system (Eq. (9) with Δ = 1 and J2 = 0.5) after a quantum quench from the Néel state. The dotted line is the fitting by Eq. (5) for the time average of \(S_2(\ell )\). The inset shows the dynamics of the 2REE at the center of the system, S2(L/2). b Same as figure a for the integrable Hamiltonian (Δ = 1,J2 = 0). Eq. (5) fits the time average well in both a and b

Full size image

We consider the Hamiltonian (10) with Δ = 1, i.e., the Heisenberg model. Again, we consider both integrable (J2 = 0) and non-integrable (J2 = 0.5) cases. An initial state of a quantum quench (t = 0) is taken as the Néel state |Néel〉≡|↑↓↑↓↑↓…〉, and the dynamics after the quench, |ψ(t)〉 = e−iHt|Néel〉, is numerically calculated by exact diagonalization. We take a time step of the evolution as dt = 0.1 and calculate the dynamics of 2RPC \(S_2(\ell ,t)\) up to t ≤ T = 300.

In Fig. 4, we show the numerical results of the 2RPC \(S_2(\ell ,t)\) after the quench (see also Supplementary Movies 1 and 2). At first, the 2REE increases linearly with time until it saturates11,40,43,44,45. After it saturates, temporal oscillations of the 2RPC are observed for a long time as long as t = T = 300 (see the inset). We consider these oscillations as being due to the finite-size effect (L = 16), and thus, we also present the time average of the Page curve, \(\bar S_2(\ell ):{\mathrm{ = }}\frac{1}{T}{\int}_0^T {\mathrm{d}}tS_2(\ell ,t)\), as an estimation of the long-time limit, \(\mathop {{{\mathrm{lim}}}}\nolimits_{t \to \infty } S_2(\ell ,t)\). As clearly seen in Fig. 4, the time average of \(S_2(\ell )\) is well-fitted by our formula (5) (dotted line) for both non-integrable and integrable cases. This illustrates the validity of our formula for pure states after a quantum quench.

Although our formula (5) is derived from a thermal state, it somehow works in the integrable cases where the states never thermalize. This success is achieved because the pure quantum states are scrambled and partially thermalized, which usually leads the states to the generalized Gibbs ensemble (GGE)41. Hence, the slope ln a quantifies the number of states in the GGE. We give a (non-rigorous) proof of the validity of the formula for the time-averaged 2RPC in Supplementary Note 2. We also note several differences between the integrable cases and the non-integrable cases. First, for the integrable cases, other physical observables, such as staggered magnetization, are not explained by the cTPQ states or thermal equilibrium states46. Second, the properties of the temporal fluctuations are different between the two cases. A fluctuation is larger in integrable cases than in non-integrable cases, as one can see in the inset of Fig. 4, where the dynamics of the 2RPC at the center of the system, S2(L/2), is plotted. We observe that the fluctuation decays algebraically with the system size L for the integrable cases, whereas it decays exponentially in the non-integrable cases47. These differences probably reflect the existence of infinitely many conserved quantities in integrable systems.

We also comment on the implication of our results for recent studies on chaos in quantum many-body systems. In refs. 5,6, the time dependence of S2(t) is related to the out-of-time-ordered correlation, which captures the essence of quantum chaos48. As we discuss in Supplementary Note 2, the convergence of \(S_2(\ell ,t)\) in time and the applicability of formula (5) thereof mean that the quantum states become scrambled in whole-spatial regions. Therefore, we expect that the time until the convergence of the 2RPC after the quench is a good diagnostic for the scrambling time49, and it would be interesting to relate those two time scales directly in future work.

Application to many-body localization transition

Thus far, we have shown that broad classes of thermal quantum states obey our analytical formula (5). Finally, we demonstrate its practical applications and advantages. The example we take here is a problem of the phase transition between an ETH phase and a MBL16,20 phase, where the volume law of entanglement provides an important diagnostic to distinguish between the two phases15,50.

MBL is defined as a disorder-induced localization in interacting systems and has recently been an active field of research because of the nontrivial interplay between the interactions and the disorder. Meanwhile, the ETH phase is the phase where ETH holds, and it appears when the strength of the disorder in a system is sufficiently weak (or absent). There is considered to be a continuous phase transition between the ETH phase and the MBL phase when the strength of disorder varies at a fixed energy density (temperature) of the system.

Since MBL essentially involves interactions, it is difficult to study them analytically, and one often resorts to numerical approaches. One of the simple diagnostics for ETH and MBL employed in such numerical studies is the n-th Rényi EE per site s n , which will exhibit a transition from being non-zero (ETH phase) to being zero (MBL phase) in the thermodynamic limit. Most previous studies have utilized sn,center = S n (L/2)/(L/2)15,50,51,52, the EE per site at the center of the system, as an estimate of s n . As we see in the inset of Fig. 2, however, the sn,center includes a systematic error of O(1/L) from the correct value of s n in the thermodynamic limit because of the deviation of the nRPC from the volume law at the center of the system. This error harms analyses of the ETH–MBL transition since the system size L accessible by numerical methods is not very large, L ~ 22.

Here, our formula (5) comes in, and extracts the 2REE per site in the thermodynamic limit in relatively small systems, as we illustrated in Fig. 2 (we give a brief justification of the validity of our formula in disordered systems in Supplementary Note 3). To substantiate the advantage of our formula (5), we study the ETH–MBL transition in the prototypical S = 1/2 spin chain for MBL20,50,53,

$$H{\mathrm{ = }}\mathop {\sum}\limits_{i = 1}^L \left( {{\mathbf{S}}_i \cdot {\mathbf{S}}_{i + 1} + h_i^zS_i^z} \right),$$
(10)

where a random magnetic field \(\{ h_i^z\} _i\) is drawn from a uniform distribution [−h, h] and a periodic boundary condition is imposed.

This model exhibits the ETH–MBL phase transition at the critical disorder strength hc. Using the von Neumann EE (\(S_{n\rightarrow 1}\)), the authors of ref. 50 estimated hc = 3.62 ± 0.2 and the critical exponent of the transition ν = 0.80 ± 0.4 at the center of the energy spectrum. However, the estimated ν violates the Harris bound ν ≥ 254. This violation can probably be understood as a finite-size effect, and we will show that the usage of formula (5) improves the situation. We note that the Harris bound for S2 is the same as that of \(S_{n\rightarrow 1}\)54.

In Fig. 5a, we perform a finite-size scaling of the 2REE per site, s2, for the first time. Those data are extracted from the fitting of the 2RPC \(S_2(\ell )\) of the eigenstates of the system, namely, s2 = ln(a) at the center of the energy spectrum. We obtain hc = 3.60 ± 0.12 and ν = 1.88 ± 0.2 (details are presented in the Methods section). Although we study smaller systems, up to L = 16, compared with the previous study50 (\(L \sim 22\)), the estimation of ν exhibits a significant improvement to satisfy the Harris bound ν ≥ 2. This highlights the usefulness of our formula (5) since a finite-size scaling of the conventional method, s2,center = S2(L/2)/(L/2), under the same conditions yields hc = 3.60 ± 0.12 and ν = 1.30 ± 0.12 (Fig. 5b).

Fig. 5
figure 5

Finite-size scaling across the ETH–MBL phase transition. a Finite-size scaling of ln(a), extracted from the fitting of the 2RPC of the eigenstates of the Hamiltonian (11), vs. L1/ν(h − hc). The estimation of the critical exponent ν is significantly improved. b Same finite-size scaling as a for s2,center = S2(L/2)/(L/2), which is a conventional estimate of the 2REE per site50

Full size image

In general, when a physical quantity in a finite-size system has an O(1/L) difference from the quantity in the thermodynamic limit, this difference deteriorates the finite-size scaling, and the critical exponent estimated from the scaling becomes completely different from that in the thermodynamic limit (for example, see the conventional ferromagnetic-spin glass transition55,56). In this sense, the critical exponent ν in our study is also sensitive to the finite-size effect of O(1/L). However, using formula (5), we can remove the finite-size effect of O(1/L), as shown in Fig. 2, resulting in an improvement in the estimation of ν.

Discussion

In conclusion, we studied the volume law of entanglement in general scrambled pure quantum states. By employing the cTPQ states, we derive analytical formulae for the nRPC for general local Hamiltonians at any temperature. In particular, we show that the 2RPC \(S_2(\ell )\) is parametrized by only two parameters (Eq. (5)), and our formula improves the finite-size scaling of the thermodynamic quantities. We numerically demonstrate that the same formula for the 2RPC works as a fitting function for general scrambled pure states other than cTPQ states, namely, excited eigenstates of general Hamiltonians and states after quantum quenches. We also propose the generalized formula that can incorporate the energy fluctuation in Supplementary Eq. (11). Finally, we employ our formula to detect the ETH–MBL phase transition. The full characterization of the 2RPC by our formula improves the estimate of the critical exponent of the transition and would resolve the controversy over the breakdown of the Harris bound at the transition.

Methods

Derivation of Equation 5

Here we present the detailed calculation of n-th Rényi EE (REE) of the cTPQ states (Eq. (5))17,18.

For any local observable \({\cal O}\) on \({\cal H}\)∑, its expectation value of a cTPQ state satisfies \(\overline {\langle\psi |{\cal O}|\psi \rangle} {\mathrm{ = tr}}({\cal O}e^{ - \beta H}){\mathrm{/tr}}(e^{ - \beta H})\), and the standard deviation from the average is exponentially small with respect to the volume of the system18 (here we denote a random average over the coefficients {z j } by \(\overline \cdots\)). In this sense, we can regard cTPQ states as a canonical example of thermal pure states that corresponds to the thermal Gibbs ensemble ρGibbs ≡ e−βH/tr(e−βH) at inverse temperature β.

The reduced density matrix of the cTPQ state ψ in A is written as

$$\rho _{\mathrm{A}}{\mathrm{ = tr}}_{\mathrm{B}}\left| \psi \right\rangle \left\langle \psi \right| = \frac{1}{{Z_\psi }}\mathop {\sum}\limits_{a_1,a_2,b_1,i_1,j_1} \pi _{12}\left| {a_1} \right\rangle \left\langle {a_2} \right|,$$
(11)

where \(\pi _{pq} \equiv z_{i_p}z_{j_p}^ \ast \left\langle {a_pb_p\left| {e^{ - \frac{1}{2}\beta H}} \right|i_p} \right\rangle \left\langle {j_p\left| {e^{ - \frac{1}{2}\beta H}} \right|a_qb_p} \right\rangle\), \(\left| {ab} \right\rangle \equiv \left| a \right\rangle \otimes \left| b \right\rangle\), and \(\left\{ {\left| a \right\rangle } \right\}_a\) and \(\left\{ {\left| b \right\rangle } \right\}_b\) are complete orthonormal bases in the subsystem A and B, respectively. The two indices \(i_p,j_q \in \left\{ {\left| j \right\rangle } \right\}_j\) run over the orthonormal basis \(\left\{ {\left| j \right\rangle } \right\}_j\). In this notation, we obtain

$${\mathrm{tr}}_{\mathrm{A}}\rho _{\mathrm{A}}^n{\mathrm{ = }}\frac{{\mathop {\sum}\limits_{(a),(b),(i),(j)} \pi _{12},\,\pi _{23}, \cdots, \pi _{n1}}}{{Z_\psi ^n}},$$
(12)

where (x) ≡ x1, x2, ⋯, x n . What we want to calculate is a random average of n-th REE, \(\overline {S_n} :{\mathrm{ = }}\frac{1}{{1 - n}}\overline {{\mathrm{ln}}\left( {{\mathrm{tr}}_{\mathrm{A}}\rho _{\mathrm{A}}^n} \right)}\). The term π12, π23⋯πn1, however, includes the product of random variables such as \(z_{i_1}z_{j_1}^ \ast,\, z_{i_2}z_{j_2}^ \ast, \cdots, z_{i_n}z_{j_n}^ \ast\) and it is difficult to calculate the average of the logarithm of it. Instead, we calculate S n by averaging the trace before taking the logarithm of it, \(\tilde S_n:{\mathrm{ = }}\frac{1}{{1 - n}}{\mathrm{ln}}\left( {\overline {{\mathrm{tr}}_{\mathrm{A}}\rho _{\mathrm{A}}^n} } \right)\). As we give a proof in Supplementary Note 4, the difference between \(\overline S _n\) and \(\tilde S_n\) is exponentially small in terms of the system size L (see also ref. 36).

By taking the random number average of \({\mathrm{tr}}_{\mathrm{A}}\rho _A^n\) with using several properties of {z i } such as \(\overline {z_i} {\mathrm{ = }}0,\,\overline {z_i^ \ast z_j} {\mathrm{ = }}\delta _{ij}\) and \(\overline {|z_i|^2|z_j|^2} {\mathrm{ = }}1 + \delta _{ij}\), we can calculate n-th REE of the cTPQ state for any Hamiltonian (below we just use S n to denote \(\tilde S_n\)). For example, the 2REE is Eq. (3), and the third REE is

$$\begin{array}{*{20}{l}} {S_3} \hfill & {\mathrm{ = }} \hfill & { - \frac{1}{2}{\mathrm ln}\left[ {\left( {{\mathrm{tr}}_{\mathrm{A}}\left( {{\mathrm{tr}}_{\mathrm{B}}e^{ - \beta H}} \right)^3 + 3{\mathrm{tr}}\,{\it{M}}} \right.} \right.} \hfill \\ {} \hfill & {} \hfill & {\left. {\left. { + {\mathrm{tr}}_{\mathrm{B}}\left( {{\mathrm{tr}}_{\mathrm{A}}e^{ - \beta H}} \right)^3 + N} \right){\mathrm{/}}\left( {{\mathrm{tr}}\,{\it{e}}^{ - \beta H}} \right)^3} \right],} \hfill \end{array}$$
(13)

where \(M{\mathrm{ = }}e^{ - \beta H}\left( {{\mathrm{tr}}_{\mathrm{B}}e^{ - \beta H} \otimes {\mathrm{tr}}_{\mathrm{A}}e^{ - \beta H}} \right)\) and \(N{\mathrm{ = }}{\sum} \left\langle {a_1b_1|e^{ - \beta H}|a_3b_2} \right\rangle\)\(\left\langle {a_2b_2\left| {e^{ - \beta H}} \right|a_1b_3} \right\rangle \left\langle {a_3b_3\left| {e^{ - \beta H}} \right|a_2b_1} \right\rangle\). After the same simplification as we do for the second REE, the term N is shown to be exponentially smaller than the other terms. Then, we obtain the final result Eq. (6). It is difficult to write down a general expression of S n for all n, but we can calculate it for given integer n ≥ 2 systematically by taking an average in Eq. (12). It would be interesting to study whether the universality that we reported for the 2REE holds for the EE (n = 1).

Finally, we also comment on the case of infinite temperature β = 0 and its relation to the previous studies. When β = 0, we can obtain a fairly simple equation for general n:

$$S_n{\mathrm{ = }}\ell {\mathrm{ln}}2 - \frac{1}{{n - 1}}{\mathrm{ln}}\left[ {\mathop {\sum}\limits_{k = 1}^n N(n,k)\left( {\frac{{2^\ell }}{{2^{L - \ell }}}} \right)^{k - 1}} \right],$$
(14)

where \(N(n,k) = \frac{1}{n}\left( {\begin{array}{*{20}{c}} n \\ k \end{array}} \right)\left( {\begin{array}{*{20}{c}} n \\ {k - 1} \end{array}} \right)\) is known as the Narayana numbers. It is possible to take an analytic continuation n → 1 and reproduce the result on the EE S1 = −trA(ρA ln ρA) calculated in ref. 13.

Numerical fitting by our formula (5)

Throughout this work, numerical fitting by our formula (5) for given data of the 2RPC \(\{ S_2(\ell )\} _{\ell = 0}^L\) is performed with the least squares method implemented in the numerical package scipy.optimize.leastsq by regarding a and K in Eq. (5) as fitting parameters.

Details on the numerical calculations on ETH–MBL transition

In the demonstration of the ETH–MBL transition, we did exact diagonalization on the Hamiltonian (11) \(H_{{\mathrm{MBL}}}{\mathrm{ = }}\mathop {\sum}\nolimits_{i = 1}^L \left( {{\mathbf{S}}_i \cdot {\mathbf{S}}_{i + 1} + h_i^zS_i^z} \right)\), where random magnetic fields \(\{ h_i^z\} _i\) are taken from a continuous uniform distribution within [−h,h] and periodic boundary condition is imposed. Since the Hamiltonian conserves the total magnetization \(S_{{\mathrm{tot}}}^z{\mathrm{ = }}\mathop {\sum}\nolimits_i S_i^z\), we fixed the magnetization sector as \(S_{{\mathrm{tot}}}^z{\mathrm{ = }}0\). We define energy density ε of the eigenstates with energy E by using the extremal values of energies, ϵ = (E−Emin)/(Emax−Emin), and focus on the case of ϵ = 0.5 throughout our analysis. For each realization of the magnetic fields, an energy eigenstate closest to ε = 0.5 is chosen and the 2RPC \(S_2(\ell )\) for that eigenstate is calculated. The number of realizations of the random magnetic fields is 1000 for all presented data. After averaging \(S_2(\ell )\) over all realizations, we perform fitting of \(S_2(\ell )\) by our formula (5) and obtain ln(a) for each value of the disorder strength h.

The scaling analysis of ln(a) and s2,center = S2(L/2)/(L/2) is performed to determine the critical disorder strength hc and the critical exponent ν in a following procedure50. First, we assume a universal scaling function of the form g(L1/ν(h − hc)) in a window of width 2w centered at hc, namely, [hc − w, hc + w]. The function g is approximated by a polynomial of degree three and the polynomial as well as hc and ν are optimized to make all data collapse into one line. In our analysis w = 2.0 is used, and (hc, ν) = (3.60 ± 0.12, 1.88 ± 0.2) is found for the case of ln(a), whereas (hc, ν) = (3.60 ± 0.12, 1.30 ± 0.13) for s2,center (the errors are estimated by a bootstrap method).

Data availability

All numerical data and computer codes used in this study are available from the corresponding author upon request.