Non-perturbative strong coupling at timelike momenta

We compute the strong coupling constant of Landau gauge QCD in the full complex momentum plane, both directly and via spectral reconstruction. In particular, we consider the Taylor coupling given by the product of ghost and gluon dressing functions. Assuming spectral representations for the latter, we first show that also the coupling obeys such a representation. The subsequent spectral reconstruction of the coupling data, obtained from 2+1 flavour lattice QCD results for the ghost and gluon, is based on a probabilistic inversion of this representation using Gaussian process regression with analytically enforced asymptotics. In contradistinction, our direct calculation relies on earlier reconstruction results for the ghost and gluon spectral functions themselves, as well as data obtained in functional QCD. Apart from its relevance for studies of resonances or scattering processes, the calculation also serves as a non-trivial benchmark of our reconstruction approach. The results show remarkable agreement, testifying to the reliability of the method.


I. INTRODUCTION
Real-time correlation functions are of great importance in the calculation of timelike observables in quantum chromodynamics (QCD). Physical scattering processes or the hadronic resonance spectrum represent prominent examples requiring first-principle input in the form of fundamental correlation functions in Minkowski spacetime. The strong coupling constant of QCD is a central ingredient in any of those, as it describes the interaction strength between the fundamental fields. One of its most salient features is asymptotic freedom, i.e., the decay towards small distances, which is well captured by perturbation theory. In contrast, the large distance or low energy behaviour, where the coupling grows large, can only be described via non-perturbative approaches such as lattice field theory or functional methods. While these approaches have proven effective for the calculation of Euclidean correlation functions, they are not welldeveloped in Minkowski spacetime.
Recently, progress has been made through the spectral functional approach [1,2], enabling a direct real-time formulation of functional methods. For applications to QCD and gravity, see [3][4][5][6]. On the lattice field theory side, direct real-time calculations are plagued by a severe sign problem. However, Minkowski correlation functions may also be obtained indirectly via spectral reconstruction of Euclidean data. This requires inverting the Källén-Lehmann (KL) spectral representation [7,8]. The applicability of Gaussian process regression (GPR) to inverse problems of this type was discussed in [9]. The method has been used for the reconstruction of ghost and gluon correlators [10], the computation of glueball masses [11], and similar problems in high energy physics [12][13][14].
In this work, we establish a spectral representation for the strong coupling constant and compute its spectral function. The calculation is facilitated by the reconstruction results for the ghost and gluon spectral functions of [10], based on propagator data from 2+1 flavour lattice QCD calculations with domain wall fermions at a physical pion mass [35,36]. In doing so, we improve on the previous reconstruction approach by incorporating known asymptotic behaviour into the GP kernel. Based on these data, we also apply GPR directly to the reconstruction of the Taylor coupling. This non-trivial benchmark of our reconstruction method yields remarkable agreement between the direct and indirect results, thereby making a strong case for the reliability of spectral reconstruction via probabilistic inversion with GPR. On the other hand, our results feature a broad range of applications in the calculation of physical observables. Knowledge of the coupling constant in the full complex plane is required, e.g., in the treatment of hadronic bound states via Bethe-Salpeter equations. Furthermore, in the calculation of physical scattering amplitudes, the strong coupling in Minkowski space is a necessary ingredient. Our main result, the spectral function of the strong coupling constant from both GPR reconstruction and via its spectral representation, is shown in Figure 2b.
This paper is organised as follows. In Section II, we derive the spectral representation for the strong coupling constant. The extension of our spectral reconstruction approach granting improved control over the asymptotics is described in Section III. Our results are presented in Section IV and we conclude in Section V. q q q q FIG. 1. qq-scattering process with a one-gluon exchange. At sufficiently large timelike exchange momenta, this process plays an important role in its respective S-matrix elements. Consequently, all internal quantities are dressed. Blue blobs represent full vertices, and the wiggly internal line is a full gluon propagator.

II. SCATTERING PROCESSES & THE TIMELIKE QCD COUPLING
Scattering processes and decays in QCD are described in terms of S-matrix elements. At low energies, the operators of the physical in and out states are complicated objects in terms of the fundamental QCD degrees of freedom. For instance, a description of the Compton scattering of protons requires the definition of the proton or, more generally, the nucleon operator in terms of its partonic constituents. Since, on the fundamental level, the partons are related to quarks and gluons, the building blocks of the respective S-matrix elements are quarkgluon and quark-photon scattering processes.
In most partonic models the fundamental scattering processes are approximated by effective models for the exchange process, such as one-gluon exchange potentials that carry the qualitative property of the gluon mass gap in QCD in terms of an effective mass. Ideally, however, they should be constructed from tree-level processes in QCD with full propagators and vertices, both of which carry on-shell, timelike, and spacelike momenta. The final S-matrix is gauge-invariant, while the tree-level components making up the individual S-matrix element contributions are not. Moreover, the S-matrix admits a spectral representation, which is not necessarily present for the gauge-fixed correlation functions.
A. Cross-section of quark-anti-quark scattering events and the S-matrix element In the present work, we undertake a first step towards such a determination of non-perturbative S-matrix building blocks in QCD. To that end, we compute the timelike strong coupling in 2+1 flavour QCD that governs the quark-anti-quark scattering process depicted in Figure 1. This diagram is at the core of many of the scattering processes used to determine the strong running coupling, It is also one of the fundamental building blocks of scattering processes in the Pomeron model [37][38][39][40]-such as the aforementioned Compton scattering of the protonwhere it is typically estimated by one-gluon exchange models. For a review, see [41]; for a recent application related to the present work, see [42].
Assuming the incoming and outgoing quarks q(p) and anti-quarksq(p) to be on-shell, qq-scattering is similar to e + e − scattering. We expect this analogy to hold for sufficiently large timelike exchange momenta p 2 > ∼ 1 GeV 2 , whereas for p 2 < ∼ 1 GeV 2 we enter the hadronic, strongly correlated regime. There, the non-trivial embedding of the scattering quarks and anti-quarks in hadrons becomes increasingly relevant, and quark-anti-quark scattering should be also considered off-shell alongside with further, more complicated processes; for a formulation in functional approaches, see [39].
Here, we concentrate on the one-gluon exchange diagram as one of the building blocks of the full S-matrix element. The associated tree-level process shown in Figure 1 consists of two full quark-gluon vertices, Γ qqA (p 1 , p 2 , p) with the on-shell momenta p 1 , p 2 for the incoming as well as Γ (3) qqA (p 3 , p 4 , −p) with on-shell −p 3 , −p 4 for the outgoing quark and anti-quark, respectively. The relative minus sign is due to the notational convention in functional computations treating all momenta as incoming. The momentum p is that of the exchange gluon with the full gluon propagator G A (p). In combination, this process can be expressed as where the (on-shell) quark wave functions Z q originate in the LSZ reduction formula. Note that the quark and gluon wave functions are defined such that the quark and gluon propagators G q (p), G A (p) are proportional to 1/Z q (p), 1/Z A (p), respectively. The scalar parts of the Euclidean propagators read where the full propagators are proportional to the identity in color space in the adjoint (gluon) and fundamental (quark) representations. The gluon propagator in the Landau gauge also carries the projection operator on the transverse subspace (see (2)), and the quark propagator is multiplied by i / p + M q (p). With (3), the standard LSZ factors carrying the pole residues are simply Z −1/2 q , as already used in (2). The S-matrix element (2) is renormalisation group (RG) invariant, as required. To see this explicitly, we reparametrise the vertices in terms of wave functions of the legs and an RG-invariant core, whereΓ (3) qqA has the transformation properties of a running coupling, and naturally occurs in the S-matrix element. Inserting (4) into the S-matrix element (2) leads us to We restrict ourselves to the limit of large transfer momentum p 2 ≡ s of the scattering event with p 1 p 3 = p 2 p 4 = s(1 − cos θ)/4 and scattering angle cos θ = p 1 p 3 /(|p 1 ||p 3 |). For small s, we approach the reliability limit of our approximations. We return to the respective discussion after deriving our results. Additionally, in a last approximation step we concentrate on the classical tensor structure γ µ T a in the full quark-gluon vertex, Here, T a is the SU (3) generator in the fundamental representation and α s (s), defined in (1), is the strong coupling of the quark-gluon scattering process in the schannel. On the equation of motion, the / p i terms vanish, and we obtain in the high energy limit. In (7), we have performed an average/sum over spins and color in the initial/final state. With (5) and (7), we arrive at with α s (p) defined in (1). Equation (8) highlights the importance of the strong coupling constant α s (s) for physical scattering processes. For the remainder of this work, we adopt the linear momentum argument p = √ s for the coupling.
In the present work, we shall compute the strong coupling α s (p) and, hence, the above S-matrix element from its spectral representation for general complex frequencies, including the timelike momenta relevant for (8). We utilise that the strong coupling can be computed from the quark-gluon vertex, the three-and four-gluon vertices, as well as the ghost-gluon vertex. The computation involves the wave functions Z q (p), Z A (p) of quarks and gluons as defined in (3) and the ghost wave function Z c (p) from The avatars of the strong couplings are then defined as the (symmetric point) dressings of the classical tensor structures, see (4) and (6). A final definition of the strong coupling in the Landau gauge is given by the propagator or Taylor coupling, that utilises Taylor's non-renormalisation theorem for the ghost-gluon vertex. This leads to the Taylor coupling, solely defined by the ghost and gluon dressing functions, All strong coupling avatars have the same universal twoloop running but differ for infrared momenta; see [43].
For an evaluation of the infrared differences between the Taylor coupling and the quark-gluon coupling, see [45]. The latter regime is not accessible within the present approximation. Hence, we use the Taylor coupling (10) for the evaluation of (8). Its corresponding spectral function ρ α is depicted in Figure 2. It allows us to compute the coupling α s (p) for complex frequencies including timelike momenta; see Figure 4. Timelike result for the strong coupling in the perturbative domain can be found, e.g., in [46,47].  (Figure 2b). We compare the spectral function computed directly via (16) (red) to that obtained via reconstruction with GPR (blue). The direct calculation uses the reconstruction results for gluon and ghost spectral functions from [10]. For the reconstruction, we use the gluon and ghost propagator data in 2+1 flavour lattice QCD from [35,36]. Both the input spectral functions and the corresponding lattice data are displayed in Figure 3. The coupling spectral functions obtained via these two complementary approaches share all qualitative features, such as peak positions and heights as well as asymptotic behaviour. The peak structure can be connected to the respective peak structure of the gluon spectral function; see Figure 3b. The error band of the reconstruction result accounts for the change in the spectral function when varying the GP kernel parameters, whereas that of the direct calculation originates from propagating the uncertainty of the input. The Euclidean lattice data for the Taylor coupling αs are displayed as gray squares in Figure 2a. We compare it to the data from its spectral representation (14) (red) as well as the reconstruction result (blue), showing that the representation holds and that the reconstruction accurately describes the lattice data.

B. Spectral representation
For the computation of (10), and hence of (8), we require the ghost and gluon propagators for timelike momenta. We assume that the propagators admit a KL representation, where we have implicitly defined the KL kernel K(p, λ). The spectral function ρ is defined via with ρ(ω) = −ρ(−ω). For the propagators of physical particles, the spectral function is the probability density for (multi-)particle excitations to be created from the vacuum in the presence of the corresponding quantum field. Consequently, in this case, the spectral function is positive semi-definite and normalisable. For propagators of 'unphysical' fields, such as gauge fields, positive semi-definiteness is no longer required and the spectral representation reduces to a statement about the analytic structure of the corresponding correlation function; see, e.g., [3,4,32,[48][49][50][51].
The ghost propagator is known to exhibit a massless particle pole in the origin, entailing a delta pole at vanishing frequency in its spectral function ρ c [3]. The gluon spectral function ρ A is continuous along the whole real frequency axis and is not expected to show distributional contributions. Taking into account the explicit forms of the spectral functions, the ghost and gluon dressing functions can be expressed as where 1/Z 0 c is the residue of the massless delta pole of ρ c , andρ c denotes the continuous part.
Given the existence of a spectral representation, the associated correlation function must obey certain symmetries and fulfill requirements about its infrared (IR) and ultraviolet (UV) asymptotics. It can be shown that the existence of spectral representations for the ghost and gluon propagators implies the existence of such a representation also for the Taylor coupling as defined in (10); see Appendix A for details. Specifically, it is given by With (14), the strong coupling spectral function is obtained from its retarded correlator via  [3] (see also Section II C), used here as input for the calculation of the coupling spectral function shown in Figure 2a via its spectral representation (14). Shaded areas represent 1σ-bands of the statistical error of the mean prediction based on the available observations and precision. Note that for the calculation of the gluon spectral function, the UV and IR asymptotic regimes are assumed to be maximally large. This leads to a small reconstruction error without accounting for systematics; see Appendix C 3 for a detailed discussion.
Now we use the definition of the Taylor coupling (10) and insert the spectral representations of ghost and gluon dressing functions (13). Then, the spectral function (15) of the coupling can be written as Since the Taylor coupling decays logarithmically in the UV, its spectral function obeys a superconvergence condition [3,51], given by In the case of the gluon propagator, this is the well-known Oehme-Zimmermann condition [52,53]. A treatment of the analytic low-frequency behaviour of continuous parts of the spectral functions has been initiated in [32]. In particular, it was shown that for correlation functions obeying a KL representation, a simple relation between the IR asymptotics of the correlator and its spectral function can be derived by differentiating with respect to the frequency. For the Taylor coupling we explicitly find Hence, if the coupling approaches zero in the origin faster than p 2 , we expect the spectral function to approach zero from below, and vice versa.

C. Lattice data
During the past two decades, lattice QCD results for Landau gauge two-point functions have advanced to an impressive quantitative level of precision; see, e.g., [54][55][56][57][58][59][60][61][62][63][64][65][66][67][68][69]. A recent review of lattice and functional results can be found in [70]. The lattice ghost dressing function and gluon propagator data used in this work have been obtained from recent calculations with 2+1 dynamical fermion flavours at the physical point [35,36]. In particular, the ensembles of gauge configurations were generated by the RBC/UKQCD collaboration in [71][72][73][74][75], leveraging the Iwasaki gauge action [76] and the domain wall fermion action [77,78] with a pion mass of 139 MeV. This choice of action (with a particular implementation of the Möbius kernel [79]) exhibits favourable chiral properties with a much smaller size in the fifth dimension than required by conventional domain wall fermions. These ensembles were utilised in [35,36] for the calculation of the ghost and gluon propagators as well as the running of the strong coupling in the Taylor (miniMOM) scheme [80][81][82], and an associated effective charge [83].
The continuum limit of the lattice data is only obtained with a proper treatment of discretisation effects. For the Landau gauge propagators this is achieved by an analysis of the physical scaling violation as described in [84], leading to continuum extrapolated propagators with the correct momentum running. The resulting gluon propagator and ghost dressing data are displayed in the insets of Figure 3. These data combined with results from functional Yang-Mills theory and QCD [43,45,85,86] have also been reconstructed in [10].
Since the lattice data for the propagators are available only on different momentum grids, the coupling as defined in (10) is computed from interpolations of the respective dressings. These are obtained by direct GPR and therefore assume no general features of the underlying correlators apart from continuity. We compute the coupling including errors for 600 logarithmically spaced points between 0.23 GeV and 2.69 GeV. For technical convenience, the coupling is extended perturbatively at large momenta in order to control the amplitude of the UV asymptotics; see Appendix B. A subset of these data is shown in Figure 2a. Here, we replace the error with the difference between the values computed as described above, and the coupling obtained from the product of the ghost and gluon spectral functions, described around (16) and in Section IV.

III. GPR RECONSTRUCTION WITH CONTROLLED ASYMPTOTICS
GPR is a popular framework for the probabilistic modelling of functions from a finite number of data points; see [87,88] for recent reviews and [89] for a textbook account. Example applications in high energy physics include the computation of parton distribution functions [90] and modelling backgrounds in detectors [91]. The method can be used to predict solutions to linear inverse problems [9], i.e., when the only available data are indirect observations of the desired function after a linear forward process. This makes the approach suitable for spectral reconstruction. Importantly, it does not in general require choosing a particular functional basis. This avoids many of the numerical artifacts like additional peak structures that are commonly encountered when employing reconstruction algorithms with predetermined families of solutions, due to the presence of unrepresentable features. We summarise the main concepts in Appendix C 1; see also [10,11] for a comprehensive introduction as well as further details and references.
As an extension to this approach, in this paper we introduce a novel technical improvement that allows us to explicitly control the asymptotic behaviour of the predictions by specifying concrete functional forms in the appropriate limits only, without restricting the expressivity of the GP model in the region of interest. When considering different design choices for GPs, one often opts for so-called universal kernels. One prominent example also used in the present work is the radial basis function (RBF) kernel (C4). The basis of kernel eigenfunctions of such universal kernels is infinite-dimensional. This allows for great flexibility in the reconstruction-universal kernels can describe any continuous function [92]. However, the GPR framework allows us to also incorporate further available prior information into the predictive distribution. In the context of spectral functions, the asymptotics in the IR and UV are often analytically tractable with perturbative or functional calculations, as well as formal relations to Euclidean data like (18). Hence, it is beneficial to introduce a bias by reducing the space of kernel eigenfunctions to the known behaviour of the target function. This can be achieved by applying Mercer's theorem [93] and constructing a kernel from the known asymptotic function φ(ω) as Since the asymptotic behaviour is only specified in the appropriate limits, the full kernel is constructed as a combination of universal and restricted kernels using smooth step functions. With this approach, it becomes possible to smoothly transition between regions with an unknown functional basis-where a generic kernel like RBF is used-and regions with a specified basis; see Appendix C 2 for further details.

IV. RESULTS
Our main result, the spectral function of the Taylor coupling (15) in QCD, is displayed in Figure 2b. It shows two variants: ρ GP α from the reconstruction of the lattice QCD data via GPR, and ρ spec α from the direct calculation based on the spectral representations of ghost and gluon propagators (16). The associated input spectral functions are shown in Figure 3. In this context, we have improved the reconstruction of the gluon propagator reported in [10] by explicitly incorporating the known IR and UV asymptotics with the method described in Section III. The error band of ρ spec α is obtained by propagating the errors of these input data. Importantly, the coupling spectral functions obtained via these two different approaches agree well within errors and share all qualitative features, such as peaks and asymptotic behaviour. In both results, we can identify two prominent peaks of similar size in positive and negative direction at roughly ∼ 0.6 GeV and ∼ 0.8 GeV, along with a smaller positive peak at ∼ 1.1 GeV. The spectral function ρ spec α (16) allows for a direct interpretation of this behaviour: it is connected to the peak structure of the gluon spectral function, which carries information about the gluon mass gap; see Figure 3b. This information is extracted reliably from the lattice data with the GPR reconstruction.
In the reconstruction of the coupling, the correct asymptotic behaviour is enforced by smooth step functions at transition points µ IR and µ UV , while fully retaining the flexibility in the enclosed region where the GP kernel remains unrestricted and universal, see Appendix C 2 for details. As mentioned above, this procedure has also been applied to the ghost and gluon spectral functions used here. It enhances significantly the stability and reliability of the prediction by connecting it to analytic results at low and high frequencies, ensuring agreement with functional and perturbative results in the relevant limits without reducing the expressivity of the GP model in the domain of interest. While the prediction shows some variation with the choice of the transition midpoints, the peak positions and heights remain remarkably stable; see Figure 5. Hence, we choose  (Figure 4b). The imaginary part explicitly shows the branch cut along the real frequency axis. The spectral function corresponds to the imaginary part of αs at the upper half plane boundary of the branch cut, divided by ω 2 . Both, the real and imaginary part, exhibit distinctive peaks which can be connected to the peak structure of the gluon spectral function; see Figure 3b. The coupling decays logarithmically for increasing |ω|.
the size of the regions dominated by the asymptotics to be as large as possible without increasing the χ 2 error of the reconstruction significantly; see Appendix C 3 and Figure 6 for details. Furthermore, changing the parameters controlling the transition to the asymptotic behaviour accounts for the majority of the variation in the spectral function, while changing the parameters of the RBF kernel produces errors at least one order of magnitude smaller; see Appendix C 3 for details. The numerical values of all kernel hyperparameters are listed in Table I. Accordingly, the size of the dynamical region carrying information about the QCD mass gap is minimised, supporting the gluonic quasi-particle picture employed in various applications such as bound state studies and transport computations. Specifically, this suggests dismissing smaller negative peaks close to the dominant quasi-particle peak-they merely reflect the asymptotic behaviour and the superconvergence condition (17). As such, they are sensitive to changes in the gauge fixing parameter and infrared closure. This suggests that they carry physically relevant information only on a subleading level.
In Figure 2a, we compare the reconstructed Euclidean Taylor couplings to the result computed from the lattice data for the ghost and gluon propagators, as described in Section II C. Using the dressing function data obtained in this way, the resulting coupling is shown to decay towards small and large momenta. In correspondence to the scale of the peaks of the spectral function-reflecting the mass gap of the theory-also the peak of the coupling itself appears at ∼ 0.6 GeV.
The blue curve in Figure 2a represents the GPR reconstruction of the Taylor coupling lattice data, corresponding to ρ GP α . The red curve represents the coupling obtained via its spectral representation (14) using the di-rectly computed spectral function ρ spec α . The calculation involves finite precision, both in the input data and in the integration. Hence we expect a small, but not negligible, relative error. The decent agreement between this result and the lattice/GPR reconstruction result provides a highly non-trivial benchmark check. The error is well within our expectations, since the result obtained from the directly computed spectral function depends on the reconstructions of the gluon and ghost propagators. If the ghost and gluon spectral functions were describing their respective propagator data to infinite precision, we would also expect perfect agreement from analytic considerations; see Appendix A. Hence, the small difference can be attributed to systematic uncertainties present in the calculation. Please note that they do not contribute to the error bands, corresponding to the purely statistical error, shown in Figure 2.
In the inset of Figure 2a, we also show the Taylor coupling divided by p 2 for small Euclidean momenta p. The derivative of this quantity is connected to the asymptotic behaviour of the spectral function in the IR by (18). We observe that in the region where lattice data are available, the slope of α s /p 2 is negative. In accordance with the analytic requirement (18), the slope of the spectral function is observed to be positive in this regime.
Finally, in Figure 4 we display the real and imaginary parts of the coupling in the full complex momentum plane. The data are obtained by evaluating the coupling spectral representation (14) with the directly calculated spectral function ρ spec α in the complex plane. The branch cut in the imaginary part, responsible for the spectral representation, is clearly visible. As expected, no further non-analyticities in the complex plane are encountered and the coupling shows the expected decay behaviour towards large frequencies.

V. CONCLUSION
In this work, we have presented results for the spectral function of the strong coupling constant in QCD obtained through a direct calculation as well as a reconstruction via GPR. Assuming spectral representations for the ghost and gluon, we have derived the spectral representation of the Taylor coupling, which is fully determined by the ghost and gluon dressing functions. With this relation, we have calculated the associated spectral function as well as the coupling itself in the full complex plane; see Figures 2 and 4. The required ghost and gluon spectral functions have been obtained using the same reconstruction method, explicitly taking into account the known asymptotic IR and UV behaviour; see Figure 3. This is facilitated by expanding the GP kernel in suitable eigenfunctions based on Mercer's theorem, extending the algorithm previously applied in [10,11]. The modification substantially improves the reliability of the approach by properly encoding the analytically tractable regimes into the prediction while preserving the expressivity and universality of the GP model in the region of interest. This extension of GPR is completely generic and may also be useful in other contexts where some analytic properties of a function to be modelled are known a priori, in particular if data scarcity is an issue.
A comparison of the results from the direct calculation and GPR reconstruction shows excellent agreement between both approaches; see Figure 2b. This independent verification provides strong support for the accuracy of the computed spectral function and also underlines the power of probabilistic inversion with GPR as a spectral reconstruction approach. In particular, the findings demonstrate that uncertainty estimates obtained within this framework are reasonable, allowing to reliably quantify the expected errors in potential downstream applications based on reconstruction results. In future work, we plan to also analyse the quark-gluon vertex coupling directly based on available lattice data; see, e.g., [94].
Our results find direct application in the calculation of non-perturbative, physical scattering processes, where the strong coupling constant needs to be known at timelike momenta. While neglecting angular dependencies, the Taylor coupling considered here carries the correct RG running and hence scale-dependence of the strong coupling constant. Furthermore, it encodes genuine nonperturbative information through the input ghost and gluon dressing functions obtained from 2+1 flavour lattice QCD. Our work hence paves the way for incorporating non-perturbative information from lattice field theory to functional methods in the calculation of timelike scattering processes.

ACKNOWLEDGMENTS
We thank J. Papavassiliou and J. Rodríguez-Quintero for discussions and collaboration in the early stages of this work. We also 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 Any product of correlation functions obeying a spectral representation allows for such a representation itself. This follows from a set of sufficient conditions for the existence of a spectral representation for an arbitrary correlation function C: (i) Holomorphicity: C is holomorphic in the upper half plane H = {z | Im z > 0}; (ii) Mirror symmetry: C(z) =C(z) and Im C(z) = 0 for Im z = 0, Re z > 0; (iii) Asymptotic decay: |z C(z)| → 0 for |z| → ∞; (iv) Spectral convergence: Heuristically, (i) and (ii) guarantee that the spectral kernel has the form 1/(z + λ 2 ) and the spectral function is defined via (12). The integration domain is restricted to λ 2 > 0 by (iii). Condition (iv) guarantees the convergence of the spectral integral. It is immediately clear that for any two correlation functions C 1 , C 2 satisfying (i)-(iii), their product C = C 1 C 2 does as well. Similarly, this also applies to (iv) (UV), stating that the spectral function ρ ∼ Im C decays fast enough for the spectral integral to converge in the UV, due to (iii). The infrared convergence condition (iv) (IR) does not need to be fulfilled; consider, e.g., C 1 = C 2 = (1/z) α with 1/2 < α < 1. Nevertheless, this can be always remedied by multiplying with an appropriate power of z. Note that this does not violate the other conditions.
The spectral representation for the strong coupling constant is then constructed as follows: by the assumption of ghost and gluon propagator obeying the KL representation, their dressing functions obey (i) and (ii). Since by its definition (10) the coupling is dimensionless, (iii) does not hold. However, division by p 2 makes (iii) and (iv) hold true. Hence, a KL representation for α s (p)/p 2 is constructed. Multiplying this representation by p 2 , we obtain the spectral representation for α s (14).

Appendix B: Asymptotic behaviour of the strong coupling
For the gluon and the ghost propagators, the leading order IR and UV asymptotics are known analytically; see [32] and references therein. In the infrared, the decoupling solution of the ghost is characterised by a constant propagator dressing Z c ≡ Z c (p = 0). On the other hand, the gluon propagator is dominated by the ghost loop polarisation diagram in the IR, since the gluon propagator itself has a mass gap and decouples in the infrared. This results in a p 2 log p 2 contribution in the IR regime; for a detailed discussion thereof, see [32].
Using the definition of the strong coupling (10), we see that it has the same, but negative, IR behaviour as the inverse gluon dressing, up to a constant contribution from the ghost dressing. From (18), we can then infer the asymptotic behaviour of the spectral function as analogously to [32]. The UV asymptotic behaviour of the strong coupling is well known from perturbative calculations and reads .
The asymptotic behaviour of the spectral function follows directly from (15) and we obtain . (B3) Appendix C: Reconstruction details

GPR basics
Here, we briefly summarise the main aspects of the GPR reconstruction procedure. For a more detailed overview, we refer to earlier works [10,11].
We assume our knowledge of the spectral function ρ(ω) before making observations of the correlator to be described by a GP prior, written as ρ(ω) ∼ GP(µ(ω), C(ω, ω)) , where µ, C denote the mean and covariance. The conditional posterior distribution for ρ(ω) given observations of the propagator G i at N G discrete Euclidean frequencies p i ≡ [p] i can be derived in closed form, This is essentially equivalent to a standard result in probability theory for the closed-form expression of a conditional multivariate normal distribution, but defined with a continuum of random variables due to being a Gaussian process, as well as additional applications of the integral transformation one seeks to invert. The equivalence becomes more concrete in practice when the GP is evaluated for a finite set of predictions; however, the choice of inference points ω is arbitrary within the given domain.
In the above expressions, µ(ω) has been set to zero, since a GP can be fully specified by its second-order statistics and the prior mean can be absorbed into C. The GP in (C3) encodes our knowledge of the spectral function after making observations of the correlator and accounting for observational noise with variance σ 2 n . The covariance C(ω, ω ) is commonly defined via a socalled kernel function with a small number of hyperparameters, which may be subject to optimisation based on the associated likelihood. A widely used parametrisation is the radial basis function (RBF) kernel, defined as where the parameter σ C controls the overall magnitude and l is a generic length scale.

Incorporating asymptotic information
With the knowledge of the IR and UV asymptotics, cf. Appendix B, an appropriate bias can be introduced. It is chosen such that the kernel is restricted to the specified functional basis as described in Section III, while retaining the flexibility of the RBF kernel in the central region. In order to achieve a smooth transition between the biased and unbiased kernels, we employ smooth step functions of the form The full kernel can then be written simply as a sum of the individual contributions, The midpoints of the transition functions θ IR/UV are specified by µ IR/UV and their steepness is controlled by IR/UV .

GP kernel hyperparameters
Since the hyperparameters of the GP kernel control the behaviour of the resulting spectral function, their choice is a pivotal step in the reconstruction. They are commonly determined via numerical optimisation by minimising (conventionally) the negative log-likelihood (NLL), where the dependence on the kernel hyperparameters σ is emphasised by an index.
The number of hyperparameters increases significantly when including the bias term that enforces the correct asymptotics (C6). Hence, the two parameters of the bare RBF kernel are chosen first by minimising (C7). The asymptotics are then introduced in the far IR/UV and shifted towards the center, all while monitoring the quality of the interpolation of the dressing by computing χ 2 at each step. We compare χ 2 instead of the NLL for different bias parameters, since the second term in (C7) constitutes a complexity penalty term. When considering an explicit functional basis, such a term is inherently in opposition to the constraint for the analytically known asymptotics and is therefore excluded. We observe that the bias kernel parameters have an open direction towards vanishing bias, e.g., for small µ IR and large µ UV . Spectral functions with µ UV > 1.5 GeV show a growing number of smaller oscillations in the UV, which are a remnant of the global length scale introduced in the RBF kernel; see Figure 5b. Accordingly, models in this parameter region can be ruled out as sensible descriptions of the underlying physics of the coupling. For µ IR < 0.25 GeV, the resulting spectral functions do not change significantly; see Figure 5a. The change of χ 2 when varying the asymptotic kernel parameters is shown in Figure 6, with the final settings used for the reconstruction indicated by crosses. These parameters are explicitly chosen to maximise the regions dominated by the asymptotics, without significantly increasing the error of the coupling reconstruction. All final parameters of the GP model used to compute the results reported in this work are listed in Table I the systematic error that arises from the choice of the model, in particular regarding different values of the hyperparameters. However, we observe that enforcing the maximally large asymptotic regimes leads to the predicted posterior covariance being comparatively small as the model is now highly restricted; see Figure 3b. Hence, the error is estimated by varying the bias parameters in a region where χ 2 is small, but the effect of different parameter choices is non-negligible, while unphysical oscillations remain largely suppressed. This region is marked in Figure 6 by horizontal bars. When considering µ UV larger than indicated in this region, a substantial amount of oscillations is introduced in the spectral function as mentioned above. For µ IR smaller than indicated in Figure 6a, the spectral function vanishes in the IR. However, it can then differ from the expected ω 2 behaviour. The largest variations in the resulting spectral functions under these changes of the hyperparameters are then used as error estimates for the reconstruction results, shown in Figure 2. Since the deviations of the predictions at the edges of the parameter space tend to be maximised in particular regions for certain parameter combinations, e.g., when µ UV and UV are both small, the error band shows a few distinct kinks.