Spectral function for overoccupied gluodynamics from classical lattice simulations

TL;DR

This paper investigates the spectral characteristics of overoccupied gluonic systems far from equilibrium using classical lattice simulations, focusing on dispersion relations and damping rates of gluonic excitations.

Contribution

It introduces a novel approach combining classical Yang-Mills simulations with linear response theory to analyze spectral functions in non-equilibrium gluodynamics.

Findings

Measured dispersion relations for gluonic excitations.

Determined damping rates for transverse and longitudinal modes.

Identified additional structures in the spectral function.

Abstract

We study the spectral properties of an overoccupied gluonic system far from equilibrium. Using classical Yang-Mills simulations and linear response theory, we determine the statistical and spectral functions. We measure dispersion relations and damping rates of transversally and longitudinally polarized excitations in the gluonic plasma, and also study further structures in the spectral function.

Equations

$$F_{jk}^{ab}(x,x^{\prime })=\frac{1}{2}⟨\left\{{\hat{A}}_j^a(x),{\hat{A}}_k^b(x^{\prime })\right\}⟩.$$

$$F_{jk}^{ab}(x,x^{\prime })={⟨A_j^b(x)A_k^b(x^{\prime })⟩}_{\mathrm{cl}}.$$

$${\rho }_{jk}^{ab}(x,x^{\prime })=i⟨\left[{\hat{A}}_j^a(x),{\hat{A}}_k^b(x^{\prime })\right]⟩.$$

$$G_{R,jk}^{ab}(t,t^{\prime },p)=\theta (t-t^{\prime }){\rho }_{jk}^{ab}(t,t^{\prime },p).$$

$${\hat{A}}_i^a(x)\to {\hat{A}}_i^a(x)+{\hat{a}}_i^a(x),$$

$$⟨{\hat{a}}_i^b(x)⟩=\int {\text{d}}^4{x}^{\prime }{G}_{R,ik}^{bc}(x,x^{\prime })j_c^k(x^{\prime }).$$

$$f(p)\sim \frac{n_0}{g^2}\theta (p_0-p).$$

$$f(p_{\max })\sim t^{-4/7}.$$

$$f(t,p)=t^{-4/7}{f}_{\mathrm{s}}(p/t^{1/7})$$

$$m^2=2N_{\mathrm{c}}\int \frac{{\text{d}}^3\mathbf{p}}{(2\pi {)}^3}\frac{f(p)}{p}.$$

$$m\sim t^{-1/7}$$

$$\frac{{\omega }_{\mathrm{pl}}}{m_{\infty }}=0.96,$$

$$\frac{{\omega }_{\mathrm{pl}}}{m_{\infty }}=\sqrt{2/3}\approx 0.82,$$

$$f(p)\approx \frac{T_{\ast }}{\omega (p)}.$$

$$T_{\ast }\equiv \frac{\frac{1}{2}{\int }_{\mathbf{p}}f(t,p)\left(f(t,p)+1\right)}{{\int }_{\mathbf{p}}\frac{f(t,p)}{\sqrt{m^2+p^2}}}\sim t^{-3/7},$$

$$\kappa \equiv {\int }_0^{\infty }\text{d}t\kappa (t),\text{where }\kappa (t)\propto \text{Tr}⟨E^i(0,\mathbf{x})E^i(t,\mathbf{x})⟩.$$

Full text

\PACS

24.85.+p, 25.75.-q, 12.38.Mh

1 Introduction

The main purpose of the program of ultrarelativistic heavy ion collisions is to create and study the properties of deconfined QCD matter in the laboratory. Our purpose here is to study strongly overoccupied color field configurations that are relevant for several of the different aspects of the collision process. In the very initial pre-equilibrium stage after the collision the dynamics is dominated by gluon saturation. The characteristic aspect of this regime is the existence of a semihard dominant transverse momentum scale, the saturation scale $Q_{\mathrm{s}}\gg\Lambda_{\mathrm{QCD}}$, generated by nonlinear interactions of the dense gluonic system. At the saturation scale the gluon field is nonperturbatively strong. This means that the gluon field strengths, or occupation numbers of gluonic states, are parametrically proportional to the inverse of the coupling $A_{\mu}A_{\mu}\sim 1/\alpha_{\mathrm{s}}$. Later in the evolution it is generally believed that the plasma reaches a state close to local thermal equilibrium. In such a thermal system a part of the degrees of freedom, namely the soft fields with momenta $p\lesssim gT$ are, similarly, nonperturbatively large. These soft fields are important for many properties of QCD matter. Thus in both cases we are in a situation where we want to understand the real time behavior of QCD systems with both a perturbative momentum scale (generically denoted here by $Q\gg\Lambda_{\mathrm{QCD}}$) and therefore weak coupling constant $\alpha_{\mathrm{s}}\ll 1$, but also gluon fields (at least for some important momentum modes) in an overoccupied state. In this regime the approximation of classical fields provides a powerful nonperturbative tool.

The standard method for understanding real time QCD dynamics in weak coupling is provided by the hard (thermal) loop (HTL) approach. Here one develops a perturbation theory based on a separation of two different momentum scales. The degrees of freedom at the hard scale $Q$ can be thought of as (classical) particles, and interact with soft ($\sim m_{D}$) modes that can be thought of as (classical) fields. In an equilibrium plasma the small coupling constant provides such a scale separation, but this approach can also be generalized to nonequilibrium systems. In addition to analytical calculations, there are also many numerical implementations (see e.g. [1, 2, 3]) of this idea, based on explicitly different descriptions of particle- and field-like degrees of freedom. Such calculations have been used to understand e.g. sphaleron transitions in thermal systems and plasma instabilities in anisotropic ones.

In a heavy ion collision the system starts from a configuration, at very early times $\tau\sim 1/Q_{\mathrm{s}}$, where there is only one scale $Q\sim m_{D}\sim Q_{\mathrm{s}}$. Equilibration is, then, the process where the two scales develop through various stages to become parametrically different, $m_{D}\ll T$. Following such a scenario, especially at its earliest stages, is problematic with a qualitatively different description of hard and soft modes. In practical terms, if the soft modes are described as fields on a lattice with lattice spacing $a$, the lattice needs to be fine enough to represent them: $m_{D}\ll 1/a$. But a particle-like physical picture of the hard modes implies that they are localized at the scale of the lattice spacing, requiring $a\gg 1/Q$ since the de Broglie wavelength of the hard particles is $\sim 1/Q$. In such a description it is therefore not possible to have $m_{D}\sim Q$.

The approach that we are following here avoids this problem. We treat all degrees of freedom on the same lattice, subject to the same lattice UV cutoff $1/a$. Thus we do not need to have a large scale separation. On the other hand, classical lattice simulations being relatively inexpensive, one can fit in a rather large separation between the hard and soft scales, even while maintaining a controllable continuum limit $Q\ll 1/a$. We treat even the hard modes as classical fields, not particles. In the HTL limit the classical treatment of the hard scales should not matter for the hard+soft interactions, since the only important thing about the hard modes are their color currents, not whether these currents are made of particles or fields. A drawback of the classical lattice approach is that the interactions between hard modes are treated incorrectly (as classical fields instead of particles as they should), although in the overoccupied regime the error is of higher order in $\alpha_{\mathrm{s}}$. Thus our system would, ultimately, thermalize to an unphysical classical equilibrium. The interactions between the hard modes are, however, much slower, and often neglected in the particle+field simulations in any case. We refer the reader to e.g. Refs. [4, 5, 6] for a discussion on the validity of the classical approximation. One can see our calculational setup as an extension of HTL setup to situations where the scale separation $m_{D}/Q$ can be varied smoothly up to large values.

This talk reviews the recent results presented in more detail in Ref. [7]. In the following we shall briefly describe our numerical method of linearized fluctuations, based on the algorithm developed in Ref. [8] and the obervables that we measure. We then introduce our test case system, the isotropic self-similar UV cascade of gluons. We will then review some of our numerical results so far and discuss interesting prospects for future projects.

2 Methods

For the time evolution of Classical Yang-Mills fields we use the standard Hamiltonian lattice [9, 10] formulation of gauge theory in real time. Here, instead of the gauge potential $A_{i}$ and the covariant derivative $D_{i}=\partial_{i}+ig[A_{i},\cdot]$, one works with link matrices connecting lattice sites $U_{i}(x)=e^{iagA_{i}(x)}$ and covariant finite differences obtained using the links. The canonical conjugate variable to the gauge potential is the chromoelectric field $E^{i}=\partial_{t}A_{i}$. The Hamiltonian setup is formulated in the temporal gauge $A_{0}=0$, where one must take care to satisfy also the Gauss’ law constraint $[D_{i},E^{i}]=0$.

Our first measurable is the “statistical function,” defined as a symmetric two point function of the classical field

[TABLE]

The statistical function measures “thermal” fluctuations in the field, and is related to the phase space density of particles in system $f(p)$. In the classical approximation the commutator is suppressed by a power of $g$ (because $A_{i}\sim 1/g$ ). Thus here $F$ is just a two point function of the classical field

[TABLE]

The statistical function at equal time is an often measured quantity in CYM calculations, but in Ref. [7] we also measure it at unequal time.

In addition to the statistical function, the other in general independent correlator in quantum field theory is the “spectral function”

[TABLE]

It is a genuinely “quantum” obervable and proportional to $\sim\hbar$ due to the the field operators commutation relations. We can nevertheless measure it in the CYM simulations by using its relation to the retarded propagator as

[TABLE]

The retarded propagator measures the linear response of the system to an external perturbation, and we measure it with the algorithm developed in Ref. [8]. We split the gauge field into a background field and a fluctuation:

[TABLE]

where the fluctuation generated by an external infinitesimal current is given bt the retarded propagator

[TABLE]

This leads to our algorithm to calculate the statistical function. We first, at time $t_{0}$, perturb the system with a current that exists only for one timestep: $j^{k}_{c}(x)\sim e^{i\mathbf{k}\cdot\mathbf{x}}\delta(t-t_{0})$. We then follow the linearized equations of motion for the fluctuation and its time derivative $a_{i}^{a}(x)=\langle\hat{a}_{i}^{a}(x)\rangle,\ e^{i}_{a}(x)$. Finally, at a later time $t>t_{0}$, we measure the correlation of the field $a_{i}^{a}(t)$ with the current $j^{i}_{a}(t_{0})$ and use this to extract the momentum space spectral function $\rho(p,t)$. A similar procedure with the electric field fluctuation can be used to obtain time derivatives of the spectral function.

3 Overoccupied cascade

Let us then move on to discuss our test case overoccupied nonequilibrium system, the isotropic self-similar cascade. This system has been extensively studied by many groups (see e.g. [11, 4, 6]), and it has been found that a kinetic theory formulation can describe the basic properties seen in the numerical studies. In this system one starts from an isotropic field configuration where only modes up to some characteristic hard momentum are occupied

[TABLE]

In practical calculations one typically replaces a strict theta function with a smoother Gaussian momentum distribution, but this detail, or the initial occupation number $n_{0}/g^{2}\sim 1/g^{2}$ or momentum scale $p_{0}$ matter little for the later self-similar time evolution of the system. This is instead fully determined by the conserved total energy density, which defines the real characteristic hard scale of the problem as $\varepsilon\sim Q^{4}/g^{2}$. In what follows we specifically define $Q=\sqrt[4]{10\pi^{3/2}g^{2}\varepsilon/\left(\sqrt{2}(N_{\mathrm{c}}^{2}-1)\right)}$, scale dimensionful quantities with $Q$ to dimensionless ones and, unless otherwise stated, plot quantities measured at $Qt=1500$.

During the time evolution of the system, the energy cascades towards UV modes in such a way that all the momentum modes up to $p_{\mathrm{max}}\sim t^{1/7}$ are occupied. The typical occupation number (e.g. at the hard scale $p_{\mathrm{max}}$) decreases with time as

[TABLE]

The scaling behavior

[TABLE]

becomes evident in the numerical result when we scale the momentum distribution of gluons by $t^{4/7}$, and plot it in terms of the scaled momentum $t^{-1/7}p$, as shown in Fig. 1 (left).

In HTL or kinetic theory one estimates the Debye or plasmon scale from the (hard) particle distribution as

[TABLE]

In the scaling regime (9) this gives an estimate

[TABLE]

for the time dependence of the soft scale. This scaling is verified numerically in Fig.1 (right). This leads us to an important aspect of this self-similar attractor system, namely that the scale separation $p_{\mathrm{max}}/m$ grows with time. By looking at different times one can therefore smoothly turn on or off the parameter that determines the validity of the HTL approximation. In the following we will be mostly working at rather late times where this scale separation is clear, in order to compare our calculation to HTL in a regime where it should be a good description of the system.

4 Spectral function and dispersion relations

We start this discussion with the statistical and spectral functions for transversely polarized quasiparticles. Relevant questions here are whether they both exhibit the same peak structure, and what are the locations and widths of the peaks. At small frequencies one would also expect to see an additional structure, the “Landau cut”. The spectral function is naturally normalized to unity (actually $\hbar$), whereas the statistical function is proportional to the number of particles in the system. In order to compare the peak structure in the two correlators it is convenient to divide the statistical function by its value at $\Delta t=t^{\prime}-t=0$, which is what we will do in the following plots. Figure 2 shows the spectral function for the transverse polarization. First of all we see that there is a very nice agreement between the two functions, which could be interpreted as a generalized fluctuation-dissipation theorem. In frequency space one can see a very nice Lorentzian shape, and extract the quasiparticle width (plasmon damping rate). One can even clearly see the Landau cut structure at small frequencies, in rough agreement with the HTL theory curve that uses only the value of $m$ extracted numerically from the data as an input parameter.

For the longitudinal polarization mode, as shown in Fig. 3 the story is very similar. There is a good agreement between the statistical and spectral functions, when the former is normalized to the equal time value. For the longitudinal polarization state the measurement is harder, since at high momenta the quasiparticle peak gets weaker and merges with the Landau cut, as could have been expected from HTL.

By extracting the location of the peak as a function of momentum, we can extract a quasiparticle dispersion relation from our data as shown in Fig. 4. Overall, within our statistics it is not possible to differentiate between the HTL functional form and a relativistic dispersion relation $\omega^{2}=m^{2}+p^{2}$. Figure 4 also shows data for the longitudinal polarization. The difference between the two polarization states is qualitatively as one would expect from HTL. Quantitatively we can characterize the dispersion relation by two different masses from the small and large momentum limits by defining the plasmon mass $\omega_{\mathrm{pl}}\equiv\omega(p\to 0)$ and the mass gap at $p\to\infty$, denoted by $m_{\infty}$. Our numerical estimate for the ratio of these scales is

[TABLE]

where HTL prediction for this quantity would be

[TABLE]

and the NLO correction has been calculated to be negative [12].

5 Further observables

Let us now discuss a little more the equal time field correlator, i.e. $f(p)$, which so far was used as a normalization factor to compare the statistical and spectral functions. In the small momentum limit the expectation from HTL theory would be to see a classical thermal distribution

[TABLE]

The effective temperature for the soft fields should be given by

[TABLE]

where for classical fields one neglects the 1 in $(f+1)$. Our results for the equal time electric field correlator are shown in Fig. 5 (left). We see that there is a significant enhancement in the infrared compared to the expectation. We do not currently have a compelling interpretation for this observation, but certainly the agreement with expectations from HTL is not as good as for the spectral function.

In the HTL theory the plasmon damping rate is also related to the effective temperature $T_{*}$. Our numerical result for the damping rate, extracted from fitting the time domain signal with an exponentially decaying function, is shown in Fig. 5 (right). Our measurements of the time dependence of $\gamma$ are consistent with the time scaling from Eq. (15). However, we are able to extract a result in a rather large range of momenta, whereas from the HTL theory we only have points at $p=0$ (shown on the figure) and $p=\infty$ [13] (not shown, but consistent within errors). Also here, the agreement of this $T_{*}$-dependent observable with HTL is not as good as for the spectral function.

6 Future prospects and conclusions

After this quick review of the first results in Ref. [7], let us briefly discuss what we see as natural next steps in this setup. Given this powerful formalism for real-time calculations the natural next step would be to use it for calculations of transport coefficients in such a nonequilibrium gluon plasma. One would be interested in e.g. the heavy quark diffusion coefficient, the shear viscosity $\eta$ or the jet quenching parameter $\hat{q}$. Probably the simplest of these is the first one. Without going into details (see e.g. [14]), for our purposes it is enough to state that this is defined as the infinite-time limit of the local unequal time electric field correlator:

[TABLE]

Here we have omitted the temporal gauge links needed to make this quantity explicitly gauge-invariant, since they become unity in temporal gauge. The correlator is in principle a straighforward quantity, but surprisingly tricky to evaluate numerically. This is due to the fact that the integrand is very strongly oscillatory, and very fine cancellations are needed to obtain the value of the integral. However, the time dependence of the integrand itself contains more physical information than just the constant. In fact, our preliminary studies strongly suggest that the time dependence contains very strong slowly oscillating modes. These could be a confirmation of the infrared enhancement seen in the gauge fixed correlator in Fig. 5. The correlator $\kappa(t)$ is manifestly gauge invariant, unlike the electric field correlators as a function of $p$ that are evaluated in Coulomb gauge. Thus the time dependence could serve as a gauge invariant confirmation of the IR enhancement.

Another natural step is to move from a 3-dimensional isotropic cascade towards the initial stage of a relativistic heavy ion collision at high energy. Here the system starts off, at leading order in the QCD coupling, as an effectively 2-dimensional boost invariant glasma [15] field. In a genuinely 2-dimensional system it seems that an HTL-type approximation cannot be applied. One way to see this is to look at the integral (10) that gives the value of the Debye mass in terms of the distribution of hard particles, e.g. in thermal equilibrium. In 3 spatial dimensions the integral is dominated by the hard modes, and thus the Debye mass is in a sense an “external” parameter for the soft fields. In 2 spatial dimensions, assuming that the distribution in the IR is close to thermal $f\sim 1/\omega$, the integral (10) receives a large, logarithmic, contribution from soft modes. Thus the Debye mass is not an external parameter anymore, but must be self-consistently determined from the soft modes. Thus it is not obvious whether a kinetic theory description exists for the 2-dimensional system, and specifically whether one would expect to see a self-similar cascade solution as a function of time. This is a situation that can perfectly well be studied in our setup of classical fields and linearized fluctuations. Concentrating at this stage on a nonexpanding 2-dimensional system (which was studied already in [16], but only at early times and with equal-time correlators), we find clear evidence for a self-similar cascade solution, whether or not this has a kinetic theory interpretation. Our preliminary numerical result for the distribution function is shown in Fig. 6.

In conclusion, we have here argued that several aspects of a heavy ion collision exhibit overoccupied $f(p)\sim 1/g^{2}$ classical gauge field configurations. This includes the initial glasma fields with the only relevant momentum scale $p\sim Q_{\mathrm{s}}$, and, later, the soft fields $p\sim gT$ in thermal system. For controlled understanding of these fields we have developed a new numerical algorithm for linearized fluctuations on top of the fully nonlinearly interacting classical gauge field. As a first test case of this setup we study an isotropic self-similar UV cascade solution. At late enough times, there exists a scale separation between the hard and soft modes. Thus one can compare the results to HTL. We have seen that the HTL approximation is rather well satisfed for the spectral function, but more uncertain for quantities related to the effective temperature of the soft modes, denoted as $T_{*}$. In particular, we are able to extract a plasmon decay rate $\gamma(p)$ as a function of momentum. In the future, we see many prospects in extending these methods to transport coefficients, with the heavy quark diffusion coefficient as a first step. We also want to study anisotropic systems relevant to the pre-equilibrium matter in heavy ion collisions, where as a first application we demostrate a self-similar scaling solution in a purely 2-dimensional gauge theory.

Acknowledgements

This work is supported by the European Research Council, grant ERC-2015-CoG-681707. T. L. has been supported by the Academy of Finland, projects No. 267321 and No. 303756. J. P. has been supported by the Jenny and Antti Wihuri Foundation and a travel grant from the Magnus Ehrnrooth foundation. K. B. and J. P. would like to thank the CERN Theory group for hospitality during part of this work. We acknowledge computational resources from CSC – IT Center for Science, Finland.