Spectral functions from the real-time functional renormalization group

We employ the functional renormalization group approach formulated on the Schwinger-Keldysh contour to calculate real-time correlation functions in scalar field theories. We provide a detailed description of the formalism, discuss suitable truncation schemes for real-time calculations as well as the numerical procedure to self-consistently solve the flow equations for the spectral function. Subsequently, we discuss the relations to other perturbative and non-perturbative approaches to calculate spectral functions, and present a detailed comparison and benchmark in $d=0+1$ dimensions.

Thus, the knowledge of spectral functions of the various strongly interacting particles is highly desirable when trying to investigate e.g. dilepton production, transport coefficients or the melting of quarkonium states in the quark-gluon plasma (QGP).
Unfortunately, extracting real-time information of strongly coupled systems is a difficult problem. The non-perturbative nature of QCD at energies below and around the phase transition prohibits the use of perturbative methods. Recently, there has been progress concerning the spectral functions of quarkonia and some transport coefficients coming from euclidean lattice simulations [1][2][3][4][5]. However, the analytic continuation of the numerical data to Minkoswki space and other problems make these investigations quite challenging and so far there are no lattice results for spectral functions of lighter hadrons.
So far our knowledge about spectral properties of thermal QCD matter comes primarily from calculations in low energy effective theories of QCD, based on a variety of different techniques including (re-summed) perturbative calculations [6][7][8] as well as non-perturbative functional approaches [9][10][11][12][13][14]. Recently, there has been great success in applying the analytically continued functional renormalization group (FRG) [15,16] to low energy effective models of QCD [17][18][19][20][21][22][23][24][25]. While many of the results from analytically continued FRG calculations have been impressive, it still is desirable to pursue non-perturbative functional calculations directly in Minkowski space. In this paper we adopt a real-time FRG approach on the Schwinger-Keldysh (SK) contour [26][27][28][29][30][31][32][33][34][35][36] to extract spectral functions in the O(N) model without the need for analytical continuation. By performing a careful perturbative analysis we show that -in the absence of spontaneous symmetry breaking -local potential approximations (LPA) are not able to generate a broadening of the spectral function. We therefore develop a truncation, based on a vertex expansion that includes momentum dependent four-point functions, which is able to capture the broadening of the spectral function as the propagators in this truncation are two-loop complete. One important feature of our method is that it is applicable for both quantum and classical-statistical field theories, such that we can compare and evaluate our results from the real-time FRG approach against non-perturbative classical-statistical real-time lattice simulations [37][38][39][40].
This paper is organized as follows: We start in section II with an introduction to dissipative classical-and quantum field theories on the SK contour and the formulation of the real-time FRG approach. After defining a d + 1 dimensional regulator scheme that respects time-ordering on the SK contour we introduce a diagrammatic notation simplifying the where S C [ϕ,φ] is the contour action on the Schwinger-Keldysh contour. Denoting the thermal distribution function of a bosonic quantum system as n qu eff (ω) = n BE (ω) + 1 2 = 2 coth βω 2 , where n BE (ω) = 1 e βω −1 is the Bose-Einstein distribution, the contour action S C [ϕ,φ] for a dissipative quantum system coupled to an external heat-bath at inverse temperature β = 1 k B T and with the rest-frame u µ = (1, 0, 0, 0) is explicitly given by where x = ∞ −∞ dx 0 d d x such that the real-time axis extends from x 0 = −∞ to x 0 = +∞ describing a time translation invariant system in thermal equilibrium [42]. While the contour action in Eq. (3) describes a dissipative quantum system with Model A type dynamics [43], the case of a non-dissipative quantum system with conservative Model C/G type dynamics 1 is obtained in the limit γ → 0 + , where the coupling to the external heat bath ultimately vanishes, but as usual in the iǫ prescription is required at intermediate steps of the calculation to ensure the correct time ordering of the propagators and convergence of the functional integral. Specifically, in the absence of interactions (λ = 0) the free propagators of the theory in momentum space   e +ip(x−y) (4) are explicitly given by G R,ab 0 (ω, p) = −1 ω 2 − E 2 p + i γ β ω δ ab , iF 0,ab (ω, p) = 0 , 1 Single component scalar theories (N = 1) classify as Model C, whereas multi-component scalar theories (N ≥ 2) feature an additional conserved current, e.g. for N = 4 one has j µ ab (x) = ǫ abcd ϕ c (x)∂ µ ϕ d (x), and therefore classify as Model G [44].
with E p = p 2 + m 2 such that in the limit γ → 0 + , the above expressions reduce to the familiar expressions for the retarded/advanced (G R/A ) and symmetric (iF ) two-point functions, who's operator definitions and basic properties are recalled in Appendix A.
Expressing the contour action in Fourier space , which as discussed in [42] guarantees the validity of the fluctuation-dissipation relations for n-point correlation functions. Specifically for two point correlation functions, the fluctuation-dissipation relation takes the form which along with the symmetry property of retarded/advanced propagators G R,ab (p) = G A,ba (−p) implies that in thermal equilibrium there is only one independent two-point correlation function. When presenting explicit numerical results, we will therefore focus our attention on the investigation of the spectral function ρ ab (ω, p), given by Besides N-component scalar quantum field theory in d spatial dimensions, we will also be interested in the corresponding classical-statistical field theories, who's dynamics can be formulated in terms of classical Langevin type field equations of motion where η a (x) represents a stochastic Gaussian white noise, with auto-correlation functions By performing the usual Martin-Siggia-Rose-Janssen-de Dominicis path-integral reformulation [45,46], the problem of calculating real-time observables in classical-statistical field theory can be formulated in an analogous fashion as a path integral in Eq. (1), where instead of Eq. (3) the classical contour action S cl C [ϕ,φ] is now given by (see e.g. [28,41]) where is the Rayleigh-Jeans distribution. By explicit comparison with Eq. (3) one finds that the classical contour action S cl C [ϕ,φ] only contains the leading O( 0 ) contributions, which as discussed extensively in the literature [41,[47][48][49][50] amounts to a change of the statistical factor between Eqns. (3) and (12), as well as the absence of the "quantum"φφφϕ interaction term in the classical-statistical field theory. We also note for completeness that the classicalstatistical theory in Eq. (12) is invariant under the symmetry transformation [42] T cl β ϕ a (ω, p) = ϕ a (−ω, p) , T cl βφ a (ω, p) = βω ϕ a (−ω, p) +φ a (−ω, p) , which again guarantees the validity of the classical fluctuation-dissipation (Kubo-Martin-Schwinger) relations for n-point correlation functions.
Due to the fact that the quantum and classical-statistical theories only differ by the presence/absence of the quantum vertex and the change of statistical factors, the real-time functional renormalization group framework allows for an efficient simultaneous discussion of both classical-statistical and quantum field theories. Since in contrast to the quantum field theory, the classical-statistical field theory can be simulated in real-time from first principles by performing real-time lattice simulations [37][38][39][40], the functional renormalization group results obtained in the classical-statistical regime can therefore be directly compared to exact numerical calculations, thus allowing for an important test of the methodology and benchmark of the quality of the underlying approximations.

B. Effective action and flow equation
Starting from the generating functional Z[J,J ] for quantum and classical-statistical field theories, the generating functional for connected correlation functions W [J,J] is given by such that connected one-and two-point correlation functions are determined by and The one-particle irreducible (1PI) effective action is obtained by a Legendre transformation of Eq. (15) with respect to the sources J andJ , for fixed values of the field expectation values φ a ,φ a , i.e.
Even though the effective action contains the full information content about the dynamics of the theory, it is notoriously hard to compute due to the functional integrations in the generating functional. The basic idea of the functional renormalization group approach is therefore to construct the effective action step-by-step, by solving a set of functional differential flow equations which successively integrate out fluctuations at different scales.
In order to construct the functional flow equations we follow standard procedure [51] and introduce a regulator term depending on the flow scale k, so that we replace the original action S[ϕ,φ] in the generating functional by a scale dependent action which includes a generic regulator term of the form Based on these modifications, the effective action Γ k [φ,φ] now depends on the scale k and is explicitly given by Based on a suitable choice of regulator functions R X k,ab (x, y), such that in the limit k → Λ the regulator suppresses all fluctuations, whereas in the limit k → 0 the all regulators vanish and all fluctuations are included, the renormalization group flow interpolates between the classical action S[φ,φ] at some ultra-violet (UV) cutoff scale k → Λ and the full effective action in the infrared, i.e.
We also note that on the Schwinger-Keldysh contour, the various regulators have to satisfy additional constraints to comply with the symmetries of the action of an equilibrium system, as will be discussed in more detail below.
By taking a renormalization group scale (k) derivative of the effective action Γ k [φ,φ] in Eq. (22), we obtain the flow equation for the effective action which upon performing a straightforward set of manipulations can be expressed as where all two-point functions in the last line are understood to be connected. By use of the relations in Eq. (17) we then arrive at the most general form for the flow equation [28] +Ṙ F k,ab (x, y)iF k,ba (y, x) +ṘF k,ab (x, y)iF k,ba (y, x) .

C. Propagators and two-point functions
The flow equation for the effective action (27) is given in terms of scale dependent propagators, which are related to the derivatives of the effective action. Denoting the second functional derivatives of the effective action as the expressions for the various propagators are then given by [27] Similarly, by taking functional derivatives of the propagators in Eq. (29) one obtains the flow equations for n-point correlation functions, which in the end have to be evaluated at the minimum of the effective action. Sinceφ = 0 and Γ φφ = 0 vanish due to discrete symmetries of the effective action, the propagators evaluated at the minimum of the effective action then simplify to 2 Using the fluctuation dissipation relation in Eq. (36) for scale dependent propagators, then implies the following relations between the different two-point functions appearing in the effective action which needs to be satisfied at any scale k.

D. Regulator functions
Even though the detailed choice of regulators is irrelevant if the functional differential flow equation is solved exactly, in practice the hierarchy of flow equations for n-point correlation functions has to be truncated at a finite order making the solution sensitive to the regulator choice. Since finding suitable regulators for real-time calculations turns out to be a rather subtle issue, we will now comment in more detail on the general conditions for the regulator functions in the real-time FRG framework and specify explicit choices below.
Clearly, the most essential property of the regulator is that it suppresses the effect of fluctuations in the real-time path integral. Expressing the regulator matrix R k,ab (ω, p) = R k (ω, p)δ ab for a space-time translation invariant system in Fourier space, as this can e.g. be achieved if the imaginary part of the bi-linear form is positive semi-definite, such that the associated term in the path integral e i∆S k [ϕ,φ] gives rise to an exponential suppression of fluctuations below the renormalization group scale.
Besides its regulating properties, it is also desirable that the introduction of the regulator does not explicitly break the symmetries of the system. Specifically, in our context of realtime dynamics in equilibrium systems, this boils down to the invariance of the regulator term ∆S k [ϕ,φ] under the symmetry transformation in eqn. (7) for quantum and eqn. (14) for classical system, which can be satisfied with Vice versa, if the regulator functions R k (ω, p) are chosen to comply with the above symmetry condition, this also guarantees the validity of the fluctuation dissipation relation for the scale (k) dependent n-point correlation functions, such that for example the fluctuationdissipation relation will automatically be satisfied at all scales.
Specifically for the real-time FRG approach, it is also highly desirable that the introduction of the regulator R k respects the time ordering properties of the retarded/advanced and symmetric propagators in coordinate space, such that e.g. the scale dependent propagator G R k (x, y) remains retarded, i.e. vanishes for space-like separations (x − y) 2 < 0, throughout the entire renormalization group evolution. Vice versa, in momentum space, this condition dictates, that the regulator term does not introduce spurious complex poles of the advanced/retarded propagators, which would result in a violation of causality. Note that there is no analogue of such a causality constraint for Euclidean FRG calculation, indicating the additional difficulties that appear in real-time QFT calculations.
Clearly the simplest possible way to comply with causality, is to employ a frequency independent (purely spatial) regulator acting as an effective mass term, such that following whereas the symmetric regulator functions R F k and RF k vanish identically in this scheme. One particular choice of the regulator function, which has been frequently employed in the literature [53] is However, a purely spatial regulator scheme has the obvious disadvantage that it can not be applied in 0 + 1 dimensions, and moreover is also not particularly suitable for higher dimensional lattice models which feature a discrete set of spatial momenta. We will therefore explore a different possibility, where inspired by the free inverse propagator the regulator takes the form [31] We emphasize that in the above expression µ k (ω, p) and γ k (ω, p) are real-valued even functions of the frequency ω, such that in addition to the real part ∝ µ k (ω, p) which corresponds to an effective mass term, the regulator also features a non-vanishing imaginary part ∝ ωγ k (ω, p), which for γ k (ω, p) > 0 corresponds to an effective damping rate. Specifically, for our 0 + 1 dimensional case study, we will choose such that both µ k (ω) and γ k (ω) diverges as k 2 in the limit k → Λ. Since the regulator diverges for sufficiently large k all fluctuations are suppressed when the renormalization group scale k approaches the UV cut-off scale Λ such that the effective action Γ k=Λ [φ,φ] = S C [φ,φ] is given by the bare action. In analogy to Euclidean FRG calculations [51], this can be easily demonstrated via a saddle point approximation of the path integral and we provide a short discussion in Appendix B for completeness.
With these, the diagramatic representation of the flow equation (27) takes the compact form where -as a novelty of our notation -a green line is shorthand notation for either blue or red and the flow equation is a sum of all allowed color permutations. Notably the introduction of this compact matrix notation is particularly useful when deriving flow equations for higher n-point functions. Since the functional differentiation of the various propagators gives rise to all possible insertions of intermediate propagators, e.g.
allows for an efficient bookkeeping with a drastically reduced number of the diagrams. Based on the diagrammatic shorthand notation, the flow equation of a generic two-point function can be compactly expressed as where also the black lines on the external legs can be either blue or red, depending on the particular two-point function under consideration.

III. EXPLICIT COMPARISON TO PERTURBATION THEORY
Before we proceed with our discussion of the real-time FRG approach, it proves insightful to analyze which set of perturbative contributions are included in the real-time functional renormalization group calculation. Generally, our strategy for this purpose will be to expand the effective action into terms proportional to powers of λ n and then write down separate flow equations for all terms ∆ (n) Γ[φ,φ] to bring them into the form of a total differential such that the integration w.r.t to the scale parameter k becomes trivial.

A. One loop contributions to propagators and vertices
Starting from the FRG flow equation for the effective action in Eq. (27) it is evident, that to one loop order only bare propagators and vertices can appear on the RHS of the flow equation. Explicit evaluation of the (scale dependent) bare propagators yields the following Since the scale (k) dependence only enters through the regulator itself, one finds the following explicit relations for the scale derivatives of the propagators which can be used to integrate the flow equations w.r.t. to k as described below. Similarly, at one loop level all vertices appearing on the RHS of the flow equation are simply given in terms of the bare vertices and take the following explicit form where we denote λ cl = − λ 3N and λ qu = − λ 12N in the following. Specifically, for the two point functions the relevant flow equations then evaluate to where the 'flavor' factor (N + 2) comes from the following contraction of O(N) indices (δ aā δ cc + δ ac δāc + δ ac δāc)δ cc = (N + 2)δ aā .
By comparison with Eqns. (48), one recognizes the RHS as total k-derivatives and the flow equations can be integrated w.r.t to the scale parameter k yielding 3 irrespective of the details of the regulator, as long as the latter ensures the suppression of UV boundary terms and does not introduce violations of causality such that the terms in the last two lines vanish.
Based on the expressions in Eq. (56), one immediately realizes that the only contribution at the one loop level is a manifestly real and local correction, which physically amounts to the familiar one loop mass shift . Hence one concludes that, in the absence of spontaneous symmetry breaking, any non-trivial modifications of the spectral shape only occur starting at the two loop level, and it is therefore important to understand how these are generated within the real-time FRG approach.
Beside the one loop correction to the two point function, we will also need the one loop corrections to the four point functions, which enters the perturbative calculation of the spectral function at the two loop level. Evidently, the one loop corrections to the four point functions can be obtained in an analogous fashion from the flow equation of the four point function  Based on the apparent symmetries of the corresponding diagrams, we can decompose the one-loop corrections to the four point functions according to where the O(N) index structure of the expression is obtained by evaluating the index contraction of bare propagators and vertices according to One where we exploited the symmetries to further compactify the expressions. We note in passing that for the quantum theory, the (tree-level) symmetry relation λ cl = 4λ qu between the local classical and quantum vertices also holds for the non-local vertex functions at the one loop level, i.e. ∆ (1) Γ φφ,φφ one loop level, in both classical and quantum theories.

B. Two loop contributions to propagators
Since the flow equation for the propagators is of one-loop form, we can obtain the two loop contribution in a similar fashion, by using one propagator or respectively one vertex at one-loop order and use bare versions for all other quantities, i.e.
where the black dot denotes the perturbative one-loop vertex and double lines denote the perturbative one-loop propagators given by By inserting the corresponding expressions into the above flow equation, one finds that the contributions to ∂ k ∆ (2) Γ φφ k (xx) fall in two topologically different categories, given by "double bubble" and "sunset" diagrams respectively. By performing a straightforward but cumbersome set of manipulations, the contributions from diagrams with a one-loop propagator, can be expressed as By combining this contribution with a corresponding set of double-bubble diagrams with a one-loop vertex, which upon further manipulations and dropping off vanishing terms can be compactly expressed in the form one finds that the sum of two contributions yields a total derivative w.r.t k, such that yielding Similar to the one loop correction ∆ (1) Γ φφ k (xx), this term is manifestly real and local providing the two loop correction to the mass shift. However, there is also the contribution from the sunset diagrams which can be compactly expressed as Clearly this contribution to the effective action is non-local and posses a non-vanishing imaginary part, which describes the collisional broadening of the spectral function. We further emphasize, that in the real-time FRG framework the sunset contribution arises entirely due to the one-loop vertex correction, indicating the importance of including nonlocal vertex structures into the truncation of the real-time FRG flow equations. By including these non-local vertex structures, as in Eq. (58), one is then able to derive the two-loop perturbative contributions to the damping rate [54], as discussed in Appendix C.

IV. NON-TRIVIAL TRUNCATIONS FOR REAL-TIME CALCULATIONS
Based on our perturbative analysis of the flow equations in the preceding section, we conclude that a two loop complete truncation scheme for the two-point function is necessary to describe the collisional broadening of the spectral function in the symmetric phase. We have also observed, that a two loop complete truncation scheme for the two point function necessarily has to include a non-local four field interaction (e.g. generated at the one loop level), indicating the the local potential approximation that is commonly used in Euclidean FRG calculations is insufficient for the purpose of real time calculations. Now in order to devise a more suitable truncation scheme, we first note that we can generally express the scale dependent effective action in a vertex expansion as Since in the symmetric phase only the n−even terms contribute, a two-loop complete expansion can be achieved by truncating the vertex expansion at the level of the four-point function (Q = 4) keeping only two-and four-point functions. Hence, the simplest possible two-loop complete expansion scheme is given by where the above truncation only takes vertices into account that can be generated at oneloop level, i.e. the (φφφφ), (φφφφ) and (φφ,φφ) vertices vanish. With regards to the nonvanishing vertex functions, we employ a generalization of the one-loop result in Eq. (58) as our ansatz with scale dependent vertex functions v X Y,k (xx). While at the one-loop level the diagonal v diag Y,k (xx) and off-diagonal v of f Y,k (xx) vertex functions are simply related by a factor of (N + 4)/2, this is not the case beyond one-loop and we generally have to distinguish between diagonal and off-diagonal vertex functions. Based on the symmetries of the effective action for an equilibrium system, the above vertex functions satisfy the following symmetry as well as the fluctuation dissipation relatioñ v X anom,k (p) = n eff (p 0 ) ṽ X cl,R,k (p) −ṽ X cl,A,k (p) .
which upon inserting the explicit expressions for the non-local vertex functions in Eq. (77-79), gives rise to the following structure of the flow equations featuring a two-loop structure of sunset diagrams in the first and second column and double bubble diagrams in the thrid column. By introducing the following short-hand notation for the one-loop integrals the flow-equations for the two-point functions then take the form where we dropped acausal contributions proportional to B R/A (0).

B. Vertex flow
Evidently, to close the system of equations we still need expressions for the vertex functions. In the following we will compare two different truncations.

One-loop vertex functions
We start by using the one-loop expressions of the vertex functions with self-consistently determined propagators. Explicitly, for the four-point functions eqs. (77-79) the pertubative one-loop expressions determined at each step of the renormalization group evolution take such that diagonal and off-diagonal vertex functions differ only by their corresponding flavor factors.
and we will use the following relation to project the flow equation onto the diagonal and off-diagonal vertex functions By performing the projection of the flow equation according to Eq. (93), the flow-equation for the projected vertex function then takes the following diagrammatic form 4 where it is important to state that Eq.
where there is no distinction between diagonal and off-diagonal index structures for the By separating the different index structures, one obtains the final result such that in contrast to Eq. (95), the flow of the vertex function in the large N limit is local in momentum space, in the sense that all vertex functions in Eq. (97) are evaluated at the same momentum. We further note that for special choices of the regulator function one can show that the right-hand side of the flow equation simplifies to a total differential [26,27,55,56] which can be solved directly by separation of variables eventually yielding the familiar result of the 2PI 1/N expansion to next-to-leading order [27,57] v diag cl (p) = which corresponds to an infinite resummation of one-loop bubble chains. Based on this analysis, we therefore conclude that the above truncation of the real-time FRG flow equations not only encompasses the correct two-loop perturbative behavior of the spectral function for generic N, but also includes all contributions up to next-to-leading order of the 2PI 1/N expansion in the large N limit. We further note that the interplay of the 2PI approach and the FRG in Euclidean time has been explored in the literature, e.g. the use of 2PI truncations in FRG calculations [58] or the use of the FRG to perform the complicated renormalization of 2PI calculations [59,60], and we expect the interplay of the approaches to be similarly useful for real-time calculations.

V. NUMERICAL IMPLEMENTATION
Due to to the nested one-loop structure of the real-time FRG flow equations, it is benefitial to employ (pseudo-)spectral methods to solve the functional differential equations numerically. We have explored two different schemes, with the first one based on a straightforward lattice discretization of frequencies, where for an arbitrary function G(ω), we store the information at a discrete set of frequencies ω i Similarly, the corresponding function G(t) in coordinate space is obtained at a discrete set with the spacing a ω = 1/a t adjusted to properly resolve the propagators at all relevant scales. Numerically, we keep track of the expansion coefficients g (t) k and g (ω) k , as well as the values of the propagators at times t = x k a t and frequencies ω = x k a ω , where x i are the Gauss Hermite points, for which G(t) and G(ω) are simply given by , where in practice we employ Gauss-Hermite quadrature, such that where w i = 1 N ψ 2 N−1 (x i ) are the corresponding quadrature weights. Based on the following relations between the expansion coefficients g (t) it is then straightforward to perform Fourier transformations, in order to efficiently calculate the right hand side of the flow equations for the two-point functions. Similarly to the FFT method, we employ a Gauss Hermite quadrature when evaluating the integrals on the right hand side of the flow equation for the four point function, and for simplicity we resort to a forward Euler scheme when solving the FRG flow equations.

VI. BENCHMARKS AND CASE STUDIES
We will benchmark the method at the example of the anharmonic oscillator, which cor- as well as on the number of field components N, which we will set to N = 1. However, it is well known that in the classical-statistical theory the coupling constant and temperature dependence are related, such that upon performing a re-scaling of the classical field equations of motion with the dependence on λ m 3 can be eliminated from the classical-statistical field theory. Of course, this is not the case in the corresponding quantum theory, such that for fixed values of the (dimensionless) thermal interaction strength λ βm 4 , the dimensionless parameter λ/m 3 effectively describes the quantum interaction strength, with λ/m 3 = 0 corresponding to the classical-statistical limit.
Before we present a series of results of real-time FRG calculations for classical and quantum systems, some sanity checks are in order, to instil confidence in our numerics. Evidently, a first important check is to compare the results obtained by different numerical methods, as shown in Fig. 1, where we compare results for the spectral functions computed using the Discrete Fourier Grid and Gauss-Hermite Representation methods described in Sec. V.
Excellent agreement between the two approaches is observed when sufficiently many discretization points are taken into account, indicating the common convergence to the correct result. Despite significant differences in the underlying implementation, we also find that within our implementation both approaches have comparable execution times on the order of 40s for calculations with one-loop vertices and 20m for calculations with self-consistent vertices, when employing the same number of N = 1024 discretization points. 5 However, we note that the convergence of the results generally appears to be somewhat better for the Discrete Fourier Grid.
Besides the convergence of the discrete representation of functions, another important sanity check for any FRG calculation is the comparison of results using different regulator schemes as the physical results in the infrared should be independent of the regulator choice.

A. Benchmarks in the classical-statistical limit
We begin with the study of classical-statistical dissipative systems, which in accordance with our discussion in Sec. II can be compared against exact numerical results from classicalstatistical simulations [37][38][39][40]. With regards to the classical-statistical simulations, we follow the methodology of previous works [39,40] and simulate the time evolution of an ensemble of N samples = 128 independent realizations by solving the discretized classical evolution equations using an Euler-Maruyama scheme with step width m∆t = 0.01. We then calculate the classical-statistical equilibrium spectral function from the un-equal time correlation function where . class.stat. denotes the average over the classical-statistical ensemble. Subsequently, we perform a discrete sine transformation to obtain the classical-statistical spectral function ρ cs (ω) in frequency space, which can then be compared directly with the classical-statistical real-time FRG calculations. So far we have investigated the spectral functions for a strongly dissipative anharmonic oscillator (γ/βm = 0.5) and we will now study the effect of reducing the dissipative coupling to the heat bath. Before we proceed, we briefly note that the effect of the dissipative coupling γ/βm is somewhat peculiar in 0+1d as, in contrast to higher dimensional theories, we expect to recover a discrete spectrum in the limit of a closed system γ/βm → 0, and the behavior could be qualitatively different in higher dimensions. Fig. 6 shows a comparison of spectral functions obtained by classical-statistical simulations and FRG calculations with one-loop vertices in the classical limit. We observe that the deviations from the classical-statistical results are increasing when we decrease the dissipative coupling γ/βm, as may be expected due to the fact that the longer lived excitations can interact with each other over a larger time scale. While for γ/βm = 0.2 the FRG calculation with one-loop vertex functions still provides a rather accurate description of the classical statistical result, the agreement becomes gradually worse with decreasing γ/βm. Especially for very small values of the dissipative coupling γ/βm < 0.05, the quasi-particle peak of the spectral function splits into a double peak, which is clearly not observed in the classical-statistical data. Similarly, also the strength of the 1 ↔ 3 resonance peak around located ω ∼ 3m generally tends to be over-estimated by the FRG calculations.
We note that the FRG calculations with scale-dependent vertex functions also become unstable for small dissipative coupling γ/βm than shown in Fig.6. In order to further investigate the instability of the self-constistent FRG method at large couplings (i.e. small dissipative couplings), we can now look at the momentum dependence of the classical, retarded vertex function v cl,R . Fig. 7 shows a comparison of the self-consistently determined data for v cl,R with the perturbative one-loop result for the same parameters as in Fig. 3.
We recognize that for small frequencies the self-consistent vertex function behaves as expected, as both real-and imaginary parts of the vertex functions are suppressed compared to the perturbative result. However, for larger frequencies we find large enhancements of the self-consistently determined vertex function over the perturbative result. Due to the rather complicated structure of the flow equation for the four point function, we are currently not sure about the exact origins of these spurious enhancements, which may be connected to the particular situation in 0+1 dimension and we hope that our procedure will work out better in higher dimensions. Besides additional studies of this behavior, it would also be useful to extract the corresponding vertex functions directly from classical-statistical simulations, which is clearly beyond the scope of this work but could potentially be achieved along the lines of [61].

B. Spectral functions in the quantum theory
Now, that we have benchmarked and assessed the range of applicability of the method at the hand of the classical-statistical theory, we can continue to investigate spectral functions in the corresponding quantum theory. A compact summary of our results is provide in Fig. 8, where we shows a comparison of spectral functions from the FRG with the one loop vertex with results from classical-statistical simulations and perturbative calculations for different values of the thermal and quantum coupling strength. We see that for small coupling all methods agree very well. When we increase the coupling we see a second peak emerging at roughly 3m due to the 1 ↔ 3 processes in the one-loop correction to the four-point function. As there is no vertex correction at the perturbative one-loop level also the oneloop spectral function fails to capture this feature. When we increase the coupling either by increasing the dimensionless combination of coupling and temperature or by driving the system more towards a strongly coupled quantum system we see that perturbation theory becomes unreliable rather quickly as there are large differences between the LO and NLO results. Specifically for large couplings, the perturbative spectral functions at the two-loop level show additional spurious peaks or may even become negative. Conversely, the FRG results remain much more well behaved throughout the observed parameter range, except perhaps for the largest combination of couplings shown in the bottom right panel. We benchmarked our methods at the 0 + 1d example of the anharmonic oscillator where we compared results from perturbation theory, our FRG calulations with both truncation schemes and results from classical-statistical simulations. Since the real-time FRG framework can be formulated in essentially the same way for classical and quantum theories, the comparison to exact results from classical-statistical simulations proved to be an important benchmark to asses the range of applicability and performance of the method.
Overall we find that the real-time FRG is able to reproduce the classical-statistical results much better than the perturbative calculations. Still, we find that a larger couplings also the FRG fails to reproduce the correct results. In case of the one-loop vertex trunctation this is connected to the perturbative origin of the vertex. In case of the fully self-consistent FRG calculations we found a spurious enhancement of the vertex functions at large momenta leading to a break down of the method for large couplings. The origin of this enhancement is still unclear but could particularly be a problem of the 0+1 dimensional theory. Another possible cause is the omission of six-point functions in our truncation.
Similar to the requirement of a scheme with two-loop complete propagators to reproduce the broadening of the spectral functions we might also need a truncation with two-loop complete four-point functions to be able to correctly renormalize the vertices.
While the formalism described in this work has been derived for N component scalar field theories in d + 1 dimensions, so far our numerical investigations have been limited to the 0+1 dimensional theory. Clearly a next important step would be to generalize our numerical investigations to higher dimensional systems, especially in 3+1d. Evidently, the comparison of real-time FRG calculations in the classical limit to classical-statistical simulations proved extremely insightful, and should also be pursued for studies in higher dimensions. We also expect that in higher dimensions it should be possible to take the limit of vanishing dissipative coupling γ/βm → 0, which would further allow to compare real-time FRG calculations in the quantum theory to results from lattice Monte-Carlo simulations and/or analytically continued FRG calculations in Euclidean space-time. Eventually, we want to generalize our framework to include fermions opening the possibility of applying our framework to low energy effective theories of QCD like e.g. the Quark-Meson model.

Appendix A: Conventions for real-time propagators
Below we summarize our conventions and also note some useful relations among the various real-time propagators. Based on the operator definitions one has along with We also note for convenience the following relations between real and imaginary parts of the various correlation functions 6 where in the quantum case we have n qu ef f = n BE + 1/2 with n BE (p) = 1/(e βp 0 + 1) being the Bose-Einstein distribution, such that In the classical case we find to be the Rayleigh-Jeans distribution. Our convention for the Fourier transformation reads We further note the following symmetry relations as well as the various relations for the real and imaginary parts Re(ρ(p)) = 0 , Re(F (p)) = +2n ef f (p 0 )ImG R (p) = −2n ef f (p 0 )ImG A (p) Im(F (p)) = 0 , along with where we used the relations θ(x) = 1 2πi ∞ −∞ dω 1 ω−iǫ e iωx along with 1 ω−iǫ = P.V. 1 ω + iπδ(ω), such that the above idenitities follow directly from the respective time orderings.
Appendix B: Evaluation of the real-time effective action in the limit k → Λ Starting from the definition of the effective action in Eq. (22), it is convenient to perform a field shift ϕ → φ + ϕ in the functional integration to separate off the contribution from the classical action. By exploiting the equations of motion to re-express the appearance of the sources J,J in terms of derivatives of the effective action, Eq. (22) can then be recast into a functional integro-differential equation for the effective action with ∆Z k [φ,φ] given by the functional where [DϕDφ] k we denotes the regulated path integral which for a Gaussian regulator can be evaluated explicitly. Expressing the functional integrations in Fourier space, one finds (up to irrelevant pre-factors) Evaluating the functional integral explicitly according to one finds that in the limit k → Λ, the relevant factors characterizing the variations w.r.t.
the sources j,j are inversely proportional to regulators, such that and the functional becomes independent of the sources j,j in the vicinity of j =j = 0 where derivatives are to be evaluated. One concludes, that in the limit k → Λ the effective action in Eq. (B1) does not receive any additional contributions from the path integral in Eq. (B2) and thus reduces to Appendix C: Perturbative contribution to the damping rate Here we will evaluate the perturbative contribution to the damping rate, which is also useful to understand the differences between classical and quantum statistical processes and in establishing the comparison to the literature that is largely based on analytic continuations of Euclidean calculations. Based on Eq. (74) the perturbative contribution is obtained by evaluating the RHS with free propagators, which take the following form in momentum space Spectral function and statistical function are then given by allowing us to proceed directly with the evaluation of diagrams. By expressing all retarded/advanced propagators in terms the spectral function using Eq. (A2), the retarded self-energy, defined as can be expressed in the form Since spectral and statistical correlation functions are purely real in coordinate space, and have well defined real and imaginary parts in momentum space, the real and imaginary parts of the can be readily evaluated, by use of relations as in (A9), which follow directly from the properties of the Heavyside step function. Since likewise the real-part can be re-constructed from Kramers-Kronig type relations, we will focus on the the imaginary part, which can be directly evaluated as Specifically in the limit γ → 0 of non-dissipative systems, the energy integrations can be performed using and the result can be compactly expressed in the form where the summation over s i = ± collects the positive and negative frequency contributions (2π) d as well as r = p − k − q to lighten the notation. By further symmetrizing the integrand of ImΣ R cl (p), the expressions can be re-cast in the form ImΣ R qu (p) = − 3π 2 (N + 2) 3 λ cl λ qu q, k s k sqsr where the pre-factors of the two terms are equal except for the different appearances of the coupling constants λ cl and λ qu . We will now concentrate on a quantum theory where we can set n ef f (E) = n(E) + 1/2 -with the Bose-Einstein distribution n(E) -and exploit the relation λ qu = λ cl /4 between the classical and quantum (tree level) vertices, it is straightforward to show that the above terms can be combined in the following way s q s k s r 4 + s r n(E k ) + 1 2 n(E q ) + 1 2 + s q n(E k ) + 1 2 n(E r ) + 1 2 + s k n(E q ) + 1 2 n(E r ) + 1 2 = (C16) n(E k ) + 1 + s k 2 n(E q ) + 1 + s q 2 n(E r ) + 1 + s r 2 − n(E k ) + 1 − s k 2 n(E q ) + 1 − s q 2 n(E r ) + 1 − s r 2 , which contain the usual quantum statistical factors for in/out-going particles in a scattering process. Collecting everything, the imaginary part of the retarded self-energy ImΣ R (p) = ImΣ R cl (p) + ImΣ R qu (p) can then be compactly expressed as in agreement with the standard result in ref. [54]. Vice versa, in the classical-statistical theory the contribution proportional to λ qu vanishes identically, and the occupancy factors n(E p ) + 1/2 are to be replaced by the Rayleigh-Jeans distribution n cl (E p ) = 1/βE p , such that ImΣ R cl (p) = − 3π 2 (N + 2) 3 λ 2 cl q, k s k sqsr [s r n cl (E k )n cl (E q ) + s q n cl (E k )n cl (E r ) + s k n cl (E q )n cl (E r )] , now yielding the classical statistical factors for in/out-going particles in a scattering process.
We have thus verified explicitly that with a suitable truncation which properly accounts for the non-local vertex structure generated at the one loop level, the real-time FRG approach correctly captures the collisional broadening of the spectral function at the two loop level.
Since for βE p ≪ 1 the statistical factors n(E p ) + 1+sp 2 ≈ 1/βE p agree approximately the classical statistical theory is expected to accurately capture the relevant contributions of excitations with energies much smaller then the temperature. However, one crucial difference is that the classical statistical theory only allows for interactions between physically occupied excitations of the system. Due to the statistical factors, the classical-statistical result behaves as ImΣ R cl (p) ∼ T 3 , such that in the limit T → 0 where classically no states are physically occupied, all contributions to the self-energy vanish identically, which is of course not the case in the corresponding quantum theory.

Appendix D: Fluctuation-dissipation relation
Below we demonstrate explicitly, that the flow equations of the two-point functions satisfy the relation A quick calculation with the one-loop forms for the vertices show that at one-loop level we have v an (p) = n eff (p 0 )(v cl,R,k (p) − v cl,A,k (p)) .
In the case of a classical theory, we have v qu,R/A (p) = 0, i.e. we can just drop the −1/4 in Eq. (D7). By plugging in the Rayleigh-Jeans distribution we again find Eq. (D8).
Going through the calculation in reverse order proofs that the existence of a general-