Replica evolution of classical field in 4+1 dimensional spacetime toward real time dynamics of quantum field

Real-time evolution of replicas of classical field is proposed as an approximate simulator of real-time quantum field dynamics at finite temperatures. We consider $N$ classical field configurations dubbed as replicas which interact with each other via the $\tau$-derivative terms and evolve with the classical equation of motion. The partition function of replicas is found to be proportional to that of quantum field in the imaginary time formalism. As the replica index $\tau$ can be regarded as the imaginary time index, the replica evolution is technically the same as the molecular dynamics part of the hybrid Monte-Carlo sampling and the replica configurations should reproduce the correct quantum equilibrium distribution after the long-time evolution. At the same time, evolution of the replica-index average of field variables is described by the classical equation of motion when the fluctuations are small. In order to examine the real-time propagation properties of replicas, we first discuss replica evolution in quantum mechanics. Statistical averages of observables are precisely obtained by the initial condition average of replica evolution, and the time evolution of the unequal-time correlation function, $\langle x(t) x(t')\rangle$, in a harmonic oscillator is also described well by the replica evolution in the range $T/\omega>0.5$. Next, we examine the statistical and dynamical properties of the $\phi^4$ theory in the 4+1 dimensional spacetime, which contains three spatial, one replica index or the imaginary time, and one real-time. We note that the Rayleigh-Jeans divergence can be removed in replica evolution with $N \geq 2$ when the mass counterterm is taken into account. We also find that the thermal mass obtained from the unequal-time correlation function at zero momentum grows as a function of the coupling as in the perturbative estimate in the small coupling region.


Introduction
Classical dynamics has been utilized to understand non-equilibrium evolution of quantum many-body systems in various fields of physics [1][2][3][4][5][6][7][8][9]. The classical equation of motion for the phase space distribution (Vlasov equation) [10] is known to provide an approximate solution of the quantum equation of motion for the density matrix (von Neumann equation) [11], provided that the classical analogue of the quantum mechanical distribution function [12] is given as the initial condition and the O( 2 ) effects are negligible. It is also known that one can numerically obtain the solution of the Vlasov equation by solving the classical equations of motion for particle ensemble representing the phase space distribution [13]. This favorable feature of classical dynamics has been invoked also in field theories [3][4][5][6][7][8][9]. For example, the classical Yang-Mills (CYM) field has been adopted to describe the initial stage of highenergy heavy-ion collisions, and have provided important insights into the non-equilibrium dynamics of the gluon field [5][6][7][8][9].
Compared with the successes in the far-from-equilibrium stages, the applicability of classical dynamics is limited when discussing equilibrium properties of quantum systems. Since the equipartition law applies to classical equilibrium, the number of high-momentum particles is overestimated and one encounters the Rayleigh-Jeans divergence. One possible way to manage the divergence is treating the hard modes above the cutoff separately. By integrating hard modes [14,15] or by introducing the mass counterterm [16], one can obtain the effective action of the classical field, soft modes below the cutoff, and utilize the action to evaluate the evolution. In CYM theory, dynamical evolution of the coupled system of classical field and particles is explicitly solved and was demonstrated to promote equilibration [6]. Still, classical field obeys classical statistics, then the cutoff momentum should be chosen to be of the order of T or smaller also in these frameworks. While the two-particle irreducible (2PI) effective action approach can treat classical field and particles on the same footing [17][18][19], the numerical cost is large and it is not yet easy to apply to realistic systems under inhomogeneous classical field.
Thus it is desirable to develop frameworks which inherit the merit of classical field dynamics but properly describe quantum statistical equilibrium after a long time evolution. Including these two features is known to be important in nuclear transport phenomena [20]. In the stochastic quantization [21], one can obtain field configurations {φ} by solving the Langevin equation, dφ x /dt = −∂S/∂φ x + ζ x , with t being the fictitious time and ζ x being the white noise, ζ x (t)ζ y (t ) = 2δ xy δ(t − t ). The field distribution approaches the quantum one, while the above Langevin equation cannot be regarded as the equation of motion to describe the real time evolution. There is a hint to incorporate the quantum statistical property into the real time evolution in the imaginary time formalism of finite temperature quantum field theory, where the field variables in 3D space are enlarged to those in 3+1D spacetime introducing the imaginary time. In the path integral representation, the thermally equilibrated quantum field distribution is described by exp(−S[φ]), where S[φ] is the 3+1D Euclidean action. In the molecular dynamics part of the hybrid Monte-Carlo (HMC) sampling [22], the Hamiltonian is set to be H = π 2 x /2 + S[φ] with π x being the canonical conjugate of the field variable φ x at a spacetime point x, and the classical equation of motion, dφ x /dt = ∂H/∂π x and dπ x /dt = −∂H/∂φ x , is solved with the initial condition of π 2 x = 1, where the time variable t is introduced in an ad hoc manner. After a long time evolution, the system reaches the equilibrium described by the classical partition function at temperature of unity, Z = Dπ Dφ exp(−H) ∝ Dφ exp(−S[φ]), then we can correctly sample the quantum field configuration in equilibrium.
In this article, we generalize and describe the "time evolution" by the molecular dynamics in the imaginary time formalism as the real time evolution of a set of classical field configurations, which are referred to as replicas. Each replica conceptually corresponds to the field configuration on each imaginary time index in the finite temperature quantum field theory, and then the replica evolution practically corresponds to the molecular dynamics part of 2/26

Replicas of classical field
Interaction btw Replicas Fig. 1 Replicas and their evolution. Replica configuration (φ τ x , π τ x ) evolves with the classical equation of motion using the Hamiltonian H. The interaction part of the replica Hamiltonian H is chosen so that the φ part of H agrees with the Euclidean action S[φ] multiplied by the lattice anisotropy ξ = a/a τ . Thus the replica partition function Z becomes proportional to the equilibrium quantum field partition function, D exp(−S[φ]).
HMC. As schematically shown in Fig. 1, the replica Hamiltonian H is given as the sum of the Hamiltonian of each replica plus the interactions between the nearest neighbor replicas, V τ,τ +1 , which is referred to as the τ -derivative term and causes thermalization of replicas. The distribution of one replica configuration regarded as the system relaxes to the quantum equilibrium distribution by the τ -derivative interactions with other replicas regarded as the heat bath. Thus, on the one hand, field variables φ reaches correct quantum statistical distribution after a long-time replica evolution. On the other hand, the replica-index average of field variables obeys the standard classical field equation of motion when fluctuations are small as shown later. Hence the replica evolution can reproduce both the quantum statistics and the classical evolution in these two limits. These features of the replica evolution are encouraging for us to consider it as a candidate which describes dynamical evolution of quantum field, while the formal justification of the replica evolution as a quantum time evolution has not been found. In order to examine the validity of replica evolution, we investigate the real-time evolution of the unequal-time two-point function at zero momentum numerically without and with the mass counterterm. We demonstrate that the thermal mass obtained from the replica evolution is consistent with the perturbative calculation results in quantum field theory.
This article is organized as follows. In Sec. 2.2, we introduce the replica evolution in quantum mechanics, and the statistical and dynamical properties of a harmonic oscillator are examined. In Sec. 3, we introduce the replica evolution in classical field, and the statistical properties of replicas are examined in the free field case. We also discuss the mass counterterm on the lattice. In Sec. 4, we show the results of real-time evolution of replicas in the φ 4 theory. By using the damped oscillator ansatz, we fit the unequal-time two-point function of replica evolution, and compare the obtained thermal mass and damping rate with the perturbative calculation results. Section 5 is devoted to the summary and perspectives. 3/26 2. Replica evolution in quantum mechanics

Statistics and dynamics of replicas in quantum mechanics
In this section, we introduce replica evolution in quantum mechanics and examine its statistical and dynamical properties. We consider the classical system described by the following Hamiltonian where the mass is set to be unity for simplicity, and D is the number of components in x and p. We now consider N replicas of canonical variables, (x τ , p τ )(τ = 0, 1, · · · , N −1), whose Hamiltonian is given by the sum of the Hamiltonians H(x τ , p τ ) over the replica indices τ and the τ -derivative terms, The periodic boundary condition is imposed in the τ direction, (x N , p N ) = (x 0 , p 0 ). The replicas are assumed to evolve in real time t according to the canonical equations of motion, The replica evolution with Eq. (3) has two distinct features: First, the ensemble of replica configurations follows the quantum statistical distribution of the spatial variables x = {x τ |τ = 0, 1, · · · , N − 1} in the long-time evolution. Second, the evolution of replica-index average agrees with the purely classical evolution, when the fluctuations among replicas are small. Let us examine the first point on the statistics of replicas. The partition function of replicas at temperature T repl = ξ is given as where β = N/ξ and τ /ξ →τ is the continuous imaginary time. Since S E [x] is the Euclidean action at temperature T = 1/β = ξ/N in the imaginary time formalism, the replica partition function Z R at temperature T repl = ξ is proportional to the quantum mechanical partition function at T = ξ/N in the large N limit. Thus observables as functions of x in quantum equilibrium are correctly obtained from the thermal average of observables O(x) T in classical equilibrium of replica configurations, 4/26 where τ in the first line is the replica index of observation, and in the second line, we replace O(x τ ) with the "replica-index average", using the translational invariance in the τ -direction. Then, a thermal expectation value of an operator O(x) T can be obtained as an expectation value of the replica-index averaged operator, In the imaginary time formalism, the thermal average of an observable in quantum mechanics is given as where the imaginary time of observationτ appears in the path integral representation. After replacing O(x(τ )) with its imaginary time average, the quantum statistical average (10) is found to be described by the replica average (7) in the large N limit. The classical equilibrium of replica configurations can be generally obtained by the longtime evolution with the canonical equation of motion, Eq. (3), due to the chaoticity of the system. Practically, it is useful to take the "replica ensemble average" instead of the long-time average, since there is no autocorrelation in the former. We prepare N conf ( 1) initial replica configurations (x (i) , p (i) )|i = 1, 2, · · · , N conf at t = 0 appropriately, solve the equation of motion, then the replica average is calculated as Let us turn to the second point. While replica ensemble simulates quantum statistical ensemble after a long time evolution, the replica-index averages of the canonical variables evolve as the classical variables when the fluctuations among replicas are small. The equation of motion for the replica-index average of the canonical variables reads where x iτ is the i-th component of x τ and (δx) 2 = τ (x τ − x) 2 /N is the fluctuations of x τ in one configuration of replicas. It should be noted that the τ -derivative terms do not operate because of the periodic boundary condition, and the second derivative terms of V

5/26
disappear from the definition of x. It is interesting to find that the first line of Eq. (13) shows the Ehrenfest's theorem where · · · denotes the replica-index average here. Then when (δx) 2 is small, these equations of motion tell us the classical nature of the replica-index average, By comparison, (δx) 2 should be a part of quantum fluctuations.

Replica evolution of harmonic oscillator
We now look further into the dynamical property of replicas of a single harmonic oscillator. By choosing V = ω 2 x 2 /2 and D = 1 in Eq. (1), the Hamiltonian is given as In order to investigate the statistical properties of replicas, it is useful to adopt the Fourier transform with respect to the replica index, where ω n = 2πn/N denotes the Matsubara frequency. With (x n ,p n ), the Hamiltonian is represented by that of the N free harmonic oscillators, The replica partition function is obtained as By using the Matsubara frequency summation formula explained in Appendix A, the logarithm of the partition function is found to be where Ω is given as The replica partition function, Eq. (19), at large N agrees with the quantum mechanical partition function at T = ξ/N , Thus we find that it is possible to obtain quantum statistical results by using a replica ensemble in equilibrium, a large number of configurations of N sets of canonical variables, (x τ , p τ ), which are prepared to be in classical equilibrium. 6/26 We next examine the time-evolution of replicas by using the time-correlation function, C(t) = x(t)x(0) . For preparation, let us recall the quantum mechanical results. The time correlation of the spatial coordinate is calculated as where x H is the operator in the Heisenberg picture, |n is the energy eigen state, c n is the expansion coefficient, |ψ = n c n |n , and we have used the relations, n−1|a|n = √ n and n+1|a † |n = √ n+1. In thermal equilibrium at temperature T , the density matrix becomes diagonal and the statistical weights are given by the Boltzmann factor, Thus the time-correlation function in thermal equilibrium is obtained as where we have used the relation n ne −nx = d/dx( n e −nx ) to obtain the second line. Since the expectation value of symmetrized product (Weyl ordering) is obtained in classical dynamics, we are interested in the time-even part part of the time-correlation function given as In replica evolution, it is easier to solve the canonical equation of motion in the Fourier transform. The equations of motion are, dx n /dt = ∂H/∂p n =p n and dp n /dt = −∂H/∂x n = −M 2 nx n , and the solution is obtained as We evaluate thermal average imposing thermal distribution to the initial field variables, x n (0) andp n (0). Since the Boltzmann weight is given as exp , the distributions ofx n andp n in equilibrium are gaussians, x n (0)x n (0) T = ξ/M 2 n δ nn and p n (0)p n (0) T = ξδ nn . Then the time-correlation function is obtained as We here adopt the replica-index average, the average of the time-correlation function over the replica indices τ , and the replica ensemble average, the average over the initial configurations.

7/26
The zero Matsubara frequency contribution in the replica formalism agrees with the quantum mechanical result in the high-temperature limit, The other Matsubara frequencies lead to higher harmonics, which do not appear in the quantum mechanical result. Yet they play an essential role to explain the value of equal-time correlation C R (0) = x 2 (0) at lower temperatures, In Fig. 2, we show the temperature dependence of the equal-time correlation function, C(0) = x 2 , which provides the amplitude of C(t). As already mentioned, C R0 (0) agrees with the quantum mechanical result at high temperatures, T ω. At lower temperatures, C R0 (0) itself deviates from the quantum mechanical result, while the other Matsubara frequency contribution lifts up the amplitude. As a result, the amplitude in the replica method increases with increasing N , and converges to quantum mechanical result in the large N limit. We confirm that the replica method reproduces the quantum mechanical result of the equal-time correlation x 2 = C Q (0) in the large N limit, which is not a total surprise since the replica method gives correct thermal quantum distribution after a long time evolution. Let us come back to the evaluation of the time-correlation functions. In Fig. 3, we show the time-correlation function of a harmonic oscillator in quantum mechanics and in the replica evolution. When the temperature T is comparable to or higher than the intrinsic frequency ω, the quantum mechanical time-correlation function is well described by the replica evolution which is dominated by the zero Matsubara frequency contribution C R0 (t). When T is smaller than ω, the amplitude of C R0 (t) is smaller than the quantum mechanical result and the contributions of higher harmonics are not negligible. As a result, C R (t) fluctuates around C R0 (t) at lower temperature as shown in the right panel of Fig. 3 for T /ω = 0.5: the result in the replica method with N = 4 at t = 0 agrees with that in the quantum mechanical result, while we find deviations at t > 0. The deviation from the quantum result can be evaluated by the fraction of sum of non-zero n amplitudes and the full amplitude,

10
which amounts to be 23.8 % (7.6%) at T /ω = 0.5 (T /ω = 1) at large N . The time-correlation function fluctuate around C R0 (t) with maximal deviation given in Eq.(31) even at very large N . Yet the global behavior of the time-correlation function is described well by C R (t) in the replica formalism in the region T /ω 0.5. In this section, we have demonstrated that the replica evolution provides configurations of x in accordance with the quantum statistical distribution, provided that the system is thermalized by some interactions with the heat bath or the system is chaotic. This should be also valid with any interaction, since Eq. (19) holds for any Hamiltonian. As for the time evolution, the quantum mechanical time-correlation function is well described by the replica evolution in the temperature region of T /ω 1 and is reasonably described at T /ω 0.5, while the higher harmonics (non-zero n) distorts the time-correlation function at T /ω 1. It should be noted that the n = 0 contribution is the same as the standard classical dynamics in the harmonic oscillator, but differences of x in different replica index modify the equation of motion even for the n = 0 modes when interactions are switched on among replicas. 9/26

Replica evolution in scalar field field
In this section, we apply the replica evolution method to the scalar φ 4 field theory on the lattice. Since the quantum field theory is the quantum mechanics of field variables on spatial points, we can utilize the method introduced in quantum mechanics also in field theories. Thus the following discussions proceeds in parallel to those in quantum mechanics.

Classical Scalar Field Theory on Lattice
We consider the φ 4 theory, where the Lagrangian is given as On a L 3 lattice, the Hamiltonian is given as where (φ x , π x ) are the canonical variables and x = ( represents the lattice spatial coordinate. Throughout this article, we take all quantities normalized by the lattice spacing a.
The classical evolution of field variables is described by the canonical equations of motion, After long-time evolution, classical field distribution relaxes to the classical statistical equilibrium, where the classical partition function is given as, Since thermally equilibrated classical field obeys the the Rayleigh-Jeans law, each momentum mode approximately carry the energy of T , and the energy density is divergent in the continuum limit, a −1 → ∞. Since high-momentum modes are not suppressed by the exponential (Boltzmann or Bose-Einstein) factor, results are sensitive to the cutoff. Thus it is necessary to choose the cutoff appropriately in order to deduce the results in quantum systems [8,9]. It is also possible to adopt the Hamiltonian with the mass counterterm to avoid the divergence of the mass [16], but we cannot avoid the classical field to relax to the classical statistical equilibrium as long as one classical field configuration evolves with the classical equation of motion.

Replica evolution
We next consider N replicas of classical field, which interact with the nearest neighbor replicas via V, given in the τ -derivative form, where (φ τ , π τ ) = {(φ τ x , π τ x ) | x i = 0, 1, . . . , L − 1} represents τ -th replica of classical field and the replica index takes the value τ = 0, 1, · · · , N − 1. We impose the periodic boundary 10/26 condition in the τ direction, (φ N x , π N x ) = (φ 0x , π 0x ). Provided that we solve the classical equation of motion with the Hamiltonian H, the partition function of replicas at temperature T repl is given as where Especially, in the case where the replica temperature is chosen to be T repl = ξ, we find In this case, we can regard S[φ] as the Euclidean action of quantum field in the imaginary time formalism, where the lattice spacing in the imaginary time direction is given as a τ = a/ξ and the parameter ξ is now interpreted as the lattice anisotropy. The prefactor (1/ξ) shows the spacetime volume of one cell in the lattice unit, a 3 a τ /a 4 = 1/ξ. The first term in the square bracket in Eq. (40) can be regarded as the squared derivative of φ with respect to the continuous imaginary time given as [∂φ/∂τ ] 2 withτ = τ /ξ. The period of the imaginary time is N/ξ = 1/T , where T is the temperature of quantum field. Now we find that the replica evolution of classical field at replica temperature T repl = ξ = N T gives equilibrium quantum field configurations at temperature T , where N is a normalization constant. The difference of the present replica evolution and the standard classical field evolution comes from the τ -derivative interaction term V. The equation of motion for φ in the replica evolution reads, where H τ = H(φ τ , π τ ) and we have erased π τ x using the equation of motion. The second term from V, which is characteristic of the replica formalism, tends to reduce the difference of the nearest neighbor replica field variables at the same spatial point, φ τ x and φ τ ±1,x , and keeps the replica ensemble in quantum equilibrium as found in the partition function Z R in Eq. (41). Without V, small difference of (φ, π) between replicas in the initial condition leads to very different field configurations after a long time evolution, since interacting classical field is generally chaotic. As a result, classical field configurations relax to classical statistical equilibrium and quantum thermal equilibrium cannot be kept. As in the cases in quantum mechanics, replica ensemble gives quantum statistical ensemble after a long time evolution, and the replica-index average of the classical field φ τ x evolves 11/26 like classical field. The equation of motion for the replica-index average of field variables reads where (δφ x ) 2 = τ (φ τ x − φ x ) 2 /N denotes the variance of φ at x in one replica configuration. Then we get the equation of motion for φ x , when the fluctuations of φ x in the replica configuration is small. The equation of motion given as Eq. (46) is the same as that for the classical field, expectation value of φ(x). When the fluctuations are not negligible, they modify the equation of motion for the classical field φ. For example, (δφ x ) 2 contributes to the mass as m 2 → m 2 + λ(δφ x ) 2 /2. As in the quantum mechanics case, the thermal average O(φ x ) T of an arbitrary observable O(φ x ) is defined as an average over the replica index τ and the thermal replica ensemble, Since the "classical field" variables are obtained by the replica-index average, fluctuations among the field configurations with different replica indices in one replica configuration should be regarded as a part of quantum fluctuations. Fluctuations in replica configurations may contain statistical and quantum fluctuations. One replica configuration would not be enough to describe a quantum state, and we need at least several replica configurations to satisfy the uncertainty principle. Further fluctuations would be considered as statistical. Thus taking both of the replica index and ensemble averages would be reasonable to take account of quantum and statistical fluctuations. While the time of the functional integration variables in Eq. (47) is (implicitly) assumed to be the same as that for the observable, these times can be different. Since the classical time evolution of canonical variables is the canonical transformation and the Hamiltonian H is a constant of motion on the classical path, the integration measure is the same, Dπ(t)Dφ(t) = Dπ in Dφ in with (π in , φ in ) being the initial field variables, and the statistical weight is also the same, exp(−H(φ(t), π(t))/ξ) = exp(−H(φ in , π in )/ξ). Thus the thermal average can be regarded as the "initial replica configuration average", provided that the initial replica ensemble is sampled according to the statistical weight and the number of samples is large enough, We adopt this prescription in the later discussions. 12/26

Partition function of Free Field
Let us discuss the equilibrium property of replicas of the free field (λ = 0), where one can obtain the partition function analytically. The Hamiltonian Eq. (36) is represented by the Fourier components, Lattice momentum and the Matsubara frequency are defined as k = (k 1 , k 2 , k 3 ) (k i = 2πm/L, m = 0, 1, . . . L − 1) and ω n = 2πn/N (n = 0, 1, . . . N − 1), respectively. Then the partition function is given by the Gaussian integral, Z (λ=0) R = k,n (ξ/ω nk ), where the integration measure is specified as DπDφ = x,τ dπ τ x dφ τ x /(2π) = k,n dπ nk dφ nk /(2π). By using the Matsubara frequency summation formula explained in Appendix A, the logarithm of the partition function for each lattice momentum k is found to be where Ω k is given as Now the partition function reads The energy expectation value for each momentum k is obtained as, where we have used the relation ξ = N T and · · · denotes the thermal expectation value at temperature T . In the low frequency limit, ω k /T 1, this energy converges to the classical value, E (λ=0) k → T . The factor in front of the parentheses converges to unity in the large N limit, ξ = N T → ∞. The first term in the parentheses is the zero point energy, and should be subtracted in field theories. The second term in the parentheses represents the thermal energy. Compared with the classical field theory, the Bose-Einstein distribution function appears and the high-momentum components are exponentially suppressed in the thermal part of energy. 13/26
By using Eqs. (56) and (57), we can evaluate spacetime dependence of the field variables and the two point functions, The thermal ensemble average for φ nk (0) and π nk (0), whose distributions are Gaussians, is taken as Then we find that the two point functions are given as In the later discussions, the following two point functions will be used and discussed, cos ω n0 t .
The first one (∆ = φ 2 ) appears in the one-loop diagram and diverges in the continuum limit. The second one (C(t)) is the unequal-time two-point function at zero momentum, referred to as the the time-correlation function in the later discussions, and is expected to show oscillatory behavior with frequency of the thermal mass. In equilibrium, the "trigger" time of the measurement, t 0 , can be taken arbitrary. 14/26

Mass renormalization
In the replica evolution with finite coupling, we need to take care of the mass renormalization as in the standard treatment of quantum field theory. We consider the contribution of the one-loop diagram and the counterterm shown in Fig. 4, then the thermal mass including the contribution from the interaction is found to be We choose the counterterm δm 2 so that it cancels the divergent contribution to the mass, λ φ 2 div /2, where φ 2 div is the divergent part of ∆ = φ 2 given in the first term in the bracket in Eq. (63). The mass term induced by the interaction, λ φ 2 /2, coincides with the factorization (Wick contraction) of the interaction term appearing in the equation of motion, λφ 3 /6 λ φ 2 φ/2.  (67) self-consistently by taking account of M dependence of ω k and Ω k . For comparison, we also show the results of the leading order estimate in the continuum and large N limit, 1/a → ∞, L → ∞ and N → ∞ [23], For m = 0, we also show the results with resummed one-loop contribution from the selfconsistent treatment [23,24], and the two-loop calculation result [24],

Numerical results of replica evolution in scalar field theory
We shall now numerically evaluate the time evolution of replicas of classical field by using the replica Hamiltonian Eq. (36) defined in the 4D spacetime including the imaginary time. In order to examine the validity of replica evolution, we discuss the time-correlation function C(t), which is the unequal-time two-point function at zero momentum. From the time correlation C(t) obtained by the replica evolution, we extract the thermal mass M and the damping rate γ, and compare them with the perturbative estimates.

Setup
We show the numerical results of time evolution of replicas on a 32 3 × 4 lattice (L = 32, N = 4) at T = 0.5 in the coupling range of 0.5 ≤ λ ≤ 10 with m = 0 and m = 0.5, as an example. We prepare the initial condition by using the Langevin equation at the replica temperature of ξ = N T = 2, where the drift constant is taken to be Γ = 0.5 and ζ τ x (t) is the white noise, ζ τ x (t) ζ τ x (t ) = δ τ,τ δ x,x δ(t − t ). Since the relation of π 2 τ x = ξ should be satisfied in equilibrium, we rescale π field at each step of the Langevin evolution. The equilibration time to prepare the initial condition is set to be t eq = 20 in the lattice unit, which is found to be reasonably long in the coupling range for λ ≥ 4. At smaller coupling, we take τ eq = 100, 60 and 40 for λ = 0.5, 1 and 2, respectively. The time step is taken to be ∆t = 0.025 and the equation of motion is solved in the leap-frog method until t = 500 after the equilibration. When we take 16/26 account of the mass renormalization, we subtract the divergent part of the induced mass by using the one-loop calculation results given in Eq. (66).
We evaluate the time-correlation function C(t) of the zero momentum component of the field variable, φ L 3 , with i being the configuration index. The equilibrium average of the correlation function is obtained as the average over the replica ensemble, as shown in Eq. (64), where the replica ensemble is prepared by the Langevin equation discussed above. The number of replica configurations is taken to be N conf = 1000. In order to reduce the statistical error, we also take average over the trigger time t 0 .
With this setup, we have solved the time evolution of replica configurations, where the ensemble average should be consistent with the equilibrium expectation value. The ensemble average of π 2 τ x is found to be (0.1 − 0.2) % larger than ξ in the present setup. The overestimate can be suppressed with smaller drift coefficient and longer equilibration time, but it takes more time for the calculation and the above deviation would be small enough.

Momentum distribution and Rayleigh-Jeans divergence
Before discussing real time evolution, let us take a look at the thermal expectation value of the momentum distribution, as a function of momentum k = (k 2 ) 1/2 . This appears in ∆ = φ 2 (Eq. (63)) and also in the energy in the form of ω 2 k |φ k | 2 = (m 2 + k 2 ) |φ k | 2 in Eq. (49). In the free field case, the momentum distribution is given as as discussed in Sec. 3.4. The first term in the bracket in Eq. (73) shows the zero point energy contribution which should be subtracted, and the second term shows the thermal part of the momentum distribution on the lattice which converges to the Bose-Einstein distribution function 1/[exp(ω k /T ) − 1] in the large N limit. In Fig. 6, we show the replica (N = 4) and classical field (N = 1) ensemble results at m = 0 and λ = 8. We here adopt the thermal mass evaluated from the time-correlation function discussed later. These results agree with the free field results on the lattice. When the divergent part is subtracted (right panel), the momentum distributions approximately show exponential decay, as the Bose-Einstein distribution does. It may be interesting to find that even in the case of classical field, an approximate exponential decay is found, while the thermal part in the free field on the lattice deviate those in the large N limit at high momenta, k > 1.5. This approximate exponential decay comes from the decomposition of the divergent and finite parts as shown in Eq. (73). But this is not the end of the story. Next, let us discuss the Rayleigh-Jeans divergence. The momentum distribution appears in energy in the form of k 2 |φ k | 2 and the number of momentum modes increase with the momentum as k 2 . Thus the energy contains the kinetic energy part of dkk 4 |φ k | 2 /2π 2 in the continuum limit. In Fig. 7, we show the momentum distribution multiplied by k 4 . We note that the classical field results (N = 1) seem to saturate to a constant value. This behavior can be understood from the decomposition in Eq.  Fig. 7 Renormalized momentum distribution multiplied by k 4 , k 4 |φ k | 2 ren , obtained in replica (N = 4, circles) and classical field (N = 1, diamonds) ensembles on the 32 3 lattice at m = 0, λ = 8, and T = 0.5. Solid curves show the thermal part of the distribution in the free field on the lattice, and the and dashed curve shows the large N limit, corresponding to the Bose-Einstein distribution multiplied by k 4 /ω k . 18/26 decomposition in Eq. (73) is not necessarily exponentially suppressed but rationally suppressed at large k. Because of the functional form, arcsinh x = log( √ 1 + x 2 + x) log(2x) at large x, the "exponential" reads exp(−Ω k /T ) exp(−2N log(ω k /N T )) = (ω k /N T ) −2N for large ω k , ω k ξ. Then for high momentum, k m and k N T , we find Then k 4 |φ k | 2 converges to a constant with N = 1. By substituting T = 0.5 and N = 1, the asymptotic value is found to be 2(N T ) 2N +1 = 0.25, which is close to the classical field results at large k. Thus we cannot fully remove the Rayleigh-Jeans divergence in the classical field in the present subtraction scheme without additional matching procedure. In contrast, the replica results (N = 4) of the momentum distribution show suppressed behavior at high momentum as the Bose-Einstein distribution does, and k 4 |φ k | 2 also decreases at large k. For the convergence of energy, k 4 |φ k | 2 needs to converge to zero faster than 1/k, then we find the constraint on N as 2(N − 1) > 1, or N > 3/2. Thus we can expect that the integral would converge to a finite value also in the continuum limit, and the Rayleigh-Jeans divergence in energy can be fully removed in the replica evolution even with N = 4.

Time-correlation function, thermal mass and damping rate
We now proceed to discuss the time-correlation function C(t) obtained from the time evolution of replica ensemble. In Fig. 8, we show the time-correlation function C(t) from replica evolution at λ = 2 and 8 with mass renormalization. We find that the time-correlation function is well described by the single damped oscillator f 1ω (t), as shown by the thick blue curves. We have obtained the thermal mass M and the damping rate γ by fitting the parameters (M, γ, A, δ) in f 1ω to C(t) obtained from the replica evolution. In Fig. 9, we show the coupling dependence of the thermal mass M obtained from the time-correlation function C(t) in replica evolution without (left) and with (right) mass renormalization. We show the fitting results using the single oscillator function f 1ω . The fitting results are close to the one loop calculation results without mass renormalization. With mass renormalization, the thermal mass is consistent with the one loop results at small coupling but considerably smaller than the one loop results in the strong coupling region. Replica evolution results at m = 0 agree with the two loop results, while the classical field results (N = 1) at m = 0 show weaker reduction from the leading order results. Thus the replica evolution is expected to give higher order interaction effects over the one loop. For more serious comparison, we need to take account of two-loop counterterms of mass and coupling, and to choose the renormalization scale µ consistent with the present lattice calculation, but these are beyond the scope of this work.
In Fig. 10, we show the damping rate γ as a function of the coupling. We compare the replica evolution results with mass counterterm in comparison with the two loop calculation 19/26  results [24], which is known to agree with the plasmon damping rate after the matching to quantum theory by substituting the leading order thermal mass estimate, M LO = λT 2 /24. In this expression, the product M γ is found to be independent of the mass. The replica evolution results with m = 0.5 seem to roughly agree with the perturbation calculation results in the small coupling region, and start to deviate from the perturbative estimate at larger coupling, λ ≥ 4. Deviations at λ ≥ 4 would be due to higher-order effects of the coupling, fluctuations, or the higher-momentum components ignored in the lattice discretization. With m = 0, the damping rates are smaller than the perturbative estimate in the strong coupling region, as in the m = 0.5 cases. At λ ≤ 4, by comparison, the damping rate tends to be larger than the perturbative estimates. Since the thermal mass is small in this region of coupling, M = 0.085 and 0.12 at λ = 0.5 and 1, respectively, we may need larger size lattice. The apparent larger damping rate at small thermal mass may be related with the fragmentation of the single particle mode into several modes. The Fourier transform at small coupling and m = 0 is found to be represented better by the superposition of several damped oscillators, each of which has a small width. If we adopt the width for each of the modes as the damping rate, the results at λ = 0.5 and 1 agree with the perturbative estimates as shown by open triangles in Fig. 10. The analysis of the Fourier transform is given in Appendix B. It should be also noted that the classical field results roughly agree with the perturbative calculation, as already noted in previous works [16,18,25].
Before closing this section, we would like to mention that C(t) is the time-symmetric part of the two-point function, C(−t) = C(t), whose transformation is referred to as the statistical function. In quantum field theory, the spectral function is more important, but we need to evaluate the commutator of unequal-time field variables by using, for example, the Poisson 21/26 bracket [25]. We show the time-odd part of the time-correlation function in a harmonic oscillator in Appendix C, but we leave evaluating the spectral function in field theories in the future work.

Summary and perspectives
We have investigated the simultaneous real-time evolution of several classical field configurations, referred to as replicas, which interact with the nearest neighbor replicas with a specific form of interaction, the τ -derivative term. Classical evolution of replica ensemble is found to provide the correct quantum field partition function in equilibrium. The average of field variables over replica indices approximately obeys the classical field equation of motion when the fluctuations among the replicas are small. The exponential suppression factor of high momentum modes appears from the sum over the replica index, which can be regarded as the imaginary time, provided that the zero point energy contribution is subtracted. We have examined the behavior of the time-correlation function (the unequaltime two-point function) at zero momentum without and with the mass counterterm. The time-correlation function in the replica evolution is expected to show oscillatory behavior with the frequency of the thermal mass M at temperatures T /M 0.5, as demonstrated in the quantum mechanics of the harmonic oscillator. The time-correlation function in the φ 4 theory seems to show reasonable behavior: The thermal mass of the zero momentum mode roughly agrees with the perturbative calculation results at small coupling. Thus the replica evolution should be useful to describe real-time behavior in equilibrium.
It would be desired to further examine the replica evolution as a candidate of the frameworks to describe non-equilibrium real-time quantum-field evolution. As long as the distribution of initial replica Hamiltonian values is properly given, distribution of field variables in replica ensemble should finally relax to the correct quantum statistical distribution, even if one starts from far-from-equilibrium configurations. In discussing non-equilibrium real-time evolution, however, it would be necessary to introduce additional time scale. It should be noted that replica evolution is conjectured to be useful in a heuristic context, but it is not derived based on some principle. Thus formal derivation or justification is desired. For example, the equivalence between the classical field theory and the Boltzmann equation [26] would be a good guide for the formal discussions. It is also interesting to discuss, for example, the O(N ) model, where there exist results of dynamical calculations using the two particle irreducible (2PI) effective action [17,18].
Once the present framework is proven to be useful in describing quantum field evolution toward equilibrium, application to the Yang-Mills field is another important subject to study. The temporal component of the vector field is usually Wick rotated in the imaginary time formalism and we cannot apply the replica evolution as it is. However, spatial components are the same in the imaginary and real time formalism and we can apply the replica evolution. Thus it is possible to examine the replica evolution of classical Yang-Mills field in the temporal gauge where the temporal components of the vector field are set to be zero. Then it is interesting to examine whether or not the quantum statistical features affect the dynamical evolution in the initial stage of high-energy heavy-ion collisions. In order to describe inhomogeneous and nonequilibrium evolution, it is necessary to take account of spacetime dependence of temperature, which may need to invoke the nonequilibrium statistical operator method [27]. 22/26 A. Matsubara frequency summation In deriving Eq. (52), we have used the Matsubara frequency summation formulae, where g(ω) is an analytic function of ω, does not have poles on the real axis, and decreases faster than 1/ω at ω → ∞, i.e. lim |ω|→∞ ωg(ω) = 0. The poles and residues of g(ω) are denoted by ω 0 and Res g(ω 0 ). Specifically, we consider the following sum where T = 1/N and ω n = 2πn/N . The derivative dS/dE is in the form of Eq. (A1), where g(ω) = 2E/(E 2 + 4 sin 2 (ω/2)) and T = 1/N . The poles and residues of g(ω) are found to be iω 0 = ±Ω = ±2arcsinh E and Res g(ω 0 ) = ±2i/ √ 1 + E 2 , so dS/dE is obtained as By integration, the sum in Eq. (A2) is found to be S = 2 N log [sinh(N Ω/2)] + const. .
The constant can be fixed as 2 log 2(1/N − 1) by considering the large E limit of S, S = 2 log E + O(1/E). By substituting E = ω k /2ξ, we obtain Eq. (52). The same formula can be used to obtain ∆ = φ 2 in Eq. (63).

B. Fourier transform of time-correlation function
In Sec. 4, the damping rate is found to be larger in the weak coupling region with mass counterterm at m = 0. Here we would like to discuss this point by using the Fourier transform 23/26 of the time-correlation function, In Fig. B1, we show ρ(ω) obtained from the replica evolution with mass counterterm at λ = 0.5 and 8 with m = 0 in comparison with the Fourier transform of the fitting function f 1ω . At λ = 0.5, the spectrum has a peak having the width of the order of 10 −2 , but the tails fall off much faster than the behavior expected from the width or the damping rate from the fitting function f 1ω , γ 2.6 × 10 −3 . One of the possible interpretations of this spectrum is to consider that there are several modes, each of which has a small damping rate but has a mass spread in the region with 10 −2 width. As an attempt, we use the two damped oscillator functions, where the two frequencies are close to each other, f 2ω (t) =A exp(−γt) {r cos[(ω + δω)t + δ] + (1 − r) cos[(ω − δω)t + δ]} .
The dotted lines show the Fourier transform of f 2ω fitted to C(t). The tail region is found to be suppressed, but not enough at λ = 0.5. Next, we consider the following fitting function, convolution of the step function and the Lorentzian, for the Fourier transform, The fitting results are shown by the solid curves, which shows both wide width of the peak and the fast fall off. The damping rate of each mode becomes smaller, γ 1.5 × 10 −4 , and roughly agrees with the perturbative estimate, as shown by open triangles in Fig. 10. At larger coupling, fitting results of the damping rate with f 1ω , f 2ω and ρ ∆ are consistent, and the peak part of the spectrum is reproduced in these functions, as shown in the right panel of Fig. B1. 24/26 C. Time-odd part of the time-correlation function in harmonic oscillator While the expectation value of symmetrized (Weyl ordered) product can be obtained in classical dynamics, we need additional care to evaluate the expectation value of the antisymmetrized product such as the commutator. For example, the time-odd part of the timecorrelation function may be obtained by using the quantum-classical correspondence for the commutator [25], where we explicitly show here and {A, B} PB is the Poisson bracket. If we ignore O( 3 ), the time-odd part of the time-correlation function would be obtained as Since the Fourier transformation and the time evolution are the canonical transformation, we can choose either (x τ , p τ ) or (x n ,p n ) in calculating the Poisson bracket and the time t 0 should be arbitrary. The zero Matsubara frequency contribution in Eq. (C2) agrees with the quantum mechanical result.
In the case of coupled oscillators such as the field theories, it is in principle possible to calculate the Poisson bracket of the unequal-time observables by using the Hessian matrix [8], while it requires to store the matrix elements of degrees of freedom squared, (2N dof ) 2 with N dof = N L 3 for one component scalar field theory on the L 3 × N lattice. We also need to multiply the matrix at each step of time, so the numerical cost is much larger than the time-correlation function discussed in this article.