Ghost spectral function from the spectral Dyson-Schwinger equation

We compute the ghost spectral function in Yang-Mills theory by solving the corresponding Dyson-Schwinger equation for a given input gluon spectral function. The results encompass both scaling and decoupling solutions for the gluon propagator input. The resulting ghost spectral function displays a particle peak at vanishing momentum and a negative scattering spectrum, whose infrared and ultraviolet tails are obtained analytically. The ghost dressing function is computed in the entire complex plane, and its salient features are identified and discussed.


I. INTRODUCTION
The complete access to the hadronic bound state and resonance structure, as well as to the non-perturbative dynamics of QCD at finite temperature and density, requires the computation of timelike correlation functions. In functional approaches, such as the functional renormalisation group (fRG) and Dyson-Schwinger equations (DSEs), the respective non-perturbative computations are carried out numerically, and hence, demand a numerical approach to timelike correlation functions.
Recently, the spectral DSE approach has been put forth [1], based on the Källén-Lehmann (KL) representation of correlation functions in terms of spectral functions. In particular, in [1] the general properties of this novel approach, including a consistent spectral renormalisation procedure, have been expounded and applied in the context of a scalar φ 4 -theory in 2+1 dimensions.
As a first step towards a full QCD treatment, we use this spectral DSE for the ghost-gluon system in Yang-Mills theory. In recent years, ghost and gluon spectral functions have been reconstructed from numerical data of Euclidean ghost and gluon propagators, see e.g. [2][3][4][5][6]. Also direct computations have been put forward, either perturbatively, e.g. [7,8], with non-perturbative analytically continued DSEs [9,10], or in a spectral approach [11]. While these direct computations unravel highly interesting structures, they are still inconclusive.
Specifically, in the gluon DSE one has to deal with a rather complicated diagrammatic representation with many non-trivial ingredients, ranging from the presence of one-and two-loop diagrams, the dependence on ghost and gluon propagators, as well as several vertices with complicated momentum dependence. Consequently, some properties of the gluon spectral functions, and in particular the potential presence and location of complex conjugate poles, are rather unstable under small variations of the vertices involved; for a detailed recent discussion, see [10].
Instead, the ghost DSE, see Figure 1, requires far less non-trivial input: apart from the ghost propagator, it depends only on the gluon propagator or, equivalently, the gluon spectral function, as well as the ghost-gluon vertex. The latter is protected by non-renormalisation, and hence shows a very mild momentum-dependence. Accordingly, in the present work we approximate this vertex by its classical counterpart.
This leaves us with a rather stable set-up: the spectral ghost DSE is solved on the basis of given input gluon spectral functions, obtained by appropriately modifying the result of [4], which was reconstructed under the assumption of a KL representation of the gluon. We also test the stability of the result under a variation of the input by tuning the whole family of scaling and decoupling solutions.
This work is organised as follows: In Section II we discuss spectral properties of Yang-Mills theory. In Section III, the spectral ghost DSE is set up, and the input gluon spectral function is discussed. We present our results for the ghost propagator and spectral function in Section IV, and discuss our findings in Section V. The appendix includes some technical details.

II. YANG-MILLS THEORY AND THE SPECTRAL REPRESENTATION
We consider 3 + 1-dimensional Yang-Mills theory with three colours, N c = 3, in the Landau gauge. The gaugefixed classical action including the ghost action reads, where ξ denotes the gauge fixing parameter and The Landau gauge is achieved for ξ = 0. Note that in (1) we have chosen a positive dispersion for the ghost. The field strength, F µν , and covariant derivative, D µ , in the adjoint representation are given by with the structure constants f abc of SU (3). Functional relations such as the flow equations in the fRG or the DSEs are one-and two-loop exact relations for the full correlation functions. For reviews see [12][13][14][15][16][17] (fRG) and [18][19][20][21][22][23][24]  functions φφ c , the propagators of the theory. They read where φ = (A µ , c,c), and the tensor T φiφj (p) carries the Lorenz and gauge group tensor structure. The scalar parts of the propagators are given by In the Landau gauge, the gluon propagator is transverse, where Π ⊥ µν (p) = δ µν − p µ p ν /p 2 denotes the standard transverse projection operator. The scalar part of the gluon propagator reads with the gluon dressing function 1/Z A (p). Similarly, for the ghost we have T ab cc = δ ab , and the scalar part reads with the ghost dressing function 1/Z c (p). In what follows we will compute (5) for time-like momenta. In fact, as we will see below, our results permit the evaluation of G c (p) for general complex momenta. Extensions of correlation functions to the complex plane are particularly interesting, in view of their relevance for the self-consistent treatment of bound-state problems, see, e.g. [25][26][27].
If the KL spectral representation [28,29] is applicable, a propagator G can be recast in terms of its spectral function ρ, The spectral function naturally arises as the set of nonanalyticities of the propagator in the complex momentum FIG. 2. Diagrammatic notation used throughout this work: Lines stand for full propagators, small black dots stand for classical vertices, and larger blue dots stand for full vertices.
plane. If (6) holds, the non-analyticities are restricted to the (linear) real momentum axis. Equation (6) leads to the following inverse relation between spectral function and the retarded propagator, where ω is a real frequency. This formulation allows us to work only with the frequency argument and set the spatial momentum to zero in practice, since the full phasespace can be restored from Lorentz invariance. Hence, for the remainder of this work, | p| will be dropped. Now we detail the above generic spectral representation (6) for the Yang-Mills propagators. We will consider the entire family (decoupling/massive and scaling) of potential infrared solutions, for detailed discussions see [30][31][32][33]. In terms of the ghost and gluon propagator inverse dressing functions Z c (p) and Z A (p), a generic decoupling behaviour is characterised by with a finite Z c = Z c (p = 0). Formally, the ghost propagator is expected to obey the KL-representation [34]. Also recent reconstructions show no signs of a violation of this property [5,6]. The ghost spectral function must exhibit a single particle peak at vanishing spectral value, with residue 1/Z c . In addition, a continuous scattering tail is expected to show up in the spectral function via the logarithmic branch cut. This leads us to the general form of the ghost spectral function, whereρ c denotes the continuous tail of the spectral function and δ(ω)/ω has to be understood as a limiting process δ(ω − m)/ω with m → 0 + . Inserting (9) in (6) leads us to a spectral representation for the dressing function, In the case where the spectral function can be normalised by solely integrating it over the whole branch cut, the normalisation is given by the value of the inverse dressing function at infinity. A detailed derivation of this is given in Appendix A. Since the inverse ghost dressing tends to zero for large momenta, the spectral function obeys dλ π λ ρ c (λ) = 0 .
Acc (p, q), see (13), data taken from [4]. The dressing is shown at the symmetric point p 2 = q 2 = (p + q) 2 for scaling and lattice-type decoupling solution, more details can be found in [4].

III. THE SPECTRAL GHOST DSE IN YANG-MILLS THEORY
In this section, we use the preparations in Section II to set up the spectral DSE for the ghost spectral function ρ c . The DSE for the ghost propagator in Yang-Mills theory is depicted in Figure 1, while the diagrammatic notation employed is summarised in Figure 2. As discussed before, its only input is the scalar part G A of the gluon propagator and the full ghost-gluon vertex. The latter consists of two tensor structures, see e.g. [33], with incoming gluon momentum p and anti-ghost momentum q, and we have dropped the momentum conserving δ-function.
The ghost-gluon vertex in (13) contains two independent vertex dressings, λ (cl) (classical tensor structure) and λ (nc) (non-classical). The non-classical dressing is proportional to the gluon momentum and hence drops out of the ghost DSE due to the transversality of the Landau gauge gluon propagator.
The classical dressing is subject to Taylor's nonrenormalisation theorem, and has a very mild momentum dependence, see e.g. [33,[47][48][49][50][51][52][53][54], see Figure 3 with the data from [33]. In Figure 3 we depict the dressing λ (nc) at the symmetric point p 2 = q 2 = (p + q) 2 for both, the scaling solution as well as the lattice-type decoupling solution. For further explanations we refer to the detailed discussion of [33]. Accordingly, we consider the approximation λ (cl) where g is the gauge coupling at the renormalisation point µ RG . Within the MOM-type scheme that we employ, the dressing functions acquire at µ RG their tree-level value, i.e., Z A (µ RG ) = Z c (µ RG ) = 1. We emphasise that our approach is by no means restricted to classical vertices only: quantum corrections may be duly accounted for, as long as the momentum loops involved can by integrated analytically. Especially, upon construction of spectral representations for higher n-point-functions, see e.g., [55][56][57][58][59][60][61][62][63][64][65][66][67], fully dressed vertices of general form can be included.
In summary, the ghost gap equation in Figure 1 can be written as with the renormalisation constantZ 3 associated with the ghost field. Similarly, the classical ghost-gluon vertex in Figure 1 contains the respective renormalisation constant, S µ Acc (p, q) = −Z 1 igf abc p µ . The (scalar) ghost selfenergy Σc c (p) with the approximation (14) is then given by, where C A = N c is the eigenvalue of the second Casimir operator of the colour group SU(N c ) in the adjoint representation. The ghost DSE of (15) can be straightforwardly rewritten in terms of the ghost dressing function, Now we write (17) in terms of spectral loops by using (6) for ghost and gluon propagators. This leads us to Σc c (p) = g 2 N c λ1,λ2 with ρ A and ρ c the gluon and ghost spectral functions, respectively. Note that the order of spectral and momentum integration have been interchanged, implicitly assuming the finiteness of the dimensionally regularised expression.

A. Spectral renormalisation
The momentum integral in (18) involves two massive propagators with spectral masses λ 1 and λ 2 . It is readily evaluated upon introduction of Feynman parameters and using dimensional regularisation in d = 4 − 2 dimensions. Calculational details and the resulting expression are given in Appendix B. After the integration over the loop momentum is carried out, we are left with two spectral integrals. They suffer from the (same) logarithmic divergence as the momentum integral, even if simply dropping the 1/ -term, that arises from the loop integral. This is a generic feature in the spectral DSE, for a thorough discussion see [1]. Moreover, the potentially divergent terms in Σc c (p) are proportional to p 2 . Such a divergence can be cured by aZ 3 , that is proportional to Σc c /p 2 , evaluated at some RG-scale µ RG . We use a standard renormalisation condition for the (inverse) dressing function, This renormalisation condition is implemented by the respective choice of the renormalisation constantZ 3 as In accordance with (14), (19) is augmented withZ 1 → 1. For a detailed discussion of self-consistent MOM-type RG-conditions for DSEs, see [69]. Eventually, this leads us to the renormalised DSE for the ghost dressing, The explicit choice of the renormalisation conditionZ c will be discussed in Section IV.

B. Iterative solution
The renormalised DSE in (20) can be evaluated analytically for general complex frequencies. For the extraction of the spectral function with (7) we choose p 0 = −i(ω+iε) with ε → 0. This leads us to where, in a slight abuse of notation, we define Σc c (ω) = Σc c (−i(ω + i0 + )). Note that the 1/(ω + i ) 2 term does not trigger a δ-function, as the self energy is proportional to (ω+i ) 2 for small ω. This allows us to use 1/(ω + i ) 2 → 1/ω 2 , as is done in (21). We also check at every iteration step that no further poles are generated, and hence (21) holds true. The explicit spectral integral expression for Σc c (ω) and its renormalised counter part can be found in Appendix B. The remaining finite spectral integrals have to be computed numerically, and the spectral function ρ c (ω) is given with (7) as .
The substitution of (21) into (22) allows us to compute the ghost spectral function as well as the propagator for Euclidean and time-like frequencies, for a given gluon spectral function ρ A . This input is discussed in Section III C. The ghost spectral function is then computed within the following iteration procedure, discussed in detail in [1], and briefly described here.
The ghost spectral function ρ (i) c , obtained after the ith iteration step, and the given input ρ A , are inserted into the spectral integral form of Σc c (p), on the righthand side of (21). Then, by means of (22), we arrive at the (i + 1)-th spectral function, ρ (i+1) c . This iteration is repeated until convergence has been reached.
The iteration commences with an initial guess for ρ c . Here we choose the canonical choice of its "classical" spectral function, i.e., a massless pole with residue one, This particular choice leads to a stable and rapid convergence of the iteration procedure. For a discussion of the numerics and their convergence properties, see Appendix C and Figure 10.

C. Gluon spectral function
As already mentioned, the input gluon spectral function is taken from [4], where a spectral reconstruction of the Yang-Mills gluon propagator fRG data from [33] has been performed. In both scaling and decoupling scenarios, the infrared behaviour of the gluon spectral function (assuming the validity of the KL representation) can be inferred from the respective infrared scaling of the gluon propagator in (4). More details can be found in [4].
For the entire family of solutions, the deep infrared limit with p → 0 is parametrised by [4], with a constant Z (IR) A . The scaling coefficient κ takes values in the range 1/2 < κ < 1, and (a) Gluon spectral functions (scaling and decoupling), see (28), based on the reconstruction of the scaling spectral function in [4] (red-dashed line). The spectral functions differ only in the infrared, shown in the inset. where the hatted dimensionless quantities in (24a) all future expressions have been rescaled with the appropriate powers of Λ QCD , e.g.,p 2 = p 2 /Λ 2 QCD . For γ G = 0, the gluon propagator in (24a) reduces to the scaling propagator.
The lattice-type propagator is obtained for a γ (lat) G that is close to the maximal one compatible with infrared QCD in the Landau gauge. The parameters (γ G ,m 2 gap ) characterise the one-parameter family of solutions. Indeed, the actual solutions in [33] are well approximated by using the functional form of the scaling solution but with the argument of (24b), and an appropriate tuning of γ G . We shall exploit this property for constructing a simple one-parameter family of gluon spectral functions, using the scaling one, ρ (dec) (ω), as our point of departure.
For completeness, we note that, for p → 0, the respective (dimensionful) ghost propagator is given by In the deep infrared, ρ (scal) (ω) is determined from (24a) and (24b), setting γ G = 0. Specifically, for ω → 0 + , we obtainρ which corresponds to the infrared tail of the full spectral function reconstructed [4], depicted in Figure 4. Similarly, in the case of the decoupling-type solutions, we arrive atρ While (26) and (27) describe the different behaviour of the scaling and decoupling spectral functions in the deep infrared, for larger spectral values the two sets of spectral functions coincide. This regime is approximately bounded from below by the first zero, ω 0 , of the scaling spectral function, shown in Figure 4, with ω 0 ≈ 0.78. A simple interpolation to the decoupling solution, based on the scaling spectral function in [4], is therefore given by with dictated by the Oehme-Zimmermann superconvergence relation for the gluon spectral function, With (28b) the total spectral weight of ρ (ω). Hence, given that the latter satisfies (29), so does the former. The scaling spectral function reconstructed in [4], satisfies (29) analytically for → 0 + in (7). In the present work we take a small ≈ 10 −7 leading to (a) Ghost spectral function: direct computation by iteration with the spectral DSE (18) and the gluon spectral function from Figure 4 for G . The inlay also indictaes the δ-function contribution in the origin, indicated by an arrow. Its amplitude is given by the value of corresponding Euclidean dressing function 1/Zc(p) at p = 0. The squares show our best fit, comp. Section IV B. for all spectral functions. For χ = 0, we get back the scaling solution with γ G = 0. The lattice gluon is achieved via (28a) for χ (lat) = 3 4 GeV 2 with Z (lat) χ = 1.86. We emphasise that fully quantitative gluon spectral functions ρ (dec) A may be achieved by means of reconstructions. While possible, this is beyond the scope of the present work. Note also that the simple analytic spectral functions ρ A (ω, χ) give semi-quantitative results for the gluon propagators, while at the same time allowing for an analytic access to the relative changes.

IV. RESULTS
With the preparation of the previous sections we now compute the ghost spectral function. The ghost DSE is solved for the three different input decoupling gluon spectral functions in Figure 4  A , we tune the mass parameter χ in (28a) such that we best agree with the lattice results from [68]. We pair each of our input G A 's with a gluon propagator from the family of self-consistent YM solutions of [33], indicated by dashed lines in the right panel of Figure 4. For G (lat) A (blue curve), we chose the solution which also matches the lattice results from [68] best. The green and yellow curves are matched with the respective solution with the same G A (0).
The renormalisation conditionZ c is now chosen such that the value of the ghost dressing function 1/Z c (p) matches that of [33] at the RG scale µ RG = 20 MeV for the respective gluon propagator. This IR renormali-sation procedure is necessary in order to compensate for the lack of self-consistency when considering the ghost DSE with fixed gluon input. The strong coupling constant is fixed to α s = 0.26.
In Figure 5a and Figure 6 we show the respective ghost spectral functions. All spectral functions shows a positive particle peak at vanishing momentum, constituted by a delta distribution. The magnitude of the corresponding residue, i.e. the particle peaks amplitude, rises with decreasing G A (0), and the residues positivity reflects the chosen positive classical dispersion of the ghost. The spectral function also has a negative scattering spectrum starting at vanishing frequency. For decreasing G A (0), one gradually approaches the scaling solution, and the spectral weight increases drastically. This also mirrors the increasing amplitude of the particle peak, which is enforced by the Oehme-Zimmermann-type superconvergence relation (11) for the ghost, for more details see Appendix A and [37]. This also leads to the known UVasymptotics for the ghost spectral function, with the ghost anomalous dimension γ c , and the UV wave function renormalisation Z (UV) c . The gluon spectral function for G (lat) A represents the lattice-type case, see Figure 4b. The respective lattice data for the ghost propagator is depicted in Figure 5b, and confirms the semi-quantitative nature of the classical vertex approximation in the ghost DSE. For smaller G A (0) → 0, the gluon propagator approaches the scaling solution. This entails, that also the ghost propagators approaches the scaling solution with 1/Z c (p) ∝ (p 2 ) −κ . The Euclidean dressing functions corresponding to the computed spectral functions for the different gluon propagator inputs are shown in Figure 5b. We show both the dressing functions obtained from the spectral Euclidean and real-time DSE and find that the spectral representation for the ghost propagator (and dressing function) holds. We also compare to the Yang-Mills results from [33], which we also used in the renormalisation condition as described above, and see that we reach very good qualitative agreement. In particular, our most scaling-like solution shows the typical scaling behaviour down to about 30 MeV. The deviations from [33] in the UV originate in the classical approximation for the full ghost-gluon vertex used here.

A. Comparison with previous works
In this section we compare our results on the ghost spectral function with that in the literature, for results with different approaches see [5][6][7]70]. The spectral function also allows us to map out the ghost propagator in the complex momentum plain, which is discussed in Section IV C including a comparison with respective results in the literature from a DSE analysis, see [10].
In [6], the ghost spectral function has been reconstructed from lattice QCD data. The results are in good agreement with our direct computation: both show a massless particle pole and a negative scattering tail. The reconstruction in [6] lacks reliability for spectral values smaller than roughly 100 MeV, as the smallest Euclidean data point used for the reconstruction is at about p = 150 MeV. In this regime, the present results from a direct spectral computation can be used as an input for future reconstructions by restricting the respective infrared completion. The same qualitative features are also found in the ghost spectral function obtained via a massive propagator expansion in [7], i.e. a massless particle pole as well as negative spectral tail.
In [5], Pade-type reconstructions of the ghost dressing spectral function from DSE and lattice data in Yang-Mills theory has been performed. These results are in contradistinction to the present result and [6,7], as the scattering tails in [5] show significant negative contributions. This corresponds to positive contributions in the propagator spectral function due to a relative sign in the definition. In addition, the UV tail of the reconstruction of DSE data is not shown in [5], and the spectral function appears to approach a constant value. A UVpositive (negative) as well as a non-vanishing tail in the propagator (dressing) spectral function violates the analytically given asymptotic fixed by Equation (7), for a detailed derivation see Appendix A.
The study of the analytic structure of the ghost propagator put forward in [70] also suggests the existence of a massless pole as well as a branch cut along the real frequency axis. As already mentioned above, Yang-Mills propagators in the whole complex momentum plane have been investigated with DSEs in [10]. The findings show good qualitative agreement with the propagators obtained from the ghost spectral function computed in the present work, but do not support a KL spectral representation of the ghost. This is discussed further in Section IV C.

B. Spectral fits
The results for the ghost spectral function with the UV asymptotics ρ (UV) (λ) in (31), the IR asymptotics ρ 0 allow for a simple fit in terms of the both asymptotics and Breit-Wigner functions for intermediate spectral values. The split into these three regimes allows for a simple parametrisation ρ (fit) c of the ghost spectral function, In (32a) we use the sigmoid function for projecting on the three regimes, 9. One of the two BSEs comprising the system that controls the scalar glueball formation. The blue (red) ellipses denote the glueball-gluon (ghost) BS amplitudes, and k± = k ± P/2, q± = q ± P/2, where q denotes the loop momentum.
where κ only carries the appropriate dimension. The intermediate regime is expanded in Breit-Wigner kernels, For our best fit, we use N = 3. The respective fit parameters are listed in Appendix C, Table I, and the fits are depicted together with the spectral functions in Figure 11. The accuracy of the fits is best evaluated within a comparison between the ghost dressing functions 1/Z c (p) obtained from the computed spectral functions and their fits. This comparison if depicted in Figure 7 for all three different input gluon propagators.

C. Results in the complex plane
We close this section with a short discussion of the potential application of the present results within bound state and resonance computations in QCD. To begin with, the behaviour of QCD correlation functions for complex-valued momenta is instrumental for the reliable computation of bound-state properties within the frameworks of the Bethe-Salpeter equations (BSEs). In this quest, the gluon and ghost propagators are of paramount importance, as may be exemplified by considering the BSEs that control the formation of glueballs in a pure Yang-Mills theory [71][72][73][74][75][76] (for lattice studies, see [77] and references therein). In fact, the present results are specifically useful for the scalar glueball: in contradistinction to its pseudo-scalar counterpart, it involves both the gluon and ghost propagators, as shown in Figure 9.
As is well-known, the need to extend the aforementioned propagators to the complex plane stems from the fact that the momentum P of the bound-state in question must satisfy P 2 = −M 2 , where M is the corresponding mass. This condition is typically implemented by introducing the rest-frame parametrization P = M (0, 0, 0, i) (see, e.g., [78]). Invariably, this complexifies the arguments of G A (q ± ) and G c (q ± ) in the BSE of Figure 9, These considerations motivate the computation of the dressing function 1/Z c (p) in the entire complex plane. To that end, we employ the KL representation of (A1), utilising the ρ c (λ) found above, and setting ip = Re ω + i Im ω. The results of this computation are shown in Figure 8.
We now compare our results for the ghost propagator in the complex plane with the spectral DSE with those from [10]. There, ghost and gluon propagators in the complex plane have been computed with complex DSEs. The gluon propagators in [10] exhibit complex conjugate singularities, and their nature and position varies greatly under small changes in the ansatz for the vertex. We emphasise, that these singularities simply indicate the limited radius of convergence of the method both for the gluon and for the ghost, for a detailed discussion see [10]. For large (angular) distances to the Euclidean axis analyticity is lost, and the method used does not produce reliable results. If reconstructed with the Schlessinger point method, the singularities observed in [10] take the form of complex conjugate poles. This has also been seen in [5], where similar reconstruction methods have been used. For further studies of the complex structure of QCD-like theories in the presence of complex conjugate poles see also the recent work [8,45].
Despite the lack of reliability for sufficiently large Minkowski frequencies, we have compared ghost dressing from [10] with the present result in this region. The imaginary part of the ghost dressing function computed there is strictly positive for time-like momenta, in qualitative agreement with our result, in particular in view of the different approximations. We have also confirmed the absence of a spectral representation of the ghost by computing the spectral function from the Minkowski dressing of [10] via Equation (7). Then, the Euclidean dressing is computed via the KL representation of Equation (6) and compared to the direct calculation. This comparison showed a significant violation of the spectral representation especially for larger Euclidean frequencies. This is to be expected, since by nature of the kernel of the spectral representation Equation (6), large Euclidean frequencies are sensitive to large spectral values, i.e. large Minkowski frequencies. These lie beyond the radius of convergence of the method used in [10], as discussed.
In summary, this analysis strongly suggests that complex conjugate poles as well as other non-analyticities in the gluon propagator beyond the real frequency axis invalidate the KL representation for both the gluon and the ghost. A more detailed discussion is deferred to future work [79]. In particular, this casts serious doubts on mixed reconstructions with a KL representation for the ghost and cc poles for the gluons. In turn, the gluon spectral function in [4] was reconstructed with the assumption of a KL representation, as outlined in Section III C. As shown in the present work, this also leads to a KL representation of the ghost. Whether this property holds true in a self-consistent solution of the coupled system, remains to be seen and is deferred to future work.

V. CONCLUSION
We have solved the Dyson-Schwinger equation for the ghost propagator in the complex plane on the basis of a given input gluon spectral functions, spanning the whole family of decoupling solutions, including the scaling limit. Our spectral DSE approach is based on the spectral DSE put forward in [1], and uses the spectral renormalisation devised there. The procedure allows for analytic solution of the momentum loop integrals by utilising the KL representation and dimensional regularisation. This facilitates the access to the full complex momentum plane, constituting the central aspect of our scheme. The present truncation uses classical vertices in the ghost gap equation, but we emphasise that the spectral DSE approach also allows for non-trivial vertex approximations, see [1].
The input data for the gluon spectral function is constructed via a decoupling-type modification of the scaling spectral function from [4]. The latter spectral function has been obtained via a reconstruction of the scaling solution fRG data of [33].
The spectral function for the ghost shows a massless pole as well as a continuous scattering tail. The classical massless mode dominates up to momenta close to the position of the maximum of the input gluon propagator. For larger momenta, the perturbative logarithmic behaviour starts to dominate, ultimately causing the dressing function to vanish in the ultraviolet. The present results and the current real-time approach with a real-time renormalisation scheme opens the door to a systematic spectral access to dynamical, time-like properties of QCD. We hope to report on a self-consistent spectral investigation of the full Yang-Mills system soon, both within the standard DSE approach and within the pinch technique. The respective results are pivotal for following studies of the resonance properties and the dynamics of QCD within the present approach. In a general manner, given a KL representation, the normalisation relation for the corresponding spectral function can be inferred from the perturbative behaviour of the propagator. Multiplying (6) by p 2 , one has In the UV, the behaviour of the dressing function Z(p) can be inferred from perturbation theory, lim p→∞ 1/Z(p) = Z −1 ∞ . For large p 2 , we can also expand the integrand, yielding defining We want to show lim p→∞ ∆(p) = 0 in order to obtain a normalisation condition for ρ via (A2) using the known perturbative asymptotics of the corresponding dressing function. In doing that, we first note that for the spectral integral in the left term of the lower line in (A2) to converge, the spectral function must obey lim ω→∞ ρ(ω) ω 2 log ω 2 → 0 .
If this requirement does not hold, ρ cannot be normalised in the above form. Based on the assumption of the existence of above representation (A1), ρ can be taken to be integrable on [0, ∞). We choose a scale Λ such that for frequencies λ > Λ, ρ is given solely by the leading UV behaviour of its corresponding propagator via (7), see also [4]. Denoting the known UV asymptotics as ρ UV , we then distinguish where ρ UV now obeys (A4). Note that by the nature of the spectral function being a tempered distribution, it can have distributional contributions such as (higher order) poles. These are allowed in our consideration as long as integrability is not violated. The parametrisation (A5) is chosen such that these contributions are contained in ρ Λ . We now split the spectral integration interval of (A3) along the split of the spectral function and conclude for finite Λ that such that for large momenta, the contribution (A6) to the spectral representation of the dressing function vanishes, Hence, in the limit of large p we are only left with spectral integral over ρ UV contributing to ∆(p) in (A3). Taking into account the known asymptotics of ρ UV from (A4) however, we find that where the lower line can already be anticipated to vanish for arbitrary constants C. However, this can also be shown rigorously by noting that upon substitution, the last line of (A8) can be reexpressed as the exponential integral function E 1 , The contribution from the lower integral boundary is finite and thus vanishes in (A8). For the upper limit we utilise the asymptotic expansion of the exponential integral and plug this back into (A8), yielding, while dropping the constant prefactor, In conclusion, recalling (A3), we arrive at which, with (A2), eventually yields the desired normalisation for the spectral function, We thus see that, in the fairly general case where ρ can be normalised via the integral in (A12), the normalisation is given by the value of the dressing function at infinity.

Appendix B: Loop momentum integration of the ghost-gluon diagram
In this appendix we detail the computation of the ghost self energy diagram Σc c . Starting at (18), we define Σc c (p) = g 2 δ ab C A λ1,λ2 with the now dimensionally regularised momentum integral The measure is now q = d d q/(2π) d .

Momentum integration
Next, we employ partial fraction decomposition and introduce Feynman parameters, i.e. utilise Upon a shift in the integration variable q → q − xp and after some manipulation, we arrive at We will not make all intermediate results explicit, such as giving the full expressions for A i and B i , which are functions of external momentum p, the spectral parameter λ 1 as well as the Feynman parameter x. Ultimately, the complete final result will be stated explicitly. The momentum integrals are now readily solved via the standard integration formulation, with m a non-negative and n a positive integer.

Feynman parameter integration
Reordering the expression in powers of the Feynman parameter x and taking the limit d → 4 − 2ε, we arrive at with γ E the Euler-Mascheroni constant. The coefficients α i and β i do not depend on x, and will be given down below. We can solve the Feynman parameter integrals analytically and simplify the first sum to obtain the final result, The coefficients α i , β i are defined as follows: and β 0 = 0 , The functions f i and g i carry the branch cuts ultimately giving rise to the spectral function and are defined by integrals over the Feynman parameter x via where we defined with ζ = λ 4 2 + λ 2 1 + p 2 2 + 2λ 2 2 (p 2 − λ 2 1 ), and 2 p 4 + 6λ 4 2 p 2 − 6p 6 log(λ 2 2 ) + 6 log(−λ 2 2 )(λ 2 2 + p 2 ) 3 − 6(λ 2 2 + p 2 ) 3 log(−λ 2 2 − p 2 ) + 13p 6 ,

Real frequencies
For a real-time expression of the ghost gluon loop, we need (B9) at real frequencies ω, i.e. I(ω, λ 1 , λ 2 ) := I(−i(ω + i0 + )). From the definitions of the functions α i (B10), β i (B11), g i (C1) and f i (B13), their respective real-time expressions are trivially obtained by the substitution p ↔ iω. The calculations were performed in Wolfram Mathematica 12.1 with the convention Im log x = π for x < 0 for the logarithmic branch cut.

Spectral integration and convergence
The spectral integrals of the form where I ren is the renormalised spectral integrand (comp. (20) or (21)), are evaluated numerically on a logarithmic momentum grid of about 200 grid points with boundary (p min , p max ) = (10 −4 , 10 2 ), identically for the Euclidean and Minkowskian axis. We use a global adaptive integration strategy with default multidimensional symmetric cubature integration rule. After spectral integration, the diagram is interpolated with splines in the Euclidean and Hermite polynomials in the Minkowski domain, both of order 3. The spectral function is then computed from the interpolants. Note that, a priori, due to (7), the domain of the ghost spectral function is given by the momentum grid. The integration domain of the spectral integral of the ghost spectral parameter has to be bounded by (p min , p max ), in order to not rely on the extrapolation of the spectral function beyond the grid points. Due to numerical oscillations at the very low end of the grid, we choose (λ min c , λ max c ) = (10 −3.5 , 10 2 ). Convergence of the integration result with respect to increase of the integration domain has been explicitly checked.
For the gluon spectral integral, the situation is different, as the spectral function is given in an algebraic form from [4]. We use the integration boundary (λ min A , λ max A ) = (10 −4 , 10 2 ).

Spectral integrands
The numerical performance of the spectral integrations presented in Appendix B is sped-up by up to two orders of magnitude by using interpolating functions of the numerical data. The interpolants are constructed by first  Table I. The change of sign around 1.2 GeV is an imprint of the oscillations of the input reconstructed gluon spectral function from [4] and is discussed in Appendix C 4.
discretising the integrand inside the three-dimensional (p, λ c , λ A ) cuboid defined by p ∈ 10 {−4,2} , λ c/A ∈ 10 {−4,4} . As for the momentum grid for the spectral integration, we use the same cuboid for the real-and imaginary-time domain. We use 60 grid points in the momentum and 160 grid points in the spectral parameter integration, both with logarithmic grid spacing. For the real-time expressions, we divide into real and imaginary part of the integrands. Both real and imaginary parts of the discretised Minkowski as well as the Euclidean expressions are then interpolated by three dimensional splines inside the cuboid. The resulting interpolating functions are then used in the spectral integration.

Spectral fits
As discussed in Section IV, we provide a ready-to-use analytic fit formula for the ghost spectral function, see (32). For our best fit we use N = 3, the fit parameters for the ghost spectral functions for all input gluon propagators are listed in Table I. We show the spectral functions and their respective fits on a log-log scale in Figure 11. For G A (0) [GeV −2 ] = 4.4 and 1.9, the spectral functions feature a change of sign between 1.2 and 1.3 GeV. These wiggles are imprints of the oscillations in the input reconstructed gluon spectral function of [4], and can be understood as numerical artefacts from the reconstruction process. However, in order to match the original Euclidean dressing function 1/Z c (p) in the UV (comp. Figure 7), it is necessary to keep the respective oscillatory behaviour in the fit.