Fermionic spectral functions with the Functional Renormalization Group

We present first results on the calculation of fermionic spectral functions from analytically continued flow equations within the Functional Renormalization Group approach. Our method is based on the same analytic continuation from imaginary to real frequencies that was developed and used previously for bosonic spectral functions. In order to demonstrate the applicability of the method also for fermionic correlations we apply it here to the real-time quark propagator in the quark-meson model and calculate the corresponding quark spectral functions in the vacuum.


I. INTRODUCTION
The spectral properties of strongly interacting matter under extreme conditions, as encountered in the early universe and compact stellar objects, are of fundamental importance for identifying the relevant degrees of freedom in the equation of state and respective transport properties. Spectral functions are real-time quantities, while the underlying equilibrium state is commonly obtained from imaginary-time (Euclidean) evaluations of the partition function. In this setting, a thermodynamically consistent computation of the spectral properties poses a major challenge since analytic continuations of the pertinent Euclidean n-point functions are required. In relativistic theories this entails a transition from Euclidean to Minkowski space-time, see for example [1][2][3][4].
In the context of the strong interaction (QCD) quark spectral functions are of particular interest. A Bayesian reconstruction method has, for example, been used in [5] and [6] to extract quark spectral functions from Euclidean data obtained from Dyson-Schwinger equations. Within the Functional Renormalization Group (FRG), which incorporates thermal as well as quantum fluctuations, Euclidean quark propagators have recently been calculated in [6] and [7]. In the present work, we focus on the calculation of real-time quark propagators. Instead of using numerical reconstruction methods, see e.g. [8,9], we perform the analytic continuation on the level of the FRG flow equations for retarded two-point correlation functions which are then solved directly in the corresponding domain of frequencies close to the real axis. Such analytic continuation methods have been put forward in [10] and [11][12][13]. In [14][15][16] the approach was extended to finite temperature, finite quark chemical potential as well as to finite spatial momenta and has been used to calculate mesonic spectral functions within the quark-meson model. The analytically continued FRG flow equations for the corresponding two-point correlation functions were thereby solved in a simple but thermodynamically consistent and symmetry preserving truncation which in the long-wavelength and static limit reduces to the leading-order derivative expansion used for the underlying effective potential. In [17] this approach was extended to calculate in-medium vector-and axial-vector meson spectral functions. For the first time we here present an FRG calculation of fermionic spectral functions obtained from analytically continued flow equations which can be solved numerically. In a first step we restrict ourselves to the vacuum and to vanishing external spatial momenta.
This work is organized as follows. In Sec. II we briefly introduce the FRG framework and its application to the quark-meson model. Our analytic continuation method and the flow equation for the real-time quark propagator and the spectral functions are discussed in Sec. III. We have solved the flow equations using both, a grid and a Taylor-expansion method as discussed in Sec. IV where we also demonstrate the particular advantages and disadvantages of either method. Results for the quark mass and the mass dressing function are presented in Sec. V while respective results for the quark propagator and the quark spectral functions are shown in Sec. VI. Various sum rules, which can be derived from the Lehmann representation of the quark propagator, are discussed in Sec. VII. We close with our summary and outlook in Sec. VIII. Further details are collected in an appendix.

II. FUNCTIONAL RENORMALIZATION GROUP AND QUARK-MESON MODEL
The Functional Renormalization Group (FRG) is a non-perturbative approach that is used for example in quantum and statistical field theory, in particular for strongly interacting systems, see e.g. [18][19][20][21][22][23][24][25] for reviews. It is formulated in (continuous) Euclidean space-time and combines Wilson's idea of the renormalization group in momentum space [26,27] with functional methods in quantum field theory.
In the following we will use the formulation pioneered by Wetterich [28] which aims at calculating the effective average action Γ k where k is the renormalization group scale. At the ultraviolet (UV) scale k = Λ, the effective average action is basically given by the bare action S of the chosen model and does not include any fluctuation effects. By lowering the scale k the effects of quantum and thermal fluctuations are gradually included until the full effective action Γ = Γ k=0 is obtained in the limit k → 0. The scale-dependence of Γ k is given by the following flow equation, also known as the Wetterich equation, where R k is a regulator function that suppresses momentum modes with momenta smaller than k, 1 and Γ (2) k is the second functional derivative with respect to the fields. Both Γ (2) k and R k can be represented as matrices in the field space of bosonic and fermionic variables, and the supertrace runs over field space as well as all internal indices also including an integration over internal momenta.
We will apply this flow equation to the quark-meson model as a low-energy effective theory for the chiral aspects of QCD with two flavors [30,31]. It includes quarks, the sigma meson and the pions as effective degrees of freedom which interact via a Yukawa-type interaction. We use the following ansatz for the effective average action of the quark-meson model in Euclidean space-time, with φ 2 = σ 2 + π 2 . This approximation, which is the leading order in a derivative expansion where the only scale-dependent object is the effective potential U k (φ 2 ), is also called the local potential approximation (LPA) [32,33]. When inserting this ansatz into the Wetterich equation, one obtains the flow equation for the effective potential, where explicit expressions for the threshold functions I k are given in App. A. At the UV scale Λ we choose the effective potential to be symmetric, and then solve the corresponding flow equation numerically, see Sec. IV. The term cσ which breaks chiral symmetry explicitly and thus plays the role of the (up/down) current quark mass in QCD, is added to the effective potential in the infrared (IR) while spontaneous chiral symmetry breaking occurs dynamically through the fermionic fluctuations which are included by solving the flow equation. The solution for the scale-dependent effective potential is then used as input for the calculation of the fermionic two-point function.
In order to obtain the flow equation for the quark twopoint function we take two functional derivatives of the Wetterich equation, Eq. (1), with respect to the fermionic fields which gives see Fig. 1 for a diagrammatic representation. Therein, q = (q 0 , q) is the internal and p = (p 0 , p) the external momentum, D = (Γ (2) k +R k ) −1 is the full regulated propagator, the vertex functions Γ (3) ψψφ are obtained from the ansatz in Eq. (2), and the remaining trace represents a summation over all internal indices as well as an integration over internal momenta. Explicit expressions for the vertex functions Γ (3) ψψφ as well as for the bosonic and fermionic regulator functions, R B and R F , are given in App. A. As in the original studies [13][14][15], here we use three-dimensional regulator functions which only regulate spatial momenta but not the energy components at the expense of some breaking of the Euclidean O(4) symmetry. This breaking was assessed and found to be negligible for external momenta well below the UV cutoff scale Λ in the Euclidean two-point functions, with still reasonably small and only quantitative effects in the time-like domain after analytic continuation [13]. While this can be avoided in principle [34,35], the three-dimensional regulators allow to perform the integration over the internal energy component or the corresponding Matsubara sum at finite temperature analytically which tremendously simplifies the analytic continuation procedure discussed in the next section.

III. ANALYTIC CONTINUATION AND SPECTRAL FUNCTIONS
In the UV, the Euclidean quark two-point function is given by with hermitian γ-matrices as can also be seen from Eq. (2). We perform the analytic continuation by treat- ing the Euclidean energy variable as a general complex argument z = ip 0 , so that with z = ω + i one obtains 2-point functions with retarded boundary conditions in the limit → 0. 2 To reproduce the standard form of the Dirac operator we furthermore introduce an overall minus sign as compared to Euclidean conventions adopted in (6), to write where we formally kept the Euclidean γ's which are usually changed to be anti-hermitian (by replacing γ → i γ) in Minkowski space-time, of course. Based on this Dirac structure, we therefore make the following ansatz for the scale-dependent quark two-point function, withˆ p ≡ p/| p|. The dressing functions can be obtained from the full two-point function as follows, The UV initial conditions for the dressing functions are thus given by The flow equation for the quark two-point function, then leads to flow equations for the individual dressing functions of the form k,ψσ (ω, p) For later convenience we also write down corresponding expressions for the inverse of the quark two-point function in Eq. (8), because we will need the imaginary parts of the retarded quark propagator for the various spectral functions, . (19) Note that this is not the regularized propagator D = (Γ (2) k +R k ) −1 used in the loops, but the (retarded) inverse of Γ (2) k,ψ in (8). As such, in the UV, it is given by as usual. At zero temperature, one generally has i.e., both are determined by a single dimensionless renormalization function Z of the invariant four-momentum k,ψ (ω, p) are then essentially the same, likewise. Here we keep them formally distinct, nevertheless, so that their flow equations can be readily extended to finite temperature and density in the future.
The quark spectral function can be obtained from the retarded propagator by taking the imaginary part, and therefore has the same Dirac structure as the propagator and the 2-point function, In the UV, the quark spectral function is given by The individual components of the spectral function can then be obtained directly from the imaginary parts of the corresponding propagator components, k,ψ (ω, p) are essentially the same as well, one then usually writes, We also note that the spectral function ρ The Lehmann representation of the retarded propagator is given by It can be used to derive various sum rules for the spectral functions as discussed in Sec. VII below. In the following we will set the spatial external momentum to zero for simplicity in this first study, p = 0, which implies G k,ψ (ω, 0) ≡ 0 as well, and drop the corresponding second argument in all momentum dependent functions. The quark propagator can then be decomposed as with Λ ± = (1 ± γ 0 )/2 and With these we define the associated quark and anti-quark spectral functions, such that ρ + k (−ω) = ρ − k (ω). These are then related to the previously defined quark spectral functions by and can be expressed in terms of the dressing functions of the two-point function as

IV. NUMERICAL IMPLEMENTATION
The flow equations for the effective potential, Eq. (3), and for the dressing functions of the quark two-point functions, Eq. (16), are solved using two different methods: the grid method and the Taylor method. Both methods use the same input parameters which are summarized in Tab. I. These parameters are chosen such as to reproduce physical vacuum values for the pseudo-particle masses and the pion decay constant in the IR, see Tab. II. We note that there are small differences between the IR values obtained from the grid method and the Taylor method which, however, will not play any role in the following since we will only focus on qualitative differences in the results from these two methods, see also [36] for a comparison of different numerical implementations.

A. Grid method
The idea of the grid method is to discretize the field variable φ 2 = σ 2 + π 2 on a grid in field space, see [31] for details. The flow equation for the effective potential then turns into a set of coupled ordinary differential equations which can be solved using standard methods. The global minimum σ 0 of the effective potential in the IR determines the expectation value of σ which is identified with the pion decay constant, σ 0 ≡ f π , at this leading-order in the derivative expansion.
The flow equation for the dressing functions of the quark two-point function can then be solved using the scale-dependent effective potential as an input. Since the flow equation for the two-point function does not couple different grid points, it is sufficient to use only one grid point φ 2 i which corresponds to the location of the global minimum at some chosen scale k. In the following we are mostly interested in the IR and therefore choose φ 2 i = σ 2 0,IR . We also note that the flow equation for the two-point function is solved down to k = 0 by using an extrapolation of the flow of the effective potential for k < k IR = 40 MeV. The same technique is used for the Taylor method.
One of the advantages of the grid method is that it does not restrict the shape of the effective potential and therefore allows for an almost arbitrary order of mesonic selfinteractions (limited only by the number of grid points). In particular, possible secondary minima of the potential can be resolved which allows to study first-order phase transitions straightforwardly.
One possible issue with the grid method arises when solving the flow equation for the retarded two-point function. As discussed in App. A, the k-integration of the flow equation for the imaginary part of the two-point function, which is closely connected to the spectral function, collapses to a sum over a few scales k i due to the appearance of Dirac delta functions. This means that for a given energy ω, the spectral function may receive only a single contribution from some intermediate scale k > k IR . Since the flow equation is always evaluated at the IR minimum φ 2 i = σ 2 0,IR , also at k > k IR , this contribution does not correspond to the actual scale-dependent minimum of U k (φ 2 ). We will show, however, that the difference in the IR between spectral functions obtained from either grid or Taylor method (where this problem does not occur) is reasonably small.

B. Taylor method
There are different versions of the Taylor method available in the literature. The general idea is to use an expansion in the form of a Taylor series around some value of the field ρ ≡ φ 2 . The chosen expansion point may be constant, see for example [37], or it may be scale dependent as discussed in the following. We first write the effective potential as where the expansion point ρ 0,k is the scale-dependent minimum and we choose K = 5. We use the same ansatz for the flow equation of the effective potential, where ∂ k U (n) k (ρ 0,k ) denotes the n-th derivative of the flow equation for the effective potential with respect to ρ, evaluated at ρ 0,k . From these two equations one can obtain flow equations for the coefficients a n and for ρ 0,k , which are given by When using that ∂ σ U k (ρ, σ)| σ= √ ρ 0,k = 0, we find a 1,k = c/(2 √ ρ 0,k ). By taking the RG-scale derivative of this relation one can express the flow equation for the scaledependent minimum in terms of the flow equation for a 1,k , which gives The same Taylor expansion can also be used for the mesonic two-point functions, see e.g. [13]. For the quark two-point function, which in the UV contains a term hσ, we will use a Taylor expansion in terms of σ 0,k instead of ρ 0,k . The quark two-point function can then be written as where we use L < 5. We can make a similar ansatz for the flow equation, where ∂ k Γ (2),(n) k (σ 0,k ) denotes the n-th derivative of the flow equation for the two-point function with respect to σ, evaluated at σ 0,k . From these two equations we can obtain flow equations for the coefficients c n,k which are given by This Taylor method has the advantage that it always uses an expansion around the scale-dependent minimum of the effective potential. In this way, the two-point function only receives contributions that correspond to the local minimum at a given scale, in contrast to the grid method where the contributions are in general not evaluated at the global minimum for intermediate RG scales k. In the next sections we will present results obtained from both methods and discuss their differences.

V. MASSES AND DRESSING FUNCTIONS
B k (ω) AND C k (ω) We will first study the flow of the Euclidean (curvature) masses as obtained from the effective potential using the grid method and the Taylor method, which are then used as input for the flow equation of the two-point function, see Fig. 2. Explicitly we have When evaluated at the scale-dependent global minimum of the effective potential, σ 0,k , these expressions represent the Euclidean curvature masses. When using the grid method as implemented in this work, however, the flow equation for the two-point function is always evaluated at the IR minimum σ 0,IR while the Taylor method uses the scale-dependent minimum σ 0,k . The scale-dependent masses can of course also be obtained when using the grid method and evaluating Eq. (42) at σ 0,k . The masses then essentially agree with the ones obtained from the Taylor method, see Fig. 2.
We note that the quark mass obtained from the grid method using a fixed value of φ 2 = ρ 0,IR is constant while when using the Taylor method the quark mass is almost zero in the UV and then significantly increases at the chiral symmetry breaking scale of k ≈ 600 MeV where also the pion and the sigma mass become different. The masses obtained from both methods agree reasonably well in the IR, i.e. up to the small deviations recorded in Tab. II which could be compensated by a small readjustment of the UV parameters.
We will now study the flow of the mass dressing function B k divided by the renormalization function Z k , see Eq. (21), of the retarded quark two-point function, as introduced in Eq. (8). At ω = 0, where they are both real, their ratio is shown in Fig. 3 where we compare the result from the grid method with those from the Taylor method for different expansion orders L. When using the grid method, we have B k=Λ = hσ 0,IR = 299 MeV while for the Taylor method we have B k=Λ = hσ 0,UV = 8.9 MeV. Despite this large difference in the UV, both methods give approximately the same result in the IR, except for the Taylor method with L = 0. Note that the Taylor method at order L = 0 here produces a result that would agree with that from the grid method, if the scale-dependent minimum σ 0,k was used in the latter in the place of the IR minimum σ 0,IR that appears for example in the derivatives of the effective potential and in the initial condition for Γ (2) which contains hσ, see Eq. (13). Such a simple modification of the grid method would neglect contributions to the flow that arise from the k dependence of the gliding minimum σ = σ 0,k on the other hand. By comparison with the Taylor results, where these contributions are contained in the second term in Eq. (41) and are thus absent at L = 0, we see that they are by no means negligible but contribute substantially to the dynamical mass generation. Already the L = 1 result is close to those from the higher orders in the Taylor expansion, however, and hence captures the main effect. We also note that the flow of the ratio B k /Z k in the static limit resembles the flow obtained for the quark mass parameter in Fig. 2 for the higher Taylor orders quite well.
In Fig. 4 we compare the flow of B k /Z k at ω = 0 with the flow of the scale-dependent quark pole mass m P ψ,k which is obtained as the (lowest) value of ω that solves B k (ω) = C k (ω) in the range where both dressing functions are still real and thus in the IR describes the stable single-particle contribution to the quark propagator here. We find that the pole mass starts to deviate from B k /Z k for k 400 MeV and that the pole mass is several MeV larger than the ω = 0 value of B k /Z k in the IR. In general, a difference between the pole mass and the mass function B k /Z k at ω = 0 is to be expected since B k /Z k has a non-trivial ω dependence which is seen in Fig. 5 where we plot the real parts of B and C (in the IR) over p 2 0 . With increasing negative values of p 2 0 = p 2 here, i.e. deeper in the timelike region, both B and C increase as long as they remain real until nonzero imaginary parts develop at their peak position as discussed in the next section. For positive p 2 0 = p 2 , i.e. in the spacelike region, C becomes pure imaginary because it is defined as C(z) = zZ(z) with z = ip 0 in the Euclidean domain, while B monotonically decreases with the Euclidean p 2 . Because the explicit flow of the two-point function vanishes when p 2 Λ 2 , B approaches gσ 0,IR = m ψ in Tab. II. In the grid method the IR value of σ 0 is the same as its UV value which means that also B k remains at its UV initial condition, Eq. (13), so that B k does not flow at all in this case. In the Taylor method one starts with a small value for B k=Λ = hσ 0,UV , as mentioned above, but the flow of σ 0,k leads to an implicit residual flow also for B k at asymptotically large momenta as seen explicitly in Eq. (41) (for L ≥ 1). This contribution ensures that B k eventually approaches the corresponding infrared value of the quark mass parameter in Tab. II with the Taylor method as well.

VI. QUARK SPECTRAL FUNCTION
We now turn to the flow of the different quark spectral functions: ρ k,ψ , ρ k,+ and ρ k,− . They all depend on the dressing functions B k (ω) and C k (ω). The infrared results, as obtained with k → 0 from either grid or Taylor method, for real and imaginary parts of B(ω) are shown in Fig. 6, the corresponding ones for C(ω) in Fig. 7.
The shapes of B(ω) and C(ω) which are also reflected in the spectral functions can be explained by considering the different particle processes that can occur within our framework. These processes can already be read off from the diagrammatic representation of the flow equation for  Taylor, L=4 C -Grid C -Taylor, L=4 FIG. 5: The real part of the functions B k and C k for k → 0 as obtained by using the grid method and the Taylor method with L = 4. We note that the quark mass obtained from the potential is the same as the value for B in the limit p0 → ∞.
the two-point function, see Fig. 1, and are given by where ψ * denotes an off-shell quark with energy ω and the other energies represent IR values. The process ψ * → ψ + π for example describes a quark with energy ω that can 'decay' into an on-shell quark with energy E ψ and a pion with energy E π . These processes allow for a clear interpretation of the shape of the real and in particular of the imaginary part of B(ω) and C(ω). The real part of B(ω) starts with a small positive slope at small energies and monotonously increases up to the first threshold at ω = E ψ + E π ≈ 420 MeV where the decay channel into a quark and a pion opens up. A second, but smaller change in the real part is visible at the second threshold at ω = E ψ + E σ ≈ 820 MeV where the process ψ * → ψ + σ becomes possible. The imaginary part of B(ω) stays at zero below the first threshold energy ω = E ψ + E π since no decay channels are available and then starts to rise quickly at the quark-pion threshold. The quark-sigma process gives rise to a small kink in the imaginary part at ω = E ψ + E σ .Up to these energies the result obtained from the grid method is in very good agreement with the results from the Taylor method for L ≥ 1. For higher Taylor orders L ≥ 3 we observe numerical difficulties at large energies. The behavior of C(ω) is is analogous: below the quark-pion threshold it stays real, with a rapidly increasing imaginary part starting there and further kinks in real and imaginary parts at the quark-sigma threshold. The leading behavior at small ω is given by C(ω) = Zω+O(ω 2 ) with Z = Z k=0 (0) ≈ 1.94. Note that both B(z) and C(z) are analytic functions at z = 0 as can be seen from their Lehmann representation and the fact that their spectral functions have no support there as we will discuss next.
We now turn to the quark spectral function ρ (C) k,ψ (ω). It starts in the UV as a simple delta function at the UV quark mass parameter with strength 1/2, cf. Eq. (24), and flows with k → 0 towards the infrared results shown in Fig. 8 as obtained from the grid method and the Taylor method. Although the UV values are very different, cf. Eq. (12), the spectral functions agree well in the IR, in particular for higher Taylor orders like L = 3. While the delta peak that is connected to the quark pole mass with the Taylor method moves from ω ≈ 9 MeV in the UV up to a value of ω ≈ 366 MeV in the IR, the pole mass obtained from the grid method changes with the flow only from ω ≈ 299 MeV in the UV to ω ≈ 375 MeV in the IR. The remaining discrepancy of about 9 MeV between the IR pole masses could in principle be compensated by a slight readjustment of the UV parameter as mentioned above. These single particle contributions at the physical mass are defined as the solutions to B(ω) = C(ω) and indicated by the vertical lines in the figures.
The difference between curvature masses, obtained from the effective potential and the physical pole masses has also been observed for mesons, see e.g. [14]. While the two need not be the same of course, but differ whenever one has frequency dependent renormalization effects as those in the ratio B(ω)/Z(ω) here, and explicitly demonstrated in Fig. 5 above, the size of the difference is generally determined by the relative distance of the closest singularity above the single-particle pole. As such the effect observed especially for the pions in the previous LPA studies was too large. To reduce this artifact one has to go beyond the leading order in the derivative expansion by least including scale-dependent wavefunction renormalization factors [38]. We should therefore be prepared that analogous improvements can also lead to sim- ilar quantitative changes here. As a next step towards a fully self-consistent calculation one should therefore extract these scale-dependent wavefunction renormalization factors and feed them back into the calculation by iteration in future as well.
We have already noted that the Taylor method at the leading order, with L = 0, misses an important qualitative effect. Higher orders on the other hand appear to converge quickly towards the grid result. We then generally have quite compelling agreement between both methods in the IR (provided the external frequency ω stays well below the cutoff scale Λ). The remaining discrepancies can in fact be considered as an indication of systematic uncertainties. The particular advantage of the Taylor-expansion method is that it provides a direct and more intuitive understanding of the evolution of the spectral function with the flow at intermediate scales k. To achieve this with the grid method one would best stop the flow at some intermediate scalek and study the behavior of the spectral functions in the fixed background with σ = σ 0,k given by the minimum of the effective potential at this intermediate scalek. In the following we focus on the physical results obtained at the end of the flow in the IR. Having established the equivalence of the methods there, we restrict to the grid results for clarity from now on.
As a first check note that ρ (C) ψ (ω) is the positive distribution that it has to be in a theory with a positive state space, and one has the inequality which is satisfied here as well. For the continuum contributions this can be seen explicitly in Fig. 9 where we show ρ antiquark spectral functions, cf. Eq. (32), which are therefore also both positive (while one has ρ ψ (ω) have the same magnitude. Consequently, the quark spectral function ρ + (ω) only exhibits one delta peak at positive energies, representing a single quark, while the antiquark spectral function ρ − (ω) only has a peak at negative energies for the single-antiquark states.
The continuum parts in ρ ± (ω) related to the various one-to-two-particle processes differ as follows: In ρ + (ω) the continuum contributions of ρ

VII. SUM RULES
The various sum rules as usual follow from expanding the Lehmann representation in Eq. (28), for small and large ω. For large ω, this leads to where we have used that the even moments of the odd function ρ (B) ψ (ω) and the odd moments of the even ρ (C) ψ (ω) both vanish. The leading order at large ω therefore corresponds to, In a renormalizable field theory, the integral of the spectral density usually diverges logarithmically and the right-hand side is then given by a formally ultraviolet divergent field renormalization constant. Here, we note that the flow of the two-point function generally vanishes for ω → ∞ and the leading-order sum rule for ρ (C) k,ψ (ω) consistent with the UV initial condition in Eq. (14) therefore reads It is satisfied exactly in the UV by definition, where the quark spectral function is given by with Y Λ = 1/2. As fluctuations are included with integrating the flow down to some scale k, the weight of these free single-(anti)particle contributions gets reduced in favor of continuum contributions from processes that are possible due to the interactions with fluctuations included above this scale. Monitoring the sum rule (53) over the flow therefore provides a valuable test of the consistency of the procedure. For this purpose, we evaluate all sum rules numerically up the the sum of the UV energies, Λ E = E ψ,Λ + E α,Λ ≈ 2324 MeV, since at this scale the FRG flow gives the first contribution to the continuum and the strength of the Dirac delta peak starts to decrease. During the flow, the quark and meson energies decrease and the threshold of the continuum moves to smaller energies until it reaches its IR value. When applied to the quark spectral function obtained from the grid method for example this then yields in the IR Herein, we have included the delta functions from the single-particle contributions to the spectral function explicitly, with scale dependent pole mass m P ψ,k and residue Y k , by using We find 2Y k=0 ≈ 0.51 which means that more than 50% of the sum rule are still provided by the singleparticle and antiparticle contributions, which only leaves less than half of the sum rule for the continuum from the interactions here, at T = µ = 0. The fact that the sum rule remains satisfied numerically to a very high degree all the way down to the infrared demonstrates the consistency of the FRG approach in that it also keeps the normalization intact, in addition to preserving Dirac and symmetry structures. At next-to-leading order, the expansion in Eq. (50) corresponds to the energy-weighted sum rules for the quark propagator, By the same argument as above, the flow for the twopoint function vanishes in this the limit which is therefore here given by the UV initial conditions in Eqs. (13) and (14) again. The energy-weighted sum rule for the quark propagator therefore here becomes Our parameters for the grid method in Table II imply for the quark mass at the UV cutoff scale m ψ,Λ = 299 MeV (the quark mass parameter m ψ is in fact scaleindependent when using the grid code). For comparison, the numerical value of the energy-weighted integral over ρ (B) (ω) in the infrared would correspond to m ψ,Λ = 314.9 MeV. This sum rule is therefore slightly violated, but only within an error of about 5%. The leading-order sum rule (52) for ρ (C) ψ (ω) and the next-to-leading-order sum rule (59) for ρ (B) ψ (ω) are determined by the purely perturbative behavior of the quark propagator. Higher moments of the spectral functions diverge (more than logarithmically) in the ultraviolet. The corresponding contributions to the propagator are suppressed by powers of 1/ω 2 relative to the leading ones at large ω and need to be obtained from the operator product expansion.
Negative moments on the other hand will converge more and more rapidly in the ultraviolet. As long as there are no contributions to the spectral function from massless excitations there will be no infrared divergences either. To obtain these moments one expands Eq. (28) for small ω, where we have again used that the even(odd) moments of ρ ψ (ω)) vanish. The non-vanishing ones then lead to the leading-order, and the next-to-leading-order sum rule, Both sides of these sum rules are now scale dependent. For the leading-order sum rule we start in the UV with 1/B k=Λ (0) = 1/m ψ,k=Λ ≈ 3.344 · 10 −3 MeV −1 and the sum rule is trivially satisfied. More importantly, however, in the IR we find 1/B k=0 (0) ≈ 1.447 · 10 −3 MeV −1 , and this then compares very well with the numerical value of the integral in Eq. (63) which gives 1.450 · 10 −3 MeV −1 .

VIII. SUMMARY AND OUTLOOK
In this work we have presented a framework to calculate fermionic spectral functions within the Functional Renormalization Group approach. Our method is based on a recently developed technique for the calculation of bosonic spectral functions that uses a well-defined analytic continuation procedure from imaginary to real energies. The fact that no numerical analytic continuation method is needed represents a distinct advantage over other approaches that have to rely on numerical reconstruction techniques like the Maximum Entropy Method.
In the present study we applied this method to the quark-meson model and calculated quark spectral functions in the vacuum. The resulting flow equations for the real-time two-point functions have been solved numerically using two different methods: the grid method and the Taylor method. Both methods produce consistent results when the corresponding flow equations are integrated all way down to the infrared, with residual discrepancies serving as indications of the systematic uncertainties. Thereby, the grid method by and large produces the more stable results, while the Taylor method provides the more direct and intuitive interpretation of the full scale-dependence of the spectral functions during the flow.
In particular, we studied the flow of the quark mass as well as the quark and anti-quark spectral functions. The different particle processes which define the shape of the spectral functions as well as the consistency with various sum rules, derived from the Lehmann representation of the quark propagator, have been discussed.
Although we have limited ourselves to the vacuum and to vanishing spatial momenta in this first work on fermionic spectral functions with the FRG, our approach can be extended to finite temperature, finite chemical potential and finite spatial momenta as already demonstrated for mesons. These extensions will also allow for the calculation of transport coefficients like the shear viscosity. Other extensions that are left to future studies include the improvement of the presently used truncation by introducing wavefunction renormalization factors or a scale-dependent Yukawa coupling. Replacing the quarks here by nucleon fields and their parity partners allows to study the corresponding baryonic spectral functions in the parity-doublet model with fluctuations beyond meanfield as in [39] in order to describe the liquid-gas transition of nuclear matter and the chiral transition at high baryon density in a unified framework. This can then furthermore be extended to include vector and axialvector mesons along the lines of [17] and study their spectral changes in dense nuclear matter.