\section{Introduction}\label{sec:Introduction}

The state-of-the-art technologies established in ultracold gases have opened the door for controlling and exploring coherent quantum dynamics of isolated many-body systems both near and far from equilibrium \cite{Bloch2012}. 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 short-range 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~\cite{Trotzky2012,Gring2012}, propagation of non-local correlations~\cite{Cheneau2012,Takasu2020a}, the Kibble--Zurek mechanism across quantum phase transitions \cite{Braun2015}, and the many-body localization (MBL) in a disordered optical lattice \cite{Choi2016,Abanin2019}. In recent years, technological developments in creating, controlling, and probing cold atoms or molecules with strong dipole-dipole interactions in an optical lattice \cite{dePaz2013,Yan2013,Hazzard2014,dePaz2016,Lepoutre2019,Fersterer2019,Patscheider2020}, Rydberg gases \cite{Schauss2015,Zeiher2016,Labuhn2016,Haenel2017,Zeiher2017,Bernien2017,Sous2018,Leseleuc2018,Guardado-Sanchez2018,Sous2019,Leseleuc2019,Keesling2019,Mizoguchi2019a}, and trapped ions \cite{Britton2012,Islam2013,Jurcevic2014,Richerme2014,Bohnet2016,Garttner2017} 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 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 large-scale many-body systems corresponding to actual experiments \cite{vidal_efficient_2003,white_real-time_2004}. However, its efficient applications are limited to one-dimensional (1D) systems. Among various candidates for approximate frameworks to tackle many-body systems, the phase-space methods, especially the truncated-Wigner approximation (TWA) \cite{Blakie2008,Polkovnikov2010} on the basis of the Wigner--Weyl correspondence, provide a realistic and widely applicable approach to quantum many-body dynamics even for higher-dimensional systems with long-range interactions \cite{Tuchman2006,Davidson_thesis,Nagao2019,Ruostekoski2005,Fujimoto2019,Davidson2017}. 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 saddle-point 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 \cite{Polkovnikov2003}. 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 \cite{Davidson2015,Nagao2020a,Wurtz2018}, 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 construct a fermionic TWA (fTWA) approach for interacting fermions, in which $so(2N)$ string variables are introduced for a fermionic mode number $N$ \cite{Davidson_thesis,Davidson2017}. The fTWA has been used to study semiclassical aspects of chaos in the Sachdev--Ye--Kitaev (SYK) model \cite{Schmitt2019,Scaffidi2019}, 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 \cite{Schachenmayer2015,Schachenmayer2015_2,Babadi2015,Pucci2016,Orioli2017,Acevedo2017,Covey2018,Mori2019,Qu2019,Morong2020a,Orioli2018,Czischek2018,Sundar2019,Signoles2019a,Lepoutre2019,Fersterer2019,Patscheider2020}. In the conventional use of the TWA for spin systems \cite{Davidson2015,Polkovnikov2010,Wurtz2018}, 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 $|\downarrow\downarrow\downarrow\downarrow\cdots\rangle$ and the staggered magnetization state $|\uparrow\downarrow\uparrow\downarrow\cdots\rangle$ 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. These properties have been investigated in Ref.~\cite{Sundar2019}. More interestingly, the DTWA can also reproduce the experimental results for Rydberg atoms \cite{Orioli2018,Signoles2019a} and dipolar atoms \cite{Lepoutre2019,Fersterer2019,Patscheider2020}, which are effectively described by spin-$\frac{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 \cite{Blakie2008,Polkovnikov2010}. This validity timescale depends on the details of the systems, such as interactions, dimensions, and initial conditions. In the zero-temperature 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~\cite{Zhu2019}. By contrast, the dependence of the validity timescale on the spatial dependence of the interaction is not well understood. In Refs.~\cite{Schachenmayer2015,Schachenmayer2015_2,Sundar2019}, the performance of the DTWA for power-law decay interaction has been discussed. Another possibility of interaction form is a step function type interaction, which approximately describes the spin-spin interaction of the Rydberg-dressed atoms \cite{Zeiher2016,Zeiher2017}. Investigating such a range-dependence of the validity timescale will be useful when the TWA is applied to analyzing such systems, in which the interaction range can be systematically controlled.

In this paper, we address the question how the validity timescale of the DTWA depends on the interaction range of the step function type interaction. 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 \cite{Pucci2016}. 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\'enyi entanglement entropy. In particular, we consider two-site R\'enyi entanglement entropy, which contains the information of one- and two-spin expectation values. In this paper, we develop a method to calculate the R\'enyi entropy within the framework of the DTWA. Our method can be applicable to not only the benchmark of the DTWA, but also the calculations of the R\'enyi 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 range of the step function type interaction. 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.~\ref{sec:model}, we explain our model and the DTWA. In Sec.~\ref{subsec:criterion}, we explain how to define the threshold time on which the DTWA is valid and compare the DTWA and tDMRG results. In Sec.~\ref{subsec:renyi_and_threshold}, we show the results of R\'{e}nyi entanglement entropy and threshold time for three different quantum-spin models. In Sec.~\ref{subsec:compare1D_and_2D}, we compare the results in 1D and 2D systems. In Sec.~\ref{sec:summary}, we summarize our results. In Appendix \ref{sec:details_of_DTWA}, we discuss the details of the derivation of the DTWA.  In Appendix \ref{sec:details_of_DTWA_sampling}, we explain the sampling scheme of the initial conditions.  In Appendix \ref{sec:derivation_of_Renyi_entropy}, we derive the expression of the R\'enyi entanglement entropy in the framework of the DTWA. In Appendix \ref{sec:comparison_with_DTWA_and_tDMRG}, we explain the details of the tDMRG calculations. In Appendix \ref{sec:new_algorithm}, we propose an algorithm to calculate the dynamics of the systems with the long-range interaction by using the tDMRG.

\section{Model and methods}\label{sec:model}
\subsection{Model}\label{subsec:model}


In this paper, we consider a family of quantum spin-$\frac{1}{2}$ systems in 1D and 2D, which is generally modeled by the Hamiltonian 
\begin{align}
\hat{H}=\frac{1}{2}\sum_{i,j, i \not=j}\sum_\mu J_{i j}^{\mu}\hat{S}_i^{\mu}\hat{S}_j^{\mu}+\bm{h}\cdot\sum_i\hat{\bm{S}}_i,\label{eq:system_Hamiltonian}
\end{align}
where $\hat{S}_i^{\mu}\; (\mu=x,y,z)$ is a spin-$\frac{1}{2}$ operator at site $i$, $J_{i j}^{\mu}$ is a  spin-exchange coupling between two distant sites, and $\bm{h}\equiv (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_{i j}^{\mu}$, we especially focus on three specific cases as follows:
\begin{align}
&J_{ij}^x=J_{ij}^y=0,\quad J_{ij}^z=-J_{i j}^{(D)},\quad (\text{Ising}),\label{eq:definition_of_Ising_interaction_1D}\\
&J_{ij}^x=J_{ij}^y=-J_{i j}^{(D)},\quad J_{ij}^z=0,\quad(\text{XY}),\label{eq:definition_of_XY_interaction_1D}\\
&J_{ij}^x=J_{ij}^y=J_{ij}^z=J_{ij}^{(D)},\quad(\text{Heisenberg}).\label{eq:definition_of_Heisenberg_interaction_1D}
\end{align}
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 the $x$ axis, i.e., $\bm{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
\begin{align}
J_{ij}^{(1)}\equiv
\begin{cases}
\vspace{0.2em}\dfrac{J}{r},\quad \text{if } |i-j|\le r,\\
0,\quad \text{otherwise}
\end{cases}
\label{eq:definition_1D_spin-spin_interaction}
\end{align}
where $J>0$ is an interaction strength and $r=1,2,\cdots$ is the interaction range. By contrast, the 2D cases ($D=2$), in which the lattice geometry is supposed to be square, are characterized by 
\begin{align}
J_{ij}^{(2)}\equiv 
\begin{cases}
\vspace{0.2em}\dfrac{J}{C_r},\quad \text{if }|\bm{R}_i-\bm{R}_j|\le R_r,\\
0,\quad\text{otherwise}
\end{cases}
\label{eq:definition_of_2D_spin-spin_interaction}
\end{align}
where $\bm{R}_i\equiv (R_{x i}, R_{y i})\equiv a(i_x, i_y)$ is the position of the $i$th lattice site, $a$ is the lattice constant, $i_x=1,2,\cdots, M_x$ and $i_y=1,2,\cdots, M_y$ are indices of the $i$th lattice sites, and $R_r$ is the distance between the $r$th neighboring sites. The total number of lattice points is given by $M=M_xM_y$. The constant $C_r$ is determined such that the following equation is satisfied:
\begin{align}
\max_i\sum_{j\not=i}J_{i j}^{(2)}=2J.\label{eq:definition_of_constant_Cr}
\end{align}
The explicit values of $C_r$ are given by $C_1=2, C_2=4, C_3=6, C_4=10, C_5=12,\cdots$. For later use, we define $N_r$ as
\begin{align}
N_{r}=
\begin{cases}
2r,\quad (\text{1D}),\\
2C_r,\quad (\text{2D}),
\end{cases}
\label{eq:definition_of_Nr}
\end{align}
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.~(\ref{eq:definition_1D_spin-spin_interaction}) and (\ref{eq:definition_of_2D_spin-spin_interaction}), are realizable in the experimental setups by means of Rydberg-dressed atoms \cite{Zeiher2016,Zeiher2017}. A system of such atoms is typically characterized by a soft-core type potential, so that interactions among atoms are almost constant in the short-distance regime and rapidly decay in the long-distance regime \cite{Henkel2010,Pupillo2010}. The couplings in Eqs.~(\ref{eq:definition_1D_spin-spin_interaction}) and (\ref{eq:definition_of_2D_spin-spin_interaction}) may describe such a situation approximately. 

In the subsequent sections we will investigate sudden-quench dynamics of the interacting spin systems \cite{Hauke2013,Schachenmayer2013,Lepori2016,Buyskikh2016} in order to characterize the limitation of the DTWA method. To be concrete, we especially consider the following direct-product wave functions as low-entangled initial states:
\begin{align}
|\psi(0)\rangle&=\prod_{i=1}^M|\leftarrow_i\rangle,\quad (\text{Ising}),\label{eq:initial_state_for_Ising_model}\\
|\psi(0)\rangle&=\prod_{i=1}^M|\downarrow_i\rangle,\quad (\text{XY}),\label{eq:initial_state_for_XY_model}\\
|\psi(0)\rangle&=\prod_{i=1}^{M/2}|\uparrow_{2i-1}\downarrow_{2i}\rangle,\quad (\text{Heisenberg}).\label{eq:initial_state_for_Heisenberg_model}
\end{align}
Here $\cket{\uparrow_i}$ and $\cket{\downarrow_i}$ denote the eigenstates of $\hat{S}_i^z$ while $\cket{\rightarrow_i}\equiv (\cket{\uparrow_i}+\cket{\downarrow_i})/\sqrt{2}$ and $\cket{\leftarrow_i}\equiv (\cket{\uparrow_i}-\cket{\downarrow_i})/\sqrt{2}$ represent the ones of $\hat{S}_i^x$. The corresponding discrete Wigner functions for these initial states are shown in Appendix~\ref{sec:details_of_DTWA_sampling}.

\subsection{Discrete phase-space approach to the R\'enyi entropy}\label{subsec:DTWA}

In this work, to characterize the performance of the DTWA, we focus on the second-order R\'enyi entanglement entropy defined by
\begin{align}
S^{(2)}_A(t)\equiv -\log\left({\rm Tr}\left\{[\hat{\rho}_A(t)]^2\right\}\right),\label{eq:definition_of_second_order_Renyi_entropy}
\end{align}
where $\hat{\rho}_A(t)\equiv {\rm Tr}_B\hat{\rho}(t)$ is the reduced density matrix associated with a subregion $A$ and $\hat{\rho}(t)$ is the density matrix of the whole system. The whole system is separated in real space into $A$ and $B$. In what follows, let us derive a discrete phase-space representation for $S^{(2)}_A(t)$.

An important ingredient to make the discrete phase-space representation is the phase-point operator ${\hat A}_{\bm{\alpha}}$ \cite{Schachenmayer2015}. 
For SU(2) spin systems, it generally takes the form
\begin{align}
\hat{A}_{\bm{\alpha}}\equiv \prod_{i=1}^M\left(\frac{1}{2}+\bm{r}_{\alpha_i}\cdot\hat{\bm{S}}_i\right), \label{eq:defintion_of_phase-point_operator}
\end{align}
where $\bm{\alpha}\equiv (\alpha_1,\alpha_2,\cdots,\alpha_M)$ with $\alpha_i \in \{(0,0), (0, 1), (1, 0), (1,1)\}$ denotes independent points in the discrete phase space and $\bm{r}_{\alpha_i}$ is a three-dimensional vector implying $\bm{r}_{(0,0)} = (+1,+1,+1)$, $\bm{r}_{(0,1)} = (-1,-1,+1)$, $\bm{r}_{(1,0)} = (+1, -1, -1)$, and $\bm{r}_{(1,1)} = (-1,+1,-1)$ \cite{Wootters1987}.
The density matrix at time $t$ can be written as \cite{Wootters1987}
\begin{align}
\hat{\rho}(t)=\sum_{\bm{\alpha}}W_{\bm{\alpha}}(0)\hat{A}_{\bm{\alpha}}(t),\label{eq:density_matrix_expressed_by_phase_point_operator}
\end{align}
where $W_{\bm{\alpha}}(0)$ is the discrete Wigner function at $t=0$ and we defined $\hat{A}_{\bm{\alpha}}(t)=e^{-i\hat{H}t/\hbar}\hat{A}_{\bm{\alpha}}e^{+i\hat{H}t/\hbar}$. The reduced density matrix $\hat{\rho}_A(t)$ is expressed by means of the discrete Wigner function at $t=0$ and $\hat{A}_{\bm{\alpha}}(t)$:
\begin{align}
\hat{\rho}_A(t)={\rm Tr}_B \left\{ \sum_{\bm{\alpha}}W_{\bm{\alpha}}(0)\hat{A}_{\bm{\alpha}}(t) \right\}.\label{eq:reduced_density_matrix_by_Discrete_Wigner_function}
\end{align}

The transformed phase-point operator $\hat{A}_{\bm{\alpha}}(t)$ contains a complete information about quantum many-body dynamics governed by ${\hat H}$. However, carrying out exact calculations for such an operator is generally impossible.  

In the descriptions of DTWA, the phase-point operator at $t$ is assumed to be factorized \cite{Schachenmayer2015}, i.e., 
\begin{align}
\hat{A}_{\bm{\alpha}}(t)\simeq \prod_{i=1}^M\left[ \frac{1}{2}+\bm{r}_i(t,\bm{\alpha})\cdot\hat{\bm{S}}_i \right]. \label{eq:ppo_dtwa}
\end{align}
The time-dependent coefficient $\bm{r}_i(t,\bm{\alpha}) \equiv 2\bm{S}_i(t)$ obeys a classical equation of motion obtained from a first-order BBGKY hierarchy truncation (see also Appendix~\ref{sec:details_of_DTWA})
\begin{align}
\hbar\frac{d}{d t}S_i^{\mu}(t)&=\epsilon_{\mu\beta\gamma}\left[h^{\beta}S_i^{\gamma}(t)+\sum_{k\not=i}J_{i k}^{\beta}S_k^{\beta}(t)S_i^{\gamma}(t)\right],\label{eq:classical_equation_of_motion_DTWA_1st}
\end{align}
where $\epsilon_{\mu\beta\gamma}$ is the Levi-Civita symbol and we used the Einstein notation for repeated Greek indices. 
Inserting Eq.~(\ref{eq:ppo_dtwa}) into Eq.~(\ref{eq:definition_of_second_order_Renyi_entropy}), we arrive at a DTWA expression of the R\'enyi entropy
\begin{align}
S_{A}^{(2)}(t) \approx -\log\left\langle\left\langle\prod_{i\in A}\left[\frac{1}{2}+2\bm{S}_i(t)\cdot\bm{S}'_i(t)\right]\right\rangle\right\rangle, \label{eq:1st_order_BBGKY_renyi_expression}
\end{align}
where $\bm{S}_i(0)=\bm{r}_{\alpha_i}/2$ and $\bm{S}'_i(0)=\bm{r}_{\alpha_i'}/2$. The doubled angular brackets mean a phase-space average weighted with two initial Wigner functions 
\begin{align}
\left\langle\left\langle f_{\bm{\alpha},\bm{\alpha}'} \right\rangle\right\rangle \equiv \sum_{\bm{\alpha},\bm{\alpha}'}W_{\bm{\alpha}}(0)W_{\bm{\alpha}'}(0)f_{\bm{\alpha},\bm{\alpha}'}.
\end{align}


For direct product states such as Eqs.~(\ref{eq:initial_state_for_Ising_model}), (\ref{eq:initial_state_for_XY_model}), and (\ref{eq:initial_state_for_Heisenberg_model}), $W_{\bm{\alpha}}(0)$ is factorized as $W_{\bm{\alpha}}(0)=\prod_{j=1}^{M}w_{\alpha_{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 nonlinear systems with a nonzero spin-exchange coupling, the R\'enyi 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 {\it complexity} in a solution of Eq.~(\ref{eq:classical_equation_of_motion_DTWA_1st}) 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.~(\ref{eq:classical_equation_of_motion_DTWA_1st}) arises in a classical trajectory of an enlarged phase space for a second-order BBGKY method \cite{Pucci2016}. The underlying idea of this method is to regard a nonseparable part of ${\hat S}^{\mu}_{i}{\hat S}^{\nu}_{j}$, which is represented by $c^{\mu\nu}_{ij}$ in a replacement ${\hat S}^{\mu}_{i}{\hat S}^{\nu}_{j} \rightarrow S_i^{\mu}S_j^{\nu}+c_{ij}^{\mu\nu}$, as an additional mechanical variable and define an approximately closed equation of motion for $S^{\mu}_{i}$ and $c_{ij}^{\mu\nu}$. We note that the first-order BBGKY truncation leads to the time-evolving equation given by Eq.~(\ref{eq:classical_equation_of_motion_DTWA_1st}). In Sec.~\ref{sec:results}, we will also exploit the second-order BBGKY method to compute the R\'enyi entropy, especially for a subregion of two sites. The detail of the BBGKY formulation will be presented in Appendix \ref{sec:details_of_DTWA}.


In this paper, we numerically solve Eq.~(\ref{eq:classical_equation_of_motion_DTWA_1st}) for the first-order BBGKY and  Eqs.~(\ref{eq:BBGKY_first_order_equation_re}) and (\ref{eq:BBGKY_second_order_re}) for the second-order BBGKY by using a fourth-order Runge-Kutta method. A time step $\Delta t$ is taken to be $\Delta t=10^{-3}\hbar/J$. We have checked that our results are converged for this $\Delta t$ in our simulation timescale. We use $M=100$ in 1D and $M_x=M_y=14$ in 2D. The computational costs for first- and second-order BBGKY equation scale as $O(M)$ and $O(M^2)$, respectively.

\section{Results}\label{sec:results}

\subsection{Criterion for the validity of the DTWA and comparison with tDMRG results}\label{subsec:criterion}

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 first- and second-order BBGKY results is small, the DTWA gives a good approximation. Here, a question arises: Which physical quantities are appropriate for comparing the first- and second-order BBGKY results? We propose that the second-order R\'{e}nyi entanglement entropy is a suitable quantity for confirming the validity of the DTWA. One advantage to use the R\'enyi entropy is that it is an unbiased quantity compared with other physical quantities such as 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 the $z$ axis, $\hat{S}_{\rm tot}^z\equiv \sum_i\hat{S}_i^z$ is conserved so that it is not appropriate for examining the validity of the DTWA.

In this paper, we calculate mean two-site R\'enyi entropies for 1D and 2D, which are defined by
\begin{align}
S^{(2)}(t)\equiv \frac{1}{M-1}\sum_{i<j,|i-j|=1}S_{ij}^{(2)}(t),\quad (\text{1D}),\label{eq:definition_mean_Renyi_1D}\\
S^{(2)}(t)\equiv \frac{1}{(M_x-1)M_y}\sum_{i<j,|R_{xi}-R_{xj}|=a}S_{ij}^{(2)}(t),\quad (\text{2D}).\label{eq:definition_mean_Renyi_2D}
\end{align}
where $S_{ij}^{(2)}(t)$ is the two-site R\'enyi entropy. Let us mention differences between our formulation of the R\'{e}nyi entropy and that of previous works \cite{Acevedo2017,Lepoutre2019,Zhu2019}. In these previous works, they calculated the single- or two-site R\'{e}nyi entropy from the expressions of the single- or two-site reduced density-matrix operator because the matrix elements of these reduced density matrices can be constructed by the expectation values of $\hat{S}_i^{\mu}$ and $\hat{S}_i^{\mu}\hat{S}_j^{\nu}$, which can be obtained by the DTWA. The advantage of our formulation is that it allows us to calculate the R\'{e}nyi entropy for multiple sites. It is easy to calculate the multiple-sites R\'{e}nyi entropy in the first-order BBGKY. We have checked the statistical convergence of the results even in the case of bipartite R\'enyi entropy. For more details, see Appendix \ref{sec:derivation_of_Renyi_entropy}.  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\'enyi entropy by using the first- and second-order BBGKY equations. In order to quantify the difference between the first- and second-order results, we define
\begin{align}
\Delta(t)\equiv \left|e^{-S_{\rm 1st}^{(2)}(t)}-e^{-S_{\rm 2nd}^{(2)}(t)}\right|/e^{-S_{\rm 1st}^{(2)}(t)},
\end{align}
where $S_{\rm 1st}^{(2)}(t)$ and $S_{\rm 2nd}^{(2)}(t)$ are the mean two-site R\'enyi entropy obtained by the first- and second-order BBGKY equation, respectively. This quantity represents a relative error of the first- and second-order results. The reason why we do not use the relative error of $S_{\rm 1st}^{(2)}(t)$ and $S_{\rm 2nd}^{(2)}(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_{\rm th}$, at which $\Delta(t)$ exceeds a small positive number $\epsilon$.
In this paper, we use $\epsilon=\frac{1}{10}$.

We note that the second-order BBGKY equation is numerically unstable as pointed out in Refs.~\cite{Orioli2017,Czischek2018}. In fact, we find the divergent behavior of the second-order BBGKY equation. For example, this behavior can be seen in $r=1$ results in Figs.~\ref{fig:compare_renyi_entropy_model_1D} and \ref{fig:renyi_entropy_1D_and_2D}. This is an intrinsic property of the second-order BBGKY equation. We have checked that this divergence behavior is not an artificial one because it does not depend on the choice of the time step $\Delta t$.

The small difference between the first- and second-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, we perform the comparison between the DTWA results and the tDMRG method in the 1D cases. See details of the tDMRG calculations in Appendix \ref{sec:comparison_with_DTWA_and_tDMRG}. Figure \ref{fig:compare_renyi_entropy_model_1D} shows the results of the mean two-site R\'enyi entropy for (a) Ising model, (b) XY model, and (c) Heisenberg model, respectively. For long-range Hamiltonians, we implement tDMRG with utilizing swap gates as detailed in Appendix~\ref{sec:new_algorithm}. We can see that the DTWA results are good agreement with the tDMRG results in a short time scale ($t\sim \hbar/J$) for all cases. In the long-range interacting cases, the DTWA results quantitatively reproduce the tDMRG results even in the long-time scale ($t\sim 10\hbar/J$). However, in the long-time scale, the second-order BBGKY results for some parameters are divergent (We can see this from the fact that there is no data point for second-order BBGKY results in the long-time regime.). This means that our criterion does not work in the long-time scale. Therefore, we expect that our criteria for the performance of the DTWA work in the short- and intermediate-time scale. In the following, we focus on these timescales.

\begin{figure}[t]
\centering
\includegraphics[width=8.5cm,clip]{compare_entropy_1d_ver1.pdf}
\caption{Comparison with the DTWA and tDMRG results. Mean two-site R\'{e}nyi entropy in 1D. (a) Ising model for $h^x=0.5J$. (b) XY model for $h^x=2.0J$. (c) Heisenberg model for $h^x=10.0J$. The black solid and dotted lines represent the tDMRG results. The open and closed symbols represent the 1st and 2nd order BBGKY results, respectively. $r$ is the interaction range.}
\label{fig:compare_renyi_entropy_model_1D}
\end{figure}%

\subsection{R\'{e}nyi entropy and threshold time}\label{subsec:renyi_and_threshold}

In this subsection, we show the results of R\'enyi entropy and threshold time for Ising, XY, and Heisenberg models.  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 models have only discrete symmetry while the Heisenberg model has continuous spin-rotation symmetry around the magnetic field.

\subsubsection{Ising model}\label{subsubsec:ising}


\begin{figure*}[t]
\centering
\includegraphics[width=17.0cm,clip]{entropy_1d_and_2d_ver1.pdf}
\caption{Time evolution after the sudden quench of the mean two-site R\'{e}nyi entropy in 1D and 2D. (a), (b) Ising model for $h^x=0.5J$. (c), (d) XY model for $h^x=2.0J$. (e), (f) Heisenberg model for $h^x=10.0J$. The open and closed symbols represent the first- and second-order BBGKY results, respectively. $r$ denotes the interaction range.}
\label{fig:renyi_entropy_1D_and_2D}
\end{figure*}%


\begin{figure}[t]
\centering
\includegraphics[width=8.0cm,clip]{threshold_time_1d_2d_ferro_Ising_mx100_xpolarized_ver1.pdf}
\caption{Threshold time as a function of number of interacting sites for the Ising model. (a) One dimension. The black solid and dashed lines represent $N_r^{0.5}$ and $N_r^{0.75}$ for guide to eye. (b) Two dimensions. The black solid and dashed lines represent $N_r^{0.4}$ and $N_r^{0.6}$ for guide to eye.}
\label{fig:T_th_for_Ising_model}
\end{figure}%

Here, we show the results of the mean two-site R\'{e}nyi entropy for the Ising model under the transverse magnetic field in Fig.~\ref{fig:renyi_entropy_1D_and_2D} (a) for 1D and Fig.~\ref{fig:renyi_entropy_1D_and_2D} (b) for 2D. In these cases, the initial condition is the fully $-x$-polarized state [see Eq.~(\ref{eq:initial_state_for_Ising_model})]. This state is the exact ground state when $h^x\to \infty$. 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 many-body systems \cite{Eisert2010}. 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 \cite{Schachenmayer2013}. In the $r\to\infty$ limit, our model becomes the Lipkin-Meshkov-Glick model \cite{Vidal2004,Latorre2005}. 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\'enyi entropy has same tendency. 

From the results of the R\'enyi entropy, we can obtain the threshold time $T_{\rm th}$ as a function of $N_r$. The results are shown in Fig.~\ref{fig:T_th_for_Ising_model} (a) for 1D and Fig.~\ref{fig:T_th_for_Ising_model} (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.~\ref{fig:T_th_for_Ising_model}. We can also see that $T_{\rm th}$ in the large-$h^x$ region is large compared to that in the 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 is mainly driven by the magnetic field. In this reason, $T_{\rm th}$ in the large-$h^x$ region is longer than that of the small-$h^x$ region.


\subsubsection{XY model}\label{subsubsec:XY}

Next, we show the results of the XY model under the magnetic field along the $-x$ direction in Fig.~\ref{fig:renyi_entropy_1D_and_2D} (c) for 1D and \ref{fig:renyi_entropy_1D_and_2D} (d) for 2D. In this case, the initial condition is the fully $-z$-polarized state [see Eq.~(\ref{eq:initial_state_for_XY_model})]. This state is an exact eigenstate when $h^x=0$.

\begin{figure}[t]
\centering
\includegraphics[width=8.0cm,clip]{threshold_time_1d_2d_ferro_xy_mx100_zpolarized_ver1.pdf}
\caption{Threshold time as a function of the number of interacting sites for the XY model. (a) One dimension. The black solid and dashed lines represent $N_r^{0.5}$ and $N_r^{0.3}$ for guide to eye. (b) Two dimensions. The black solid and dashed lines represent $N_r^{0.5}$ and $N_r^{0.3}$ for guide to eye.}
\label{fig:T_th_for_XY_model}
\end{figure}%

From these results, we obtain the threshold time as a function of $N_r$. The results are shown in Fig.~\ref{fig:T_th_for_XY_model} (a) for 1D and Fig.~\ref{fig:T_th_for_XY_model} (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.~\ref{fig:T_th_for_XY_model}). 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.

\subsubsection{Heisenberg model}\label{subsubsec:Heisenberg}

We show the results of the Heisenberg model under the magnetic field along the $-x$ direction in Fig.~\ref{fig:renyi_entropy_1D_and_2D} (e) for 1D and Fig.~\ref{fig:renyi_entropy_1D_and_2D} (f) for 2D. In this case, the initial condition is the Ne\'el state [see Eq.~(\ref{eq:initial_state_for_Heisenberg_model})]. 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. 

\begin{figure}[t]
\centering
\includegraphics[width=8.0cm,clip]{threshold_time_1d_2d_antiferro_heisenberg_mx100_neel_ver1.pdf}
\caption{Threshold time as a function of the number of interacting sites for the Heisenberg model. (a) One dimension. The black solid and dashed lines represent $N_r^{0.5}$ and $N_r^{0.8}$ for guide to eye. (b) Two dimensions. The black solid and dashed lines represent $N_r^{0.5}$ and $N_r^{1.0}$ for guide to eye.}
\label{fig:T_th_for_heisenberg_model}
\end{figure}%

The threshold time as a function of $N_r$ is in Fig.~\ref{fig:T_th_for_heisenberg_model}. We can see that the threshold time increases as power law of $N_r$ (see solid and dotted lines in Fig.~\ref{fig:T_th_for_heisenberg_model}). 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_{\rm th}$. For example, the oscillation of $h^x=10J$ shown in Fig.~\ref{fig:T_th_for_heisenberg_model} (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 $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-nearest-neighbor 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_{\rm th}$ as a function of the interaction range appears. 


\subsection{Comparison with one- and two-dimensional results}\label{subsec:compare1D_and_2D}

\begin{figure}[t]
\centering
\includegraphics[width=8.5cm,clip]{compare_threshold_total_ver1.pdf}
\caption{Comparison with the threshold time for Ising model in 1D and 2D. (a) Ising model. (b) XY model. (c) Heisenberg model. Open and closed symbols represent 1D and 2D results, respectively.}
\label{fig:T_th_for_1D_2D}
\end{figure}%

Here, we compare the results of the threshold time in 1D and 2D. The results are shown in Fig.~\ref{fig:T_th_for_1D_2D} (a) for the Ising case, Fig.~\ref{fig:T_th_for_1D_2D} (b) for the XY case, and Fig.~\ref{fig:T_th_for_1D_2D} (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_{\rm 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 spins 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.~\ref{fig:T_th_for_1D_2D}.). 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_{\rm c}=0.5J$ for 1D \cite{Sachdev_text_book} and $h_{\rm c}^x\simeq 3.044J$ for 2D \cite{Jongh1998,Rieger1998} 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 in the 1D case. Therefore, the 1D and 2D results deviate clearly. 




\section{Summary}\label{sec:summary}

In summary, we investigated the timescale on which the DTWA can quantitatively describe quantum dynamics of spin-$\frac{1}{2}$ models with a step function type spin-spin interaction in the short- and intermediate timescales. In order to corroborate this, we developed a formulation of the R\'enyi entropy within the DTWA framework. Using this formulation, we evaluate the R\'enyi entropy after a sudden quench in the Ising, XY, and Heisenberg models under a uniform magnetic field in 1D and 2D. The R\'enyi entropy is calculated by the DTWA and its extension including the second-order 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\'enyi entropy in the DTWA and second-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 also an 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 long-time evolution. Therefore, it is important to confirm the validity of the DTWA in the sweep dynamics case. 

\begin{acknowledgments}

The tDMRG calculations in this work are performed with ITensor library \cite{Fishman2007}. This work was supported by MEXT Quantum Leap Flagship Program (MEXT Q-LEAP) Grant Number JPMXS0118069021 (M.K. and I.D.), JST CREST No. JPMJCR1673 (I.D.), and JSPS KAKENHI Grant Number JP20K14389 (M.K.), JP20K14377 (S.G.), JP18K03492 (I.D.), and JP18H05228 (I.D.).

\end{acknowledgments}


