SU(3) truncated Wigner approximation for strongly interacting Bose gases

We develop and utilize the SU(3) truncated Wigner approximation (TWA) in order to analyze far-from-equilibrium quantum dynamics of strongly interacting Bose gases in an optical lattice. Specifically, we explicitly represent the corresponding Bose--Hubbard model at an arbitrary filling factor with restricted local Hilbert spaces in terms of SU(3) matrices. Moreover, we introduce a discrete Wigner sampling technique for the SU(3) TWA and examine its performance as well as that of the SU(3) TWA with the Gaussian approximation for the continuous Wigner function. We directly compare outputs of these two approaches with exact computations regarding dynamics of the Bose--Hubbard model at unit filling with a small size and that of a fully-connected spin-1 model with a large size. We show that both approaches can quantitatively capture quantum dynamics on a timescale of $\hbar/(Jz)$, where $J$ and $z$ denote the hopping energy and the coordination number. We apply the two kinds of SU(3) TWA to dynamical spreading of a two-point correlation function of the Bose--Hubbard model on a square lattice with a large system size, which has been measured in recent experiments. Noticeable deviations between the theories and experiments indicate that proper inclusion of effects of the spatial inhomogeneity, which is not straightforward in our formulation of the SU(3) TWA, may be necessary.


I. INTRODUCTION
Quantum simulators built with synthetic quantum platforms that are highly controllable have been applied for studying quantum many-body physics in and out of equilibrium. Examples of such quantum simulators include ultracold gases in optical lattices [1][2][3][4][5], Rydberg atoms in optical tweezer arrays [6], trapped ions [7], and superconducting circuits [8,9]. Of particular interest is far-from-equilibrium quantum dynamics of isolated many-body systems described by the tight-binding Hubbard-type models, which can be simulated with ultracold gases in optical lattices. The quantitative accuracy of such analog quantum simulators for non-equilibrium lattice systems has been examined through direct comparisons with outputs from exact computational methods for some special cases, such as the exact diagonalization for small systems [10] and the matrix-product-state (MPS) approaches for onedimensional (1D) systems [3,4]. With the high accuracy confirmed, results obtained from optical-lattice quantum simulators have been exploited in order to test approximate computational methods for quantum many-body dynamics in higher dimensions. For instance, it has been shown in Ref. [5] that the non-equilibrium dynamical mean-field theory can quantitatively capture dynamics of the three-dimensional (3D) Hubbard model subjected to a periodic driving. Moreover, in Ref. [11], the Gross-Pitaevskii truncated-Wigner approximation (GPTWA), * knagao@physnet.uni-hamburg.de which is a semiclassical phase-space method on the basis of the GP mean-field theory [12,13], has been directly compared with experimental data regarding dynamics of the 3D Bose-Hubbard model in a weakly interacting regime after a quantum quench. It has been shown that the outputs of GPTWA with no free parameter are in good agreement with experimental data for early-time regions.
In recent years, some experimental works have explored quantum quench dynamics of strongly interacting ultracold gases in two-dimensional (2D) and 3D optical lattices [14,15]. In Ref. [15], an experimental group at Kyoto University has studied sudden-quench dynamics of equal-time single-particle correlation functions for a strongly interacting 174 Yb gas loaded into a deep 2D lattice. In contrast to 1D systems, it is generally hard to numerically simulate time evolution of correlation functions in 2D and 3D even on a short timescale. It has been found in Ref. [15] that the ordinary GPTWA cannot fully capture characteristic properties of the correlation propagation after sudden quenches, e.g., peak and dip properties observed in the correlation signals and saturated values of the correlation at relatively long times. This can be attributed to the fact that in the strongly interacting regime the adequate classical limit of the system is not condensates of coherent bosons described by the GP theory.
In Ref. [16], Davidson and Polkovnikov have introduced a promising phase-space approach for analyzing strongly interacting Bose-Hubbard systems. This method is called the SU(3) TWA [hereafter SU(3)TWA]. For sufficiently large local interactions, the Bose-Hubbard model reduces to an effective pseudospin-1 model acting on a projected Hilbert space [17,18]. In the SU(2) TWA method, which is typically discussed and used in the context of experiments of large-spin systems and arrays of trapped Rydberg atoms [19,20], this effective model is treated as a Hamiltonian consisting of the SU(2) spin operators for S = 1 [13]. However, for the SU(3) TWA, the model is translated into a Hamiltonian consisting of SU(3) matrices, which gives an alternative phase-space representation of the system with extra five dimensions in addition to the three dimensions of the SU(2) phase space. Since the local interaction terms of the effective model can be linearized in the SU(3) matrices, the local particle and hole fluctuations, which produce key effects on the dynamical properties of the strongly interacting regime, are accurately captured at the level of the semiclassical approximation [16]. The TWA method based on the GP trajectories is not suitable to formulate those fluctuations in the strongly interacting limit, just as the Bogoliubov approximation for weakly interacting dilute Bose gases fails to describe the quantum phase transitions to the Mott-insulator phases at low temperatures [21]. We therefore expect that, the SU(3) TWA may simulate the dynamics in the strongly interacting regime of the experiment [15], beyond the capability of the GPTWA, and also the SU(2) TWA.
In their original work, the performance of the SU(3)TWA was tested by applying it to a fully connected spin-1 model, which has an all-to-all spin-exchange (or hopping) term and can be numerically diagonalized even at a large size. However, its quantitative accuracy in realistic cases, where the hopping reaches only nearest neighbors and the system size is large, has not been examined so far. Furthermore, an effective model that they used to describe Bose-Hubbard systems is valid only for highfilling cases. Therefore, their formulation is not directly applicable to unit-filling Bose-Hubbard systems, which are typically considered in the context of the quantumsimulation studies. We note that a numerical calculation of the SU(3)TWA for a unit-filling experimental setup has been presented in Ref. [22]; however, its explicit formalism has not been provided so far.
The goal of this paper is to examine the performance of the SU(3)TWA in simulating quench dynamics of strongly interacting Bose gases in a 2D optical lattice [15]. We extend the previous formalism, which was applied to an effective pseudospin-1 model for the Bose-Hubbard model with large filling factors and strong interactions [23,24], to the unit-filling case [17,18] corresponding to the experimental setup. As a technique to evaluate the phase-space integration emerging in the SU(3)TWA, we will employ two different approaches, i.e., the Gaussian approximation for the (continuous) Wigner function [16] and the discrete TWA (DTWA) approach [19,20,25]. In particular, the DTWA approach is thought to be better than the Gaussian approach. Indeed, the numerical sampling of the DTWA can be readily carried out without approximation of the probability distribution functions (see also Refs. [19,25]). In this paper, we also study the performance of a DTWA sampling for the SU(3)TWA via large-scale numerical simulations for a fully connected spin-1 model. A numerical simulation on the basis of the DTWA scheme will be compared with the experimental data as well as that of the Gaussian approximation.
The remainder of this paper is organized as follows: In Sec. II, we introduce an effective pseudospin-1 model for the Bose-Hubbard Hamiltonian in a strongly interacting regime and a fully connected spin-1 model, respectively. In Sec. III, we formulate the SU(3)TWA for the effective model. In Sec. IV, we study the Gaussian approximation and the DTWA approach for SU(3) phase-space variables. In Sec. V, using the SU(3)TWA, we calculate quench dynamics of equal-time single-particle correlation functions for a strongly interacting Bose gas in a 2D optical lattice. There, we compare some semiclassical results with actual experimental data obtained in Ref. [15]. In Sec. VI, we conclude this paper and present outlooks for future studies.

II. MODELS
In this paper, we study time evolution of a strongly interacting Bose gas loaded into an optical lattice. To describe this system, we consider the Bose-Hubbard Hamiltonian on a certain lattice structure [26,27] whereâ † j andâ j are the creation and annihilation operators of bosons at site j. The angular brackets j, k indicate a nearest-neighbor link on the lattice. The real parameters J and U denote the hopping amplitude and interaction strength, respectively. A ratio of the parameters, U/J, can be widely controlled by tuning the opticallattice depth [15] or utilizing a Feshbach-resonance technique [14].
In a strongly interacting regime of Eq. (1), fluctuations of occupation per site are sufficiently suppressed from the mean fillingn. Therefore, only a subset of local Fock states near the mean filling is relevant to strongly interacting dynamics governed by Eq. (1). If the interaction is sufficiently strong, i.e., U/(nJ) 1, one can safely assume that only three Fock states, i.e., |n − 1 j , |n j , |n + 1 j are relevant to time evolution of the interacting bosons. In a projected Hilbert space spanned by such a local basis, the Bose-Hubbard Hamiltonian (1) is approximated as an effective pseudospin-1 model [17,18], which is given bŷ where δν − = 1 + 1/n − 1 andŜ ± j =Ŝ x j ± iŜ y j . The pseudospin operatorŜ µ j (µ = x, y, z) satisfies the SU(2) Lie algebra The three-leg tensor µνγ is the fully antisymmetric structure constant satisfying xyz = − yxz = yzx = · · · = 1. Hereinafter, the repeated greek indices indicate the contraction of tensors. It should be noticed that if one takes the high-filling limit, i.e.,n 1, the effective model is simplified [23,24] aŝ where B can be interpreted as a magnetic field applied along the z-axis. In the previous work [16], the SU(3)TWA was applied to this high-filling model defined on a cubic lattice. However, in order to analyze experimental systems with a setup ofn = 1 as realized in Ref. [15], it is required to use Eq. (2) rather than the high-filling model. In Sec. III, we will explain how one generalizes the SU(3)TWA to Eq. (2). In Sec. IV, we present detailed investigations on Monte Carlo integration methods employed for SU(3)TWA simulations. To examine quantitative validity of our numerical approaches, especially a DTWA approach for SU(3) phase-space variables, we will revisit a fully connected spin-1 model, which is a model studied in Ref. [16]. The Hamiltonian of the fully connected model is given bŷ The spin-exchange coupling term describes all-to-all connections between distant spin operators. Hence, each lattice point has a coordination number z = M − 1. As M increases, the valid timescale of the SU(3)TWA for this model becomes longer for a certain U/(zJ) [16]. Furthermore, due to a characteristic property described in Appendix A, exact quantum dynamics of this model can be easily simulated by using classical computers even for a considerably large M . Accordingly, the fully connected model is suitable for examining the performance of the sampling methods. See also Appendix A for details about how to implement exact numerical simulations of this model.

III. SU(3) TRUNCATED-WIGNER APPROXIMATION
The first step for building the SU(3)TWA for spin-1 models is to rewrite their Hamiltonian by means of eight numbers of SU(3) matrices [16]. Let us consider a set of SU(3) generators {X µ } (µ = 1, · · · , 8) obeying the SU(3) Lie algebra Here f µνγ is a fully antisymmetric structure constant accompanied by the SU (3) group. If we take the Jordan-Schwinger mapping into account, each generator can be written in the bi-linear form of the SU(3) Schwinger bosons [17,18,23] To reproduce the original Hilbert space, the particle number must be preserved per site by a constraint nb † nbn = 1. The value of f µνγ depends on the detail of T mn µ . Our choice for T µ will be shown later in Eq. (11), and the corresponding f µνγ will be given by Eq. (13). The SU(3) matrices T µ form a complete set of 3 × 3 matrices, so that an arbitrary local operator acting on the three-state Hilbert space is expressed as a linear combination of these matrices. Using this property, one can linearize local interaction terms in spin-1 models, such as U 2 j (Ŝ z j ) 2 , in terms of SU(3) matrices. Specifically for the effective model (2), if the interaction U is sufficiently large compared tonJ characterizing the hopping term, then the Hamiltonian is regarded as being almost linear in SU(3) matrices. Therefore, the SU(3)TWA for this model is expected to be valid during a long timescale. Furthermore, if the hopping term is negligible, the SU(3)TWA becomes exact at all times because there exists no truncation error stemming from higher-order derivatives of the time-evolving equation for the Wigner function [13].
Let us generalize the SU(3)TWA formalism to the arbitrary filling model (2). First, we express the effective Hamiltonian by means of the local SU(3) generators denoted byX A key point is that the local interaction term of the SU(2) spin operators is translated into a linear combination of such SU(3) generators as Then, we make a Wigner-Weyl transform of the Hamiltonian and obtain a classical Hamiltonian for the SU(3) phase-space variables where δν + = 1 + 1/n + 1. The SU(3)TWA states that within a semiclassical approximation the time evolution of the expectation value of an operatorΩ, i.e., Ω (t) , can be represented in terms of saddle-point trajectories of SU(3) variables, which are governed by H W and weighted with a Wigner quasi-probability distribution function where 0,µ is the integration measure and Ω W is a Weyl symbol ofΩ. The classical trajectory X cl (t) obeys Hamilton's equation associated with the SU(3) Lie algebra This equation of motion is integrated under an initial condition X (j) 0,µ is distributed according to W (X 0 ). The width of the Wigner function gives quantum-fluctuation corrections to saddlepoint or mean-field results, which formally correspond to the time-dependent Gutzwiller approximation with a single-site cluster consisting of three levels.
If we take the high-filling limit for the classical Hamiltonian (8), all the terms involving X disappear. Therefore, these additional variables are responsible for different consequences between the high and low filling descriptions. It should be noted that a constant term has been eliminated from Eq. (8) because it does not affect Eq. (10). The above formalism will be used in Sec. V to analyze the experimental setup in Ref. [15].
In this paper, we will utilize the following representation for the SU(3) matrices, according to the notations by Davidson and Polkovnikov [16]: These matrices are normalized as It is confirmed that, in this specific representation, nonzero values of f µνγ are given by Of course, this is not the unique choice. Instead of this representation, one can also use the Gell-Mann matrices, which are more familiar in high-energy physics [28].

IV. MONTE CARLO INTEGRATIONS
In this section, we study Monte Carlo integration methods for evaluating the phase-space integration of the initial Wigner function. In Ref. [16], an approximate Gaussian-Wigner function has been used to perform numerical simulations. This Gaussian approximation is a simple and efficient prescription for resolving a kind of minus-sign problem in TWA simulations, which means that the exact Wigner function defined by means of the Schwinger-boson coherent states typically takes negative values. In Sec. IV, to simulate the experimental setup, we will indeed employ the Gaussian approach.
As an alternative sampling scheme that allows us to avoid the appearance of negative-valued Wigner function, we also use a DTWA approach [19]. This approach is formulated on the basis of the discrete-Wigner representation of a finite Hilbert space quantum system. The concept of the discrete-Wigner representation has been invented by Wootters in Ref. [29]. In this section, by extending the previous DTWA method for SU(2) spin systems [19], we develop a DTWA approach suited for the SU(3)TWA. To this end, we will introduce phasepoint operators for the SU(3) generators, each of which is represented as a 3 × 3 matrix.

A. Gaussian approximation
In the Gaussian approximation for exact Wigner functions, an appropriate Gauss distribution is used to approximately express initial density matrices within a class of positive-definite functions [16]. To be specific, let us consider a fully polarized state along the x axis, i.e., ρ 1 = |S x = 1 S x = 1|. Its matrix form is given by To obtain the corresponding Gaussian-Wigner function, we make the following ansatz with free parameters R = (R µν ), m = (m µ ), and σ = (σ µ ): This distribution defines the first and second order moments of the SU(3) phase-space variables The free parameters are determined such that the Gaussian-Wigner function exactly reproduces the firstand second-moments of the density matrix (14), i.e., The angular brackets in the right-hand side mean the quantum-mechanical average withρ 1 . To determine R in practice, we diagonalize an 8 × 8 matrix corresponding to a connected and symmetrized correlation function with respect to the density matrix The eight-dimensional matrix R is constructed from the eigenvectors, which are obtained when C µν is diagonalized. Each eigenvalue gives the squared covariance σ 2 µ . The mean value m µ is the rotation of the vector ( X µ ), i.e., R µν X ν = m µ . The direct calculation leads to the following result: With these parameters, the Gauss distribution (15) randomly generates the phase-space variables reproducing the exact low-order moments of the state in Eq. (14).
In the projected Hilbert space for the effective pseudospin-1 models, the deep Mott-insulator state, which is approximately realized in a sufficiently deep optical lattice, is expressed as a direct product state of ρ 2 = |S z = 0 S z = 0|. The matrix form ofρ 2 is given by The corresponding parameters of the Gauss distribution function are calculated as

B. SU(3) discrete-Wigner representation
Let us consider a discrete-Wigner representation for a finite-level system, whose Hilbert space is spanned by three basis vectors {|0 , |1 , |2 }. The key building blocks for this representation are the so-called phase-point oper-atorsÂ α , which are 3 × 3 matrices acting on the Hilbert space. The integer index α = (a 1 , a 2 ) (a 1 , a 2 = 0, 1, 2) expresses a point in the discrete phase space Γ, which now contains nine points. The phase-point operators are also called the Stratonovich-Weyl kernels [30].
The phase-point operators are important because they define a Wigner-Weyl transform of quantum-mechanical operators. In the discrete-Wigner representation, the Weyl symbol of an operatorΩ is defined as its projection onto a point α ∈ Γ: Specifically, such a projection of a given density matrix ρ leads to the discrete-Wigner function The pre-factor 1/3 is needed to ensure the unity normalization of the Wigner function α∈Γ w α = 1, see also below. By analogy with continuous cases, where the coordinate and momentum operators (x,p) define a continuous phase-point operator, the discrete phase-point operator should have the following properties [29]: 1. Hermiticity:Â † α =Â α for any α ∈ Γ. Then, the phase space functions are real as long as the corresponding operators are Hermitian.
Such discrete phase-point operators can also be made for general cases where the Hilbert space is in N dimensions (N ≥ 2 is a primal number) [29]. Furthermore, it is possible to construct a discrete number-phase representation for Bose systems, whose Hilbert space is spanned by generators of the Heisenberg-Weyl group, and it provides a DTWA-like semiclassical approximation for their quantum dynamics if the allowed occupancy of particles is sufficiently large [31].
As an inverse transformation of Eqs. (28) and (29), the operatorsΩ andρ are linearly expanded inÂ α such that Then, the expectation value ofΩ forρ reads as The summation in the last expression is taken over the whole Γ. In the second equality, we have used the trace orthogonality ofÂ α . The concrete forms of Ω α and w α are specified after one determinesÂ α for all α = (a 1 , a 2 ) such that they satisfy the required conditions as presented above. If we adopt Wootters's representation of the phase-point operators [29], we have It is convenient to expandÂ α in the generators of the SU(3) Lie algebra, i.e., Its projection coefficient x µ (α) = Tr[Â αXµ ] is the discrete Weyl symbol ofX µ . After direct calculations, we obtain the following discrete phase-space variables for A (0) α : x 4 (α) = 2δ a1,1 cos 4πa 2 3 , as a combined eight-dimensional vector on each phase point, α = (0, 1), (1,2), and (2, 0) correspond to the following configurations, respectively: Two classical spins x(α) and x(α ) at different points α = α are not orthogonal to each other. Indeed, these have a finite inner product even for α = α In the DTWA simulation, such discretized spins are randomly distributed according to w α and give a set of initial conditions for the classical trajectories. The discussions of the DTWA for the SU(3) systems will be presented in Sec. IV C.
To clarify the sampling weight of DTWA simulations, which will be used in the following sections, let us calculate the discrete Wigner function for the Mott insulator state [Eq. (24)] by using A (0) α . It results in a positivedefinite distribution function This result means that in the Mott-insulator state three configurations at α = (1, 0), (1, 1), (1, 2) are realized with equal probability 1 3 while other ones have the zero probability. Therefore, we can directly evaluate the average with the Wigner function in numerics without further approximation of the distribution function. However, the positivity of Eq. (36) is not a general property. For example, the x-polarized state in Eq. (14) yields oscillatory terms in the distribution While the first term with δ a1,1 is always positive, the second term with δ a1,0 and δ a1,2 takes negative values due to the oscillating contributions.
As mentioned in previous works [29,32], the definition of the phase-point operators is not unique. In general, there exists a non-singular (or regular) transformation, A α →Ŝ −1Â αŜ , which retains the required properties of the phase-point operators [29]. This type of ambiguity will be utilized in Appendix C to construct a reasonable set of phase-point operators for given density matrices.

C. SU(3)DTWA
Here we formulate the DTWA for the SU(3) phasespace variables. Throughout this paper, we refer to this approach as the SU (3)DTWA.
Let us consider real-time dynamics of a many-body spin-1 system described by a HamiltonianĤ. The initial density matrixρ 0 =ρ(t = 0) can be expressed as an expansion in a tensor product of local phase-point operatorsρ where w α ≡ w α1,··· ,α M is a many-body discrete-Wigner function defined in the M -body phase space Γ M ≡ Γ 1 ⊗ · · · ⊗ Γ M . Note that M typically represents a total number of sites for lattice systems. Each local operator A αj acts on the site j. Such an expansion is expected to exist for any states because a set ofÂ αj forms a local operator basis. An operatorΩ that we are interested in has also an expansion given bŷ Then, the expectation value ofΩ at time t > 0, i.e., The propagation function U W (β, α; t) connecting two Weyl symbols Ω β and w α is defined by whereÂ α = M j=1Â αj andÛ (t) = e − i Ĥ t is the unitary time-evolution operator. This propagator contains complete information of quantum many-body dynamics governed byĤ. However, the unitary transformation given byÛ (t)Â αÛ † (t) changes the tensor product into complicated operator strings in the Hilbert space, so that the exact evaluation of U W (β, α; t) is generally impossible.
The TWA for quantum dynamics is nothing else but an appropriate semiclassical approximation for the phasespace propagator U W (β, α; t) [33]. In the treatment discussed in Ref. [19], one makes the following directproduct ansatz for the many-body phase-point operators at time t > 0: The time dependence of x The classical Hamiltonian H W can be derived by replac-ingX where x Inserting Eq. (44), we finally arrive at the SU(3)DTWA representation of Ω (t) : If we putΩ =X (j = k) and perform the summation over β ∈ Γ M , we have the formulas In typical cases, initial density matrices are factorized with respect to the single-body index j. Then, the discrete-Wigner function reads as Therefore, we obtain As learned from these expressions, the only difference of the SU(3)DTWA from the standard SU(3)TWA comes from their probability distributions for the phase-space variables. In other words, the classical dynamics in the SU(3)DTWA still happen in the continuous phase space. Compared to the Gaussian approximation, the DTWA method features a numerical advantage that it allows to sample spin configurations with positive probabilities for typical product states, which give rise to negative probabilities in the exact continuous representation [19,20]. In the literature such as Ref. [19], examples are presented, demonstrating that the DTWA improves revival properties of the quantum dynamics, which the Gaussian approximation fails to capture. The direct comparison between the two methods will be presented in Sec. V B.
We mention that our description, which explicitly uses the phase-point operators and, therefore, explicitly defines a discrete-Wigner function for a density matrix, is distinct from a similar discrete-sampling approach for general SU(N ) systems developed in Ref. [25]. The latter approach has not introduced any phase-point operators explicitly, but instead has utilized a quantumtomography-like methodology to define probability distributions for each phase-space variable. This state-ofthe-art sampling technique, which is also called the generalized DTWA (GDTWA) [25], has been already applied to actual experimental setups of large-spin systems such as 52 Cr gases [34] and 167 Er gases [35], and the performance has been evaluated against the experimental data. In Sec. V, we compare this sampling scheme to our schemes, specifically for the 2D Bose-Hubbard model with a small size.
To implement the tomography technique for the SU(3) TWA, we decompose each SU(3) matrix T µ in its diagonalized basis, i.e., denote the eigenvectors of T µ associated with the eigenvalues λ (s) µ . Note that, generally speaking, the matrices T µ cannot be simultaneously diagonalized. We compute an expectation value of T µ with a density matrix ρ to obtain Following Ref. [25], the coefficients p (s) µ are regarded as the probabilities for the discrete spins d µ ∈ {λ where s λ We determine the probabilities p Therefore, in the TWA simulations, d 1 , d 2 , d 6 , and d 7 randomly choose either 1 or −1 with an equal probability, while d 3 = d 4 = d 5 = 0 and d 8 = 2/ √ 3 for all samples. Note that the fluctuations of each variable are statistically independent of those of the other ones. More detailed discussions of the tomography technique are found in Ref. [25]. In Appendix D, we add a supplemental discussion on the relationship between this tomography method and our DTWA scheme, associated with the reproducibility of a second-order moment for a pure state.

D. Fully connected spin-1 model
To compare the SU(3)DTWA with the Gaussian SU(3)TWA, we study the fully connected spin-1 model (4). To be specific, we calculate sudden-quench dynamics of several physical quantities by using the SU(3)TWA with the Gaussian-Wigner function and the SU(3)DTWA, respectively, and compare these semiclassical results with the exact ones.
In Fig. 2, we numerically simulate the time evolution of the fully connected spin-1 model of Eq. (4) after sudden quenches from the x-polarized direct-product state  Fig. 2(b)]. In the lower panels in Fig. 2, we also show the time evolution of M −1 j (Ŝ z j ) 2 (t) starting from the same initial state. It should be noticed that the slight recurrence of the oscillation observed in Fig. 2(d) at late times after t ≈ 60 /U are not captured within the semiclassical approximation as expected in typical TWA simulations [13,16].
In Fig. 2, we also simulate the same dynamics by using the SU(3)DTWA approach. For all the panels, the SU(3)DTWA results (red dotted lines) reasonably reproduce the same dynamics as those of the Gaussian SU(3)TWA. As explained in Appendix C, for the DTWA results in Fig. 2, we have prepared a statistical mixture of random initial conditions characterized by multiple sets of phase-point operators. A similar technique has been used in Ref. [32]. We emphasize that if we only use the Wootters representation for samplings, it will fail to correctly produce the dynamics [see also Fig. 7(a)].
In Fig. 3, we compute the expectation value M −1 j (Ŝ z j ) 2 (t) for another initial state  discrete (red dotted) SU(3)TWA results reproduce the first and second peaks of the exact expectation value (blue solid) within t < 30 /U . As U/J decreases, the timescale, during which the exact quantum dynamics are reasonably captured by the semiclassical expressions, is shortened. This tendency can be attributed to the nonlinearity of the system that gives rise to a significant error in the exact time evolution of the many-body Wigner function. In Fig. 3(b) corresponding to U = 125J, both semiclassical approaches only recover the first peak within t < 10 /U , however, they fail to describe the second peak, especially, its amplitude.
After leaving from the early-time stage, the SU(3)TWA clearly deviates from the exact dynamics. In particular, it is clearly seen in Fig. 3 that the Gaussian SU(3)TWA tends to saturate into a steady value but not to make a recurrence of the oscillation, both for U = 250J and 125J. Interestingly, especially in Fig. 3(b), while the SU(3)DTWA also fails to describe the exact dynamics for t > 10 /U , but it exhibits an oscillatory behavior rather than saturation. However, it should be emphasized that the discrete Monte Carlo sampling does not affect the quantitative timescale itself, during which the quantum dynamics are almost accurately captured within the semiclassical expressions. This seems to be reasonable because the classical equations of motion for the continuous and discrete cases are the same.
To close this section, we have demonstrated that the SU(3)DTWA is nearly as accurate as the Gaussian approximation with respect to simulating the quench dynamics. In the next section, we apply these techniques to analyses of the experimental results for 2D Bose-Hubbard systems [15].

V. APPLICATION TO THE 2D BOSE-HUBBARD SYSTEM
In this section, we apply the SU(3)TWA approaches for studying far-from-equilibrium dynamics of the Bose-Hubbard model on a square lattice at unit filling.  We specifically analyze dynamics of equal-time singleparticle correlation functions after a quench from a Mottinsulating state to a parameter region near the quantum critical point [15]. Theoretical studies on dynamics of equal-time correlation functions have been reported in Refs. [11,[36][37][38][39][40][41][42][43][44][45].

A. Experimental setup
First we briefly summarize the details of the experimental setup in Ref. [15]. Takasu and his coworkers have measured sudden-quench dynamics of the single-particle correlation functions inside the 2D Mott-insulator phase in the following steps: 1. They prepared a unit-filling Mott insulator of an ultracold 174 Yb gas in an optical square lattice with s = V 0 /E R = 15. The energy scales V 0 and E R denote the optical lattice depth and the recoil energy of this system, respectively. The prepared system is well described by the direct product Fock state for bosons wheren j |n j = n j |n j .
2. The lattice depth was abruptly decreased from s = 15 to s = 9. The time to ramp down the lattice depth is approximately 0.1 ms. The lattice depth after the quench implies U/J = 19.6.
3. After the quench, the resulting dynamics was observed by measuring the time-of-flight interference pattern that can be converted to the equal-time single-particle correlation functions, where r j = (x j , y j ) indicates each site on the square lattice with units of the lattice constant d lat = 266 nm. The real-space summation is performed under the conditions |x j − x j | = ∆ x and |y j −y j | = ∆ y , and we write ∆ = (∆ x , ∆ y ). Recall that M is the total number of lattice points.
In this work, as a simplified setup, we neglect harmonic trap potentials in numerical simulations. We simply assume that all the atoms participate in a uniform Mottinsulator state before the quench. Effects due to spatial inhomogeneity of the gases will be discussed in Sec. V D.
To close this subsection, here we note that the qualitative behaviors of the dynamics of the spatial correlation functions measured after quantum quenches can change depending on the initial states that we take. For instance, for the coherent state as the initial states, which describes a coherent condensation of bosons at the non-interacting limit, sudden changes of the interaction, from zero to weak interactions, result in observing fine oscillations in time of the density-density equal-time correlation function, reflecting the coherent motion of the Bogoliubov quasiparticles [11]. By contrast, if we choose the Mottinsulator states as the initial conditions, and propagate the states with the Hamiltonian with the same interactions (i.e., quenches from infinite to weak interactions), we observe propagation of a peak signal without fine oscillations in the same correlation function [11]. Its propagation velocity is well explained by the single-particle excitation spectrum of the Hartree-Fock approximation. Reliable TWA results on this kind of initial-state dependence of the quench dynamics can be found in our previous study for the 2D Bose-Hubbard model with a large filling factor [11].

B. Small size case
Before proceeding to our main results corresponding to the experimental setup, let us consider the quench dynamics for a small-size 2D Bose-Hubbard system, say 9 sites, in order to compare outputs of the SU(3)TWA approaches with those of the exact numerical calculation. For simplicity, we focus on the sudden-quench limit, in which the ramp-down time is neglected. Figure 4 shows the numerical results for the suddenquench dynamics of K ∆ (t) for a small-size Bose-Hubbard system. The simulation setup has M = 3 2 = 9 sites and we adopt periodic boundary conditions [46]. In Fig. 4, the full-quantum dynamics of the Bose-Hubbard system is evaluated by integrating the time-dependent Schrödinger equation of the Hamiltonian (1) (gray dotted line). The maximum occupation of the local site is n max = 2, hence, the three lowest states, i.e., |0 , |1 , |2 , are allowed in this simulation. We observe that the correlation functions at ∆ = (1, 0) and (1, 1) form a first-peak region in the time range of 0 < tJ/ < 0.5. At later times, tJ/ > 0.5, the time evolution of correlations exhibits an almost undamped oscillation reflecting its small size.
In Fig. 4, we also simulate the same dynamics by using the SU(3)TWA for the effective-model Hamiltonian (2) according to the Gaussian and discrete-Wigner approaches of Monte Carlo samplings. The unit-filling Mott-insulator state is given by Eq. (53), i.e., |Ψ ini ≈ |Ψ Mott . We observe that both the Gaussian SU(3)TWA (red circle) and SU(3)DTWA (blue triangle) quantitatively capture the first-peak region in the range of 0 < tJ/ < 0.5, especially its initial growth, its time point of the center of the region, and its correlation intensity. However, the later-time dynamics for tJ/ > 0.5 can not be well captured within the SU(3) semiclassical representation. Indeed, the semiclassical results exhibit almost saturated behaviors rather than the temporal oscillation with a large amplitude. It should be emphasized that the difference between two semiclassical results in the later-time dynamics comes from our choice of the initial distribution for the phase-space variables. Interestingly, it is clearly seen that, around t = 0.5 /J in Fig. 4, the SU(3)DTWA gives a slightly better result, i.e., shows deeper dips of correlations. For this comparison, the SU(3)DTWA can be seen as a better description than the Gaussian SU(3)TWA. Figure 4 also displays the simulation result on the basis of the tomography technique as presented in Sec. IV C. We numerically find that it is closer to the Gaussian result, rather than the DTWA one. This coincidence to the Gaussian simulation indicates that the tomography technique also provides a reasonable sampling scheme for the initial condition. Since there is no considerable deviation from the Gaussian result, in the following discussions, we do not use the tomography technique.
It is interesting and helpful to calculate the quench dynamics by using the GPTWA for the strongly interacting Bose-Hubbard system as a reference. In order to carry out an efficient simulation, we have used an approximate Gaussian distribution representing the Fock states [11]. The details of the GPTWA will be briefly reviewed in Appendix B. In Fig. 4, the GPTWA simulation (green square) fails to describe the correlation intensity in the first-peak region while it reproduces well a very early growth of the correlation function at ∆ = (1, 0) within tJ/ < 0.2. Therefore, for the purpose of simulating the strongly interacting dynamics, the SU(3)TWA certainly provides a better description than the GPTWA.

C. Comparison to the experimental results
We calculate the quench dynamics for a larger-size system corresponding to the experimental setup. Figure 5 shows the correlation function K ∆ (t) for M = 20 2 = 400 with periodic boundary conditions. First, we prepare the system in the unit-filling Mott-insulator state (t < 0), and then abruptly decrease the lattice depth until t = 0. For t > 0, the system evolves in time at U = 19.6J. While the dynamics of the effective pseudospin-1 model (2) is computed in the SU(3)TWA simulations, that of the Bose-Hubbard model (1) with no truncation of the local Hilbert space is computed in the GPTWA. We note that the GPTWA result in Fig. 5 is a reproduction from Ref. [15].
In Fig. 5(a), we observe that all the semiclassical results explain well the growth of the nearest-neighbor correlation at ∆ = (1, 0) in the early-time stage within t < 0.1 /J. In addition, these reasonably describe the correlation offset at t = 0. The experimental data show a peak in the time domain of 0 < tJ/ < 0.2. At longer times, the measured correlation gradually satu-rates to a steady value. In the comparison performed in Fig. 5(a), the experimental result is seemingly closer to the GPTWA rather than the SU(3)TWA. In particular, the peak position and the correlation intensity in the time window indicated by Fig. 5(a) are relatively close to the ones simulated by the GPTWA. This is in contrast to the small-size case in Sec. V B, where the SU(3)TWA is closer to the exact dynamics and can provide a reasonable first-peak region at short times. Notice that the correlation intensity of the experiment is typically lesser than both SU(3)TWA and GPTWA results.
Next, we focus on longer distances, say ∆ = (1, 1) [ Fig. 5(b)] and ∆ = (2, 0) [ Fig. 5(c)]. The experimental data are seen to achieve a peak during 0.1 < tJ/ < 0. , it is hard to locate the center of the first-peak region in the experimental data because of significant noises. Within the error bars, we expect that there exists a peak region somewhere in the range of t < 0.5 /J. For these long distances, correlation intensities in the GPTWA are seen to be suppressed because it cannot capture strong quantum fluctuations in the parameter regime. In particular, no clear peak region is observed in the simulation even at short times. Therefore, its agreement to the experiment is worse. By contrast, the SU(3)TWA, which is expected to describe local quantum fluctuations in the regime more accurately, can produce a reasonably strong correlation, which is comparable to the experiment, and clear peak regions in the range of t < 0.5 /J. Hence, in this case, these SU (3) simulations are closer to the experiment.
Finally, let us mention that, in the experiment, not only the nearest-neighbor correlation but also the longerdistance ones exhibit a finite and non-negligible offset at t = 0. However, according to the SU(3)TWA and the GPTWA, such an offset for longer distances should be almost zero. We will discuss this point in details in the next section.

D. Discussions
In the direct comparisons for the large system, we observed that the experimental results of the spatial correlation function are closer to the GPTWA especially at short distances while the SU(3)TWA looks better at long distances. As learned from the numerical simulations for the small size, the SU(3)TWA should work better more than the GPTWA in the strongly interacting parameter regime. Moreover, we also recognized that the nonzero offsets of the correlations at distances except for nearest neighbors are not consistent to all the semiclassical results. We argue that the above unexpected observations could be attributed to some contributions present in the actual experiment, which are not precisely taken into account in our SU(3)TWA simulations.
First, we discuss occupations of bosons allowed in the SU(3)TWA for the Bose-Hubbard systems. The formalism of SU(3)TWA for bosons is constructed under assumptions that the local Hilbert space is truncated up to three states. If the dimension of the reduced state space is extended from three to more, it will improve the simulated result, more or less, quantitatively. To perform this extension, one needs to increase the local phase space furthermore. For instance, if five states are relevant locally, SU(5) matrices should be chosen as a phase-space variable. However, we may expect that higher occupations give no significant effect, at least in our current case, in which the strength of the interaction is large enough to suppress them. In order to justify this expectation, in Appendix E, we will clarify the degree to which occupations greater than 2 affect the quench dynamics of the correlation function in the parameter regime of the experiment by utilizing an exact numerical calculation for a small size.
Second, we make a comment on effects of an inhomogeneous trap potential. In the experimental setup in Ref. [15], the prepared initial state actually contains a strongly correlated superfluid component with incommensurate fillings due to a harmonic trap while the region of the Mott insulator with unit filling is much larger. Such a contribution is not dominant over the whole gas, but not completely negligible. In the Supplemental Material of Ref. [15], an MPS calculation has been performed for a 1D trapped Bose gas in the presence of narrow superfluid regions in the system. A numerical result shows that a finite offset appears at the end point of quenches at several distances in addition to the nearest neighbor. This strongly indicates that the presence of superfluid contributions, more or less, affects the time evolution of the correlation function. To our current techniques for the SU(3)TWA, it is difficult to initialize a system into such inhomogeneous states as prepared in the MPS simulation. In future works we will develop an efficient technique to treat this kind of initialization problem.

VI. CONCLUSIONS AND OUTLOOKS
In conclusion, we have analyzed far-from-equilibrium dynamics of strongly interacting Bose gases in an optical lattice by using the SU(3)TWA on the basis of different Monte Carlo sampling schemes. In the middle of this paper (Sec. IV), the SU(3)DTWA approach has been developed as a sampling scheme, and applied to the fully connected spin-1 model with a large size in order to examine  this approach. We demonstrated that the SU(3)DTWA is nearly as accurate as the Gaussian SU(3)TWA in simulating time evolution after sudden quantum quenches.
In the main part of this paper (Sec. V), we have applied the SU(3)TWA to sudden quench dynamics of a strongly interacting Bose gas in the 2D optical lattice. The semiclassical methods on the basis of the GPTWA, the SU(3)DTWA, and the Gaussian SU(3)TWA have been compared with exact numerical calculations for the 2D Bose-Hubbard model with a small size. We recognized that the SU(3)DTWA and the Gaussian SU(3)TWA can provide better descriptions than the GPTWA in a strongly interacting regime. The numerical results on the basis of those semiclassical methods have also been compared with the recent experiment at Kyoto University. We found that at short distances, the experiment is closer to the GPTWA while, at relatively-long distances, it is reasonably close to the SU(3)DTWA and the Gaussian SU(3)TWA. We argued that this observation can be attributed to parts of the actual experimental realization including an inhomogeneous trap potential, which are not precisely taken into account in our numerical simulations.
Beyond the scope of this work, it would be interesting to develop a cluster TWA approach [47] for the strongly interacting Bose-Hubbard systems. For applications of this strategy in higher dimensions than 1D, a reasonable reduction scheme of dimensions of cluster phase-space variables may be required to make simulations realistic and efficient. Subject to a fixed M , an arbitrary state of this system is spanned by a Fock vector labeled by two non-zero integers ν 1 ≥ 0 and ν 2 ≥ 0, where 0 ≤ ν 1 +ν 2 ≤ M . This basis state is a simultaneous eigenstate forΠ 3 andΠ 8 , therefore, Π 3 |ν 1 , ν 2 = (ν 1 − ν 2 )|ν 1 , ν 2 , The rest of the operators, e.g.,Π 1 , behave as a ladder operator connecting different Fock stateŝ One can evaluate the matrix element of the Hamiltonian between |ν 1 , ν 2 and |ν 1 , ν 2 , i.e., ν 1 , ν 2 |Ĥ fc |ν 1 , ν 2 . Its dimension algebraically increases with M , so that one can implement the exact numerical analysis on computers even at large M . In Fig. 6, we compute the time evolution of the expectation value j (Ŝ z j ) 2 (t) = 2M/3− Π 8 (t) / √ 3 for M = 8 using two different approaches: the dashed line is a numerical integration of the time-dependent Schrödinger equation for the Hamiltonian matrix expressed in terms of the collective spinΠ α , whereas the solid line is for the spin Hamiltonian in terms ofŜ α j . The initial state for this simulation is the zero-magnetization direct-product state |Ψ(t = 0) = j |S z j = 0 . The perfect agreement of two results manifests that the collective-spin expression gives a more efficient way to have the same result than the straightforward approach.

Appendix B: Gross-Pitaevskii truncated-Wigner approximation for the Bose-Hubbard Hamiltonian
For the TWA in the coherent-state phase space, the classical time evolution of the Bose-Hubbard Hamiltonian is governed by the discrete GP equation associated with the Heisenberg-Weyl group: The classical function H W (α, α * ) = (Ĥ BH ) W is the Weyl symbol ofĤ BH given by If we write α cl (t) as a solution of the GP equation with conditions α cl (t = 0) = α 0 , the expectation value of an operatorΩ, i.e., Ω (t) = Tr[Ωρ(t)] = Tr[Ω(t)ρ(t = 0)] is reduced to the following phase-space integration form (for details, see [11][12][13]): Here dαdα * = π −M M j=1 dRe[α j ]dIm[α j ] is the measure of the phase-space integration. The weight function over the phase space is the Wigner function defined by means of the coherent state basis We note that the GPTWA typically provides quantitative descriptions of real-time dynamics of the Bose-Hubbard systems when they have a sufficiently small interaction or sufficiently large filling factor. In recent years, this type of semiclassical method has been applied to multiple dynamical problems of lattice bosons, e.g., see Refs. [11,15,[49][50][51] for details.

Appendix C: Details of the SU(3)DTWA simulation
First let us present the numerical sampling when the x-polarized state is chosen as our initial state. Generally speaking, the discrete-Wigner function representing such a superposed state exhibits negativity. To carry out an efficient numerical simulation, we take the following steps: As the first step, we prepare a polarized down-spin state along z-axis at t = −π ≡ t 0 If we use the Wootters representation for the phase-point operator, the corresponding discrete-Wigner function is positive. Therefore, it is easy to sample randomized spins from the distribution.
Then, we shine a global pulse such that it evolves |Ψ 0 into the desired target state, i.e., the polarized state in the x-axis |Ψ 0 = j |S x j = 1 . Such a spin-flip process is designed via a unitary time evolution described byÛ p (t) = e − i Ĥ p(t−t0 ) with a Hamiltonian If the unitary operation of the pulse is applied from t = t 0 to t = 0, each spin state is locally flipped such that |S z j = −1 → −|S x j = 1 . The minus sign at the final state gives no effect on the expectation value. Notice that the time evolution governed byÛ p is exactly simulated by the SU(3)TWA becauseĤ p is linear in the phasespace variables. At t = 0, the prepared random values of the spin configurations are expected to obey the Wigner distribution of |Ψ 0 = j |S x j = 1 .  Fig. 7(a) is calculated by using the Wootters representation A (0) α and following the above procedure. In what follows, we write S 0 as a statistical ensemble of the discretized phasespace variables sampled from the discrete-Wigner function for A (0) α . It is clear that the SU(3)DTWA with S 0 fails to reproduce the exact dynamics even though the Gaussian approach can do so. This consequence seems to be related to the fact that realizable configurations in S 0 are quite restricted compared with those belonging to the Gaussian distribution.
To resolve this problem, we utilize a prescription in which we define a few other sets of the phase-point operator and make a statistical mixture of them as done in Ref. [32]. As a non-trivial example, we can construct the following two sets of phase-point operators instead of and result in X The discrete sampling scheme based on the ensemble S 1 (S 2 ), however, fails to reproduce the moment in the phase-space representation. In fact, the phase-space average X 2 1 → x 2 1 produces 0 (1) as checked via direct computations. Hence, there is an underestimation (overestimation) of the quantum correlation in the classical ensemble generated by the naive phase-point-operator method. Note that, if in the beginning the squared op-eratorX 2 1 is linearized in the SU(3) matrices, and after that it is transformed to the phase-space quantities, the point-operator method accurately reproduces the moment. The statistical mixture S 1 ∪ S 2 that we have made in the previous appendix adequately averages the fluctuations of the classical variable belonging to each ensemble, and, as the consequence, produces the exact value of the moment as the phase-space average [namely, in this case, (0 + 1)/2 = 1/2]. This observation catches an underlying reason of the success of the DTWA simulation for the state |S x = 1 , which is prepared after the unitary evolution of |S z = −1 (see also Appendix C). We expect that the tomography scheme will also give the adequate sampling for the simulation, but it is not explicitly implemented in this paper. Thorough analyses about the connection between our DTWA scheme and the tomography method will be addressed elsewhere, which are beyond the central purpose of this work. To visualize how the three-state truncation works in the parameter regime of the experiment, we numerically integrated the time-dependent Schrödinger equation for the 2D Bose-Hubbard Hamiltonian with M = 2 2 = 4 and some values of n max . Recall that n max means the maximum occupation of each site. The initial state of the following simulation is the unit-filling and homogeneous Mott-insulator state.
In Fig. 8, we show exact numerical results for the quench dynamics of the single-particle correlation function â † jâ j . The interaction during the time evolution is set to U/J = 20, which is close to the actual value of the experiment, i.e., U/J = 19.6. The results for n max = 3 (green solid line) and n max = 4 (red dotted line) agree with each other, indicating that 4particle occupations are completely suppressed at least until t = 50 /U = 2.5 /J. Although the result for n max = 2 (blue dashed line), which corresponds to the assumptions of the SU(3)TWA, fails to perfectly reproduce the result for n max = 3, it captures very well the short-time evolution of the peak region of the correlations within t ≤ /J. Indeed, the peak region at early times agrees well with the one for n max = 3 and the intensity of the correlation is close to the exact one. As the system evolves in time, the deviation between the results for n max = 2 and 3 gradually gets significant.