Shear viscosity of classical fields in scalar theory

TL;DR

This paper calculates the shear viscosity of massless classical scalar fields in $^4$ theory using the Green-Kubo formula, revealing a slow decay in the correlation function that dominates viscosity and may be a universal feature of dense matter.

Contribution

It introduces a scaling function approach to shear viscosity in classical scalar fields and identifies a slow decay process influencing viscosity, connecting to molecular dynamics findings.

Findings

Two distinct damping time scales in correlation functions.

Slow decay dominates shear viscosity in the massless classical field.

Presence of damped oscillations in early correlation behavior.

Abstract

We investigate the shear viscosity of massless classical scalar fields in the $\phi^4$ theory on a lattice by using the Green-Kubo formula. Based on the scaling property of the classical field, the shear viscosity is represented using a scaling function. Equilibrium expectation value of the time-correlation function of the energy-momentum tensor is evaluated as the ensemble average of the classical field configurations, whose time evolution is obtained by solving the classical equation of motion starting from the initial condition in thermal equilibrium. It is found that there are two distinct damping time scales in the time-correlation function, which is found to show damped oscillation behavior in the early stage around a slow monotonous decay with an exponential form, and the slow decay part is found to dominate the shear viscosity in the massless classical field theory. This kind of slow decay is also known to exist in the molecular dynamics simulation, then it may be a generic feature of dense matter.

Tables

Table 1. Table 1: Parameters in the exponential decay part of C 12 ​ ( t ) subscript 𝐶 12 𝑡 C_{12}(t) ( A 𝐴 A and Γ Γ \Gamma ) and the scaling function( f η ​ ( λ ​ T ) = λ 2 ​ T ​ η subscript 𝑓 𝜂 𝜆 𝑇 superscript 𝜆 2 𝑇 𝜂 f_{\eta}(\lambda T)=\lambda^{2}T\eta ) at λ = 0.5 , 1 , 3 , 10 𝜆 0.5 1 3 10 \lambda=0.5,1,3,10 and 30 30 30 on the 16 3 , 32 3 superscript 16 3 superscript 32 3 16^{3},32^{3} and 64 3 superscript 64 3 64^{3} lattices. We show the results at T = 1 𝑇 1 T=1 .

$\lambda$	$L$	$A$	$\Gamma /{\lambda }^2$	$f_{\eta }\left(\lambda T\right)$
0.5	16	14$\times {10}^{-3}$	17$\times {10}^{-4}$	8.4
	32	13$\times {10}^{-3}$	8.1$\times {10}^{-4}$	16
	64	11$\times {10}^{-3}$	6.7$\times {10}^{-4}$	17
1	16	11$\times {10}^{-3}$	9.2$\times {10}^{-4}$	12
	32	12$\times {10}^{-3}$	5.2$\times {10}^{-4}$	22
	64	11$\times {10}^{-3}$	6.4$\times {10}^{-4}$	17
3	16	9.6$\times {10}^{-3}$	3.8$\times {10}^{-4}$	25
	32	10$\times {10}^{-3}$	4.4$\times {10}^{-4}$	23
	64	10$\times {10}^{-3}$	4.5$\times {10}^{-4}$	23
10	16	6.6$\times {10}^{-3}$	2.4$\times {10}^{-4}$	27
	32	8.2$\times {10}^{-3}$	2.6$\times {10}^{-4}$	31
	64	7.3$\times {10}^{-3}$	2.5$\times {10}^{-4}$	29
30	16	4.8$\times {10}^{-3}$	1.1$\times {10}^{-4}$	43
	32	4.5$\times {10}^{-3}$	0.82$\times {10}^{-4}$	54
	64	4.7$\times {10}^{-3}$	0.96$\times {10}^{-4}$	49

Equations

$$\mathcal{L}=\frac{1}{2}{\partial }_{\mu }\varphi {\partial }^{\mu }\varphi -\frac{1}{2}{m}^2{\varphi }^2-\frac{\lambda }{24}{\varphi }^4.$$

$$H=\frac{1}{2}\sum _{\mathbf{x}}\left[{\pi }^2(x)+{\left({\partial }_i\varphi (x)\right)}^2+m^2\varphi (x{)}^2+\frac{\lambda }{12}{\varphi }^4(x)\right],$$

$$\dot{\varphi }(x)=\frac{\partial H}{\partial \pi (x)},\dot{\pi }(x)=-\frac{\partial H}{\partial \varphi (x)}.$$

$$T_{ij}$$

$$\eta =\underset{\omega \to 0}{\lim }\frac{1}{T}{\int }_0^{\infty }dt{e}^{i\omega t}{C}_{12}(t),$$

$$C_{12}(t)=V⟨{\tau }_{12}(t){\tau }_{12}(0){⟩}_{\mathrm{eq}},$$

$${\tau }_{12}(t)\equiv \frac{1}{V}\int d^{3}x{\tau }_{12}(\mathbf{x},t),$$

$$C_{12}(t)=V⟨{\tau }_{12}(t){\tau }_{12}(0){⟩}_{\mathrm{eq}}\stackrel{t\to \infty }{\to }V⟨{\tau }_{12}(0){⟩}_{\mathrm{eq}}^2=0.$$

$$⟨\mathcal{O}{⟩}_{\mathrm{eq}}\simeq$$

$$\dot{\pi }(x)$$

$$⟨R(x)R(x^{\prime })⟩=$$

$$\mathcal{Z}=\int \mathcal{D}\pi \mathcal{D}\varphi \exp (-H/T)=$$

$$=$$

$$⟨\pi (\mathbf{x}{)}^2{⟩}_{\mathrm{eq}}=T,⟨|\pi (\mathbf{k}){|}^2{⟩}_{\mathrm{eq}}=T.$$

$$\varphi \to \varphi /\sqrt{\alpha },\pi \to \pi /\sqrt{\alpha },\lambda \to \alpha \lambda .$$

$$T_{\lambda }=⟨{\pi }_x^2{⟩}_{\lambda }=\alpha ⟨{\pi }_x^2{⟩}_{\alpha \lambda }=\alpha T_{\alpha \lambda }.$$

$$C_{12}(t;\lambda ,T)=$$

$$\eta (\lambda ,T)=$$

$${\lambda }^2{C}_{12}(t;\lambda ,T)=$$

$${\lambda }^{2}T\eta (\lambda ,T)=$$

$${\omega }_{\mathbf{k}}=2\sqrt{{\sin }^2(k_1/2)+{\sin }^2(k_2/2)+{\sin }^2(k_3/2)}$$

$$T_{\mathrm{fit}}=1.002\pm \mathrm{0.007.}$$

$$C_{12}^{(\exp )}(t)=A\exp (-\Gamma t).$$

$$\rho (\omega )=$$

$$=$$

$$\simeq$$

$$\rho (\omega )=$$

$$C_{12}^{(\mathrm{sub})}(t)=$$

$${\rho }^{(\exp )}(\omega )=$$

$$\eta =\underset{\omega \to 0}{\lim }\frac{\rho (\omega )}{T}.$$

$$\eta =\underset{\omega \to 0}{\lim }\frac{\rho (\omega )}{T}\simeq \underset{\omega \to 0}{\lim }\frac{{\rho }^{(\exp )}(\omega )}{T}=\frac{A}{T\Gamma }.$$

$$\eta (\lambda ,T)=\frac{f_{\eta }(\lambda T)}{{\lambda }^{2}T}.$$

$$\eta /T{|}_{\lambda T=1,m^2=0}>20\times \eta /T{|}_{\lambda T=1,m^2=-0.5}.$$

$$A$$

Full text

\subjectindex

A51, D28

1 Introduction

Classical fields have been proven useful in describing various physical systems which have actually quantum nature, such as the condensate of Bose gases GP , inflation in the early universe inflation , nuclear structure Serot:1984ey , and gluon fields in heavy-ion collisions McLerran:1993ni ; Gelis . They can be used to obtain the ground states, for example, of Bose gases by using the Gross-Pitaevskii equation GP and those of nuclei in the relativistic mean field theory Serot:1984ey .

Classical field theory can be also utilized to understand the dynamical evolution of the system inflation ; McLerran:1993ni ; Gelis . In dynamical evolution of a quantum system, particles can be created by a coupling with the classical field, and its creation and their mutual collisions produce entropy and lead to eventual thermalization. Furthermore, it has been shown that the dynamics of the classical field itself is responsible for entropy production in a semiclassical approximation and thus play a significant role in the realization of an approximate local thermal equilibrium so that the hydrodynamical description of the time-evolution afterward becomes valid in the case of Yang-Mills field CYMtherm .

The classical field dynamics and its possible relevance to entropy production is considered as one of the key ingredients for a hydrodynamical description of a quantum system, then it would be intriguing to explore the properties of nonequilibrium processes as well as equilibrium statistical physics of classical fields. As far as we are aware of, it seems that systematic investigations of nonequilibrium properties such as the transport coefficients of classical fields are limited Homor ; Aarts ; Aarts:2001yn . Taking account of the possible role of classical fields in thermalization, it is natural to expect that the classical field has its own viscosity Homor . The Green-Kubo formula, shown later in Eq. (5), tells us that the shear viscosity can be obtained from the integral of the time-correlation function of the energy-momentum tensor, $C_{12}(t)$ to $t=\infty$. Since the Green-Kubo formula is obtained from the linear response theory and should be valid independent of the theoretical framework, it is possible to obtain the shear viscosity provided by classical fields from the thermal average of $C_{12}(t)$.

In this work, we discuss the shear viscosity of the classical field in the massless $\phi^{4}$ theory using the Green-Kubo formula. In order to make thermalized initial states, we utilize the Langevin equation composed of streaming, diffusion and stochastic force (latter two are collectively called the Langevin terms) to generate configurations in thermal equilibrium; the achievement of thermal equilibrium is checked by examining the equipartition of energy. Once the thermal equilibrium is realized, the Langevin terms are switched off and the time evolution of the system is solely governed by the intrinsic dynamics of the $\phi^{4}$ theory. Then we calculate $C_{12}(t)$ and obtain the shear viscosity of the $\phi^{4}$ field. It should be noted that the $\phi^{4}$ theory has the scaling property, which leads to several interesting consequences. For example, $C_{12}(t)$ and the shear viscosity are represented by some scaling functions as shown later.

It is found that $C_{12}(t)$ in massless $\phi^{4}$ theory has a slow exponential decay part, which is not found in the massive case Homor . Due to the existence of such a slow exponential decay, the shear viscosity of the massless classical scalar field is much larger than that in the massive case Homor . In reality thermal fluctuations induce a thermal mass also in the massless $\phi^{4}$ theory with vanishing bare mass, while the thermal mass squared is proportional to the coupling and small in the weak coupling region, where the exponential tail of $C_{12}(t)$ is found to grow. Furthermore, the massless $\phi^{4}$ theory would be more relevant to the classical Yang-Mills theory, which also has vanishing bare mass and the scaling property, and has been utilized to describe the initial stage of high-energy heavy-ion collisions. The study of the classical massless $\phi^{4}$ theory actually provides insight into the classical Yang-Mills field dynamics fluctuation .

One of the purposes of the present work is a detailed analysis of the time correlation of the classical energy-momentum tensor $C_{12}(t)$, which turns out to exhibit rich phenomena beyond naïve expectation such as a simple damped oscillation, especially in the case of the massless $\phi^{4}$ theory. In fact, such an analysis is necessary to extract the viscosity in a physically meaningful way. It contains two main components; a fast damped oscillation seen in the early stage and a slow monotonous decay with an exponential form. We also make a detailed analysis of the spectral function which contain interesting physics of the excitation modes of the system. The decay rates of these two components are different by two or three orders in magnitude, and the slow decay part dominates the integral of $C_{12}(t)$. We show that subtraction of the exponential term is necessary to make the integration to $t=\infty$ in the Green-Kubo formula well-defined. The shear viscosity is then obtained from the integral of $C_{12}(t)$, and the scaling function for the shear viscosity is determined.

One may recall that there is a subtlety in discussing thermalization and equilibrium using classical fields. Classical field theory can describe the low momentum modes of quantum systems. In the equilibrium of classical fields, an approximate equipartition of energy holds and each mode has the energy comparable to the temperature. Thus the distribution of massless modes obeys the the Rayleigh-Jeans law, which agrees with the Bose-Einstein distribution at low momenta, and the energy density of the classical field is divergent in the continuum limit. There are several ways to avoid the Rayleigh-Jeans divergence. One of the ways is to introduce the counter terms in the action. The renormalizability of classical scalar fields was studied Aarts , where it is found that the static thermal expectation values in the classical scalar theory are renormalizable and the classical fields show surprisingly slow thermalization. However, the behavior of higher-dimensional observables such as $C_{12}(t)$ is dominated by hard particles of order temperature and is sensitive to the ultraviolet cutoff in the classical field theory. A sophisticated way to circumvent this problem is to take account of the coupling with particles Aarts:2001yn ; Hatta:2011ky , to include quantum statistical modification of the classical variables Ohnishi:1994wn , or to introduce stochastic contribution from the hard modes Greiner:1996dx as expected from the Mori-Zwanzig projection operator formalism MoriZwanzig . It is, however, beyond the scope of the present work to treat the classical field dynamics as an infrared effective dynamics of a quantum field theory. Instead we shall confine ourselves to investigating the transport properties within the classical field theory by just introducing a finite lattice spacing. Nevertheless we expect that our study will constitute a basic ingredient for an understanding of the results given in full quantum simulations in future.

This paper is organized as follows. In Sec. 2, we first explain the lattice formulation of the classical field in the $\phi^{4}$ theory, the classical field ensemble, and the Green-Kubo formula for the shear viscosity. We also discuss the Langevin terms, which promote thermalization, and the equipartition nature in the equilibrium of the classical fields. In Sec. 3, we show the numerical results of the time-correlation function of the energy-momentum tensor and its Fourier spectrum. The shear viscosity is obtained from the $\omega\to 0$ limit of the Fourier spectrum. We also discuss the coupling dependence of the shear viscosity. In Sec. 4, we summarize our work.

2 Classical Scalar Field and Its Shear Viscosity

2.1 Classical Scalar Field Theory on Lattice

We consider the $\phi^{4}$ theory, where the Lagrangian is given as

[TABLE]

On a $L^{3}$ lattice, the Hamiltonian is given by

[TABLE]

where ($\phi(x),\pi(x)=\dot{\phi}(x)$) are the canonical variables and $\partial_{i}$ denotes the forward difference operator in the $i$-th direction, i.e. $\partial_{i}\phi(x)=\phi(x+\hat{e}^{i})-\phi(x)$. We take all quantities normalized by the lattice spacing $a$ throughout this article.

We discuss the evolution of the scalar field in the classical field approximation, where the field variables $(\phi(x),\pi(x))$ are regarded as $c$-numbers. Instead of taking account of fluctuations of the field operators in the evolution in each classical field configuration, statistical fluctuations are taken into account by using the classical field ensemble as described later. Then the above Hamiltonian is regarded as a classical Hamiltonian, which gives the classical equation of motion,

[TABLE]

We define the off-diagonal matrix elements of the energy-momentum tensor of scalar fields on the lattice as

[TABLE]

Here we adopt the central difference, $\partial^{\rm c}_{i}\phi(x)=\left[\phi(x+\hat{e}^{i})-\phi(x-\hat{e}^{i})\right]/2$. With this prescription, $\partial^{\rm c}_{i}\phi(x)$ is located in the same space-time point as $\phi(x)$ and the matrix element becomes more symmetric in the space directions. Then the central difference $\partial^{\rm c}_{i}$ is found to give a better definition of the energy-momentum tensor than the forward difference operator $\partial_{i}$.

2.2 Green-Kubo Relation

We evaluate the shear viscosity of the classical field by using the Green-Kubo formula which is a powerful relation in the linear response theory zubarev1996statistical .

[TABLE]

where $T$ and $V$ denote temperature and volume, respectively, $\tau_{12}(t)$ is the space-averaged off-diagonal matrix element of the energy-momentum tensor, and $\langle\cdots\rangle_{\mathrm{eq}}$ represents the expectation value in equilibrium.

In the long time-difference limit $t\to\infty$, the correlation disappear and $C_{12}(t)$ is expected to vanish.

[TABLE]

The final equality is expected from the isotropy of thermal systems. Thus when the correlation disappears faster than $1/t$, the integral in Eq. (5) converges to a finite value. In actual calculations, the correlation has a long-time tail and it is not easy to obtain the converged integral. As demonstrated later, we find that it is possible to evaluate the integral to infinity safely by subtracting the slowly decreasing exponential function and adding the exponential contribution separately.

2.3 Equilibrium in Classical Field Ensemble

In the Green-Kubo formula, we need the equilibrium expectation value of the time-correlation function of the energy-momentum tensor. We here obtain the equilibrium expectation value by using configurations in the classical field ensemble. We prepare $N_{\mathrm{conf}}$ initial classical field configurations in equilibrium. Then we regard the average of observables in this ensemble as the expectation values in classical equilibrium,

[TABLE]

where $\langle\mathcal{O}\rangle_{\rm eq}$ is the expectation value of the observable $\hat{\mathcal{O}}$, $N_{\mathrm{conf}}$ is the number of classical field configurations, and $\mathcal{O}_{i}$ is the observable in the $i$-th configuration.

In order to obtain many equilibrium classical field configurations efficiently, we introduce the artificial Langevin terms in the thermalization stage. In the $\phi^{4}$ theory, it takes a long time for the system to equilibrate especially at low $T$ or $\lambda$, then a large computational resource is required to prepare a large number of equilibrium configurations. By introducing the Langevin terms, thermalization is hastened. We add the Langevin terms, combination of the diffusion and stochastic forces, to the right hand side of the equation of motion for $\pi(x)$,

[TABLE]

where $\gamma$ is the diffusion coefficient and $R(x)$ represents the fluctuation, which is assumed to obey a Gaussian distribution function, namely the white noise. The strength of the fluctuation is determined by the Einstein relation,

[TABLE]

where $T$ is the temperature.

When equilibrium is reached, classical field configurations are distributed according to the Boltzmann weight. The partition function of the classical field on the lattice is given as

[TABLE]

where $\bm{k}$ denotes the lattice momentum $\bm{k}=2\pi(n_{x},n_{y},n_{z})/L$ ($n_{i}=0,1,\cdots L-1$), $\pi(\bm{k})=\sum_{\bm{x}}\pi(\bm{x})e^{-i\bm{k}\cdot\bm{x}}/\sqrt{L^{3}}$ is the Fourier transform of $\pi(\bm{x})$, and $H_{\phi}$ is a part of the Hamiltonian containing only $\phi$.

Note that $\pi(\bm{k})$ is complex in general and $\pi(\bm{k})$ and $\pi(-\bm{k})$ satisfies $\pi(-\bm{k})=\pi(\bm{k})^{*}$, then the integral should be regarded to be $d\pi(\bm{k})d\pi(-\bm{k})=2d\,\mathrm{Re}\pi(\bm{k})\,d\,\mathrm{Im}\pi(\bm{k})$. Only in those cases with all the momentum components being [math] or $\pi$, $\pi(\bm{k})$ is real.

Since the Hamiltonian is quadratic in $\pi(\bm{x})$ as well as in $\pi(\bm{k})$, the expectation values in equilibrium should be given as

[TABLE]

Thus we can measure the temperature of classical equilibrium system from the expectation values of $\pi(\bm{x})^{2}$ or $|\pi(\bm{k})|^{2}$, where the expectation value does not depend on the position and the momentum. In section 3.2, we check this relation.

2.4 Scaling property

In actual calculations, we concentrate on the time evolution of classical fields, where equation of motion has the scaling property: The equation of motion, Eqs. (3), is invariant under the transformation,

[TABLE]

Then the phase space trajectories are the same with this transformation, which is referred to as the $\lambda$ deformation hereafter.

The $\lambda$ deformation results in $H_{\lambda}\to H_{\alpha\lambda}/\alpha$, with $H_{\lambda}$ being the Hamiltonian with the coupling $\lambda$, and the temperature defined above can be regarded as

[TABLE]

By using the $\lambda$ deformation, the time-correlation function and the shear viscosity at a given $(\lambda,T)$ are found to be related with those having the same $\lambda T$ value as,

[TABLE]

Then it is enough to vary $\lambda$ and to keep $T$ fixed in order to evaluate $C_{12}(t)$ and $\eta$ at various combinations of $\lambda$ and $T$. These relations tantamount to the existence of the scaling functions which determine $C_{12}(t;\lambda,T)$ and $\eta(\lambda,T)$,

[TABLE]

It should be remembered that these scaling functions depend also on $m^{2}$.

3 Numerical results

3.1 Calculation Setup

In actual calculations, we concentrate on the case with $m=0$. The coupling is chosen to be $\lambda=0.5,1,3,5,10$ and $30$. The calculations are performed with the three different lattice sizes $L=16,32$ and $64$, and we shall mainly show the results with $L=32$, although significant/insignificant lattice-size dependence will be discussed in some important cases. Equilibrium classical field configurations are generated by introducing the Langevin terms with $\gamma=0.1$. Time evolution with the Langevin terms is calculated with a fixed lattice temperature $T=1$ for a duration of $t_{\mathrm{eq}}=1/\gamma$, where we have confirmed that classical equilibrium is reached. After preparing $N_{\mathrm{conf}}=1000$ classical field equilibrium configurations for each value of $\lambda$, we perform the time evolution calculation without the Langevin terms until $t=10^{4}/\lambda^{2}$. The time evolution is calculated in the leap-flog scheme with the time step of $\Delta t=0.01$.

3.2 Examination of Equilibrium

We first check that classical equilibrium is reached correctly by using the classical equilibrium relation in Eq. (12) at various frequencies. Thermalization process with the Langevin terms is discussed in Appendices A and B.

In Fig. 1, we show $\langle|\pi(t,\bm{k})|^{2}\rangle_{\rm eq}$ as a function of the frequency

[TABLE]

on a $32^{3}$ lattice after the equilibration process with the Langevin terms. We find that $\langle|\pi(t,\bm{k})|^{2}\rangle_{\rm eq}$ is almost constant and is close to the given temperature, $T=1$. By fitting a constant $T_{\mathrm{fit}}$ to $\langle|\pi(t,\bm{k})|^{2}\rangle_{\rm eq}$, we obtain,

[TABLE]

Thus the classical equilibrium is confirmed to be reached, where the measured temperature is consistent with the target one within the relative error of 0.2 %.

3.3 Time-correlation function of energy-momentum tensor

We shall now discuss the time-correlation function $C_{12}(t)$ of the energy-momentum tensor in the $\phi^{4}$ theory. In Fig. 2, we show the time-correlation function of the energy-momentum tensor for several values of $\lambda$ in the short time period, $t<5$. In the time region of $t<5$, the correlation function seems to show damped oscillation around a constant value. The oscillation period is around $\tau_{\mathrm{short}}=1.2-1.4$. Then the frequency is evaluated as $\omega_{\mathrm{short}}=2\pi/\tau_{\mathrm{short}}\sim(4.5-5.2)$.

In contrast, longer time results show that the background decreases with a larger rate for larger $\lambda$, as shown in Fig. 3. We find that the correlation is found to oscillate around a slow monotonous decay. This slow decay may be approximated by an exponential function,

[TABLE]

We show the fitted results using one exponential function by black dotted lines in Fig. 3, whose parameters are summarized in Table 1, and their $\lambda$ dependence is shown in Fig. 4. The fitting parameter $A$ is not very sensitive to $\lambda$, while the $\lambda$ dependence of $\Gamma$ is strong. When we add a constant or another exponential term, $\chi^{2}$ is slightly improved but the correlation among parameters makes parameter uncertainties much larger. Since $C_{12}(t)$ should approach zero in the long-time limit, $t\to\infty$, we assume that the decay is assumed to be given by the single-exponential function in the later discussions.

3.4 Fourier spectrum

Next, let us discuss the Fourier spectrum $\rho(\omega)$ of the time-correlation function of the energy-momentum tensor, $C_{12}(t)$,

[TABLE]

where $\Delta t=2t_{\mathrm{max}}/N$. The last line in Eq. (22) is nothing but the discrete Fourier transformation, and we obtain the value of $\rho(\omega)$ at $\omega=n\pi/t_{\mathrm{max}}$. The zero frequency limit $\lim_{\omega\to 0}\rho(\omega)$ is equal to the integral of $C_{12}(t)$ upto $t=\infty$, and is related with the shear viscosity as given in Eq. (26). The relation between the Fourier spectrum and the spectral density is described in Appendix C.

Integration up to $t=\infty$ turns out to be rather involved because of the exponential but slow decay of $C_{12}(t)$ at large $t$. As an example, we show the integral of $C_{12}(t)$ to $t$ in Fig. 5. At $\lambda=1$, the decay rate $\Gamma\simeq 5\times 10^{-4}$ is small, and the integral is not yet saturated at $\lambda t=1000$ as shown by the lower dashed curve. At $\lambda=10$, the integral is convergent but is still increasing at $\lambda t=1000$. In order to avoid these problems, we first subtract the exponential decay part from $C_{12}(t)$, perform the integration, and add the Fourier transformation of the exponential decay part, which is identically equal to $A/\Gamma$. With this separation, the integral to $t\to\infty$ is obtained efficiently for not only large but also small $t$ as shown by the open circles in Fig. 5.

The Fourier transformation is performed in a similar manner. We make the Fourier transformation of the subtracted time-correlation function $C_{12}^{\mathrm{(sub)}}(t)$, and the Fourier transformation of the exponential decay part is added after the discrete Fourier transformation in the third line of Eq. (22),

[TABLE]

The contribution of the exponential function Eq. (21) to the Fourier spectrum takes a form of a Lorentzian,

[TABLE]

In Fig. 6, we show the Fourier spectra at $\lambda=1,3$ and $10$ on a $32^{3}$ lattice. We find two structures, a peak at $\omega=0$ and a bump at $\omega\simeq 4$. At small $\omega$, the Fourier spectrum is dominated by the Lorentzian function from the slow decay part, shown by the dotted curves in the upper panel of Fig. 6. With decreasing $\lambda$, the decay rate $\Gamma$ becomes smaller and hence the peak becomes sharper. The peak height $A/\Gamma$ at the zero frequency has milder $\lambda$ dependence.

At high frequencies, a broad bump appears around $\omega=(2-6)$. As already mentioned, the bump around $\omega=4$ is attributed to the oscillatory behavior in the early stage. Since this oscillatory behavior exists irrespective of the magnitude of $\lambda$, a remnant of the oscillation may partly constitute the high frequency bumps. We fit the oscillation curves observed in $C_{12}(t)$ in the early stage by three damped oscillators, and the Fourier transform of them are shown by dotted curves in the lower panel of Fig. 6. As a reference, the spectra of the original $C_{12}(t)$ are shown by solid curves, which coincide well with the dotted curves, the spectra of three damped oscillators.

This broad bump also appears in the free case, $\lambda=0$. Then the frequency of $C_{12}(t)$ should be given as the sum of frequencies of two momentum modes, $\omega=\omega_{\bm{k}_{1}}+\omega_{\bm{k}_{2}}$, in the free case. The single mode frequency is in the range of $0\leq\omega_{\bm{k}}\leq 2\sqrt{3}$, then the sum of two frequencies would be in the range $0\leq\omega\leq 4\sqrt{3}\simeq 7$ without the interaction, which agrees with the spread of the bump. The frequency difference among various two-momentum modes causes the decoherence and damping of the time-correlation function in the early stage. In Appendices A and B, we demonstrate that the two-momentum modes explain the damped oscillatory behavior in the early stage in free scalar theory with Langevin terms.

3.5 Shear viscosity

By using the Fourier spectrum $\rho(\omega)$, the shear viscosity is obtained from the low-frequency limit,

[TABLE]

As already discussed, the Fourier transform of the exponential decay part of $C_{12}(t)$ appears as a Lorentzian and dominates the Fourier spectrum at low frequencies. The contribution of the subtracted part is less than 1 % of the exponential decay part on a $32^{3}$ lattice. Therefore, it is sufficient to consider the Lorentzian part in calculating the shear viscosity with a few percent accuracy. In this approximation, the shear viscosity is found to be

[TABLE]

where $A$ and $\Gamma$ are the fitting parameters introduced in Eq. (21).

We now discuss the scaling function $f_{\eta}(\lambda T)$ given in Eq. (18). While the actual numerical calculations are performed at $T=1$, we can evaluate the shear viscosity at various $(\lambda,T)$ by using $f_{\eta}$ as,

[TABLE]

In Fig. 7, we show the scaling function $f_{\eta}(\lambda T)$ as a function of $\lambda T$ on the $16^{3},32^{3}$ and $64^{3}$ lattices. While $f_{\eta}$ measured on a $16^{3}$ lattice show decreasing behavior at small $\lambda T$, results on $32^{3}$ and $64^{3}$ lattices seem to converge to a constant at small $\lambda T$, $f_{\eta}(\lambda T)\simeq 17$, which implies that $\eta\propto 1/\lambda^{2}$ for small $\lambda T$. Incidentally, it may be interesting to note that the shear viscosity in a rarefied gas is proportional to the mean free path and hence depends on the coupling constant as $\lambda^{-2}$. At large $\lambda T$, $f_{\eta}$ increases slowly, and lattice size dependence is found to be small.

Here we compare the present results with those in the previous work Homor which examined the shear viscosity near the critical point of the second order phase transition in the classical $\phi^{4}$ theory with negative bare mass squared ($m^{2}=-0.5$). With the spontaneous symmetry breaking and the thermally induced mass, the results in Ref. Homor can be regarded as those in the massive cases except for the region close to the critical temperature. The shear viscosity in Homor is proportional to the lattice temperature and insensitive to the coupling constant. Furthermore, the obtained value of it is quite small compared to that in the present work. At $\lambda T$=1, for example, the shear viscosity divided by the lattice temperature, $\eta/T$, at $m^{2}=-0.5$ is smaller than $\eta/T$ at $m^{2}=0$ by a factor more than 20,

[TABLE]

Thus it is found that the shear viscosity of the massless classical scalar field is more sensitive to the coupling constant and significantly larger than the shear viscosity in the massive case Homor . The differences from the previous results are caused by the existence of the slow exponential-decay part in the time-correlation function, which is sensitive to the coupling constant. While one may consider that thermal fluctuations induce the thermal mass and the massless property would be lost at finite temperature, the thermal mass squared induced by thermal fluctuations is proportional to the coupling constant and is expected to be small in the weak coupling region, where the long-time tail of exponential form appears in the time-correlation function.

4 Summary

In this article, we have studied the shear viscosity of classical scalar fields on a lattice by using the classical field equation of motion, the Green-Kubo formula, and the classical field equilibrium configurations. Using the Green-Kubo formula, we can obtain the shear viscosity from the time-correlation function of energy-momentum tensor in equilibrium. We have adopted the classical field ensemble method, where the classical field equilibrium configurations are generated by evolving the system with the artificial Langevin terms. We have confirmed the equilibrium is reached from the frequency dependence of the expectation value of the canonical momentum squared, $\langle|\pi(\bm{k},t)|^{2}\rangle$.

We have found that the time-correlation function shows oscillatory behavior around a slow monotonous decay with an exponential form. The Fourier spectrum at low frequency is dominated by the slow decay, then the shear viscosity is also well evaluated by the contribution from the exponential decay.

The Fourier spectrum is found to have two components, a peak at $\omega\simeq 0$ and a bump at $\omega\simeq 4$. The latter appears also in the free field theories. By comparison, the former comes from the slow exponential decay part and suggests that there exists a long-lived excitation mode, whose pole is on the imaginary axis, $\omega=-i\Gamma$. Since the damping rate is different from that of the high frequency peak by two or three orders, this long-lived mode should have the collective nature, such as the hydrodynamic modes.

We have determined the scaling function for the shear viscosity of massless classical field, $f_{\eta}(\lambda T)=\lambda^{2}T\eta(\lambda,T)$, which is invariant under the scaling transformation. The scaling function $f_{\eta}$ is found to take almost constant values at small $\lambda T$, which implies that $\eta$ behaves as $1/\lambda^{2}$ as a function of the coupling constant and this behavior happens to be in agreement with the behavior in the rarefied gas, and increases slowly at large $\lambda T$. These characteristics of the value and $\lambda$-dependence of $\eta$ are quite different from those obtained in the previous work Homor that investigated the shear viscosity near the critical point of the second order phase transition in the classical $\phi^{4}$ theory with negative bare mass squared. The differences are due to the existence of the exponential-decay part in the time-correlation function.

While we have discussed the shear viscosity in the massless case ($m=0$), the dependence on the mass were discussed in Ref. Homor . It is interesting to examine the role of the slow decay also in finite mass and negative mass squared cases. It will be also valuable to discuss the shear viscosity of classical Yang-Mills field, which is considered to describe the initial stage of high-energy heavy-ion collisions and has the scaling property as in the case of the massless classical field theory. Another important direction to study is to include the high momentum contributions as well as the quantum effects. Equipartition property of classical fields leads to the Rayleigh-Jeans divergence and nonrenormalizability in the large cutoff limit, $a^{-1}\to\infty$, and classical treatment fails to describe high momentum (quantum) dynamics. In order to evaluate the shear viscosity of quantum systems quantitatively and to compare it with the perturbative calculation result Jeon , it is desired to take account of low and high momentum modes with a proper matching procedure in a nonperturbative way Aarts:2001yn ; Hatta:2011ky ; Ohnishi:1994wn .

Acknowledgments

This work is supported in part by the Grants-in-Aid for Scientific Research from JSPS (Nos. 15H03663, 16K05350, and 19H01898), and by the Yukawa International Program for Quark-hadron Sciences (YIPQS).

Appendix A Time correlation in free scalar theory with Langevin terms

It would be instructive to discuss the thermalization stage described by the Langevin equation. We show the time-correlation function of the energy-momentum tensor, $C_{12}(t)=V\langle\tau_{12}(t)\tau_{12}(0)\rangle$, in the free scalar theory with the Langevin terms in Fig. 8. We set the parameters of the Langevin terms as ($\gamma=0.1,T=1.0$). The red line shows $C_{12}(t)$, and the black line shows the fitting function. We find that $C_{12}(t)$ dumps exponentially while it slightly oscillates around the exponential decay part. We fit $C_{12}(t)$ in the exponential function, $C_{12}(t)=Ae^{-\Gamma t}$, where fitting parameter $A$ is the amplitude of the correlation function and $\Gamma$ is the decay rate, which is the reciprocal of the relaxation time, $\tau=\Gamma^{-1}$. After the fitting, we obtain the set of fitting parameter,

[TABLE]

The decay rate is consistent with the diffusion coefficient,

[TABLE]

Then the exponential damp of the correlation function is found to come from the diffusion term added in the Langevin equation, Eq. (9).

Appendix B Analytic treatment of the time-correlation function with Langevin terms

Let us consider the time-correlation function of the free scalar theory with the Langevin terms in a model treatment. For preparation, let us recall the simplest Langevin equation,

[TABLE]

where $R$ is a white noise, $\langle R(t)R(t^{\prime})\rangle=2\gamma T\delta(t-t^{\prime})$. We can absorb the diffusion term by considering the combination, $\exp(\gamma t)p(t)$, then it is possible to integrate the equation of motion,

[TABLE]

Now let us apply the similar idea to the equation of motion of the free scalar theory together with the Langevin terms, which reduces to the Langevin equation for each momentum mode,

[TABLE]

where $R_{\bm{k}}$ is a complex white noise, $\langle R_{\bm{k}}^{*}(t)R_{{\bm{k}}^{\prime}}(t^{\prime})\rangle=2\gamma T\delta_{{\bm{k}}{\bm{k}}^{\prime}}\delta(t-t^{\prime})$. By combining these equations, we obtain the equation of one variable $X_{\bm{k}}=\pi_{\bm{k}}+i\Omega_{\bm{k}}\phi_{\bm{k}}$, as in the simple Langevin equation,

[TABLE]

We assume that $\gamma$ is small enough, $\gamma<\omega_{\bm{k}}/2$, and $\overline{\omega}_{\bm{k}}$ is real. Since the first equation (35) reads

[TABLE]

then we can integrate the right hand side and get,

[TABLE]

In order to obtain $\pi_{\bm{k}}$ and $\phi_{\bm{k}}$ separately, we invoke the two solutions, $\overline{\Omega}_{\bm{k}}=\pm\overline{\omega}_{\bm{k}}-i\gamma/2$,

[TABLE]

By using these two, we find that $\phi_{\bm{k}}(t)$ is given as the sum of the diffusion and fluctuation parts,

[TABLE]

where $\overline{\gamma}=\gamma/2\overline{\omega}_{\bm{k}}$.

The squared averages are calculated as

[TABLE]

We now confirmed that the equilibrium value $\langle\phi_{\bm{k}}^{*}\phi_{\bm{k}}\rangle_{\mathrm{eq}}=T/\omega_{\bm{k}}^{2}$ is obtained at $t\to\infty$ from the fluctuation part, Eq. (45). Then in Eq. (46) we have utilized this equilibrium value and that for $\pi_{\bm{k}}$, $\langle\pi_{\bm{k}}^{*}(0)\pi_{\bm{k}}(0)\rangle_{\mathrm{eq}}=T$.

For the later use, let us consider the average of four field product, $\langle(\phi_{\bm{k}}^{*}\phi_{\bm{k}})^{2}\rangle$. We concentrate the fluctuation part, and we find that

[TABLE]

The factor 2 comes from the two ways of contractions, and the result looks natural. Nevertheless, readers may care the product of four white noise terms. In order to be concrete, we express the white noise in a finite time step case, $\int_{t_{i}}^{t_{i}+\Delta t}dtR_{\bm{k}}(t)=\sqrt{\gamma T\Delta t}\xi_{i}$ and $\xi_{i}=\xi^{R}_{i}+i\xi^{I}_{i}$, where $\xi^{R}_{i}$ and $\xi^{I}_{i}$ are random numbers in the $i$-th time step and obey the normal Gaussian distribution,

[TABLE]

Then some of the quartic average are given as

[TABLE]

It is interesting to find that there is no additional term appearing in the case of $i=j=k=\ell$. By using these relations, the quartic average of $\phi$ is obtained as

[TABLE]

where $f(t^{\prime})=e^{-\gamma(t-t^{\prime})/2}\sin^{2}\overline{\omega}_{\bm{k}}(t-t^{\prime})/\overline{\omega}_{\bm{k}}$.

We shall now evaluate the time-correlation function. The spatial average of the energy-momentum tensor ($xy$ element) is given as

[TABLE]

The time-correlation function can be evaluated as

[TABLE]

where the connected average $C_{\phi}$ is given as

[TABLE]

where we have used the relation $\phi_{\bm{k}}^{*}(t)=\phi_{-{\bm{k}}}(t)$. The equilibrium average of four field product factorizes as in the case of the Gaussian noise in Eq. (51). Then, the contribution from the unconnected part, $\langle\phi_{\bm{k}}^{*}(t)\phi_{\bm{k}}(t)\rangle_{\mathrm{eq}}\langle\phi^{*}_{{\bm{k}}^{\prime}}(0)\phi_{{\bm{k}}^{\prime}}(0)\rangle_{\mathrm{eq}}$, in the third line of Eq. (54) disappears, since the equilibrium distribution is isotropic and the equilibrium average of $\tau_{12}$ in Eq. (53) vanishes. Then, by using the relation $\phi_{-{\bm{k}}}(t)=\phi_{\bm{k}}^{*}(t)$, the contributions from ${\bm{k}}^{\prime}=-{\bm{k}}$ as well as those from ${\bm{k}}^{\prime}={\bm{k}}$ are found to survive, and we obtain the last line of Eq. (54).

We find that only the diffusion part of $\phi_{\bm{k}}(t)$ contribute to $C_{\phi}$, since the fluctuation part is independent of $\phi_{\bm{k}}(0)$.

[TABLE]

We have used the equilibrium average of the four field product, $\langle(\phi_{\bm{k}}^{*}\phi_{\bm{k}})^{2}\rangle_{\mathrm{eq}}=2(T/\omega_{\bm{k}}^{2})^{2}$. This value is different from that for the real variable obeying the Gaussian distribution. The time-correlation function is now found to be

[TABLE]

Let us evaluate $C_{12}(t)$ a little more in the Einstein-Debye treatment of the lattice momentum, $1/V\times\sum_{\bm{k}}\to\int^{\Lambda}d{\bm{k}}/(2\pi)^{3}$. We here ignore the diffusion effect on the oscillating frequency, $\overline{\omega}_{\bm{k}}\to\omega_{\bm{k}}\simeq k$. Then the time-correlation function is evaluated as

[TABLE]

where $x=2\Lambda t$ and $\gamma_{c}=\gamma/2\Lambda$.

In Fig. 9, we show $C_{12}(t)$ in free scalar theory with Langevin terms. The qualitative behavior in the numerical results is well explained. We adopt the cutoff $\Lambda=(6\pi^{2})^{1/3}$, which gives $\int^{\Lambda}d{\bm{k}}/(2\pi)^{3}=1$. We also show the results using the lattice momentum in Eq. (57). We consider two cases of the lattice momentum. Case (I): We define the lattice momentum from the forward difference, $k_{i}=2\sin(\pi n_{i}/L)(n_{i}=0,1,\ldots L-1)$. Case (II): We use $k_{i}=\sin(2\pi n_{i}/L)$ expected from the central difference. The latter results show smaller values of $C_{12}(t)$, and roughly agrees with the numerical calculation also in strength.

The agreement of the numerical and analytical results of $C_{12}(t)$ justifies the present treatment, taking the ensemble average of classically evolved field configurations, to obtain thermal expectation values including the shear viscosity in the classical scalar field theory. In the free scalar theory with Langevin terms, $C_{12}(t)$ is found to show the exponential decay in the late-time region as in the case with interaction. This behavior suggests that the interaction term in the $\phi^{4}$ theory may be approximately taken into account by the Langevin terms in the long time region. It is also suggested that the damped oscillatory behavior of $C_{12}(t)$ in the early stage and the broad bump in the spectral function can be understood as the contribution of the two-momentum modes, the contribution from the pile-up of the modes having ${\bm{k}}$ and $-{\bm{k}}$. The frequency of the two-momentum mode depends on ${\bm{k}}$, then various two-momentum modes causes the decoherence and damping of the time-correlation function in the early stage and the broad bump in the spectral function.

Appendix C Spectral density

We here briefly mention the relation of the Fourier spectrum and the spectral density. The spectral density $\rho_{\mathrm{den}}(\omega)$ is defined as the spectral function (Fourier spectrum) of the normal Green’s function, $\langle\left[\tau_{12}(t),\tau_{12}(0)\right]\rangle_{\rm eq}$. In general, the Fourier spectrum of the product and the spectral density are related for $\omega\neq 0$,

[TABLE]

In the low frequency limit, this relation goes to the simple relation,

[TABLE]

Thus we can obtain the shear viscosity from the spectral density in the low frequency limit, $\omega\to 0$ kapusta2006finite ,

[TABLE]

In the main part of this article, we have used the Fourier spectrum $\rho(\omega)=VI(\omega)$ instead of the spectral density $\rho_{\mathrm{den}}(\omega)$ for simplicity.