{'fraction_non_alphanumeric': 0.0388200957751268, 'fraction_numerical': 0.01079595364257176, 'mean_word_length': 4.605463786531131, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Time-dependent many-variable variational Monte Carlo method for nonequilibrium strongly correlated electron systems

We develop a time-dependent variational Monte Carlo (t-VMC) method for quantum dynamics of strongly correlated electrons. The t-VMC method has been recently applied to bosonic systems and quantum spin systems. Here, we propose a time-dependent trial wave function with many variational parameters, which is suitable for nonequilibrium strongly correlated electron systems. As the trial state, we adopt the generalized pair-product wave function with correlation factors and quantum-number projections. This trial wave function has been proven to accurately describe ground states of strongly correlated electron systems. To show the accuracy and efficiency of our trial wave function in nonequilibrium states as well, we present our benchmark results for relaxation dynamics during and after interaction quench protocols of fermionic Hubbard models. We find that our trial wave function well reproduces the exact results for the time evolution of physical quantities such as energy, momentum distribution, spin structure factor, and superconducting correlations. These results show that the t-VMC with our trial wave function offers an efficient and accurate way to study challenging problems of nonequilibrium dynamics in strongly correlated electron systems.


I. INTRODUCTION
Quantum systems with strong many-body correlations in equilibrium show intriguing properties such as the metal-insulator transition 1 and high-temperature superconductivity 2,3 . Recently, because of potential routes to realizing intriguing phenomena that are not attainable in the equilibrium, strongly correlated electron systems driven out of equilibrium have attracted much attention. In fact, owing to the development of experimental techniques, we have been able to control or realize unprecedented phases and their phase transitions by applying strong and short pulse of external fields such as intensive laser pumping [4][5][6][7][8][9][10][11][12][13] .
For satisfactory theoretical understanding of quantum dynamics of many-body systems, we have to solve manybody time-dependent Schrödinger equation i d dt |ψ(t) = H(t) |ψ(t) , where |ψ(t) and H(t) represent the wave function and the Hamiltonian at time t, respectively. The formal solution of the time-dependent Schrödinger equation is given by |ψ(t) = T exp(−i t 0 H(s)ds) |ψ(0) . Here, T represents the time ordering. However, such an approach is tractable only for small many-body systems, because the Hilbert space grows exponentially as the system size increases. Furthermore, reduction to an effective single-particle problems such as the time-dependent Hartree-Fock method 14 does not give us accurate results. To treat larger systems accurately, there exist several numerical methods such as time-dependent density matrix renormalization group method(DMRG) [15][16][17] , and nonequilibrium dynamical mean field theory(DMFT) 18 and quantum Monte Carlo(QMC) method 19 . However, DMRG and DMFT have difficulties in treating large systems in two or three spatial dimensions when one wishes to treat spatial correlations and fluctuations accurately. In order to include non-local correlations, the dynamical cluster approximation(DCA) has been proposed as an extension of the DMFT 20 . However, it requires high computational costs for a large cluster size. Although the QMC method can treat finite systems exactly, its applications are very limited due to the notorious negative-sign problem. Recently, to overcome the above difficulties, Carleo et al. developed the time-dependent variational Monte Carlo (t-VMC) method and optimized Jastrow factors in bosonic systems 21,22 . This method can be formulated based on the time dependent variational principle (TDVP) 14,23,24 . The t-VMC method has been applied not only to bosonic systems 21,22 but also to spin models 25 . However, to the best of our knowledge, its application to correlated electron systems has not been successful yet. This limitation may be ascribed to the difficulty of constructing an accurate trial wave function for such systems. Therefore, the proposal of accurate trial wave functions suitable for nonequilibrium electron systems in the t-VMC method is desirable.
The purpose of this study is to propose an accurate and efficient trial wave function for strongly correlated electron systems out of equilibrium in the t-VMC framework. To achieve this purpose, we focus on highly accurate trial wave functions for ground states of correlated electron systems. In such systems, many studies have attempted to construct an accurate and efficient trial wave function in the variational Monte Carlo (VMC) framework [26][27][28][29][30][31] . In Ref. [31], Tahara and one of the authors have reduced the biases by using the quantumnumber projections and introducing many variational parameters to one-body part. They have adopted a generalized pair-product wave function as a one-body part because it can flexibly describe different competing phases such as correlated metals, antiferromagnetic states and superconducting states. This improved trial wave function has proven to be highly accurate for ground states of strongly correlated electron systems [31][32][33] . In this paper, we show that this trial wave function is an accurate and efficient one even for nonequilibrium strongly correlated electron systems in the t-VMC framework.
The organization of this paper is as follows. In Sec. I I, we introduce the TDVP which enables us to obtain an optimal time-dependent trial wave function and formulate the t-VMC method by using the TDVP. Section II I describes a trial wave function with a large number of variational parameters for nonequilibrium strongly correlated electron systems. In Sec. IV, we show the accuracy and efficiency of our trial wave function by presenting several benchmark results. Finally, we summarize our work in Sec. V.

II. TIME-DEPENDENT VARIATIONAL PRINCIPLE
The time-dependent variational principle (TDVP) proposed by McLachlan is a variational principle for timedependent wave functions 23 . In this principle, we consider a distance between i d dt |ψ α and H |ψ α where α = {α k |k = 1, · · · , N p } represent time-dependent variational parameters. By definition, the distance satisfies the inequality where the equality holds if |ψ α is the solution of the time-dependent Schrödinger equation. Here, the norm |Ψ is defined as the square root of an inner product of a wave function |Ψ , i.e. |Ψ = Ψ|Ψ . If we could optimize the variational parameters at each time step such that the equality holds, we obtain the exact solution of |ψ α . If a trial wave function well approximates the exact solution of the time-dependent Schrödinger equation, the value of the lower bound should be small. Based on this idea, we optimize variational parameters at each time-step such that the distance is minimized. Originally, the TDVP was applied in the field of quantum chemistry 34,35 . Recently, a similar principle has been applied to the matrix product state for quantum spin models 24,36 , the bosonic Jastrow-type wave function for the Bose-Hubbard model 21,22 and the Gutzwiller approximation for strongly correlated electron systems [37][38][39] .
Although exact time evolution is unitary, and thus, the norm ψ α |ψ α is conserved, it is not necessary conserved in TDVP [Eq. (1)]. To remove the restriction on the norm, we use a TDVP for norm-independent dynamics 24 , The details of the TDVP for norm-independent dynamics is described in Appendix A. Based on this TDVP, we can derive the differential equation of the time-dependent variational parameters. Namely, by solving the minimization problem on the distance (2), we obtain the time evolution of the variational parameters 21,22,24,36 : where a matrix S and a vector g are described as respectively. In the t-VMC method, we estimate an expectation value A = ψα|A|ψα ψα|ψα by the Markov-chain Monte Carlo method. The derivative operators O k and O † k are defined by using real space configurations of electrons {x} as respectively. Here, The differential equation [Eq. (3)] is called TDVP equation 24,36 . This TDVP equation can also be derived by minimizing the time-dependent action 24,40 . If we use the time-dependent variational principle for imaginary time evolution t = −iτ and solve TDVP equation by using Euler method, we obtain the stochastic reconfiguration scheme proposed by Sollera 26 . The TDVP equation has a symplectic property 24,40,41 . This property leads to the energy conservation if the Hamiltonian is timeindependent and we could calculate the derivative of parametersα k exactly.
In this study, in order to solve the TDVP equation, we use the fourth-order Runge-Kutta method which provides us with a stable and efficient way to perform the time integration. Note that the Runge-Kutta method is not a symplectic integral method. Furthermore, there are stochastic errors in the Monte Carlo calculation of quantities such as O † k H and O † k O l . These cause the breaking of the symplectic property of the TDVP equation. Nevertheless, we observed that the energy is conserved with high accuracy as we show in Sec. IV.

III. TIME-DEPENDENT VARIATIONAL WAVE FUNCTION
In the t-VMC method, the choice of trial wave functions is important. As a trial wave function, we adopt the form of which has been used for equilibrium systems in Refs. [31][32][33]. Here, L represents quantum-number projections which recovers the symmetries the wave function should have throughout the time evolution, and P(t) represents correlation factors. For the one-body part |φ(t) , we employ the pair-product wave function. In addition, we include backflow correlations in the pair-product wave function for lattice model 29,42 . The time-dependent variational parameters are included in the correlation factors as well as in the one-body part with the backflow correlations. Note that these variational parameters should be treated as complex numbers because the variational parameters evolve as complex numbers in the present method. In this section, we describe each component in detail.

A. One-body Part
In the conventional VMC, the Slater determinant with small variational parameters is used as the one-body part. In order to improve the conventional one-body part, we assume the form of the pair-product wave function with many variational parameters 31 : where N is the number of electrons, N s is the system size, c † iσ (c iσ ) is a creation (annihilation) operator of an electron with spin σ on the site i, and the pairing amplitude f ij is treated as variational parameters. This pairproduct wave function is a general form of a Hartree-Fock-Bogoliubov-type wave function which allows antiferromagnetic and superconducting orders 31,43,44 . Thus, it takes an advantage of flexibly describing the paramagnetic metals, antiferromagnetic ordered states and superconducting states with any type of frequency independent gap on equal footing.

B. Correlation Factors
By operating the correlation factors on the pairproduct wave function, we can include many-body correlation effects beyond the mean-field level. In this study, we use the Gutzwiller factor P G 45 and the Jastrow factor P J 46 , i.e., P = P G P J . The Gutzwiller factor which was introduced by Gutzwiller 45 has the form of where n iσ = c † iσ c iσ and g is the variational parameter. The Gutzwiller factor punishes the double occupation of electrons on the same site in real space configurations.
In the limit g → ∞, the Gutzwiller wave function P G |φ corresponds to a state which contains no double occupation. The Gutzwiller factor is a simple way to improve a mean-field wave function such as a Slater determinant and the pair-product wave function. However, it was numerically proven that the Gutzwiller wave function cannot describe the nonmagnetic Mott transition in any finite dimensional systems 47 . The main reason for this is that the Gutzwiller factor only includes the on-site correlation. Although some doubly-occupied (doublon) and empty (holon) sites exit in the Mott insulator where charge fluctuations are allowed, the doublon and holon have to be bound in realizing an insulating behavior.
In order to describe the Mott transition, the Jastrow factor is introduced 27 : where are the variational parameters, and r i (r j ) represents the position vector of the site i(j). The Jastrow factor can be represented by using a doublon number operator and a holon number operator. The electron number operator n i is written as is the holon number operator. By using this relation, Jastrow factor becomes Equation (12) shows that the Jastrow factor includes off-site repulsive correlations between doublon-doublon (holon-holon) pairs and attractive correlations between doublon-holon pairs if v ij > 0. This doublon-holon attractive correlations in the Jastrow factor play a crucial role for the Mott transition 27,48 .

C. Quantum-Number Projections
In general ground states of finite quantum systems, the symmetries of the Hamiltonian must be preserved even if the symmetry-breaking occurs in the thermodynamic limit. Furthermore, such symmetries must be preserved even after time evolutions as long as external fields do not break the original symmetries. However, conventional trial wave functions often break symmetries of the Hamiltonian.
The quantum-number projection enables us to recover the symmetries of the trial wave function 49 . Here, we introduce two quantum-number projections: the spin projection L S and the momentum projection L K . The spin projection L S is the projection onto the state with a total spin S and z-component of spin S z = 0. This projection has a form of the integration over spin space: where Ω = (α, β, γ) is the Euler angle, P S (cos β) is the Sth Legendre polynomial, and R(Ω) = e iαS z e iβS y e iγS z is the rotational operator. When the one-body part |φ and the real space configurations {x} satisfy the condition of S z = 0, we can omit the integrations over γ and α as follows: We can efficiently estimate the integration over β by the Gauss-Legendre quadrature 50 . The momentum projection L K is the projection onto the state with a total momentum K. This projection has a form: where T R is the translation operator with the translation vector R.

D. Backflow Correlations
One efficient way to include correlation effects in a trial wave function is to include backflow correlations 29,42,[51][52][53] . Recently, Tocchio et al. proposed a way of introducing backflow correlations into a Slater determinant for lattice models and found that the backflow correlations substantially improve ground-state energy of frustrated electronic systems in the region above intermediate strength of coupling 29,42 . In a way similar to the approach by Tocchio et al., the backflow correlations for lattice models can be implemented in the pairproduct wave functions with the momentum projection as where |φ b and PfX b (x −R ) represent the pair-product wave function with backflow correlations and the Pfaffian of the skew-symmetric matrix Here, x R represents the real space configuration which is created by shifting a configuration x by a translation vector R. The site of the n-th (m-th) electron is represented by i n (i m ). T R (i n ) is the site characterized by the position vector r in + R. For simplicity, we do not consider the momentum projection in the following equations. The pairing amplitude with backflow correlations is defined by where {η µν τ τ ′ } represent variational parameters. τ (τ ′ ) is taken over τ (τ ′ ) that satisfies the following condition : We usually choose r max as the range of hopping terms in the Hamiltonian. Here, we drop the electron indices n and m for simplicity of notation, and we define Θ µσ i,i+τ (x) as (22) where δ i,i+τ represents the Kronecker delta, h iσ = 1−n iσ and A x = x|A|x / x|x . We impose that η µν τ τ ′ has the inversion symmetry and is independent of spin indices, namely η µν is replaced with 1 when Θ 1σ i,i+τ (x) = 0 for any τ and σ. Note that the introduction of backflow correlations make computational costs heavy because we need to recalculate the element of pairing wave function whenever we generate a candidate of the next sample and calculate expectation values of off-diagonal operators.
Finally, we explain the difficulties in operating the spin projection on the pair-product wave function with backflow correlations. From the definition of the pairing amplitudes with backflow correlations, the trial wave function with the spin projection and backflow correlations is described as Here, we need to take the summation over real space configurations x ′ because the skey-symmetric matrix depends on x ′ . However, we cannot usually take this summation because the number of x ′ grows exponentially as the system size increases. Thus, it is difficult to operate the spin projection on a trial wave function with backflow correlations.

IV. RESULTS
In this section, we show the accuracy and efficiency of the t-VMC method. For benchmark tests, we con-sider the quench dynamics in fermionic Hubbard models. We compare the t-VMC results with numerically exact results obtained by calculating the formal solution of the time-dependent Schrödinger equation and timedependent DMRG.

A. Model and Setting
For benchmark tests, we consider the fermionic Hubbard model which is defined as where t hop represents the hopping amplitude between nearest-neighbor sites and U (t) represents the timedependent onsite interaction. In this paper, we set t hop = 1.0. In two-dimensional cases, we treat the square lattice. We consider the linear-ramp quench where the strength of interaction U (t) is linearly changed from U i to U f during time t q : where U i and U f represent the strength of interaction before and after the quench protocol, respectively. For equilibrium systems, there exist many theoretical studies on the Hubbard models 33,54-57 . In the present models, the ground state at half-filling is believed to be the antiferromagnetic insulator for any nonzero U/t hop due to the Fermi surface nesting. Away from half-filling, this model shows rich phases such as correlated metals, antiferromagnetic metal and d-wave superconducting states 33,[55][56][57] . According to the VMC calculations, such superconducting states appear above U/t hop 6.0 in an interval of the doping concentration 33 .
In all the t-VMC calculations, we impose a boundary condition so that the closed shell condition is satisfied for U (t) = 0.0. We choose the discrete time step ∆t = 1.0×10 −2 /U (t) for U (t) > t hop and ∆t = 1.0×10 −2 /t hop for U (t) < t hop in the Runge-Kutta procedure.

B. Comparison with Exact Results
In this section, we mainly show the t-VMC results compared with "exact" results. Here, the "exact" results mean those obtained by directly calculating the formal solution of the time-dependent Schrödinger equation: where the wave function is completely expanded in the full Hilbert space. The errors in this calculation arise only from the time discretization with the amplitude of O(∆t M ). To preserve the total energy and the norm of the wave function with good accuracy, we have to choose a small time-grid ∆t and a high order M . In this study, we choose ∆t = 1.0 × 10 −3 and M = 3 for systems at half filling and ∆t = 1.0 × 10 −2 and M = 5 for doped systems. Although the exact results can be obtained irrespective of the interaction strength, the system size is severely restricted to small systems because of the exponentially growing computational cost with the increase of the system size. In order to check whether our trial wave function is accurate even for larger systems, we here compare the t-VMC results with time-dependent DMRG results which are expected to be highly accurate in the thermodynamic limit.
In this section, we mainly show the results obtained by using a trial wave function with all quantumnumber projections but without backflow correlations |ψ(t) = L K=0 L S=0 P(t) |φ(t) . Although we also tried a trial wave function with backflow correlations |ψ(t) = L K=0 P(t) |φ b (t) in our benchmarks, we did not observe any clear improvement at least for small systems. Figure  1 shows an example of the time evolution of the double occupancy for the one-dimensional Hubbard model at half filling and (U i , U f , t q , N s )=(0.0, 12.0, 5.0, 10). This result indicates no effect of the backflow correlations. Note that backflow correlations may improve the description of dynamics in large systems with geometrical frustrations in a strong coupling region because these correlations are important for the description of ground states for such systems 29,42 .

One-and two-dimensional Hubbard model at half-filling
In this section, we calculate the time evolution of several physical quantities. The physical quantities we have measured are the double occupancy d(t), the momentum distribution n(t; k) and the spin structure factor S s (t; q) defined as respectively. Here, k and q are wavenumbers in the Brillouin zone. In addition, S i = 1/2 σ,σ ′ c † iσ σ σ,σ ′ c iσ ′ , where σ are the Pauli matrices. Figures 2 and 3 show the time evolution of several quantities in one dimension for (U i , t q , N s )=(0.0, 5.0, 16) and in two dimensions for (U i , t q , N s )=(0.0, 5.0, 4×4), respectively. The dimension of the Hilbert space in these systems is Ns C N × Ns C N ≈ 10 8 . However, our trial wave function has only several hundreds parameters. Nevertheless, the t-VMC results well reproduce the exact results. These results show that our trial wave function offers a highly accurate and efficient description of quantum dynamics for strongly correlated electron systems.
Here, we briefly comment on a dynamical transition in the Hubbard models at half-filling. Several works showed that the dynamics of the jump ∆n(t) is different depending on the strength of interaction [58][59][60] . For weak interactions, the jump ∆n(t) decreases gradually from ∆n(t) = 1.0 to a constant. On the other hand, for strong interactions, the jump ∆n(t) exhibits a collapse-andrevival oscillation. In Fig.2(c), we observe clear collapseand-revival oscillations in one-dimensional system especially at U f = 8.0. However, Fig.3(c) shows that only weak oscillations are detected even at large U f in the twodimensional system. The difference appears to show a qualitative difference between one-and two-dimensional systems. However, it is difficult to conclude whether a dynamical transition to the collapse-and-revival oscillation happens since the system sizes are too small to measure physical properties right at the Fermi surface, especially in two-dimensional case. To study nonequilibrium properties such as the dynamical transition, we need to treat larger system size. We leave its analysis for a future study.
Next we show the dependence on time-dependent trial wave functions. In Fig. 4(a), we present the t-VMC re-   occupancy d(t), (c) jump of momentum distribution at Fermi energy ∆n(t) = n(t; (π, 0)) − n(t; (π, π/Ns)), and (d) spin structure factor Ss(t; q = (π, π)) for the linear-ramp quenches in Hubbard model on square lattice at half filling. The parameter set is chosen as (Ui, tq, Ns)=(0.0, 5.0, 4×4). Symbols and curves represent the t-VMC results and the exact results, respectively. sults for several different time-dependent trial wave functions for the quench from U i = 0.0 to U f = 4.0. We note that |φ instead of |φ(t) represents the one-body part fixed through the time evolution at |φ(t = 0) optimized in the ground state before quenching. In two of the trial wave functions, we optimized all variational parameters with the momentum projection and the results agree with the exact ones. Other results are obtained without optimizing some part of variational parameters or without operating the momentum projection. As seen in this figure, these results do not reproduce the exact results with substantial discrepancies. Especially, we find that the result obtained by optimizing only the Gutzwiller factor clearly disagrees with the exact one. The main reason for this disagreement is that Gutzwiller-type wave function cannot describe insulating states because of a lack of long-range off-site correlations which the Jastrow factor includes. In fact, the result obtained by optimizing only the Gutzwiller-Jastrow factor shows qualitative agreement with the exact one. This tendency is similar to that in equilibrium systems where a Gutzwiller-type wave function fails in reproducing the physical properties of the Hubbard models in finite dimensions 27,47 . One might think that variational parameters in one-body part are not so important to describe the dynamics qualitatively because optimizing the one-body part affects only the amplitude of the oscillation after the quench in Fig.4  (a). However, as we see clearly in Sec.IV B 3, variational parameters in one-body part play an important role even in a qualitative description of nonequilibrium states. In Fig. 4(a), we do not see any improvements by operating the spin projection to L K=0 P(t) |φ(t) . However, for the quench to strong interaction U f = 8.0 in Fig.  4(b), the difference between the two wave functions with the momentum projection is more obvious than that for U f = 4.0. These results suggest that it is better to operate both of the quantum-number projections on the trial wave function especially in the strong coupling region. These results show that in order to obtain accurate results by the t-VMC method, we should operate the quantum-number projection and optimize all the variational parameters. This accuracy sensitive to trial wave functions is similar to that for ground states (see Ref. [31][32][33] and Appendix.B).

Size dependence
In the previous subsection, we have shown the benchmark results for small systems at half-filling. In order to check the accuracy of our trial wave function for larger systems, we here compare the t-VMC results with highly accurate DMRG result in the thermodynamic limit. Since the VMC method offers the results for finitesize systems, we need to investigate the system-size dependence of our results.
In Fig. 5, we present the t-VMC results of ∆n(t) for sudden quench protocol (t q = 0.0) from U i = 0.0 to U f = 1.0 in the one-dimensional Hubbard model at half-filling. In the t-VMC calculations, we operated only the momentum projection to reduce numerical cost since spin projection is not so important in the region of small interaction (See Fig. 4). To check the system-size dependence, we showed three t-VMC results with different system sizes. For comparisons, we also show the results obtained by DMRG, DMFT and DCA 20 . As an impurity solver of the DMFT and DCA calculations, the iterative perturbation theory (IPT) was employed. In the DCA calculation, the reciprocal wavevector K satisfies the following condition: K = 2nπ/N c , where n represents integer and N c represents a cluster size. Here, N c = 64 . As seen in this figure, the jump ∆n(t) obtained by DMRG relaxes slowly with an oscillation. This feature is observed in both the DCA and t-VMC results although the DMFT result does not show the oscillation clearly after t > 1.5. However, the DCA result shows a clear deviation from the DMRG result at long time even when large clus-ter size is used. Tsuji and his coworkers have reported that this deviation comes from the quantum corrections from higher-order diagrams neglected in the IPT 20 . On the other hand, our t-VMC results have a strong sizedependence but approach the result of DMRG at long time as the system size increases. These results imply that our trial wave function in the t-VMC method have successfully included the quantum fluctuation beyond the DCA result.

Superconducting correlation in hole-doped Hubbard model in two dimensions
The doped Hubbard model on the square lattice is one of the simplest models for studying the high-T c superconductivity in copper oxides. In this model for equilibrium, many theoretical studies have proposed the existence of d-wave superconducting states in the hole-doped region with strong on-site interaction 33,[55][56][57] . Recently, some experiments have reported that nonequilibrium superconducting states in copper oxides have been realized even at room temperature 12,13 . In order to theoretically determine whether time-evolved states realize superconducting state, it is necessary to calculate time evolutions of pairing correlation functions.
As a benchmark for the d x 2 −y 2 -wave pairing correlation functions P d (t; r), we consider the hole-doped Hubbard models in two dimensions. The pairing correlation function is defined as Here, ∆ d (r i ) represents the d x 2 −y 2 -wave superconducting order parameter defined as where f d (r) = δ ry,0 (δ rx,1 + δ rx,−1 ) − δ rx,0 (δ ry,1 + δ ry,−1 ) (30) is the form factor which describes the d x 2 −y 2 -wave symmetry and r = (r x , r y ). In Figs  We again note that |φ instead of |φ(t) represents the one-body part fixed through the time evolution at |φ(t = 0) optimized in the ground state before quenching. As seen in Fig. 6(e), only the correlation function at the largest distance obtained by using P G (t)P J (t) |φ shows a large deviation from the exact result, i.e., its value is one order of magnitude lower than the other ones. This result implies that only optimizing the correlation factors is insufficient in reproducing pairing correlations in time evolution. In fact, by optimizing the amplitudes of singlet pairings f ij in one-body part in time evolution, the t-VMC result at the largest distance well reproduces the exact one. It is important to obtain the long-range part of the max|P d (t; r)| accurately because it enables us to detect the emergence of the superconducting phase in large systems. Therefore, to describe different nonequilibrium states flexibly, it is crucial to optimize the one-body part.
From these benchmark results, even the superconducting correlation can be well reproduced by using our trial wave function, which may inspire studies along this line in the t-VMC method. One of the intriguing studies is on the influence of strong laser pulse on superconductivity in correlated electron systems. Some recent works have shown that the hopping amplitude is reduced by applying strong laser and the relative strength of interaction to the transfer become effectively larger than before [61][62][63] . By using this effect, the d-wave superconductivity may grow because, for ground states, it grows as the on-site interaction increases. However, in order to confirm whether nonequilibrium states show a true long-range order or not, calculations for larger systems are required. Studies on nonequilibrium superconducting states in larger systems will be reported elsewhere.

V. SUMMARY AND OUTLOOK
In summary, we have developed a time-dependent variational Monte Carlo method (t-VMC) for strongly correlated electron systems and have shown the benchmark results for the fermionic Hubbard model out of equilibrium. By comparing our t-VMC results with the exact results, we found that our trial wave function well reproduces exact time evolutions in both one and two dimensions. These results show that our trial wave function offers an accurate and efficient description of nonequilibrium states in strongly correlated electron systems.
One of the advantages of the VMC method is its wide applicability. In fact, the VMC method can be applied to complicated ab initio effective models derived by the downfolding scheme and contributed to identifying mechanism of physical properties in real materials 64,65 . Applications of the t-VMC method to such models are intriguing future issues. Furthermore, the VMC method can be applied not only to purely electronic systems but also to electron-phonon coupled systems 66 . Therefore, it would be possible to study the phenomena of relaxation process and photoinduced phase transitions through phonon modes in real materials such as copper oxides.
where α and α represents variational parameters and its complex conjugates, respectively. Although the norm of the wave function is preserved under the exact time evolution, the evolution on the manifold M may break the norm-conservation. To remove the restriction on the norm, the norm-dependent Lagrangian L(α, α) should be normalized. Thus, the modified LagrangianL(α, α) for norm-independent dynamics is introduced: The variation of the corresponding action δS with respect to variations of parameters ψ α | → ψ α | + δψ α | is given by Stationarity of the modified action δS = 0 leads to the variational equation on the manifold M Thus, the condition of minimizing the modified action for norm-independent dynamics is equivalent to Eq. (2) in the full Hilbert space. Based on Eq.(2) or Eq.(A1), we can easily derive the Euler-Lagrange equation described asα where a matrix S and a vector g are described as S kl = ∂ ∂α k ∂ ∂α l ln ψ α |ψ α , g k = ∂ ∂α k ψ α |H|ψ α ψ α |ψ α .

Appendix B: Dependence on Trial Wave Functions for Ground States of Hubbard Models
In this appendix, to gain insight into the dependence of physical properties on trial wave functions in the nonequilibrium states, we compare them with the benchmarks for the ground states of the Hubbard models. Tables I and II show how physical properties depend on trial wave functions in the ground state of the Hubbard models at and away from half-filling, respectively. To show the accuracy of our VMC results, we also show the results obtained by the exact diagonalization (ED) method. In these Tables, |φ F and |φ opt represent the Fermi sea state and optimized pair-product wave function, respectively. In all the cases in Table I, there are large discrepancies, especially in the jump ∆n, between the results obtained by using the Gutzwiller-type wave function (GWF) P G |φ F and those of the ED method. The main reason for this is that GWF cannot describe insulating states as we described in Sec.III. By operating the Jastrow factor, the VMC results at half-filling are qualitatively consistent with those of the ED method. For both half-filled and hole doped models, the energies obtained by using our best trial wave function L K=0 L S=0 P G P J |φ opt agree with those of the ED method, and the relative errors are less than 0.5%. Figure  7 shows the pairing correlations obtained by using different trial wave functions in the doped Hubbard model on the square lattice for U/t hop = 8.0, n = 12/16. As shown in Fig.7, L K=0 L S=0 P G P J |φ opt has the best accuracy of the d x 2 −y 2 -wave superconducting correlations for all the distances. These trends on trial wave functions are similar to those in the t-VMC method.
TABLE I. Comparison of physical quantities obtained by the exact diagonalization (ED) method with those by using different trial wave function for one-dimensional Hubbard model at half-filling. E/Ns, ∆n and Ss(π) represents energy per site, jump of the momentum distribution near the fermi energy, and the spin structure factor, respectively. The numbers in parentheses denote the statistical errors in the last digits.