Spectral function of fermions in a highly occupied non-Abelian plasma

This is a self-archived version of an original article. This version may differ from the original in pagination and typographic details. Author(s): Title: Year: Version: Copyright: Rights: Rights url: Please cite the original version: CC BY 4.0 https://creativecommons.org/licenses/by/4.0/ Spectral function of fermions in a highly occupied non-Abelian plasma © 2022 The Author(s). Published by Elsevier B.V. Published version Boguslavski, K.; Lappi, T.; Mace, M.; Schlichting, S. Boguslavski, K., Lappi, T., Mace, M., & Schlichting, S. (2022). Spectral function of fermions in a highly occupied non-Abelian plasma. Physics Letters B, 827, Article 136963. https://doi.org/10.1016/j.physletb.2022.136963 2022 Physics Letters B 827 (2022) 136963 Contents lists available at ScienceDirect Physics Letters B www.elsevier.com/locate/physletb Spectral function of fermions in a highly occupied non-Abelian plasma K. Boguslavski a,∗, T. Lappi b,c, M. Mace b,c, S. Schlichting d aInstitute for Theoretical Physics, Technische Universität Wien, Wiedner Hauptstr. 8-10, 1040, Vienna, Austria bDepartment of Physics, University of Jyväskylä, P.O. Box 35, 40014, University of Jyväskylä, Finland cHelsinki Institute of Physics, P.O. Box 64, 00014, University of Helsinki, Finland dFakultät für Physik, Universität Bielefeld, 33615, Bielefeld, Germany a r t i c l e i n f o a b s t r a c t Article history: Received 10 July 2021 Received in revised form 3 February 2022 Accepted 12 February 2022 Available online 24 February 2022 Editor: J.-P. Blaizot Keywords: Quark-gluon plasma Heavy-ion collisions Nonequilibrium QFT Thermal QFT Spectral function Hard-thermal loop We develop a method to obtain fermion spectral functions non-perturbatively in a non-Abelian gauge theory with high occupation numbers of gauge fields. After recovering the free field case, we extract the spectral function of fermions in a highly occupied non-Abelian plasma close to its non-thermal fixed point, i.e., in a self-similar regime of the non-equilibrium dynamics. We find good agreement with hard loop perturbation theory for medium-induced masses, dispersion relations and quasiparticle residues. We also extract the full momentum dependence of the damping rate of the collective excitations. ©2022 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). Funded by SCOAP3. 1. Introduction Nonperturbatively strong color fields play an important part in the physics of ultrarelativistic heavy ion collisions and the early universe. The pre-equilibrium initial stages of the collision are characterized by highly occupied glasma field configurations [1–4] resulting from the collision of two dense gluonic systems. Similarly states with large occupation numbers of scalars or gauge bosons can emerge from instabilities in the reheating of the early universe [5–7]. Also in an equilibrated quark-gluon plasma, the infrared sector is dominated by gluon fields with high occupation numbers. Even if the energy density of the system is dominated by bosons, it is important to understand the interactions of fermions with these strong bosonic fields. In heavy-ion collisions, jets formed by light energetic quarks are created in the earliest stage of the collision. They propagate through the dense gluonic system, which has an effect on their energy loss [8]. Electromagnetic observables provide another experimental window into the earliest stage of heavy-ion collisions [9,10], and originate from quarks which are the only carriers of electric charge in the medium. In order to develop a microscopic description of such observables, one must understand the interactions of quarks with a system of overoccupied gluon fields. *Corresponding author. E-mail address: kirill.bogusla[email protected] (K. Boguslavski). Since the dynamics of bosonic states with occupation numbers of the order of the inverse self-coupling ∼1/αs, can be naturally described in terms of classical fields [11–13]to leading order in αs, classical statistical field simulations are commonly used in heavyion physics [14,15] and cosmology [16]. We will here use this classical field picture to study the microscopic properties of Dirac fermions (quarks) interacting with a over-occupied non-abelian gauge field. The interaction of a fermion with the background field is encoded in its spectral function. The purpose of this paper is to compute this spectral function. Based on earlier calculations of spectral functions for gluons [17,18] and the real time lattice fermion code developed in [19,20], we will first develop a numerical method to calculate the spectral function in an out-of-equilibrium overoccupied background gauge field configuration. We are studying fermion interactions with a strongly overoccupied gluon field, thus the dynamics is dominated by gluons and the physical situation is very different from systems at large baryon density (see, e.g., recent work in [21]). Our classical-statistical method is similar to the ones used to extract spectral functions in scalar theories both far from equilibrium [22–24]and for a thermal system [25–27], and can also be applied to study the dynamics of fermionic excitations in the presence of scalar or abelian gauge fields. We will compute the spectral function in momentum space, in both the time and frequency domains. From this spectral function we can extract medium-induced masses, dispersion relations https://doi.org/10.1016/j.physletb.2022.136963 0370-2693/©2022 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). Funded by SCOAP3. K. Boguslavski, T. Lappi, M. Mace et al. Physics Letters B 827 (2022) 136963 and quasiparticle residues for the different spinor structures of the spectral function. These quantities will then be compared to predictions of hard-thermal loop (HTL) perturbation theory [28–31]. We can also extract the damping rate of fermionic quasiparticles for different momenta, which is a much more nontrivial quantity to obtain in perturbation theory [28,32,33,29]. This paper is structured as follows. We first describe the numerical method for extracting the fermion spectral function in Sec. 2, and test it for the analytically solvable case of free fermions. We then move to a nontrivial background field in Sec. 3, where we first briefly describe the overoccupied universal UV-cascade gluon field configuration that we are using, and then present our numerical results and compare them to the expectation from HTL perturbation theory. We briefly conclude in Sec. 4. For completeness, the HTL formulas for the spectral function from the literature are provided in Appendix A. 2. Spectral functions from classical-statistical lattice simulations 2.1. Classical-statistical simulations We consider a non-abelian SU(Nc)gauge theory discretized on a lattice with N3 ssites and lattice spacing as. We use Nc=2in this work. The gauge fields are expressed in terms of lattice gauge links U(t, x) ≈exp igasAj(t,x)and electric field variables Ej(t, x) ≈ ga2 s∂tAj(t, x)in temporal axial gauge (A0=0), where g=√4παs denotes the gauge coupling. The evolution equations for the gauge field sector then result from the lattice Hamiltonian HYM=1 g2as x,i Tr[Ei(t,x)2]+1 2 j ReTr[1−Uij(t,x)],(1) where Uij(t, x) =Ui(t, x)Uj(t, x +ˆ ı)U† i(t, x +ˆ j)U† j(t, x)are the usual plaquette variables and ˆ ı, ˆ jdenote the unit lattice vectors in the i, j =1, 2, 3spatial directions. While the gauge fields are treated as classical fields, fermions are described in terms of quantum mechanical field operator ˆ ψ(t, x), whose evolution is governed by the Hamiltonian ˆ HW=1 2 xˆ ψ†(t,x), γ0−i/Ds+mˆ ψ(t,x)(2) in the presence of the classical background gauge fields. With regards to the lattice discretization of fermions, we follow previous works [19,20] and discretize the Hamiltonian with a tree-level improved Wilson Dirac operator −i/Dsˆ ψ(t,x)=(3) 1 2 n,i Cn−iγi−nrWU+nˆ ı(t,x)ˆ ψ(t,x+nˆ ı) +2nrWˆ ψ(t,x)−−iγi+nrWU−nˆ ı(t,x)ˆ ψ(t,x−nˆ ı). Here i =1, 2, 3is the spatial Lorentz index and rW=1is the Wilson parameter. The parallel transporters accross multiple lattice sites are given by products of individual link matrices, and are denoted by U+nˆ ı(t, x) =n−1 k=0Ui(t, x +kˆ ı)and U−nˆ ı(t, x) = n k=1U† i(t, x −kˆ ı). For a leading order tree-level improvement [19, 20]we set the coefficients Cnas C1=4/3, C2=−1/6 and Cn>2= 0. Since the equation of motion for the fermion field operator iγ0∂0ˆ ψ(t,x)=−iγjDs W,j+mˆ ψ(t,x)(4) is linear in the fermion field ˆ ψ(t, x), it can be conveniently solved in terms of a mode function expansion [34,35]. This means that we expand the operator ˆ ψ(t, x), in terms of creation and annihilation operators of particles (b) and anti-particles (d) with definite momenta pat a reference time t ˆ ψ(t,x)=1 √V λ,pˆ bλ,p(t)φu λ,p(t,x)+ˆ d† λ,p(t)φv λ,p(t,x), (5) where λ =1, ···, 2Nccollectively labels the spin and color indices. The operator structure is determined by the action of the ˆ bλ,p(t) and ˆ d† λ,p(t)at a fixed reference time t, where the creation and annihilation operators satisfy the usual equal-time anti-commutation relations ˆ bλ,p(t), ˆ b† λ,p(t)=ˆ dλ,p(t), ˆ d† λ,p(t)=Vδp,pδλλ.(6) With these fixed equal-time commutation relations, the time evolution of the field operator ˆ ψ(t, x)is, on the other hand, encoded in the set of 4NcN3 swave-functions (or N3 scolored spinors) φu/v λ,p(t, x). They describe the propagation of a state that is given by a plane wave at the reference time t, i.e., satisfying the initial condition φu λ,p(t,x)t=t=uλ(p)e+ip·x(7) φv λ,p(t,x)t=t=vλ(p)e−ip·x(8) Each of these wave-functions satisfies the Dirac equation in the classical background field (4). Due to the large phase-space occupancy of gluons, the fermionic sector is suppressed relative to the gauge fields by one power of αsin weak coupling. Working at leading order accuracy, we can therefore neglect the backreaction of fermions on the dynamical gauge fields, just as we are neglecting gluonic quantum corrections in the gluon field dynamics. Neglecting the backreaction also makes our calculation computationally significantly less demanding, as we will explain further below. Employing a leap-frog type scheme for the discretized time evolution then results in the following update rules for the gauge field and fermion sectors [20] Uj(t,x)=eiat/asEj(t−at/2,x)Uj(t−at,x), (9) Ej(t+at/2,x)−Ej(t−at/2,x)=(10) −at as j=iUij(t,x)+Ui(−j)(t,x)ah , φu/v λp(t+at,x)−φu/v λp(t−at,x)=(11) −2iatγ0−iγjDs W,j[U]+mφu/v λp(t,x), with [.]ah denoting the anti-Hermitian part of a matrix. This can be solved iteratively to calculate the time evolution.1 2.2. Spectral function of fermions Generally, the spectral function is defined as the expectation value of the unequal time anti-commutator of fermion fields ραβ(x,y)=ˆ ψα(x), ˆ ¯ ψβ(y)ψ(12) where ˆ ¯ ψ=ˆ ψ†γ0and α, βare Lorentz indices, which should not be confused with the indices λ, λdenoting the spin and color 1If not stated otherwise, we employ at/as=0.005 for the discretized time evolution. 2 K. Boguslavski, T. Lappi, M. Mace et al. Physics Letters B 827 (2022) 136963 states, that we will write as subscripts. Here .ψdenotes expectation values of fermionic operators in the presence of gauge fields, and .denotes the classical-statistical average over gauge field configurations, as is usually performed for observables in the classical-statistical framework [22]. Since the spectral function has a 4 ×4matrix structure, it is useful to decompose it into scalar (S), pseudo-scalar (P), vector (V), axial-vector (A) and tensor (T) components according to ρ=ρS+iγ5ρP+γμρμ V+γμγ5ρμ A+1 2σμνρμν T.(13) Due to rotational symmetry the vector spectral function is proportional to the momentum, and we can express its spatial components in terms of a scalar function ρVas ρi V(p)=pi Ep ρV,(14) where Epis the free dispersion relation that will be discussed in Sec. 2.3. The temporal and spatial components of the spectral function can then be extracted as ρ0 V=1 4Tr(ργ0), ρV=−Eppj 4p2Tr(ργ j). (15) On a discrete lattice, rotational symmetry is broken, and such a relation does not hold exactly. For momentum modes far from the UV cutoff that we will consider it should, however, be satisfied. On the lattice we will determine the function ρVby replacing piby the effective lattice momentum ˆ picorresponding to the discretization of the derivative operator (Sec. 2.3). By inserting the mode function expansion of the fermion fields in Eq. (5)into the definition of the spectral function Eq. (12), one obtains ραβ(x,y)= λ,pφu,α λ,p(x0,x)φu,γ λ,p(y0,y)∗(16) +φv,α λ,p(x0,x)φv,γ λ,p(y0,y)∗γγβ 0, where we used the equal-time anti-commutation relations in Eq. (6)to evaluate the anti-commutators. By changing to central and difference coordinates X =(x +y)/2 and x=(x −y)in the spatial direction, we perform a spatial average over the position X and a Fourier transform w.r.t. to difference coordinates xaccording to ραβ(x0,y0,p)=(17) 1 Vd3Xd3xe−ip(x−y)ραβ(x0,x,y0,y). Since we are studying a system where expectation values are translationally invariant, the momentum space spectral function does not depend on the central coordinate X. We thus arrive at a compact expression for the spectral function in terms of the mode functions as ραβ(x0,y0,p)=1 V λ,q˜ φu,α λ,q(x0,p)˜ φu,γ λ,q(y0,p)∗(18) +˜ φv,α λ,q(x0,p)˜ φv,γ λ,q(y0,p)∗γγβ 0, where ˜ φλ,q(x0, p) =d3xφλ,q(x0, x)e−ip·xdenotes the spatial Fourier transform of the wave-functions. We note that the wave functions ˜ φλ,q(x0, p)depend on two momentum arguments: q which is the wavenumber at the initial time t, and pwhich is the momentum where the wave function is evaluated. By choosing the reference time for the mode function expansion in Eq. (5) as t=y0, we can simplify the momentum structure of the spectral function. The initial condition, Eqs. (7) and (8), corresponds to ˜ φλ,q(t, p) ∝δ(3)(p −q)in momentum space. This can be used to evaluate the sum over the momenta q, leading to an expression for the spectral that is particularly convenient for numerical evaluation ραβ(x0,y0,p)(19) =1 V λ˜ φu,α λ,p(x0,p)u†,γ λ(p)+˜ φv,α λ,−p(x0,p)v†,γ λ(−p)γγβ 0. In general, the knowledge of the full set of 4NcN3 swavefunctions is required to construct the time evolution of the fermion field operator ˆ ψ(t, x)according to Eq. (5). However, the spectral function in Eq. (19)can be expressed in terms of the 4Nc components of a single momentum mode φu/v λ,p(t, x). Since each momentum mode pcan be computed completely independently, the calculation of the fermion spectral function is computationally significantly less demanding than simulations including the backreaction of dynamical fermions. This makes it possible to use significantly larger lattices, leading to a better resolution of different momentum scales. Our algorithm to calculate the fermion spectral function can be summarized as follows: 1. Generate a configuration of lattice gauge links Uand electric fields E, and evolve it by the classical Yang Mills equations (9), (10)up to the reference time tfrom which the spectral function is measured. 2. Select a subset of Nmodes momentum (p)modes, for which the spectral function is computed, and initialize the Nφ=4Nc× Nmodes fermion wave-functions φu/v λ,p(t)at the reference time t according to Eqs. (7), (8). 3. Solve the Dirac equation in (11)for all Nφmodes along with the classical Yang Mills equations (9), (10)to compute the time evolution for t>t. 4. Calculate the spectral function ρ(t, t, p)for the Nmodes momenta pand t>tby projecting out the appropriate plane wave component according to Eq. (19). This algorithm is completely analogous to the one for the gluon spectral function developed in [17,18]. One initializes a fluctuation in a specific momentum mode, evolves forward in time in coordinate space, and projects back to the momentum state after the evolution. In the case of gluons, this projection involves a projection to the appropriate polarization state, while for fermions one uses the free spinors u, vto project out the appropriate helicity and positive and negative energy states. Before we proceed with the calculation of fermion spectral functions based on the above algorithm, some further comments on the gauge dependence are in order. Evidently, the fermion spectral function defined in Eq. (12)is a gauge dependent quantity, whose non-perturbative calculation requires the implementation of a suitable gauge fixing procedure. While the temporal axial gauge condition A0=0is naturally implemented in the Hamiltonian lattice gauge theory formulation, this leaves the residual gauge freedom to perform time independent gauge transformations. We eliminate the residual gauge freedom by fixing Coulomb gauge ∂jAj(t, x) =0at the time t=twhen the calculation of the spectral function is initialized, i.e., between the first and the second step in the above algorithm. We note that this procedure is analogous to the linear response framework employed in Refs. [36,17,18] to extract the gluon spectral function, where similarly, one fixed 3 K. Boguslavski, T. Lappi, M. Mace et al. Physics Letters B 827 (2022) 136963 Coulomb gauge and subsequently studies the response of the gauge fields to plane wave perturbations in order to extract the spectral function. 2.3. Benchmark for free fermions We first illustrate and check the method by calculating the spectral function for free fermions. This is achieved by setting the background gauge links Ui(t, x) =1and electric fields Ei(t, y) = 0. The free fermion spectral function is simply given by [31] ρfree(ω,p)=2πsgn(ω)pμγμ+mpδ(ω2−E2 p). (20) In the lattice discretization the three-momentum part of the fourmomentum pμin this expression must correspond to the discretization of the derivatives that we are using. Thus we have pμ=(ω, ˆ p), where the effective quasi-particle momentum is ˆ pi=− n Cn as sin n2πni Ns,(21) with the discrete lattice momentum mode index ni=0, ···, Ns− 1[20]. The Wilson term generates an effective mass that makes the doubler modes more massive (and breaks chiral symmetry), so that mp=m+2rW i,n nCn as sin2nπni Ns.(22) In terms of the effective momentum ˆ pand mass mpthe energy of the single free fermion satisfies the usual relativistic dispersion relation Ep=ˆ p2+m2 p. It is convenient to express the spectral function in terms of particle and anti-particle excitations as ρfree(ω,p)=2π 2Ep+(p)δ(ω−Ep)+−(−p)δ(ω+Ep), (23) where ±(p)denote the usual projections of the Dirac components2 ±(p)=γ0Ep−γjˆ pj±mp.(24) The free spectral function ρfree(t, p)in the time domain is then obtained as ρfree(t,p)=1 2Ep+(p)e−iEpt+−(−p)eiEpt =γ0cos(Ept)+iγjˆ pj Ep−mp Epsin(Ept). (25) Comparing this to the general spinor decomposition of the spectral function in Eq. (13)we see that the only non-vanishing components are the scalar Imρfree S=−mp Epsin(Ept), the temporal part of the vector Reρ0,free V=cos(Ept)and the spatial part of the vector spectral function Imρfree V=− sin(Ept). We show our numerical results for the free spectral function in Fig. 1, where we present the time evolution of the components ImρS, Reρ0 Vand ImρVcalculated on a 643lattice with momentum asp =(0.098, 0.195, 0.29)and mass parameter mas=0.003125 corresponding to nearly massless fermions. An excellent agreement between continuous lines, depicting the analytical expressions in 2Note that the particle and anti-particle projections are given by +(p) = 1 Ncλuλ,p¯ uλ,pand −(p) =1 Ncλvλ,p¯ vλ,p. Fig. 1. Components of the spectral function ρS, ρ0 Vand ρVof free (Wilson) fermions for a fixed momentum asp =(0.098, 0.195, 0.29)as a function of t. Solid curves correspond to the analytical results in (25). Eq. (25), and points, corresponding to the numerical lattice data, is observed for all components,3validating our procedure to calculate spectral functions. 3. Nonperturbatively computed spectral functions We now turn to the investigation of quark spectral functions in a non-equilibrium plasma. We consider a highly occupied plasma of gluons, as described by the initial phase-space distribution of gluons g2fg(t=0,p)=n0 Q pe−p2 2Q2,(26) with p =|p|and where n0/g21is the initial occupancy and Qis the characteristic energy scale. Such initial conditions can be represented by a classical-statistical ensemble of fluctuating gauge fields, which we implement numerically as in Ref. [17]. Such overoccupied gluonic systems have been studied in several recent works [37–43,14,17]; they encounter a rapid memory loss about the details of the initial conditions, and subsequently experience a self-similar scaling behavior where the dynamics of the phasespace distribution g2fg(t,p)=(Qt)αfs(Qt)βp/Q),(27) can be described in terms of a scaling function fsand universal scaling exponents α=−4/7, β=−1/7[37–43,14]. Since the scaling behavior in Eq. (27)can be realized for a variety of different initial conditions [44,39,45,14,46], this non-thermal fixed point state represents a generic non-equilibrium state of a highly occupied plasma, and we will calculate the quark spectral function in this self-similar scaling regime. Here we will start from a moderate occupancy of n0=0.24and first consider quark spectral functions at a fixed reference time Qt=1500, which is well within the self-similar regime. We will then investigate the Qt dependence of the spectral functions. In the scaling solution the time dependence of the hard scale and screening scale is known, and the dependence of the fermion spectral function on the reference time Qt can be used to understand its structure in terms of these 3We have also checked that numerical results for the vanishing components vanish to machine precision of 10−16 for this test case. 4We note that in order to allow for a direct comparison, our initial conditions and our choice for the extraction time Qt are the same as in Ref. [17], where the gluon spectral function was extracted. 4 K. Boguslavski, T. Lappi, M. Mace et al. Physics Letters B 827 (2022) 136963 Fig. 2. The components of the spectral function ρ0 Vand ρVas functions of t≥0. microscopic scales of the gluon field configurations. Note that the dependence of the spectral function on the relative time t−thappens at a much shorter timescale ∼1/mgthan the dependence of the universal cascade solution on Qt, which is a consequence of the self-similar dynamics. Thus measurements of the spectral function at different Qt effectively study different quasi-static systems characterized by different scale separations between the hard and soft scales. If not stated otherwise, our simulations are performed on Ns=256 lattices for nearly massless fermions m =0.003125 Q with lattice spacing Qa s=0.75. 3.1. Spectral functions in relative time Starting from the initial conditions in Eq. (26), we evolve the classical Yang-Mills simulations up to the time Qt=1500, where we start the calculation of the quark spectral function. Based on the algorithm presented in Sec. 2, we then directly obtain the different components of the spectral function ρ(t+t, t, p)in the time domain. Due to the underlying symmetries, and since we consider massless fermions, we will focus on the non-vanishing vector components Reρ0 Vand ImρVof the spectral function.5They are depicted in Fig. 2for a range of momenta p/Q=0, ..., 0.21. Based on the results in Fig. 2one observes that the spectral function in the time domain features a damped oscillatory behavior, with Reρ0 V(t+t,t,p)≈e−γ(t,p)tcos(ω(t,p)t), ImρV(t+t,t,p)≈−e−γ(t,p)tsin(ω(t,p)t). (28) Clearly, the main differences to the free spectral function discussed in Sec. 2.3 concern the finite damping rate γ(p)as well as the non-trivial dispersion relation ω(p), which is nonzero even at p = 0due to the (non-)thermal mass induced by the medium. 5Numerically, we find that the pseudoscalar, axial vector and tensor components, as well as Imρ0 Vand ReρV, are suppressed by at least 2 orders of magnitude. Fig. 3. The spectral function ρ+(t,ω,p)at Qt=1500. 3.2. Spectral functions in the frequency domain Next, in order to obtain the corresponding spectral functions in the frequency domain, we perform a Fourier transform with respect to the time difference t=t−taccording to ρ0 V(t,ω,p)=2 ∞  0 dtcos(ωt)Reρ0 V(t+t,t,p), ρV(t,ω,p)=−2i ∞  0 dtsin(ωt)ImρV(t+t,t,p), (29) where we assumed that Reρ0 V(t+t,t,p)and ImρV(t+t,t,p) are even / odd functions in tfor fixed reference time t. We note that in practice, the integrals in Eq. (29)are approximated with Qtmax ∼400 −500 for the upper integration limit. We use zero padding, which implies that we interpret the Fourier transform as a usual integral with a continuous argument ωthat we evaluate using standard integration techniques at more intermediate frequencies than provided by a discrete Fourier transform. We have checked that using a Hann windowing function in the Fourier transformation similarly to Ref. [18]did not change the results. We provide a compact summary of our results in Fig. 3, where we present a three dimensional view of the behavior of the quark spectral function ρ+(t,ω,p)=ρ0 V(t,ω,p)+ρV(t,ω,p)(30) as a function of frequency ωand momentum p, noting that based on Eq. (29)the corresponding spectral function for anti-quarks ρ−(t, ω, p) =ρ0 V(t, ω, p) −ρV(t, ω, p)can be directly obtained as ρ−(t, ω, p) =ρ+(t, −ω, p). Starting from a symmetric spectral function at zero spatial momentum p =0, one observes that the spectral function becomes asymmetric along the frequency direction for p >0, with a dominant peak at a positive frequency ω+(p)and a rapidly decreasing peak at negative frequency ω−(p). While the positive frequency peak corresponds to the usual quasiparticle excitation of a quark, the excitations at ω−are referred to as ‘antiquark holes’ or ‘plasminos’ and arise from collective excitations, which emerge in thermal equilibrium [29,31]or in a non-equilibrium state as in this work. 3.3. Comparison to HTL perturbation theory The properties of the gluon spectral function in the same field configurations that we are studying have been extensively compared to HTL perturbation theory in Ref. [17]. We will here perform a similar comparison for the quark spectral function. The general structure of the HTL spectral function is given by 5 K. Boguslavski, T. Lappi, M. Mace et al. Physics Letters B 827 (2022) 136963 Fig. 4. The spectral function ρ+as a function of ω. The HTL predicted Landau damping part (green dashed) and a fit to the quasiparticle peaks (black dash-dotted) are shown separately, together with the HTL curve resulting from their sum as in (34)(blue continuous). As in all figures, error bars are shown for the data curves. They are computed as the standard error of the mean and are here of the order of the line width. ρHTL +(ω,p)=2πβ+(ω/p,p)(31) +2π[Z+(p)δ(ω−ω+(p)) +Z−(p)δ(ω+ω−(p))], where ω±(p)and Z±(p)denote the positions and residues of the quasi-particle and plasmino poles, while β+describes the contribution from the Landau cut. We provide the detailed expressions in App. A, noting that in HTL perturbation theory, these quantities are uniquely determined in terms of the momentum pand the quark screening mass mf. In leading order HTL perturbation theory this is given by m2 f=CFd3p (2π)3 g2fg(p) p,(32) with CF=(N2 c−1)/(2Nc). Within our numerical simulations we determine the quark screening mass from the relation m2 f= CF 2Nc m2 g. Here, following [17], the gluon asymptotic mass m2 gis obtained from the self-consistent solution of m2 g=2Nc (N2 c−1)d3p (2π)3 g2Tr ET(p)E∗ T(p) p2+m2 g ,(33) with the transverse field correlator ET(p)E∗ T(p) =(δij −pipj/p2)× Ei(p)E∗ j(p). Once the mass parameter m2 fis determined, HTL perturbation theory gives us a prediction for the fermion spectral function without any further parameters. In order to compare our results to HTL perturbation theory, we fit our numerical data to the following functional form ρHTL+γ +(ω,p)=2πβ+(ω/p,p)(34) +2Z+(p)γ+(p) (ω−ω+(p))2+γ2 +(p)+2Z−(p)γ−(p) (ω+ω−(p))2+γ2 −(p), where ‘HTL+γ’ refers to the HTL form supplemented with a finite width γ. While the leading order HTL spectral function (see Appendix A) features stable quasiparticles represented by delta peaks, this parametrization allows for a finite width of the peaks, which are taken to have a Lorentzian form. The free parameters in the fit are the locations of the quasiparticle peaks ω±(p), their residues Z±(p)and the widths γ±(p)for each value of the momentum p. The Landau cut contribution β+(ω/p, p)is taken to be the one from HTL perturbation theory (see Appendix Afor the explicit functional form). We demonstrate the quality of this fit, (henceforth denoted as ‘HTL+γ’ referring to Eq. (34)) in Fig. 4, where full fits are shown as blue continuous lines. The individual quasi-particle, plasmino and Landau damping contributions are also shown separately as black dash-dotted and green dashed lines. Overall, one observes an excellent agreement between our data and the HTL+γfits. Small deviations occur only for p =0, where the approximation of a width much smaller than the dispersion is not valid, and in the vicinity of ω≃−p, where the Landau cut is smeared due to interactions. We use the HTL+γfits to extract the dispersion relations ω±(p), residues Z±(p), and damping rates γ±(p)separately for each momentum pfrom our numerical lattice data, averaging over the direction p/pwhere available. Error bars are obtained as the sum of the fitting error and the standard error of the mean. The extracted values of ω±(p), Z±(p)and γ±(p)are shown in Fig. 5 as functions of momentum, together with the HTL predictions for ωHTL ±(p)and ZHTL ±(p)depicted as dashed or dotted lines. Beside the extracted values for Ns=256, Qa s=0.75, we also show results 6 K. Boguslavski, T. Lappi, M. Mace et al. Physics Letters B 827 (2022) 136963 Fig. 5. Extracted values for the dispersion relations ω±(p), residues Z±(p), and damping rates γ±(p)using (34). Results from a smaller lattice with N=64, Qa s= 1are shown for comparison. HTL predictions for ωHTL ±(p)and ZHTL ±(p)are added as gray dashed or dotted lines. The red arrows show the position of the fermion mass mf. obtained with larger lattice spacing N=64, Qa s=1, to indicate that, apart from possibly the width γ, the results are not very sensitive to discretization artifacts. We find that our results for the dispersion relations and residues agree remarkably well with the predictions from HTL perturbation theory, for which red arrows indicate the position of the fermion screening mass mfcomputed within HTL. Even the expected non-monotonic behavior of the ω−(p)dispersion is clearly visible in our data, and one also observes that, as expected from HTL, the plasmino excitation gets strongly suppressed for momenta p mf. Clearly, the most significant deviation from leading order HTL perturbation theory is the emergence of a finite decay width of quasi-particles and plasminos γ±(p)depicted in the lower panel of Fig. 5. While perturbative calculations of the fermion damping rate suffer from an infrared sensitivity to the soft gauge field propagator [29],6our non-perturbative calculation can yield first principles insights into the magnitude and momentum dependence of the damping rate. Generally, we find that γ+(p)is smaller, but of comparable size to the quark screening mass mf. One also observes from Fig. 5that the fermion damping rate γ+(p)decreases monotonically as a function of momentum, which is qualitatively different from gluonic spectral functions in non-equilibrium overoccupied plasmas, where a monotonically increasing damping rate has been observed [17,18]. 3.4. Time evolution So far we have studied the behavior of the quark spectral function at a fixed reference time Qt=1500, in the self-similar evolution of a highly occupied gluon plasma. By focusing on the behavior of the quark spectral function at vanishing spatial momentum p =0, we will now investigate the non-equilibrium evolution of the quark screening mass and damping rate, where at different evolution times Qt different separations of hard and soft scales in the system can be accessed [39,14,42,43]. We present our results for the zero momentum spectral function ρ+(t, ω, p=0)in the top panel of Fig. 6, where we show the frequency dependence of the spectral function at different times Qt=245, 735, 1960. When plotting all dimensionful quantities in terms of Q, the qualitative features of the quark spectral function in Fig. 6a) still remains essentially the same at all times. By expressing all dimensionful scales in units of the mass mF(t), this statement can be made quantitative, as shown in Fig. 6b), where all curves fall on top of each other to good accuracy, indicating that mf(t)is the only relevant scale. The time dependence of the fermion mass mf(t) ≡ω±(t, p=0) and damping rate γ(t, p=0)are depicted in Fig. 6c). We find that the time dependence of the effective quark mass exhibits an approximate mf(t)/Q∝(Qt)−1/7scaling behavior, as can be expected by evaluating the perturbative expression in Eq. (32)for the self-similar scaling behavior of the gluon distribution in Eq. (27). Direct comparison of the perturbative expression in Eq. (32), which is shown in terms of a black dashed line in Fig. 6, indicates that the extracted value of mf(t)can be described rather accurately with deviations up to a 10% level. Similarly to the effective quark mass, the quark damping rate γ(t, p=0)also decreases as a function of time, as visible in Fig. 6c). More precisely, its time-evolution is approximately the same γ(t, p=0) ∼mf(t)in the plotted time range. This is in contrast to perturbative HTL expectations [28,29], where the damping rate is expected to be proportional to the effective temperature that scales as γHTL(t, p=0) ∝g2T∗(t) ∼Q(Qt)−3/7in the selfsimilar regime. The latter would imply that the associated damping rate would decrease more rapidly in time than the thermal mass, resulting in increasingly sharp quasi-particle peaks at late times. Such behavior has indeed been observed for the gluon spectral function in Ref. [17]. In contrast, we find that due to the similar decrease of quark mass and damping rate, the spectral functions 6The analytical expression for γ(p=0)has been calculated in thermal equilibrium in Ref. [28]. However, this calculation does not directly give a precise estimate in the overoccupied regime. We are not aware of an extension of this calculation to our non-equilibrium system. It is interesting to note that in thermal equilibrium the gluon [47]and fermion [28]damping rates are similar in magnitude whereas for our system the quark damping rate is an order of magnitude larger than the gluon damping rate extracted in Ref. [17]. 7 K. Boguslavski, T. Lappi, M. Mace et al. Physics Letters B 827 (2022) 136963 Fig. 6. Spectral function ρ+=ρ0 Vat p =0for different times as a function of frequency, with all dimensionful quantities rescaled a) by Qand b) by mF(t). c) Evolution of the mass mF(t)and width γ(t, p=0)for the zero mode p =0as function of time Qt, shown on a log-log scale. Open circles / crosses correspond to extractions in the frequency / time domain using Eq. (34)/ Eq. (28), respectively, which are in excellent agreement with each other. The black dot-dashed line shows the HTL fermion mass using Eq. (32). in Fig. 6b) do not experience a significant sharpening of the quasiparticle peaks over the course of the evolution, contrary to the perturbative expectation. 4. Conclusions and outlook In this work, we have presented a novel method to perform non-perturbative real time calculations of fermion spectral functions in highly occupied plasmas. Based on a classical-statistical description of bosonic quantum fields, the fermion spectral function can be calculated by solving linearized evolution equations for fermions in the background of dynamical bosonic fields. Since only an individual momentum mode needs to be simulated at the same time, obtaining the spectral function is comutationally much less demanding than a full simulation of the fermion sector [19,20]. Based on this approach, we studied the behavior of the quark spectral function in the vicinity of a so-called non-thermal fixed point where the non-equilibrium plasma exhibits a self-similar scaling behavior. We observe Landau damping and clear quasiparticle peaks for which we extracted dispersion relations, decay widths and residues as function of the momentum. Generally the dispersion relation and residues are well reproduced by leading order HTL perturbation theory, with a single parameter – the quark screening mass mf– which we extract consistently within the HTL framework. Beyond the familiar structures of leading order HTL perturbation theory, we find that the non-perturbative spectral functions also exhibit a finite decay width γ+(t, p), and we extracted its time and momentum dependence from our simulations. Unexpectedly, the damping rate of the zero momentum γ+(t, p=0)decreases much slower than in HTL perturbation theory and remains of the same order as the mass γ+(t, p=0) ∼ mf(t), a feature that has been observed also in lower dimensional gluon spectral functions [18]. Beyond the results presented in this paper, the methodology to perform non-perturbative calculations of fermion spectral functions provides an interesting new tool to benchmark and perhaps improve perturbative calculations in the presence of strong gauge or scalar fields. Some possible extensions could include, e.g., the analysis of quark spectral functions in an expanding QCD plasma, or the investigation of the behavior of highly-energetic or heavyflavor quarks, which we intend to pursue in the future. Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgements We would like to thank A. Kurkela, A. Piñeiro Orioli, J. Peuron, S. Sharma and L. v. Smekal for collaboration on related topics. This work has been supported by the European Research Council under grant no. ERC-2015-CoG-681707, by the EU Horizon 2020 research and innovation programme, STRONG-2020 project (grant agreement No 824093) by the Academy of Finland, project 321840, by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) through the CRC-TR 211 “Strong-interaction matter under extreme conditions” Project number 315477589, and by the Austrian Science Fund (FWF) under project P34455-N. We gratefully acknowledge the National Energy Research Scientific Computing Center, a U.S. Department of Energy Office of Science User Facility supported under Contract No. DE-AC02-05CH11231, the Vienna Scientific Cluster (VSC), Austria, (project 71444), and CSC -IT Center for Science, Finland, for providing computational resources. The content of this article does not reflect the official opinion of the European Union and responsibility for the information and views expressed therein lies entirely with the authors. Appendix A. Spectral functions in HTL perturbation theory We recall here the results of perturbative calculations in the HTL framework that serve as a point of reference for interpreting our numerical results. We consider here the fermionic spectral function as computed within the hard-loop (HTL) framework at leading order. The fermionic HTL self-energy is known (see [48,31], and [30,49]in the case of a more general non-thermal state) and reads (ω,p)=m2 fd 4π γμˆ Kμ P·ˆ K(A.1) with ˆ K=(1, k/k), P=(ω+i, p)and metric signature (1, −1). The fermion mass mfis given by m2 f=(d−1)g2CF 4ddp (2π)d 2fg(p)+Nf(fq(p)+¯ fq(p)) p(A.2) with CF=(N2 c−1)/(2Nc), Nffermionic fields with distributions fq(p)(and ¯ fq(p)) for particles (anti-particles) and the distribution function of gauge fields fg(p). We are here working in the limit fg1, where the fermionic contributions can be neglected and we can connect the formula to the asymptotic mass of gluons as 8