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

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.


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,2,3,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,6,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.
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,12,13] to leading order in α s , classical statistical field simulations are commonly used in heavy-ion 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,23,24] and for a thermal system [25,26,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 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,29,30,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

Classical-statistical simulations
We consider a non-abelian SU(N c ) gauge theory discretized on a lattice with N 3 s sites and lattice spacing a s . We use N c = 2 in this work. The gauge fields are expressed in terms of lattice gauge links U (t , x) ≈ exp (iga s A j (t , x)) and electric field variables s denotes the gauge coupling. The evolution equations for the gauge field sector then result from the lattice Hamiltonian where x) are the usual plaquette variables andî, denote the unit lattice vectors in the i, j = 1, 2, 3 spatial 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 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 Here i = 1, 2, 3 is the spatial Lorentz index and r W = 1 is 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 . For a leading order tree-level improvement [19,20] we set the coefficients C n as C 1 = 4/3, C 2 = −1/6 and C n>2 = 0.
Since the equation of motion for the fermion field operator 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 p at a reference time t where λ = 1, · · · , 2N c collectively labels the spin and color indices. The operator structure is determined by the action of theb λ,p (t) andd † λ,p (t) at a fixed reference time t, where the creation and annihilation operators satisfy the usual equal-time anti-commutation relations 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 4N c N 3 s wave-functions (or N 3 . 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 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 α s in 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]

Spectral function of fermions
Generally, the spectral function is defined as the expectation value of the unequal time anti-commutator of fermion fields whereψ =ψ † γ 0 and α, β are Lorentz indices, which should not be confused with the indices λ, λ denoting the spin and color 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 × 4 matrix structure, it is useful to decompose it into scalar (S), pseudo-scalar (P), vector (V), axial-vector (A) and tensor (T) components according to 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 ρ V as where E p is 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 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 ρ V by replacing p i by the effective lattice momentump i corresponding 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 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 ∆x according to 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 x)e −ip·x denotes the spatial Fourier transform of the wave-functions. We note that the wave functionsφ λ,q (x 0 , p) depend on two momentum arguments: q which is the wavenumber at the initial time t, and p which is the momentum where the wave function is evaluated. By choosing the reference time for the mode function expansion in Eq. (5) as t = y 0 , 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 In general, the knowledge of the full set of 4N c N 3 s wavefunctions 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 4N c components of a single momentum mode φ u/v λ,p (t , x). Since each momentum mode p can 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 U and electric fields E, and evolve it by the classical Yang Mills equations (9,10) up to the reference time t from which the spectral function is measured.
2. Select a subset of N modes momentum (p) modes, for which the spectral function is computed, and initialize the N φ = 4N c × N modes fermion wavefunctions φ 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 N modes momenta p and t > t by 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, v to 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 A 0 = 0 is 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 ∂ j A j (t, x) = 0 at the time t = t when 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 Coulomb gauge and subsequently studies the response of the gauge fields to plane wave perturbations in order to extract the spectral function.

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 U i (t , x) = 1 and electric fields E i (t , y) = 0. The free fermion spectral function is simply given by [31] In the lattice discretization the three-momentum part of the four-momentum 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 iŝ with the discrete lattice momentum mode index n i = 0, · · · , N s − 1 [20]. The Wilson term generates an effective mass that makes the doubler modes more massive (and breaks chiral symmetry), so that In terms of the effective momentump and mass m p the energy of the single free fermion satisfies the usual relativistic dispersion relation E p = p 2 + m 2 p . It is convenient to express the spectral function in terms of particle and anti-particle excitations as where Λ ± (p) denote the usual projections of the Dirac components 2 The free spectral function ρ free (∆t, p) in the time domain is then obtained as 2 Note that the particle and anti-particle projections are given by 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 V and Imρ V calculated on a 64 3 lattice with momentum a s p = (0.098, 0.195, 0.29) and mass parameter ma s = 0.003125 corresponding to nearly massless fermions. An excellent agreement between continuous lines, depicting the analytical expressions in Eq. (25), and points, corresponding to the numerical lattice data, is observed for all components, 3 validating our procedure to calculate spectral functions.

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 with p = |p| and where n 0 /g 2 1 is the initial occupancy and Q is 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,38,39,40,41,14,42,43,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 phase-space distribution can be described in terms of a scaling function f s and universal scaling exponents α = −4/7, β = −1/7 [37,38,39,40,41,14,42,43]. 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 n 0 = 0.2 4 and 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 microscopic scales of the gluon field configurations. Note that the dependence of the spectral function on the relative time t − t happens at a much shorter timescale ∼ 1/m g than 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 N s = 256 lattices for nearly massless fermions m = 0.003125 Q with lattice spacing Qa s = 0.75.

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 V and Imρ V of the spectral function. 5 They are depicted in Fig. 2 for a range of momenta 4 We 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. 5 Numerically, we find that the pseudoscalar, axial vector and tensor components, as well as Imρ 0 V and Reρ V , are suppressed by at least 2 orders of magnitude. Reρ 0 V (t + ∆t, t, p) ≈ e −γ(t,p)∆t cos(ω(t, p)∆t) , Imρ V (t + ∆t, t, p) ≈ −e −γ(t,p)∆t sin(ω(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 = 0 due to the (non-)thermal mass induced by the medium.

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 −t according to where we assumed that Reρ 0 V (t + ∆t, t, p) and Imρ V (t + ∆t, t, p) are even / odd functions in ∆t for fixed reference time t. We note that in practice, the integrals in Eq. (29) are approximated with Q∆t max ∼ 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 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.

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 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. Appendix A, noting that in HTL perturbation theory, these quantities are uniquely determined in terms of the momentum p and the quark screening mass m f . In leading order HTL perturbation theory this is given by with C F = (N 2 c − 1)/(2N c ). Within our numerical simulations we determine the quark screening mass from the Here, following [17], the gluon asymptotic mass m 2 g is obtained from the self-consistent solution of with the transverse field correlator E T (p)E * T (p) = (δ ij − p i p j /p 2 )E i (p)E * j (p). Once the mass parameter m 2 f is 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) 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 A for 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 p from our numerical lattice data, averaging over the direction p/p where 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 Z HTL ± (p) depicted as dashed or dotted lines. Beside the extracted values for N s = 256, Qa s = 0.75, we also show results 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 m f computed 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 m f .
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 in- frared sensitivity to the soft gauge field propagator [29], 6 our 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 m f . One also observes from Fig. 5 that 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]. 6 The 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].

Time evolution
So far we have studied the behavior of the quark spectral function at a fixed reference time Qt = 1500, in the selfsimilar 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 m F (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 m f (t) is the only relevant scale.
The time dependence of the fermion mass m f (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 m f (t)/Q ∝ (Qt) −1/7 scaling 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 m f (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) ∼ m f (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) ∝ g 2 T * (t) ∼ Q(Qt) −3/7 in the self-similar 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 in Fig. 6b) do not experience a significant sharpening of the quasiparticle peaks over the course of the evolution, contrary to the perturbative expectation.

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 nonthermal fixed point where the non-equilibrium plasma exhibits a self-similar scaling behavior. We observe Landau damping and clear quasi-particle 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 parameterthe quark screening mass m f -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) ∼ m f (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.