Reconstructing QCD Spectral Functions with Gaussian Processes

We reconstruct ghost and gluon spectral functions in 2+1 flavor QCD with Gaussian process regression. This framework allows us to largely suppress spurious oscillations and other common reconstruction artifacts by specifying generic magnitude and length scale parameters in the kernel function. The Euclidean propagator data are taken from lattice simulations with domain wall fermions at the physical point. For the infrared and ultraviolet extensions of the lattice propagators as well as the low-frequency asymptotics of the ghost spectral function, we utilize results from functional computations in Yang-Mills theory and QCD. This further reduces the systematic error significantly. Our numerical results are compared against a direct real-time functional computation of the ghost and an earlier reconstruction of the gluon in Yang-Mills theory. The systematic approach presented in this work offers a promising route towards unveiling real-time properties of QCD.

Introduction. The resolution of many open questions in quantum chromodynamics (QCD) requires the knowledge of time-like observables and hence the computation of real-time correlation functions. Applications range from the hadronic resonance spectrum over scattering processes to transport and non-equilibrium phenomena in heavy-ion collisions. For example, the computation of the glueball spectrum via Bethe-Salpeter equations relies on the time-like propagators for gluon and ghost, both of which are reconstructed in the present work. Likewise, QCD transport coefficients used in hydrodynamic simulations can be computed diagrammatically from the real-time gluon propagator. Similarly, phenomenological QCD transport models with their underlying assumption of a quasi-particle nature of the gluon can hugely benefit in multiple ways from the present results. First of all, a reliable computation of the gluon spectral function may offer much-needed support for the quasi-particle assumption of these models, as well as give access to its limitations. Secondly, the QCD gluon spectral function itself can feature as a direct input and pivotal building block in these models. Together with further time-like correlation functions, this offers a path for a systematic quantitative improvement of phenomenological transport approaches towards first-principle transport in QCD.
By now, Euclidean correlation functions in QCD are accessible within first-principle approaches such as lattice simulations or functional equations. In contradistinction, accessing real-time properties remains a notoriously hard task. Minkowski correlation functions may be obtained from Euclidean data via spectral reconstruction, exploiting the Källén-Lehmann (KL) representation [1,2]. This requires computing the spectral function via an inverse integral transform. In the present work, we approach the * urban@thphys.uni-heidelberg.de; corresponding author. problem with Gaussian process regression (GPR). The applicability of GPR to inverse problems of this type has been discussed in [3]. Specifically, it was shown how GPs can be used to obtain probabilistic models of functions for which only weighted averages are available.
We apply GPR to the reconstruction of ghost and gluon spectral functions based on recent results from 2+1 flavor lattice QCD with domain wall fermions at a pion mass of 139 MeV [4,5]. Furthermore, we improve the systematic error control by incorporating additional data in the infrared (IR) and ultraviolet (UV) regimes from functional renormalization group (fRG) and Dyson-Schwinger (DSE) computations in Yang-Mills theory and QCD [6][7][8][9][10][11][12], mostly obtained within the fQCD collaboration [13].
Spectral representation. The KL spectral representation of the two-point correlation function in momentum space reads G(p 0 ) = ∞ 0 dω π ω ρ(ω) ω 2 + p 2 0 = ∞ 0 dω K(p 0 , ω) ρ(ω) , (1) with the KL kernel K(p 0 , ω) and ρ(−ω) = −ρ(ω). In the vacuum, the spatial momentum dependence of the propagator can be obtained via a Lorentz boost, simply by p 2 0 → p 2 with p 2 = p 2 0 + p 2 . With Equation (1), the spectral function is obtained from the retarded propagator via ρ(ω) = 2 Im G(−i(ω + i0 + )) . ( For asymptotic states, the spectral function is the probability density for (multi-)particle excitations created from the vacuum in the presence of the corresponding quantum field. Consequently, in this case the spectral function is positive semi-definite. For propagators of 'unphysical' fields, such as gauge fields, the spectral representation may still hold. However, the spectral function can then also have negative parts, and the existence of a spectral representation simply constrains the allowed complex structure of correlation functions; see e.g. [8,12,[14][15][16].
In this work, we reconstruct ghost and gluon spectral functions of 2+1 flavor QCD under the assumption that both admit a KL representation. It can be shown that the total spectral weight vanishes, respectively for both the ghost and gluon spectral functions, ρ c and ρ A . For the gluon, this is the wellknown Oehme-Zimmermann superconvergence (OZS) condition [17,18]; for recent discussions with general fields, see [8,12,16]. These works also include a treatment of the analytic low-frequency behavior of continuous parts of the spectral functions, initiated in [8].
A general spectral function ρ consists of a continuous partρ and a sum of particle and resonance peaks (proportional to the δ-function and its derivatives). In this work, we assume that the gluon spectral function only consists of a continuous part ρ A =ρ A satisfying Equation (3). This is the generic structure suggested by all functional equations describing the gluon propagator due to the ghost being massless. While derivatives of δ-functions are formally also allowed, we exclude these structures from our ansatz due to the absence of a generic mechanism generating the required roots of the inverse gluon propagator on the real momentum axis. In turn, due to the 1/p 2 behavior of the Euclidean lattice ghost propagator in the IR, the associated spectral function exhibits a particle peak at vanishing frequency in addition to its continuous part, i.e.
Reconstruction with GPR. Starting from early developments in the context of geostatistics in the 1950s [47], today GPR is widely employed in a variety of settings for the probabilistic modeling of functions from a finite number of observations; see [48,49] for reviews and [50] for a modern textbook account. Recently, the method has been applied to the reconstruction of parton distribution functions from lattice QCD [51]. In this section, we summarize the main ingredients for spectral reconstruction with GPR based on the developments reported in [3]. A short introduction to GPR for function prediction as well as further details and references are provided in Appendix A.
We assume our knowledge of the spectral function ρ(ω) to be described by a GP, written as where µ(ω), C(ω, ω ) denote the mean and covariance functions. Importantly, in this approach we do not restrict the space of possible solutions by choosing a specific functional basis, which often leads to spurious artifacts in the reconstruction in order to compensate for unrepresentable features. Instead, the GP defines a distribution over families of functions with rather generic properties, specified via the kernel parametrization described below. The KL integral in Equation (1) is a linear transformation that preserves Gaussian statistics. Hence, given Equation (5) one may obtain statistical predictions G i at N G specified momenta p i as Here, N denotes a multivariate normal distribution, to be distinguished from distributions over function space denoted by GP. Statistical uncertainties associated with individual prediction pointsμ i may be computed from the diagonal of the covariance matrix asσ i = C ii . Conversely, the framework also enables inference in the opposite direction. The inherent analytic tractability associated with Gaussian statistics allows formulating the conditional distribution for ρ(ω) given observations G i in closed form. The full expression may then be derived as dηdη K(p i , η)C(η, ω)  Figure 2. The results agree within the given statistical uncertainties as shown in the bottom panels, where the posterior GPs for the correlators are evaluated at the fixed momenta provided by the lattice data, which is then subtracted leaving the error bars intact. The total mean squared errors amount to ∼5e-6 for the ghost and ∼4e-5 for the gluon.
The GP in Equation (7) encodes our knowledge of the spectral function after making observations of the propagator and accounting for observational noise with variance σ 2 n . The corresponding expressions for the dressing function instead of the propagator can be immediately obtained by inserting an additional factor of p 2 i at every occurrence of the KL kernel K(p i , ω) in Equations (6) and (7).
The flexibility of the approach makes it possible to also incorporate further available prior information in various forms into the predictive distribution in the same manner, yielding similar though somewhat more complicated expressions. This may include e.g. direct observations of ρ and its derivatives, assumptions about the asymptotic behavior, or global normalization constraints.
In order for GPs to be useful for modeling, the covariance C(ω, ω ) may be defined via a so-called kernel function. It is commonly parametrized using a small number of hyperparameters, which may be subjected to optimization based on the associated likelihood. The mean function µ(ω) is often set to zero, since its contribution can be fully absorbed by the kernel. Typically, the latter is the sole focus of the optimization procedure. However, a custom mean function may still be useful in certain situations in order to incorporate prior beliefs about the functional form of the expected solution, which can improve the calculation by providing a better starting point for the optimization routine.
A frequently used kernel parametrization is the radial basis function (RBF) kernel, also called squared exponential. It is defined as where the parameter σ C controls the overall magnitude and l is a generic length scale. The RBF kernel has been established as the standard choice for many applications due to a number of attractive features, such as universality [52] and every function in its prior being infinitely differentiable. It is also used for our first results on spectral reconstruction with GPR presented in this work. Nevertheless, designing custom kernels for specific problems has been shown to greatly increase the usefulness of the approach in various settings and is also promising here. In particular, it may be interesting to construct kernel functions that can be integrated analytically against the KL kernel, such that the frequency integrals in Equations (6) and (7) may be carried out analytically instead of numerically. To this end, one could potentially employ functions of Breit-Wigner type as done for the spectral function itself in [8]. In contradistinction, we may use them to instead define a suitable GP kernel, thereby still avoiding the restriction to a specific functional basis as previously mentioned. We comment on this and other possible improvements to our reconstruction approach in the conclusion. Furthermore, we emphasize that the present approach in principle does not require us to choose a specific set of nodes ω i . In fact, instead of computing a discrete set of point predictions or coefficients of a predefined functional basis, the prediction for ρ is obtained as a function of ω, albeit only implicitly via the kernel formulation. In particular, the GP also allows computing all of the derivatives of the prediction analytically at any pointincluding the associated statistical uncertainties-by differentiating the expressions in Equation (7) with respect to ω (as well as ω for the covariance). A finite set of nodes ω i is chosen only at inference time in order to evaluate the GP, however, the choice is completely arbitrary within the given domain. This property is one of the most attractive features of GPR for spectral reconstruction and probabilistic function prediction in general.
Input Data. In the past two decades, increasing interest in the momentum behavior of the fundamental two-point Green's functions in QCD as well as further correlation functions of higher order has triggered respective lattice calculations in particular of Yang-Mills and QCD propagators; see e.g. [53][54][55][56][57][58][59][60][61][62][63][64][65][66][67]. The lattice data for the ghost dressing function and gluon propagator employed in this work are shown in Figure 1. They are obtained from recent simulations of 2+1 flavor QCD at the physical point [4,5]; see Appendix B 1 for further details and references. Additional input data and benchmarks are provided by one-parameter families of solutions from functional computations in Yang-Mills theory and QCD [6,8,11,12], which are matched to the continuum-extrapolated lattice data as shown in Reconstruction Results. The GPR for the reconstruction of the ghost spectral function is performed using the aforementioned standard RBF kernel. We extend the lattice input data for the dressing function into the deep IR and simultaneously fix the low-frequency asymptotics of the spectral function using a direct real-time result in Yang-Mills theory obtained via the spectral ghost DSE [12] (see also Appendix B 2). This is achieved by treating the spectral DSE result as an additional observation. Our procedure uniquely determines the non-zero value of ρ c for ω → 0 + , but also increases the reliability of the solution in the most interesting central region with respect to the kernel hyperparameters. Using just the lattice data without the extension by the spectral DSE result leads to a much higher variance in the solution space, with widely different asymptotic behaviors of solution candidates in the IR. The kernel hyperparameters are chosen by optimizing the associated likelihood of observations with an additional Gaussian hyperprior, which we achieve through a fine-grained grid scan; see Appendix C for details. The reconstructed spectral function in Figure 2a accurately reproduces the dressing function data within the uncertainties displayed in Figure 1a, with a total mean squared error of ∼5e-6.
The features of our prediction are strikingly similar to the aforementioned Yang-Mills result shown in Figure 4a in Appendix B, even though only the IR limit is incorporated into the reconstruction. This is expected heuristically, since the ghost only interacts with the quarks indirectly via the gluon vertices, and the effects of introducing dynamical quarks must hence be of higher order. The similarity is particularly notable considering that the methods are conceptually very different.
For the reconstruction of the gluon spectral function, the lattice input data are extended into the UV using an earlier fRG computation [8], which is quantitatively reliable in this regime. We discuss this in more detail in the next paragraph and in Appendix B 2. As for the ghost, this extension leads to greatly enhanced stability of the reconstruction with respect to the kernel hyperparameters. In particular, it ensures convergence to zero for ω → ∞, whereas with just the lattice data we often observe convergence to a non-zero constant and in some cases even pathological divergences. A modified frequency scale is used in the RBF kernel in order to suppress spurious oscillations in the IR and UV tails. The hyperparameters are again obtained via optimization of the likelihood with Gaussian hyperpriors while approximately enforcing the OZS condition; see Appendix C for details. The reconstruction shown in Figure 2b accurately reproduces the lattice data within the given uncertainties, as shown in Figure 1b, with a total mean squared error of ∼4e-5. While also being fully consistent, deviations from the lattice propagator are somewhat stronger than for the ghost dressing function and seem to become more pronounced in the IR. This is likely caused by the comparably large uncertainties of the lattice data at small momenta.
The peak structure of the spectral function appears similar to an earlier reconstruction of the Yang-Mills propagator in the fRG framework [8], shown in Figure 4b in Appendix B. We emphasize that the UV extension is done with the Yang-Mills data of [6] instead of the full 2+1 flavor results from [11]. This is detailed in Appendix B 2 and facilitates the comparison with the Yang-Mills reconstruction [8]. In particular, the positions of the leading positive peaks approximately coincide, with ω ≈ 0.818 for the present result and ω ≈ 0.835 for the fRG reconstruction. This reflects the approximate coincidence of the peaks of the Euclidean gluon dressing functions shown in Figure 3a in Appendix B. We also note that a small peak to the right of the second local minimum is present in both reconstructions. This feature may be a generic reconstruction artifact since it is not necessitated by theoretical considerations, but is observed in both results from conceptually very distinct methods. However, the comparably large uncertainties in this region also include plausible solutions without additional zero-crossings.
Significant differences between the two reconstructions are observed mainly in the overall peak height and width. Generally, the QCD result for the gluon is expected to differ more strongly from the pure gauge theory than the ghost due to the direct coupling to quarks. However, differences may also be attributed in part to the limited availability and precision of data and the resulting difficulty in resolving highly peaked structures. We find that generating narrower peaks with greater amplitudes by allowing the kernel's magnitude parameter σ C to increase and the length scale l to decrease leads to stronger oscillations in the solution. This is a common feature of conceptually similar reconstruction approaches, such as linear regression with a Tikhonov regularizer (also called ridge regression), which has been applied e.g. in [31]. Introducing such a regularization scheme, which is equivalent to assuming a Gaussian prior, leads to a favoring of solutions that are closer to zero. This additional bias can introduce the unwanted oscillations. Within the GPR approach, the kernel hyperparameters provide more detailed control over the regularization and can be tuned to deliberately suppress such unphysical features. However, this may result in reconstructions that are naturally flatter, which must be taken into account when interpreting and utilizing the result. This demonstrates one of the key advantages of GPR, namely the possibility to dynamically adjust the resolution depending on the available amount and quality of the input data, while still matching the observations as accurately as possible.
Although the obtained spectral functions reproduce the lattice data to high accuracy, the asymptotic behaviors of the mean predictions in the deep IR and UV differ from the analytic results derived in [8]. In particular, different scaling exponents are observed and the gluon spectral function shows the opposite sign in the UV. Nevertheless, the analytically expected behavior is still plausibly contained within the computed errors, which are comparably large in these regimes. This indicates that not enough prior information is available to the GP from just the data in order to accurately resolve the tails of the spectral functions, which may come as no surprise. While this issue does not affect the reconstruction in the region of interest, it may be problematic for precision computations that use these results as inputs. In order to directly enforce the correct asymptotics, potential approaches are the incorporation of the analytically known behaviors into the prior means of the GPs or finding more suitable choices for the kernel functions. Furthermore, exploiting the available analytic results to provide additional prior information about the derivative structure may be particularly helpful in stabilizing the tail behavior. To achieve this, one may again write down the joint distribution of the predicted spectral function at any frequency and its associated derivatives to arbitrary order in closed form and derive the conditional posterior distribution similar to Equation (7).
Conclusion. In this work, we apply Gaussian process regression to the reconstruction of ghost and gluon spectral functions in 2+1 flavor QCD at the physical point. These spectral functions are the pivotal building blocks of diagrammatic representations for bound state equations such as Bethe-Salpeter and Faddeev equations, see e.g. [68][69][70], as well as transport coefficients, see e.g. [23,71].
Importantly, the gluon spectral function has a pronounced quasi-particle peak, the position of which is related to the mass gap in QCD. This extends previous vacuum and finite-temperature results in Yang-Mills theory [8,23] to physical QCD. Our findings provide non-trivial QCD support to the phenomenological use of quasi-particle gluon spectral functions for transport computations; see [72] for a recent review. Moreover, the present results can be directly employed as first-principle QCD inputs in order to systematically improve the respective phenomenological approaches towards a firstprinciple treatment of QCD transport processes.
These promising phenomenological applications of the present results also highlight the necessity of further improving the reconstruction approach itself, for which a number of potential directions can be envisaged. This includes the aforementioned possibility of designing custom kernels for the problem at hand, potentially with analytic integrability against the KL kernel. Constructing suitable, expressive kernels may also be automated and improved through the use of hyperkernels [73] or techniques such as deep kernel learning [74]. To account for some variability in the kernel hyperparameters, one may replace the maximum likelihood approach by an integral over parameter space using a suitable hyperprior which encodes any prior assumptions. Alternatively, optimal hyperparameters may also be selected based on a data-driven machine learning approach, using datasets consisting of pairs of correlators and associated spectral functions.
Furthermore, the flexibility of the GPR framework allows the incorporation of various supplementary constraints derived from theoretical arguments, such as information about derivatives, known asymptotic behaviors, or normalization conditions. This is expected to further improve the accuracy and reliability of the reconstruction, in particular for the IR and UV tails of the spectral functions that are otherwise difficult to resolve. This will be the subject of future work, accompanied by direct functional computations of further spectral properties along the lines of [12,75].
The immediate next steps in our endeavor towards unveiling real-time properties of QCD are the application and extension of the present numerical method to quark propagators as well as correlation functions computed at finite temperature. This will enable quantitative studies of hitherto theoretically inaccessible non-equilibrium dynamics of QCD in the transport phase of heavy-ion collisions within a first-principle approach.
Acknowledgements. We thank A. Cyrol, F. Gao, L. Kades, J.Y. Lin, J. Papavassiliou and A. Rothkopf for discussions. We thank P. Boucaud and F. De Soto for an earlier collaboration and for their help in the preparation of the lattice data. We are indebted to the RBC/UKQCD collaboration, especially to P. Boyle, N. Christ, Z. Dong, N. Garron, C. Jung, B. Mawhinney and O. Witzel, for access to the lattices used in this work. We also thank the members of the fQCD collaboration [13] for discussions and collaborations on related subjects. JRQ acknowledges the support of MICINN PID2019-107844GB-C22 grant. Numerical computations have used resources of CINES, GENCI IDRIS (project id 52271) and of the This appendix serves as a brief introduction to GPR for function prediction using a finite number of direct or indirect observations, based primarily on [3]. We adopt the notation used in the main text for consistency, however, the general formalism presented here is also applicable outside of the specific context of spectral reconstruction for quantum field theory. For a modern, comprehensive textbook treatment of the topic, we refer the interested reader to [50]. For a brief, pedagogical introduction to GPR with simple code examples, we recommend [76]. In the context of inverse theory, [77] provides a recent review.
We first discuss GPR for the case where direct observations are available for the function to be modeled. We assume our knowledge of the function ρ(ω) to be encoded in a GP with mean and covariance functions µ(ω), C(ω, ω ), denoted by where the covariance is assumed to be symmetric, i.e. C(ω, ω ) = C(ω , ω). As per the definition of a GP, any finite set of function evaluations at N sample points ω i follows a multivariate normal distribution, Similarly, we can write down the joint distribution of a set of observationsρ i at pointsω i and the value of the function at an arbitrary point ω as (A3) where boldface type denotes vector and matrix quantities. Here, we have definedμ ≡ µ(ω i ),Ĉ i (ω) ≡ C(ω i , ω), andĈ ij ≡ C(ω i ,ω j ). σ 2 n defines the point-wise variance of additional measurement noise which may be present in the observationsρ. Due to the analytic tractability of multivariate Gaussians, the conditional distribution of function values ρ(ω) given observationsρ may then be derived as The covariance is parametrized by a suitable kernel function, whereby one may encode any prior beliefs about the types of solutions one expects by choosing an appropriate form for the problem at hand. For an introduction to constructing GP kernels of various types as well as strategies to apply and combine them, we recommend the kernel cookbook [78]. A kernel's hyperparameters, denoted here byα, may be subjected to optimization by maximizing the associated likelihood, where we have writtenĈα to emphasize the dependence on the hyperparameters. Instead of directly maximizing p(ρ|α) as a function ofα, one conventionally minimizes the negative log likelihood (NLL), (A6) Since simply finding and employing the maximum likelihood configuration of hyperparameters may ignore relevant additional structures in the distribution, one can also integrate outα using suitable hyperpriors to account for some variability.
Based on the formulation of GPR for direct observationsρ at pointsω, one can derive the expressions for inference from indirect observationsĜ at pointsp as discussed in the main text by applying the forward process of the associated linear inverse problem, in our case the KL integral defined in Equation (1). This involves all terms related to the observations that depend on the discrete set of pointsω, which are promoted back to the continuous domain and subsequently integrated out to yield the nodesp instead.

Appendix B: Input Data
Combining the data from lattice simulations and functional computations as described in the main text requires matching the scales through renormalization. In this work, we always rescale the functional methods results to match the lattice data in the appropriate regime.

Lattice Simulations
The lattice data employed in this work were obtained from configurations generated by the RBC/UKQCD collaboration-first introduced in [79-83]-with 2+1 dynamical quark flavors using the Iwasaki [84] and domain wall fermion [85,86] actions, respectively for the gauge and quark sectors, at the physical point (a pion mass amounting to 139 MeV) by the particular implementation of the Möbius kernel [87]. These developments were then exploited in [4,5] in order to calculate the gluon and ghost propagators as well as the strong coupling in a particular scheme [88][89][90], and an effective charge stemming from it [91]. A description of this calculation is given, for instance, in [62].
In computing propagators that properly feature the physical running with momenta, data should be thoroughly cured from lattice regularization artifacts. In particular, as explained in [4], our results are obtained after a careful scrutiny of discretization artifacts, thereby accounting for the continuum-limit extrapolation, following [92]. As a noteworthy remark, a recent work [67] has revealed the key role played by the procedure of [92] for an adequate removal of discretization artifacts in achieving a consistent description of Yang-Mills two-and threepoint correlators, involving both lattice and DSE results.
The resulting ghost dressing function and gluon propagator data are displayed in Figures 1a and 1b, respectively. They are compared against their counterparts obtained from evaluating Equation (1) for the reconstructed spectral functions shown in Figure 2, as well as the results from functional methods described in the following section. The dressing functions of all input datasets are compared in Figure 3 to further illustrate their similarities and differences.

Functional Methods
We briefly summarize results from functional computations in Yang-Mills theory and QCD that are employed in this work to provide additional prior information for the reconstruction. For reviews on the application of functional methods in this context, see e.g. [93][94][95][96].
We use the real-time Yang-Mills results from [12] to extend the lattice QCD data of the ghost dressing function into the deep IR, as shown in Figure 3a. The approach also provides direct access to the associated spectral function, which we employ to fix the low-frequency asymptotic behavior of the reconstruction. It is obtained via the spectral ghost DSE, building upon the technique of spectral renormalization [75]. Making use of Equation (1) for the ghost and gluon propagator, the momentum integrals appearing in the loop diagrams of the ghost propagator DSE can be solved analytically. This preserves the full analytic momentum dependence and allows evaluating the equation on the real momentum axis. The spectral function can then be directly extracted from the real- time propagator DSE via Equation (2); see Figure 4a for a comparison to the reconstruction result of the present work. As input gluon spectral function, the reconstruction result of [8] based on the scaling solution obtained via the fRG in [6] is used. Assuming a spectral representation for the gluon propagator, in both scaling and decoupling scenario the IR behavior of the gluon spectral function follows directly from the propagator [8]. This is utilized to modify the given scaling spectral function such that we obtain a decoupling-type gluon propagator matching the value of the given lattice propagator well within the given uncertainties. Due to its mild momentum dependence, the ghost-gluon vertex is assumed to be classical. The lattice QCD data for the gluon propagator are extended towards the UV using earlier results from functional computations in Yang-Mills theory [6]. Differences to the 2+1 flavor QCD result for the gluon propagator reported in [11], being based on [7], are comparably small in the relevant momentum range. A stronger deviation can be observed in the dressing functions, as shown in Figure 3b. Despite these differences, the reconstruction still produces remarkably reliable results, cf. Figure 1b. Nevertheless, we aim to replace the Yang-Mills UV extension by the 2+1 flavor QCD data from [11] in order to further optimize the accuracy of the result and mitigate any potential issues. For related results and further correlation functions see [9,10,97,98]. More specifically, the fRG results in [6] are derived within an advanced approximation where the momentum dependence of all vertices is approximated at the symmetric point, for respective DSE results see [99]. For our purposes, this data set provides the optimal trade-off for momentum range versus accuracy. Due to the high numerical precision, the results are particularly well-suited as an input for spectral reconstruction. The Yang-Mills data have already been employed for this purpose in [8] and we use this earlier reconstruction for comparison; see Figure 4b. In summary, the extension of the 2+1 flavor lattice data with the high precision Yang-Mills data up to momenta p 2 = 10 2 GeV 2 allows a more direct comparison (in terms of scales) with the Yang-Mills reconstruction in [8], while only modifying the large frequency tail of the gluon spectral function for frequencies ω 5 GeV, see Figure 4.

Appendix C: Implementation
In this section, we comment on certain points of the implementation in more detail. We first address numerical aspects of the optimization and a discussion of the required computational effort. Subsequently, we provide further information about data usage, kernel design choices and theoretical constraints for the particular reconstructions reported in this work.

Hyperparameter Optimization and Computational Cost
To find optimal values for the kernel's hyperparameters, we perform a fine-grained grid scan of the NLL with additional hyperpriors where necessary. Alternatively, the NLL may also be minimized with a gradientbased ansatz using a standard optimizer such as L-BFGS. However, mapping out the posterior distribution in more detail tends to be highly instructive for the problem at hand. It is also less prone to numerical problems such as unstable directions and violation of positive definiteness of the covariance, as these can be identified early on, and should hence be preferred when feasible. This is also where the bulk of the computational effort goes, as it involves calculating for each individual grid point the comparably expensive inverse and determinant of the covariance matrix, which naively scales like O(N 3 ). For very large datasets where their direct evaluation becomes infeasible, one may resort to cheaper linear solvers for the inverse and stochastic approximations of the determinant, but this is unlikely to become necessary in this particular context. Cost may also be mitigated by scanning the parameter space hierarchically, starting at low resolution and zooming into the interesting regions. The whole procedure is trivially parallelizable, as each grid point can be treated independently. At the scale of the present work, each instance was handled by a standard CPU node with low performance requirements. Some first tests were also conducted on a single machine, where mapping out the parameter space for each reconstruction with medium resolution took a few hours at most. In comparison to finding the optimal hyperparameters, the subsequent inference step is negligibly cheap. Of course, the total computational effort for the reconstruction is dwarfed by the requirements of the large-scale lattice simulations described in Appendix B 1, which are orders of magnitude more expensive.

a. Ghost
In the case of the ghost spectral function, we treat the low-frequency asymptotics extracted from the direct DSE computation in Yang-Mills theory as an additional observation for the GP. This is only possible for the ghost, as a similarly direct determination of the Yang-Mills gluon spectral function is currently not available. The procedure is implemented by including the value of ρ at ω = 0 in the construction of the joint distribution of observations and predictions. In particular, one needs to compute additional expressions for the covariances of the point ρ(0) and the correlator data. This requires some programming headache, but carries no further conceptual difficulty.
As stated in the main text, we use the standard RBF kernel and identify optimal hyperparameters via a highresolution grid scan. We note an unstable direction in the magnitude parameter σ C , which is cured by subjecting it to a zero-mean Gaussian hyperprior. As an illustrative example, the heatmap for the NLL including this additional regularization term for σ C is shown in Figure 5.

b. Gluon
In the case of the gluon spectral function, no real-time result in Yang-Mills theory is available to fix the asymptotics. However, as an additional theoretical constraint we require the solution to respect the aforementioned OZS condition defined in Equation (3). While one might expect this to further complicate the reconstruction, it actually helps in narrowing down the space of plausible solutions. The condition can simply be enforced approximately by treating it as an additional indirect observation and checking it a posteriori. The associated transformation is here just the convolution with ω instead of the KL integral. We confirm that the OZS condition is fulfilled with a relative accuracy of ∼1%, computed by evaluating the ratio of the left-hand side of Equation (3) and the same expression using the modulus of the integrand, i.e. ∞ 0 dω |ωρ A (ω)|. As mentioned in the main text, we find it helpful to modify the standard RBF kernel by non-linearly rescaling the frequency as ω →ω = ω 4 (1+ω 4 ) −1 before computing the squared distance. This leads to a strongly improved asymptotic stability of the reconstructed spectral function, in particular at large frequencies, compared to just using ω itself. The procedure may be interpreted either as a non-stationary modification of the kernel or as a preprocessing step for the data to the same effect.