Bipartite reweight-annealing algorithm to extract large-scale data of entanglement entropy and its derivative in high precision

Introduction.- With the rapid development of quantum information, its intersection with condensed matter physics has been attracting more and more attention in recent decades [1, 2]. One important topic is to probe the intrinsic physics of many-body systems using the entanglement entropy (EE) [3, 4, 5, 6]. For instance, among its many intriguing features, it offers a direct connection to the conformal field theory (CFT) and a categoical description of the problem under consideration [7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26]. Using EE to identify novel phase and critical phenomena represents a cutting-edge area in the field of quantum many-body numerics. A particularly recent issue is the dispute at the deconfined quantum critical point (DQCP) [27, 28, 29]. The EE at the DQCP, e.g., in the JQ model [30, 31], exhibits totally different behaviours compared with that in normal criticality within the Landau-Ginzburg-Wilson paradigm [20, 32, 33, 34, 35, 36]. According to the prediction from the unitary CFT [37, 38], the EE with a cornered cutting at the criticality should obey the behaviours $s=al-b\mathrm{ln}l+c$, where $s$ is the EE and $l$ is the length of the entangled boundary, in which the coefficient $b$ can not be negative. However, some recent QMC studies show that $b$ is negative. which seemingly suggests the DQCP in the JQ model is not a unitary CFT, probably instead of a weakly first order phase transition [22, 35, 34]. In contrast, another recent work indicates that the sign of $b$ is dependent on the cutting form of the entangled region. For a tilted cutting, $b$ is positive and consistent with the emergent SO(5) symmetry at the DQCP [33]. All in all, the relationship between entanglement entropy and condensed matter physics has been growing increasingly closer recently.

How to obtain high-precision EE via quantum Monte Carlo (QMC) [39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49] with less computation cost and low technical barrier is still a huge challenge in large-scale quantum many-body computation. Although many algorithms have been developed to extract the EE [50, 51, 52, 53, 54, 55, 56, 57, 58, 59], some of which can also achieve high precision, the details of the algorithm becomes more and more complex. The key point of extracting the Rényi ratio $R_{A}^{(2)}=Z_{A}^{(2)}/Z^{2}$ (here we choose the Rényi order to be two without loss of generality) is to calculate the overlap of these two different space-time manifolds, $Z_{A}^{(2)}$ and $Z^{2}$, as shown in Fig.1 (a) and (b), and then get the ratio directly [50, 54]. Once the ratio $R_{A}^{(2)}\rightarrow 0$, the calculation becomes extremely difficult. Due to the area law of EE, the $R_{A}^{(2)}\propto e^{-al}$ decays to zero quickly in large systems, where $l$ is the perimeter of the entangled region. To overcome this difficulty, in the early days, the incremental method of the entangled region was involved [50, 51]. Its key point is one by one adding the lattice site into the entangled region and product the ratio of the process to obtain the final ratio. The process can be written as $R_{A}^{(2)}=Z_{A}^{(2)}/Z^{2}=\prod_{i=0}^{N_{A}-1}Z_{A_{i+1}}^{(2)}/Z_{A_{i}}^{(2)}$, where the $i$ means the number of lattice sites of the entangled region, i.e., $Z_{A_{0}}^{(2)}=Z^{2}$ and $Z_{A_{N_{A}}}^{(2)}=Z_{A}^{(2)}$. Through calculating each intermediate-process ratio $Z_{A_{i+1}}^{(2)}/Z_{A_{i}}^{(2)}$, the high precision $R_{A}^{(2)}$ could be extracted. The shortage of this method is the number of lattice site should be integer which means the splitting of the process must be finite, thus some ratios $Z_{A_{i+1}}^{(2)}/Z_{A_{i}}^{(2)}$ may still be closed to zero even after splitting. We have to note that, the replica manifold are changing during the calculation due to the intermediate processes in this method. It makes the technical barrier of QMC higher.

To solve the finite splitting problem mentioned above, the continually incremental algorithm of QMC has been developed [57, 22, 53]. A virtual process described by a general function $Z_{A}^{2}(\lambda)$ is involved, where $Z_{A}^{2}(\lambda=1)=Z_{A}^{2}$ and $Z_{A}^{2}(\lambda=0)=Z^{2}$. The problem becomes to calculate the ratio $Z_{A}^{2}(\lambda=1)/Z_{A}^{2}(\lambda=0)$ which is equal to $\prod_{i=0}^{N_{A}-1}Z_{A}^{(2)}(\lambda_{i+1})/Z_{A}^{(2)}(\lambda_{i})$. Here the $\lambda$ is continuous from 0 to 1, thus the interval $[0,1]$ can be divided in to any pieces according to the calculation requirement. However, not only the $\lambda$ needs to be changed during the sampling of $Z_{A}^{2}(\lambda)$, but also the entangled region still has to be varied under each fixed $\lambda$, which actually strongly takes huge computational cost and highly improves the technical requirements of the method. Moreover, due to the virtually non-physical intermediate-processes, the results of these intermediate processes $Z_{A}^{2}(\lambda\neq 1,0)$ can not be effectively used which leads a lot of wastes.

In this paper, we propose a simple method without changing the space-time manifold during simulation, and the intermediate-process values are physical and valuable. The high-precision EE and its derivative can be obtained in less cost and low-technical barrier from now on.

Method.- The EE of a subsystem $A$ coupled with an environment $B$ is defined by the reduced density matrix (RDM). Specifically, The $n$th order Rényi entropy is defined as $S^{(n)}=\frac{1}{1-n}\ln[\mathrm{Tr}(\rho_{A}^{n})]=\frac{1}{1-n}\ln R_{A}^{(n)}$, where $R_{A}^{(n)}=Z_{A}^{(n)}/Z^{n}$ and $\rho_{A}=\mathrm{Tr}_{B}\rho$. From the above equations, we know that $Z_{A}^{(n)}\propto\mathrm{Tr}(\rho_{A}^{n})$ while $Z^{n}$ is the proportional factor.

The normalization factor $Z^{n}$ sometimes is not important, for example, when we only concern the dynamical information of the entanglement Hamiltonian (e.g., entanglement spectrum) [26, 60, 61, 62, 63, 64]. In these cases, only the manifold of $Z_{A}^{(n)}$ needs be simulated. However, when we consider the EE calculation, the factor becomes very important for the detailed value. In fact, the hardest difficulty of the calculation of EE comes from simulating the ratio $R_{A}^{(n)}=Z_{A}^{(n)}/Z^{n}$. That’s the reason why the EE algorithms always have to change the simulating manifold between $Z_{A}^{(n)}$ and $Z^{n}$.

Different from the traditional method calculating the ratio $R_{A}^{(n)}$ directly, we start the calculation of $Z_{A}^{(n)}$ and $Z^{n}$ respectively. Here you may feel the sense why we do not need to change the manifold during the simulation. Given a certain distributional function, such as $Z_{A}^{(n)}(J)$, where $J$ is a general parameter (e.g., temperature, the parameter of Hamiltonian, …). For convenience, we define $Z_{A}^{(1)}\equiv Z$ without losing generality. For two parameter points $J^{\prime}$ and $J$, the ratio can be simulated via QMC sampling:

where the $\langle...\rangle_{Z_{A}^{(n)}(J)}$ means the QMC samplings have been done under the manifold $Z_{A}^{(n)}$ at parameter $J$. $W(J^{\prime})$ and $W(J)$ denote the corresponding weights with different parameters $J^{\prime}$ and $J$ for the same configuration sampled by QMC. It means that we simulate the system at parameter $J$, and then obtain a amount of configurations with weight $W(J)$, at the same time, we can calculate the related weight $W(J^{\prime})$ by treating the parameter as $J^{\prime}$ for the same configuration. The ratio of ${W(J^{\prime})}/{W(J)}$ can be calculated in each sampling to gain the final average as Eq.(1).

In principle, the ratio ${Z_{A}^{(n)}(J^{\prime})}/{Z_{A}^{(n)}(J)}$ of any $J^{\prime}$ and $J$ can be solved in the way above. However, we still need to consider how to keep the importance sampling in our QMC simulation. Obviously, if the $J^{\prime}$ and $J$ are close enough, the ratio is close to 1 and is easy for calculation. Otherwise, the ratio would be close to zero or infinite which can not be well-sampled by QMC. It can be easily understood in the reweighting frame as shown in Fig.1 (c), if we want to use a well-known distribution $Z_{A}^{(n)}(J)=\sum W(J)$ to calculate another distribution $Z_{A}^{(n)}(J^{\prime})=\sum W(J^{\prime})$ via resetting the weight of samplings, the weight before/after resetting of the same configuration should be close to each other. Then, it is still an importance sampling in this sense when the $J^{\prime}$ and $J$ are close enough [65, 66, 67, 68]. In fact, this spirit is also the key point of the incremental trick in previous EE algorithms. Therefore, we introduce the continuously incremental trick here to fix the problem:

where $J_{0}=J$ and $J_{N}=J^{\prime}$, other $J_{i}$ is between the two in turn. Thus the QMC can be kept in importance sampling through this reweight-annealing spirit [66, 68].

In this way, we are able to obtain any ratio ${Z_{A}^{(n)}(J^{\prime})}/{Z_{A}^{(n)}(J)}$ in realistic simulation. But it still can not give out the target solution of $Z_{A}^{(n)}(J^{\prime})/Z^{n}(J^{\prime})$. The antidote comes from some well-known point. Considering we have calculated the values of ${Z_{A}^{(n)}(J^{\prime})}/{Z_{A}^{(n)}(J)}$ and ${Z(J^{\prime})}/{Z(J)}$ from the method above, the problem $Z_{A}^{(n)}(J^{\prime})/Z^{n}(J^{\prime})=?$ will be solved if we know a reference point $Z_{A}^{(n)}(J)/Z^{n}(J)$. A simplest reference point is that the $Z_{A}^{(n)}(J)/Z^{n}(J)=1$ when the ground state is a product state $|A\rangle\otimes|B\rangle$. A product state is easy to be gained, for example, we can add an external field in a spin Hamiltonian to polarize all the spins, or a simplest way is decreasing the coupling between $A$ and $B$ to $0$. Of course, other known reference points are also allowed, e.g., the status at infinite temperature or obtaining one known point through other numerical method. In fact, a lot of many-body Hamiltonians have a known ground state at the limit of parameter.

Now the whole scheme for extracting the Rényi EE is complete and realizable, it is not limited to the detailed QMC methods and many-body models. To further understand it and test its performance, we will show a $J_{1}-J_{2}$ spin model [69, 70] via using the stochastic series expansion (SSE) QMC [39, 40, 41, 42, 43, 71] as the example in the followings.

$J_{1}-J_{2}$ model.- We simulate a 2D spin-1/2 antiferromagnetic (AFM) Heisenberg $J_{1}-J_{2}$ model as an example to obtain its EE. The Hamiltonian is

where $\langle ij\rangle$ means the nearest neighbor bond, and $J_{1}$ and $J_{2}$ are the coupling strengths of thin and thick bonds as shown in Fig.2 and it’s ground-state phase diagram (see Fig.2 (c)) has been determined accurately by the quantum Monte Carlo method [69]. The inverse temperature $\beta=2L$ suffices to successfully achieve the desired data quality with high efficiency. In the following simulation, we fix the $J_{2}=1$ and tune the $J_{1}$ from $0_{+}$ to $1$. It’s worth noting that the ground state is a dimer product state when $J_{1}\rightarrow 0$, where the $Z_{A}^{(n)}/Z^{n}=1$ once the dimers are not cut by the entangled edge.

In the SSE frame, the Eq.(1) is equal to

where the $n_{J_{1}}$ is the number of $J_{1}$ operators in the SSE sampling, no matter which space-time manifold $Z_{A}^{(n)}$ is simulated. The details of this equation can be found in the Supplementary Materials (SM).

In the realistic simulation, we need to calculate ${Z_{A}^{(2)}(J_{1}^{\prime})}/{Z_{A}^{(2)}(J_{1}=0_{+})}$ and ${Z(J_{1}^{\prime})}/{Z(J_{1}=0_{+})}$ respectively. And then obtain the final ratio ${[Z_{A}^{(2)}}/Z^{2}]|_{J_{1}^{\prime}}$ based on $[{Z_{A}^{(2)}}/Z^{2}]|_{J_{1}=0_{+}}=1$. Although the parameter is continuously tunable which is same as the non-equilibrium method  [57, 22, 53], the incremental path of our method is physical and meaningful. It is a real parameter path of the cencerned Hamiltonian. In another word, under similar computational cost, the previous method obtains a point while ours gains a curve of EE data. Even ignoring the low-technical-barrier advantage of our method, the efficiency of our calculation is also much higher in this sense.

Cornerless cutting.- Firstly, we try to calculate the EE with cornerless cutting as shown in Fig.2 (a). According to the previous works [20, 22, 23], only the entangled edge without dimers gives the correct results consistent with CFT prediction. In the Fig. 3 ($\mathrm{a}_{1}$), we shows several curves of EE data with different $J_{1}$. The area law fitting data are shown in Table. 1. According to the theoretical prediction [72], $b=N_{G}/2=1$ in the Néel phase of the AFM $J_{1}-J_{2}$ Heisenberg model, where $N_{G}$ means the number of the Goldstone modes. The Table. 1 shows consistent results that the $b\sim 1$ at $J_{1}=1.0,0.9,0.8,0.6$ and $b$ becomes smaller while $J_{1}$ gets closer to the $O(3)$ critical point $J_{c}$. The theoretical calculation [73] pointed out that the $b=0$ at a Wilson-Fisher $O(N)$ quantum criticality of $d\geq 2$ systems. Our result at $J_{c}$ in the Table also supports this prediction. We note that the recent work [74] found the finite size effect in the spin-1/2 AFM Heisenberg model is strong which obviously affects the fitting of $b=1$, and a good fitting needs some more corrections considering the finite size effect. However, we find the simple fitting is not bad in our results. The reason may be the total system we chose is a rectangle and the region $A$ is a square, while they chose a square total system and rectangle $A$. Other QMC works with similar cutting choice as ours also obtains the $b\sim 1$ through the non-equilibrium algorithm [57, 22]. Actually, the finite size effect is more serious in a square total system than a rectangle one, which can also be reflected in Fig.  3 ($a_{3}$) and ($b_{3}$). In the figures we found the derivative of EE obviously converges to the critical point faster in a rectangle case than a square one. We will discuss these data more in the following section.

Another advantage of our method is naturally obtaining the EE of different parameter values, as shown in Fig. 3 ($\mathrm{a}_{2}$). It makes QMC possible to probe the phase transition via EE in 2D and higher dimensional systems as same as the density matrix renormalization group (DMRG) usually does in 1D [75, 76, 77, 78, 79]. In the Fig. 3 ($\mathrm{a}_{2}$), the convexity of the function changes at the critical point, which can also be seen in the derivative of EE more obviously. In the following section, we will introduce a much simpler way to calculate the derivative of EE without incremental process and show the peak of the derivative locates at the phase transition point. It’s worth noting that sometimes the original function of EE directly probe the phase transition while sometimes the derivative does. This issue will be carefully discussed in our coming work ¹¹1Zhe Wang, Zheng Yan, et al. In preparation..

Cornered cutting.- For the cornered cutting case, the $b\sim 0.08$ at the 2+1d O(3) quantum criticality is also known according to previous theoretical and numerical works [81, 82, 83, 20, 22]. In the Fig. 3 ($b_{1}$), the fitting obtains the consistent result $b=0.08(1)$ at $J_{c}$. Similar to the cornerless case, the function EE of $J_{1}$ also displays a change of the convexity at the quantum criticality as shown in Fig. 3 ($b_{2}$).

EE derivative.- It can be proved in the SM that the derivative of the $n$th Rényi EE can be measured in the form:

where the $J$ is a general parameter and $n$ is the Rényi index, the first average is based on the distribution of $Z_{A}^{(n)}$ and the second is $Z$. Taking the $J_{1}-J_{2}$ model as an example with fixed $J_{2}=1$ and $n=2$, set $J_{1}$ as the tunable parameter here and note $H$ is a linear function of $J_{1}$, the Eq.(5) becomes ${dS^{(2)}}/{dJ_{1}}=2\beta\langle{H_{J_{1}}}/{J_{1}}\rangle_{Z_{A}^{(2)}}-2\beta\langle{H_{J_{1}}}/{J_{1}}\rangle_{Z}$, where the ${H_{J_{1}}}$ means the $J_{1}$ term of the $H$. In the SSE frame, it is similar as measuring energy value, which is very simple. The details can be found in the SM.

This conclusion inspires us that we do not need calculating dense data of EE to obtain the derivative. Instead, simulate the average, $2\beta\langle{H_{J_{1}}}/{J_{1}}\rangle_{Z_{A}^{(2)}}-2\beta\langle{H_{J_{1}}}/{J_{1}}\rangle_{Z}$, at the $J_{1}$ we concern directly is enough. We found a similar spirit has been used in calculating the derivative of the Rényi negativity versus the inverse temperature [84]. Through this method, we calculate the derivative of EE curves in the Fig. 3 ($\mathrm{a}_{2}$) and Fig. 3 ($b_{2}$). As shown in Fig. 3 ($\mathrm{a}_{3}$) and Fig. 3 ($b_{3}$), the peaks of the EE derivative really locate at the phase transition point.

In fact, this measurement of the EE derivative also points out another way for calculating the EE value through an integral

where the ${dS^{(n)}}/{dJ}$ can be gained through Eq. (5) and the EE $S^{(n)}(J_{0})$ at the reference point should be known, such as a product state. We demonstrate the equivalence of the two ways (Eq. (4) and Eq. 6) by taking the $J_{1}-J_{2}$ model as an example, as shown in the Fig. 4, either in cornerless or cornered case.

We note Jarzynski’s equality [85] can also be used in our methods, similar to the previous non-equilibrium algorithms [53, 57]. But we found that there is almost no acceleration effect for the non-equilibrium version compared with the equilibrium QMC, the details in which will be discussed in our coming work ²²2Zhe Wang, Zhiyan Wang, Yi-Ming Ding, Zheng Yan, et al. In preparation..

Conclusion.- Overall, we develop a practical and unbiased scheme with low-technical barrier to extract the high-precision EE and its derivative from QMC simulation. The space-time manifold needn’t be changed during the simulation and the measurement is an easy diagonal observable. All the intermediate-process measurement-quantities are meaningful and physical, which is unwasteful and useful. It makes QMC possible to probe the novel phases and phase transition by scanning EE in a large parameter region in 2D and higher dimensional systems as same as the DMRG usually does in 1D.

We display a simulating-path from a dimer product state to a Néel phase by taking the $J_{1}-J_{2}$ AFM Heisenberg model as an example. A diverging peak of EE derivative arises at the phase transition point. The universal coefficient of the sub-leading term of the EE has been successfully extracted in high precision, either at O(3) criticality or continuous-symmetry-breaking phase. Our method is not only limited to boson QMC, but can also be used to other QMC approaches, such as the fermion QMC ³³3Weilun Jiang, Zheng Yan, et al. In preparation., for highly entangled quantum matter.

Acknowledgement We thank the helpful discussions with Jiarui Zhao, Dong-Xu Liu, Yin Tang, Zehui Deng, Bin-Bin Chen, Yuan Da Liao, and Wei Zhu. ZY acknowledges the collaborations with Yan-Cheng Wang, Z. Y. Meng and Meng Cheng in other related works. The work is supported by the start-up funding of Westlake University. The authors thank the high-performance computing center of Westlake University and the Beijng PARATERA Tech Co.,Ltd. for providing HPC resources.