Performance evaluation of the discrete truncated Wigner approximation for quench dynamics of quantum spin systems with long-range interactions

The discrete truncated Wigner approximation (DTWA) is a powerful tool for analyzing dynamics of quantum-spin systems. Since the DTWA includes the leading order quantum corrections to a mean-field approximation, it is naturally expected that the DTWA becomes more accurate when the range of interactions of the system increases. However, quantitative corroboration of this expectation is still lacking mainly because it is generally difficult in a large system to evaluate a timescale on which the DTWA is quantitatively valid. In order to investigate how the validity timescale depends on the interaction range, we analyze dynamics of quantum spin models subjected to a sudden quench of a magnetic field by means of both DTWA and its extension including the second-order correction, which is derived from the Bogoliubov-Born-Green-Kirkwood-Yvon equation. We also develop a new formulation for calculating the second-order R\'enyi entropy within the framework of the DTWA. By comparing the time evolution of the R\'enyi entropy computed by the DTWA with that by the extension including the correction, we find that both in the one- and two-dimensional systems the validity timescale increases algebraically with the interaction range.


I. INTRODUCTION
The state-of-the-art technologies established in quantum optics have opened the door for controlling and exploring coherent quantum dynamics of isolated manybody systems both near and far from equilibrium [1]. In the context of condensed-matter and solid-state physics, charge-neutral atoms loaded onto an optical lattice have been extensively studied as an analog quantum simulator for the tight-binding Hubbard-type models with shortrange interactions. Owing to its controllability and cleanness, one can gain access to fundamental questions about dynamical properties of Hubbard-type systems. The recent topics explored in experiments include thermalization dynamics of an isolated quantum system [2,3], propagation of non-local correlations [4,5], the Kibble-Zurek mechanism across quantum phase transitions [6], and the many-body localization (MBL) in a disordered optical lattice [7,8]. In recent years, technological developments in creating, controlling, and probing cold atoms or molecules with strong dipole-dipole interactions in an optical lattice [9][10][11][12][13][14][15], Rydberg gases [16][17][18][19][20][21][22][23][24][25], and trapped ions [26][27][28][29][30][31] have enabled quantum simulation of various quantum-spin systems with long-range interactions. In particular, Rydberg gases can be manipulated by means of the optical tweezer techniques, so that these offer an intriguing playground to explore novel quantum magnetism and non-equilibrium dynamics of localized spins caused by variable-range interactions.
While these experimental advances have stimulated * Electronic address: kunimi@ims.ac.jp theoretical studies of quantum many-body dynamics of systems with various interaction ranges, they are still limited due to the lack of available computational tools. As a quasi-exact numerical method, the time-dependent density-matrix renormalization group (tDMRG) has been typically utilized for simulating largescale many-body systems corresponding to actual experiments [32,33]. However, its efficient applications are limited to one-dimensional (1D) systems. Among various candidates for approximate frameworks to tackle manybody systems, the phase space methods, especially the truncated-Wigner approximation (TWA) [34,35] on the basis of the Wigner-Weyl correspondence, provide a realistic and widely-applicable approach to quantum manybody dynamics even for higher-dimensional systems with long-range interactions [36][37][38][39][40][41]. Employing the TWA, quantum dynamics are reduced to a semiclassical problem of simulating randomly distributed classical trajectories in a phase space, each of which obeys a saddlepoint or mean-field equation of motion for a given quantum system. The TWA gives quantitative descriptions of quantum dynamics even at long times if the system is in a certain classical limit or nearly non-interacting limit. More precisely, the TWA is asymptotically exact at short times [42]. It implies that, within the TWA, there exists a threshold timescale separating semiclassically simulatable and non-simulatable regimes of quantum dynamics depending on the choice of the phase space. As demonstrated in some works [43][44][45], by increasing the number of phase-space variables, one can improve the validity of TWA descriptions for strongly-correlated lattice systems composed of bosons or spins. Such an increased phase-space approach is referred to as the SU(N ) or cluster TWA. Furthermore, it is also possible to con-struct a fermionic TWA (fTWA) approach for interacting fermions, in which so(2N ) string variables are introduced for a fermionic mode number N [37,41]. The fTWA has been used to study semiclassical aspects of chaos in the Sachdev-Ye-Kitaev (SYK) model [46,47], which consists of all-to-all (infinitely long-range) two-body hoppings. For describing dynamics of quantum spin systems, the discrete TWA (DTWA) has been widely applied in various contexts [13][14][15][48][49][50][51][52][53][54][55][56][57][58][59][60][61]. In the conventional use of the TWA for spin systems [35,43,45], to be efficient, its Monte-Carlo sampling part employs a Gaussian approximation for the continuous Wigner distribution function. On the other hand, the DTWA utilizes a discrete Wigner function for sampling phase-space variables instead of the continuous Wigner distribution. Since the discrete Wigner representation is defined for the basis of local-spin eigenstates rather than coherent states, it can express typical initial states such as the all down spin state | ↓↓↓↓ · · · and the staggered magnetization state | ↑↓↑↓ · · · as a positive-valued distribution. Thanks to this advantage, the DTWA accurately describes all the initial moments of these states and can capture some revival properties of quantum dynamics beyond the Gaussian approximation. More interestingly, the DTWA can also reproduce the experimental results for Rydberg atoms [51,52] and dipolar atoms [13][14][15], which are effectively described by spin-1/2 and large-S models, respectively.
Although the DTWA is a powerful tool to analyze quantum spin systems, there are some problems. One is a timescale on which the DTWA is valid. Generally, the TWA framework gives quantitatively valid results in a short time regime [34,35]. This validity timescale depends on the details of the systems, such as interactions, dimensions, and initial conditions. In the zerotemperature ground states or thermal equilibrium cases, it is well established that the mean-field approximation gives the exact results in large-S limit or infinite dimensions (or all-to-all coupling). As for the quantum dynamics, it has been shown that the validity timescale becomes longer when the size of the spin S increases [62]. By contrast, the dependence of the validity timescale on the spatial range of interactions has not been investigated. Investigating such a range-dependence of the validity timescale will be useful when the TWA is applied to analyzing quantum spin dynamics of systems consisting of Rydberg-dressed atoms [17,19], in which the interaction range can be controlled.
In this paper, we address the question how the validity timescale of the DTWA depends on the interaction range. Our approach is to use higher-order corrections of the DTWA, which are derived by using the Bogoliubov-Born-Green-Kirkwood-Yvon (BBGKY) hierarchy equation [54]. When the difference between the DTWA and its higher-order corrections is small, we can expect that the DTWA is a good approximation.
In order to compare the DTWA with its higher order corrections, we focus on the second order Rényi entan-glement entropy. This quantity plays an important role in many situations, such as thermalization, MBL, etc. Recently, the Rényi entropy has been experimentally observed [63][64][65]. In this paper, we develop a new method to calculate the Rényi entropy within the framework of the DTWA. Our new method can be applicable to not only the benchmark of the DTWA but also the calculations of the Rényi entropy in higher dimensions, which are difficult to access by other methods.
From the comparison between the DTWA and its extension including higher-order corrections, we show that the validity timescale of the DTWA becomes longer as increasing the interaction range. We confirm this property for three different kinds of quantum-spin models, namely Ising, XY, and Heisenberg model under a uniform magnetic field in 1D and 2D. This result means that the DTWA becomes better as the classical limit is approached.
This paper is organized as follows: In Sec. II, we explain our model and the DTWA. In Sec. III A, we explain how to define the threshold time on which the DTWA is valid. In Sec. III B, we show the results of Rényi entanglement entropy and threshold time for three different quantum-spin models. In Sec. III C, we compare the results in 1D and 2D systems. In Sec. IV, we summarize our results. In Appendix A, we discuss the details of the derivation of the DTWA. In Appendix B, we explain the sampling scheme of the initial conditions. In Appendix C, we derive the expression of the Rényi entanglement entropy in the framework of the DTWA. In Appendix D, we compare the DTWA and tDMRG results in one dimension. In Appendix E, we propose a new algorithm to calculate the dynamics of the systems with the longrange interaction by using the tDMRG.

A. Model
In this paper, we consider a family of quantum spin-1/2 systems in 1D and 2D, which is generally modeled by the Hamiltonian whereŜ µ i (µ = x, y, z) is a spin-1/2 operator at site i, J µ ij is a spin-exchange coupling between two distant sites, and h ≡ (h x , h y , h z ) is a uniform magnetic field, respectively. Throughout this paper, we impose open boundary conditions and write M as an even integer expressing the total number of lattice sites. As a concrete form of the coupling J µ ij , we especially focus on three specific cases as follows: From the top to bottom, let us refer to these as a ferromagnetic Ising, ferromagnetic XY, and antiferromagnetic Heisenberg models, respectively. Furthermore, we assume that a magnetic field is applied along x-axis, i.e., h = (h x , 0, 0). The details of J (D) ij depend on the spatial dimension of the lattice D. For the 1D cases (D = 1), it has the properties where J > 0 is an interaction strength and r = 1, 2, · · · is the interaction range. By contrast, the 2D cases (D = 2), in which the lattice geometry is supposed to be square, are characterized by where is the position of i-th lattice site, a is the lattice constant, i x = 1, 2, · · · , M x and i y = 1, 2, · · · , M y are indices of the i-th lattice sites, and R r is the distance between r-th neighboring sites. The total number of lattice points is given by M = M x M y . The constant C r is determined such that the following equation is satisfied: The explicit values of C r are given by C 1 = 2, C 2 = 4, C 3 = 6, C 4 = 10, C 5 = 12, · · · . For later use, we define N r as where N r approximately denotes the number of connections per spin quantifying how many spins are connected to each spin. Approaching the boundaries from the center of the system, the actual number of connections decreases from N r due to the finite system size and the open boundary. It is worth noting that these types of long-range interaction, Eqs. (5) and (6), are realizable in the experimental setups by means of Rydberg-dressed atoms [17,19]. A system of such atoms is typically characterized by a softcore type potential, so that interactions among atoms are almost constant in the short-distance regime and rapidly decay in the long-distance regime [66,67]. The couplings in Eqs. (5) and (6) may describe such a situation approximately.
In the subsequent sections we will investigate suddenquench dynamics of the interacting spin systems in order to characterize the limitation of the DTWA method. To be concrete, we especially consider the following directproduct wave functions as low-entangled initial states: Here |↑ i and |↓ i denote the eigenstates ofŜ z i while The corresponding discrete Wigner functions for these initial states are shown in Appendix B.
In the descriptions of DTWA, the phase-point operator at t is assumed to be factorized [48], i.e., The time-dependent coefficient r i (t, α) ≡ 2S i (t) obeys a classical equation of motion obtained from a first-order BBGKY hierarchy truncation (see also Appendix A) where µβγ is the Levi-Civita symbol and we used the Einstein notation for repeated Greek indices. Inserting Eq. (16) into Eq. (12), we arrive at a DTWA expression of the Rényi entropy where S i (0) = r αi /2 and S i (0) = r α i /2. The doubled bracket means a phase-space average weighted with two initial Wigner functions For direct product states such as Eqs. (9), (10), and (11), W α (0) is factorized as W α (0) = M j=1 w αj (0). It means that each local spin variable can fluctuate independently and the entropy for a subsystem results in zero at t = 0. The subsystem entropy remains zero during the time evolution if the Hamiltonian is entirely decoupled into local parts and each part is linear in SU(2) matrices. For non-linear systems with a nonzero spinexchange coupling, the Rényi entropy becomes nonzero as a consequence of the many-body time evolution. The semiclassical expression for the subsystem entropy states that the amount of entanglement across the boundary of two subregions is related to the degree of complexity in a solution of Eq. (17) that is, if it is possible to write down, provided as a complicated non-linear function of initial conditions.
A higher order correction beyond the DTWA description based on Eq. (17) arises in a classical trajectory of an enlarged phase space for a second order BBGKY method [54]. The underlying idea of this method is to regard a non-separable part ofŜ µ iŜ ν j , which is represented by c µν ij in a replacementŜ µ iŜ ν j → S µ i S ν j + c µν ij , as an additional mechanical variable and define an approximately closed equation of motion for S µ i and c µν ij . We note that the first order BBGKY truncation leads to the time evolving equation given by Eq. (17). In Sec. III, we will also exploit the second order BBGKY method to compute the Rényi entropy, especially for a subregion of two sites. The detail of the BBGKY formulation will be presented in Appendix A.
In this paper, we numerically solve Eq. (17) for the 1st order BBGKY and Eqs. (A12) and (A13) for the 2nd order BBGKY by using a 4th order Runge-Kutta method. A time step ∆t is taken to be ∆t = 10 −3 /J. We use M = 100 in 1D and M x = M y = 14 in 2D.

A. Criterion for the validity of the DTWA
In this subsection, we introduce a criterion for giving an estimation of the timescale within which the DTWA is quantitatively valid. Our approach is based on the assumption that when the difference between the 1st order and 2nd order BBGKY results is small, the DTWA gives a good approximation. Here, a question arises; Which physical quantities are appropriate for comparing the 1st order and 2nd order BBGKY results? We propose that the 2nd order Rényi entanglement entropy is a suitable quantity for confirming the validity of the DTWA. One advantage to use the Rényi entropy is that it is an unbiased quantity compared with other physical quantities such as the spin expectation values and spin-spin correlations. The latter quantities strongly depend on the dynamics and symmetry of the systems. For example, if the system has spin rotational symmetry along z-axis, S z tot ≡ iŜ z i is conserved so that it is not appropriate for examining the validity of the DTWA.
In this paper, we calculate mean two-site Rényi entropies for 1D and 2D, which are defined by where S (2) ij (t) is the two-site Rényi entropy. Let us mention differences between our formulation of the Rényi entropy and that of previous works [13,55,62]. In these previous works, they calculated the single-site or two-site Rényi entropy from the expressions of the single-or twosite reduced density matrix operator because the matrix elements of these reduced density matrices can be constructed by the expectation values ofŜ µ i andŜ µ iŜ ν j , which can be obtained by the DTWA. The advantage of our formulation is that it allows us to calculate the Rényi entropy for multiple sites. It is easy to calculate the multiple-sites Rényi entropy in the 1st order BBGKY. For more details, see Appendix C. Here, we focus on the case where the two sites are nearest neighbors in 1D. In 2D, we consider the neighboring sites of the x-direction only.
We calculate the two-site Rényi entropy by using the 1st and 2nd order BBGKY equations. In order to quantify the difference between the 1st and 2nd order results, we define where S 1st (t) and S (2) 2nd (t) are the mean two-site Rényi entropy obtained by the 1st and 2nd order BBGKY equation, respectively. This quantity represents a relative error of the 1st and 2nd order results. The reason why we do not use the relative error of S (2) 1st (t) and S (2) 2nd (t) is to avoid the divergence of the relative error because S (2) (0) = 0 in our initial conditions. Comparing these quantities, we can define a threshold time T th , at which ∆(t) exceeds a small positive number . In this paper, we use = 1/10.
We note that the 2nd order BBGKY equation is numerically unstable as pointed out in Refs. [56,57]. In fact, we find the divergent behavior of the 2nd order BBGKY equation. For example, this behavior can be seen in r = 1 results in Fig. 1. This is an intrinsic property of the 2nd order BBGKY equation. We have checked that this divergence behavior is not artificial one because it does not depend on the choice of the time step ∆t.
The small difference between the 1st and 2nd order results is a necessary condition for the DTWA to be good approximation. This criterion is based on the assumption that the 1st and 2nd order results can approximate the exact results. Even if the difference is small, our criterion is meaningless if the DTWA cannot reproduce the exact results. To corroborate that our criterion indeed works, in Appendix D we perform the comparison between the DTWA results and the tDMRG method in the 1D cases.

B. Rényi entropy and threshold time
In this subsection, we show the results of Rényi entropy and threshold time for Ising, XY, and Heisenberg model. We consider three different spin models in order to indicate that the statement that the validity timescale of the DTWA increases with increasing the interaction range holds regardless of the integrability and symmetry. In 1D, the Ising model with transverse magnetic field is integrable while the XY and Heisenberg models under a uniform magnetic field are nonintegrable. The Ising and XY has only discrete symmetry while the Heisenberg model has continuous spin rotation symmetry around the magnetic field.

Ising model
Here, we show the results of the mean two-site Rényi entropy for the Ising model under the transverse magnetic field in Fig. 1 (a) for 1D and Fig. 1 (b) for 2D. In these cases, the initial condition is the fully −x-polarized state [see Eq. (9)]. This state is the exact ground state when h x → ∞. We can see the growth of the entanglement in an early stage of the dynamics. This behavior is a typical one in the quench dynamics of manybody systems [69]. In the long-range interacting systems, the growth of the entanglement is slow compared to the short-range interacting systems. This tendency is consistent with the previous work [70]. In r → ∞ limit, our model becomes the Lipkin-Meshkov-Glick model [71,72]. In this model, the dynamics is constructed by a small number of quantum states. The bipartite entanglement entropy is bounded by the logarithm of the system size. Although this property is related to the bipartite entanglement, we can naively expect that the two-site Rényi entropy has same tendency.
From the results of the Rényi entropy, we can obtain the threshold time T th as a function of N r . The results are shown in Fig. 2 (a) for 1D and Fig. 2 (b) for 2D. We can see that the threshold time increases as a power law of N r (see solid and dotted lines in Fig. 2. We can also see that T th in the large h x region is large compared to that in small h x region. This behavior can be understood from the fact that the DTWA yields exact results when interaction terms are absent. In the large h x region, the dynamics are mainly driven by the magnetic field. In this reason, T th in the large h x region is longer than that of the small h x region.

XY model
Next, we show the results of the XY model under the magnetic field along −x direction in Fig. 1 (c) for 1D and 1 (d) for 2D. In this case, the initial condition is the fully −z-polarized state [see Eq. (10)]. This state is an exact eigenstate when h x = 0.
From these results, we obtain the threshold time as a function of N r . The results are shown in Fig. 3 (a) for 1D and Fig. 3 (b) for 2D. The results are similar to those of the Ising model. We can see that the threshold time increases as power law of N r (see solid and dotted lines in Fig. 3). The dependence of the magnetic field is also similar to that of the Ising model. The large h x case is better than the small h x case.

Heisenberg model
We show the results of the Heisenberg model under the magnetic field along −x direction in Fig. 1 (e) for 1D and Fig. 1 (f) for 2D. In this case, the initial condition is the Neel state [see Eq. (11)]. In contrast to the previous cases, we consider the antiferromagnetic interaction and nonuniform initial condition. The reason is as follows. If we consider the ferromagnetic Heisenberg model with a fully-polarized initial condition under the uniform magnetic field, the resultant dynamics is the Larmor precession motion, which is not affected by the interaction. During this dynamics, the entanglement is exactly zero. Therefore, we need to consider another situation.
The threshold time as a function of N r is in Fig. 4. We can see that the threshold time increases as power law of N r (see solid and dotted lines in Fig. 4). The dependence of the magnetic field is also similar to that of the Ising and XY models. In contrast to the previous two models, we can see an oscillation behavior of T th . For example, the oscillation of h x = 10J shown in Fig. 4 (a) can be seen clearly. This is due to the relation between the interaction and the initial condition. To explain this behavior, we consider r = 1 and r = 2 cases in an early time regime. In the r = 1 case, each pair of spins coupled via the Heisenberg interaction are aligned antiparallelly in the initial state such that the initial state has relatively low energy. However, in the r = 2 case, the next-nearestneighbor interactions couple parallelly aligned pairs of spins in the initial state such that the initial state has much higher energy than the r = 1 case. In other words, the injected energy by the quench alternates when r increases one by one. Therefore, the oscillating behavior of T th as a function of the interaction range appears.

C. Comparison with one-and two-dimensional results
Here, we compare the results of the threshold time in 1D and 2D. The results are shown in Fig. 5 (a) for the Ising case, Fig. 5 (b) for the XY case, and Fig. 5 (c) for the Heisenberg case. From these results, we can see that the threshold time increases as power law of N r in all the models and spatial dimensions. We can also find that T th results are almost overlapped in 1D and 2D at the same h x except the Ising model for h x = 4.0J, which will be discussed later. These results suggest that the threshold time depends on the number of interacting spins. We conclude that when the number of interacting spin is increased, the validity timescale of the DTWA becomes longer.
Here, we remark on the result of the Ising model for h x = 4.0J (see blue triangle symbols in Fig. 5.). Unlike the other results, the 1D and 2D results are clearly deviated. This deviation may be attributed to the distance of the parameter from the quantum critical point. The critical field strength of the transverse field Ising model for r = 1 is h x c = 0.5J for 1D [73] and h x c 3.044J for 2D [74,75] in our notation. This means that the distance of the magnetic field from the quantum critical point is different in 1D and 2D. In 2D, h x = 4.0J is closer to the critical value than 1D case. Therefore, the 1D and 2D results deviate clearly.

IV. SUMMARY
In summary, we investigated the timescale on which the DTWA can quantitatively describe quantum dynamics of spin-1/2 models. In order to corroborate this, we developed a new formulation of the Rényi entropy within the DTWA framework. Using this new formulation, we evaluate the Rényi entropy after a sudden quench in the Ising, XY, and Heisenberg models under a uniform magnetic field in 1D and 2D. The Rényi entropy is calculated by the DTWA and its extension including the secondorder correction that is derived from the BBGKY hierarchy equation. Comparing these results, we determined the threshold time, on which the relative error of the exponential of the Rényi entropy in the DTWA and 2nd order BBGKY results exceeds 10%. We found that the threshold time increases as a power law function of the interaction range (or the number of interacting spins per site). This result suggests that the accuracy of the DTWA becomes better as the classical limit (in this case, all-to-all coupling) is approached. This behavior is consistent with the properties of the equilibrium cases.
In this paper, we focused on the sudden quench dynamics of the quantum spin systems. The sweep dynamics, in which a parameter of the system varies slowly, is an also important problem. The adiabatic sweep of the parameter across the phase transition point leads the Kibble-Zurek mechanism. The DTWA can be applicable to this phenomenon. The adiabatic sweep is necessary for longtime evolution. Therefore, it is important to confirm the validity of the DTWA in the sweep dynamics case. In this appendix, we discuss the derivation of the DTWA for more details. See also Ref. [54].
A starting point for deriving the BBGKY hierarchy equation is the von Neumann equation for the density matrix operator: whereρ(t) ≡ e −iĤt/ ρ(0)e +iĤt/ . By using the phasepoint operator, the density matrix operatorρ(t) can be written as [68]ρ Substituting Eq. (A2) into Eq. (A1), we obtain the equa- tion forÂ α (t): Here, we introduce the partial trace of the phase-point operatorÂ where Tr i and Tr ij are the trace over the Hilbert space except the site i and sites i and j (i = j), respectively. To derive the BBGKY hierarchy equation, we use the following cluster expansion A i (t, α) andB ij (t, α) can be expanded aŝ where r i (t, α) and c µν ij (t, α) are expansion coefficients and determined by solving the classical equation of mo-tion, which will be discussed below, and we use the Einstein's notation for Greek indices in Eq. (A9). We note that the relation c µν ij (t) = c νµ ji (t) holds. SettingB ij (t, α) = 0 in Eq. (A6), we can obtain the 1st order BBGKY hierarchy equation, which corresponds to the conventional DTWA. The Wigner-Weyl symbol of the spin operatorŜ µ i becomes The equation of motion for S µ i (t) can be obtained by the time derivative of Eq. (A4): where we used Eqs. (A3) and (A6) and setB ij (t, α) = 0. We note that c µν ij (t) = 0 in the 1st order BBGKY hierarchy equation becauseB ij (t, α) = 0.

(B4)
The above discrete Winger functions are semi-positive definite and normalized; αi w αi (0) = 1. Therefore, we can regard the above discrete Wigner functions as a probability distribution function and use the Monte Carlo sampling for the initial states.
As pointed out in Ref. [54], there is ambiguity of the sampling scheme of the initial states because we have degrees of freedom of the definition of the phase-point operator (or definition of r α ). For example, if we use two different phase-point operators, we can decompose the density matrix operator asρ(0) = α [W α (0)Â α /2 + W α (0)Â α /2], where W α (0) andÂ α are the discrete Wigner function and phase-point operator for different definitions. This means that we have infinite number of possible choices of the sampling. According to Ref. [54], the suitable choice of the phase-point operator depends on the model and initial condition and they proposed some better choices rather than using Eqs. (B2), (B3), and (B4). To implement the modified sampling, we introduce the following quantities For the Ising model and XY model, we sample r αi from the following set with equal probability [54]: In this appendix, we derive the expression of the Rényi entropy in the DTWA. Here, we consider the subsystem A and its complement B. The total system is given by sum of the A and B. The reduced density matrix for subsystem A is defined bŷ where Tr B denotes the trace over the subsystem B. The second order Rényi entropy is defined by Using the discrete Wigner function [see Eq. (A2)], we can write Eq. (C1) aŝ In the 1st order BBGKY hierarchy equation, we approximate the phase-point operator asÂ α (t) M i=1Â i (t, α). Substituting this expression into Eq. (15), we obtainρ where we used Tr iÂαi (t) = 1. To obtain the Rényi entropy, we calculate where the initial conditions of S i (t) and S i (t) are sampled from W α (0) and W α (0), respectively. This expression implies that we can obtain the second order Rényi entropy by using the replica method. The procedure is as follows: First, we prepare two independent copies of the initial states and calculate the time evolution for two copies independently. Then, we calculate the ensemble average of i∈A [1/2 + 2S i (t) · S i (t)]. The second order Rényi entropy in the 1st order BBGKY is given by Next, we derive the expression of the Rényi entropy in the 2nd order BBGKY hierarchy. On the contrary to the 1st order BBGKY hierarchy, we restrict the subsystem size to two sites. This is due to a practical reason. In the 2nd order BBGKY, we need to approximate the phasepoint operatorÂ α (t) by using the cluster expansion. If we consider a large subsystem, we need expressions for a higher order cluster expansion ofÂ α (t), which is difficult to write down. Therefore, we only consider two site Rényi entropy S (2) ij (t). The reduced density matrix operator in the 2nd order BBGKY becomeŝ where we used Eq. (A6). Using this expression, we obtain where the initial conditions for c µν ij (t) and c µν ij (t) are sampled from W α (0) and W α (0), respectively, and we also use the Einstein's notation for Greek indices. From Eq. (C8), we can obtain the two-site Rényi entropy in the 2nd order BBGKY hierarchy.
If the Hamiltonian consists of two-site operators, bond terms, the time evolution of MPS can be performed by operating Trotter gates to MPS [32,33,77,78]. Even though the Hamiltonian has long-range interactions, one can perform the time evolution of MPS with utilizing the swap gates [79]. It should be noted that the swap gates to be operated are not unique, and that even the number of required swap operations can be different. Less swap operations require less computational resources. One may come up with a good choice of swap gates if the types of bond terms are limited likewise Bauernfeind et al. [80]. If the Hamiltonian consists of many types of bond terms likewise the long-range models such as Eq. (1), finding out a good choice is quite a exhausting task. In this appendix, we present an algorithm which automatically produces an efficient (maybe not best) choice of swap operations.
The Hamiltonian consisting of two-site operators can be expressed asĤ and one can compute the Trotter gates exp(−iδtĤ i,j ) from bond termsĤ i,j . At first step, we group the pair indices of bond terms [i, j] so that bond terms in the same group commute each other. We also tries to group bond terms with the same distance j − i and order groups in ascending order of the distance. The grouping can be accomplished by Algorithm 1 [81]. Next, we determine the arrangement of site indices where the Trotter gates in each group are operated. As shown in Algorithm 2, swap operations required from one arrangement to another arrangement can be obtained by a sort algorithm implemented only by adjacent swap operations such as the bubble sort or the gnome sort. Therefore, it is suffice to determine only the arrangement of site indices. In order to find an arrangement requiring less swap operations, we reorder bond terms in a group in ascending order in the sense of the present arrangement not to disturb the present arrangement so much. This ordering determines the arrangement of site indices in bond terms and we have to insert remaining indices between bonds. Similarly to the case of the bond terms, ordering based on the present arrangement will not disturb the present arrangement so much. Here, we have a simple alternative: insert remaining indices based on the first index of bond pairs or the second one. From the alternative, we can obtain two candidates of an arrangement. The procedure for obtaining two candidates is summarized in Algorithm 3, and one example is given in Fig 7. From the two candidates, we select one candidate with considering possibilities in the next group. With using Algorithm 3, one can obtain two candidates in the next group for each candidate. Based on consequent four candidates, we choose one arrangement for the present group which is contained in the best candidates. By iterating this process over groups, one can obtain a sequence of arrangement and swap operators required for performing time evolutions of long-range Hamiltonians.
Algorithm 4 produces the efficient gates in Bauernfeind et al. [80] when bonds = [(1, 2), (1, 3), . . . , (1, M )], and thus we consider that gates from the algorithm are efficient. From Trotter gates calculated from bond gates and the ordered list of swap operations given by Algorithm 4, one can obtain the list of gates corresponding to the first order decomposition of the time-evolution operator i<j exp(−iδtĤ i,j ). The second order decomposition is obtained by successive operations of gates in the reversed list. Furthermore, higher order decompositions can be obtained from compositions of the second order decompositions [76].