Spectral representation of lattice gluon and ghost propagators at zero temperature

We consider the analytic continuation of Euclidean propagator data obtained from 4D simulations to Minkowski space. In order to perform this continuation, the common approach is to first extract the K\"all\'en-Lehmann spectral density of the field. Once this is known, it can be extended to Minkowski space to yield the Minkowski propagator. However, obtaining the K\"all\'en-Lehmann spectral density from propagator data is a well known ill-posed numerical problem. To regularize this problem we implement an appropriate version of Tikhonov regularization supplemented with the Morozov discrepancy principle. We will then apply this to various toy model data to demonstrate the conditions of validity for this method, and finally to zero temperature gluon and ghost lattice QCD data. We carefully explain how to deal with the IR singularity of the massless ghost propagator. We also uncover the numerically different performance when using two ---mathematically equivalent--- versions of the K\"all\'en-Lehmann spectral integral.


Introduction
A plethora of practical Quantum Field Theory calculational tools, both analytical and numerical, have been developed in a Euclidean setting, despite living in a Minkowski spacetime. In particular, non-perturbative approaches to Quantum Field Theory tend to rely solely on the Euclidean formulation of the theory due to its technical advantages. This is also the case for the discretized lattice approach to Yang-Mills theory, or most of its continuum studies using the Dyson-Schwinger equations, although there are also examples of solving the Dyson-Schwinger equations directly in Minkowski space [1,2,3]. The Euclidean formulation accesses only space-like momenta and, therefore, phenomena that are associated with time-like momenta cannot be investigated within Euclidean Quantum Field Theory. Moreover, whereas in perturbation theory the analytical continuation of the Euclidean correlation functions into the entire complex momenta Argand plane relies on the usual Wick rotation, it is not clear that the same rule can be applied for the nonperturbative regime. Oftentimes the analytical structure of the perturbative correlation functions are known, allowing to go from the Minkowski to Euclidean spacetime and vice-versa in an unam-biguous way. However, for the non-perturbative regime, the analytical structure of the correlation functions is difficult to determine per se and the analytical continuation of the Euclidean space correlation functions into the complex plane becomes a much harder problem. Additionally, with numerical data in particular, we only have a discrete set of data at our disposal which we wish to analytically continue to the complex plane, while it is well-known that such a continuation is only unique when departing from a function known over an open subset of C.
For two-point correlation functions, a possible strategy is to compute the Källén-Lehmann spectral density ρ(ω) from the Euclidean data and reintegrate ρ(ω) to determine the propagator for Minkowski p 2 ≤ 0. The spectral density encodes information about the spectra and, at finite temperature, it can be related to transport coefficients and thermodynamical properties of e.g. the quark-gluon plasma [4]. In this paper we discuss the computation of the spectral density from Euclidean data, focusing on the inversion of the propagator data and on how to improve the method set out in [5].
One might wonder why it is relevant to investigate the spectral representations of gauge variant degrees of freedom like gluons, ghosts or quarks in the first place, since in general their propagators will be gauge dependent and therefore so will the corresponding spectral densities. But because the physical bound state spectrum is gauge invariant despite being constructed from gauge variant gluon, ghost and quark propagators [6], the spectral function of a bound state propagator must therefore be nontrivially influenced by the analytic structure of the underlying constituents, see for instance [6,7,8,9,10,11] for discussions and examples. This explains the relevance of various studies devoted to the spectral properties of a priori unphysical degrees of freedom. The Landau gauge in particular is a standard choice to study the physical bound state spectrum of QCD, see [6,7,8,9,11], exactly because there are ample estimates, both from the continuum as from the lattice side, about the necessary Green functions that serve as input for the Bethe-Salpeter equations. A beautiful example is the seminal paper [12], where the Bethe-Salpeter equations for glueballs are considered, making explicit use of the spectral properties encoded in the spectral functions of gluons and ghosts, reported earlier in [13]. Our research output gives a direct lattice verification of how trustworthy these latter gluon and ghost spectral functions actually are.
The story continues at finite temperature, as at some sufficiently high temperature above deconfinement, we expect to achieve some kind of physical quasi-particle behaviour for the spectral functions of deconfined quarks and gluons [14,15]. From the previous perspective, it may also be interesting to note that recent lattice studies went beyond Landau gauge, and saw only a very mild dependence on the gauge parameter for at least gluon and ghost propagators [16,17,18]. This suggests that the gluon and ghost spectral densities will themselves only mildly depend on the gauge parameter, so a study of the lattice Landau gauge spectral functions can teach us some generic properties of confined degrees of freedom, whilst allowing to confront other analytical estimates which are based on certain essential underlying assumptions and approximations which do not plague lattice computations.
The Euclidean two-point function G(p 4 ) can be expressed in terms of its Källén-Lehmann spectral density ρ(ω) as (see e.g. [19,20,21]) with p 4 the imaginary frequency and ω 0 an IR-cutoff, potentially zero. Henceforth, the representation in Eq.
(2) will be referred to as the p 2 -formalism, the strategy followed earlier in [5].
From the antisymmetry property for a spectral density corresponding to a propagator of scalar degrees of freedom [21], ρ(−ω) = −ρ(ω), it follows that Eq. (1) can equivalently be written as This representation will be referred to as the ip-formalism. When presented with a finite data set for G(p 4 ), one can rewrite Eq. (2) or (3) as a matrix equation where K is the respective integral kernel represented as a matrix, and try to solve it for ρ(ω).
Writing K in terms of its singular value decomposition ij u i S ij v † j , where S ij is a rectangular diagonal matrix, the least-squares solution to Eq. (4) is then given by with s i being the singular values of the matrix K. However, since the singular values of the matrix K span a very large range, the matrix is ill-conditioned, meaning that even a small error in the input values can cause a huge variation in the output. Therefore, the inversion problem represented in Eq. (4) has to be regularized in order to be able to compute a solution for ρ(ω). A modification has to be made such that the condition number is reduced and the resulting problems limited.
Tikhonov regularization solves this inversion problem by adding a term proportional to the norm of the solution to the residual sum of squares 1 : In terms of the singular values of K the solution is then given by By comparing with Eq. (5), we find that the Tikhonov parameter α 2 dampens the effect of the smallest singular values s i .
For Tikhonov regularization to work, an appropriate choice for the parameter α 2 has to be made, but there is no unique way to select the value of α 2 . In this paper we will use the Morozov discrepancy principle [22], which seems a good choice to invert lattice data since such data always feature a statistical error σ i for every data point G i . The Morozov discrepancy principle states that α 2 should be selected such that where i σ 2 i is the total variance in the data. The α 2 obeying this constraint is guaranteed to be unique [22] and this choice for the regularization parameter means that the quality of the reconstruction matches the quality of the original data set. It can be shown that in the limit σ i → 0 and N → ∞, the Morozov solution converges to the exact solution [22]. It is also worth noting that Tikhonov regularization using the Morozov criterion can alternatively be understood from a Bayesian approach as "historic MEM", with a default model m = 0. This default model choice is well motivated, as the UV asymptotics of the ghost and gluon spectral functions predict that they will tend to zero, corresponding with a default model m = 0 in the UV. But the Tikhonov functional Eq. (6) is special among different possible choices of prior distributions, as it can be solved analytically and thus numerical difficulties associated with solving non-Gaussian priors can be avoided. This paper reports on applying the procedure outlined above to gluon and ghost two-point functions obtained from lattice QCD. However, before applying the Tikhonov procedure to such lattice data, two different implementations of Tikhonov regularisation are outlined in Section 2, corresponding to solving either Eq. (2) or Eq. (3). Then, both of these methods were applied to three different toy models, the results of which are detailed in Section 3. The first of these toy spectral densities is an everywhere positive distribution, chosen to represent an observable physical particle. Indeed, via the optical theorem the spectral function can be related to an observable probability [20], implying its positive definiteness. By contrast, the second and third toy spectral densities display positivity violations, mimicking the expected behaviour for unphysical (confined) particles such as gluons and ghosts. Lastly, the two methods were applied to lattice gluon and ghost T = 0 data sets, the results of which are given in Section 4.

Survey of the method
The lattice data for the propagator come with known statistical uncertainties and, furthermore, given that the different momenta are computed from the same set of gauge configurations, the different momenta are statistically correlated. These two effects can be taken into account in the variational principle behind the Tikhonov regularization scheme, replacing (6) by the new minimizing functional 2 where Σ is the covariance matrix. For lattice data, the covariance matrix can be computed from the different gauge configurations as where N Conf is the number of gauge configurations used to compute the propagator, G k (p i ) is the propagator for gauge configuration k at momentum p i , and G(p j ) is the lattice estimation for the propagator. However, for both the gluon and ghost lattice data we have checked that the covariant matrix is an almost diagonal matrix and, therefore, herein we shall only consider a diagonal covariance matrix Σ ij = σ 2 i δ ij (no sum), where σ 2 i is the variance of G(p i ): Solving for the minimum of function (9) involves computing ∂J α /∂ρ = 0, which results in Defining c := Kρ − G, then gives However, since c itself depends on ρ, substituting Eq. (13) into the definition of c, the following linear system is found: Solving this linear system for c at a given value of α, the spectral function ρ can be reconstructed using Eq. (13). From the above definitions, it follows that the reconstructed propagator written in terms of c reads

The p 2 -formalism
Starting from Eq. (2), the Tikhonov functional (9) becomes 3 Upon repeating the functional equivalent of the steps taken above, the following expressions for M andρ are found:

The ip-formalism
Repeating these steps but starting from Eq. (3) instead, yields resulting in It is worth noting that as ω 0 → 0, 3 Note that we are not taking into account the correlation between the different momenta, just the variances.
This is proportional to the M found when inverting a Laplace transform [23], whereas the spectral representation in Eq. (1) can be viewed as a double Laplace transform. This will become relevant when discussing the observed difference in reconstruction quality between the two methods. Although both representations are mathematically equivalent, the associated inversion procedures perform differently at the numerical level.

Construction of toy models
In principle, any of the above methods can be used to compute the spectral function and from it rebuild the propagators. However, from the numerical point of view, given the different characteristics of the matrix M = K K T , the two procedures can behave quite differently. Therefore, before applying the inversions to the reconstruction of the propagators, we investigate their performance on three toy models.
A Breit-Wigner type model with M = 3 and γ = 1 (dimensionless). This toy model for the spectral function was investigated in [24]. A similar type of toy model, albeit with a wider peak, was also used in the paper [5].
A "Bessel" model without IR-cutoff, with J n (ω) being the Bessel functions of the first kind, was constructed to obey the same sum rule as gluons and ghosts are supposed to obey, namely ∞ ω0 ρ(ω)ω dω = 0.
This model is extended to ω < 0 by demanding ρ(−ω) = −ρ(ω). In Appendix A we have recollected the argument why the spectral functions of the gluons and ghosts, assuming the associated propagators have a Källén-Lehmann spectral representation to begin with, must obey the sum rule (25); see [25,26,6,27] for further reference. In short the sum rule (25) can be obtained using the large momentum behavior of the propagator, whereto, thanks to asymptotic freedom, perturbation theory applies. Obviously Eq. (25) implies that the spectral function can no longer be positive-definite. Consequently, one has to resort to inversion strategies that can accommodate for such spectral functions [28,29,5,30,31,24,32], which excludes e.g. the popular (standard) Maximum Entropy Method [33,34,4]. Alternative methods include Padé rational function approximation [35] or machine learning-based methodologies [36], of which it still needs to be established if these also perform well for unphysical Green functions. Yet another recipe for inversion was proposed in [37] based on analytical insights, but to our knowledge it has never been tested in practice, most likely due to the reason mentioned in [38]: the required precision to obtain sensible results is unrealistic in numerical computations. Besides these numerical approaches, analytical estimates of spectral functions can also be made. Oftentimes these are performed in conjunction with numerical tools; see [13,32,39,4,40,41] and references therein for examples.
A third model, not only obeying Eq. (25) but also featuring an IR-cutoff, was constructed with the spectral function where Again the spectral function was extended to ω < 0 using ρ(−ω) = −ρ(ω).

Data building and analysis for the toy models
In order to mimic the conditions of a lattice simulation, we proceed as follows. For a given spectral density (23), (24), (26), the "propagator" G orig is computed using either Eq. (2) or Eq. (3). From this G orig , N bootstrap data sets G are generated satisfying a Gaussian distribution with mean value G orig and variance ( G orig ) 2 , i.e. the G are distributed according to a probability distribution G ∼ N (G orig , ( G orig ) 2 ), where is the noise level (in percentage) of the samples. Furthermore, for each data set N res momenta are uniformly sampled in the interval p ∈ [−10, 10] for the ip-formalism, of which the range p ∈ [0 , 10] are squared for the p 2 -formalism. The p are dimensionless, as there is no physical scale here. The choice for this particular range is motivated by the fact that all relevant features of the toy models lie within this range.
For each of the G bootstrap samples, the inversion is performed using the two formalisms discussed previously. For each we have used N bootstrap = 1000 samples. In this way the distribution of the optimal α parameter as a function of and N res could be studied. Typically, the distribution of the optimal α is Gaussian or almost Gaussian. The exception occurs for large enough N res where the α distributions show a tail that touches the point α = 0. For these small optimal α values, the initial ill-defined inversion reappears, the inversion fails and the original propagator data cannot be reconstructed. To prevent these pathological cases we require the optimal α distribution to be compatible with a normal law, i.e. that they can be fitted by a Gaussian law. From the practical point of view, this is enough to prevent the small optimal α values.
The quality of the reconstructed spectral function ρ re , defined as the spectral function returned by the inversion methods derived in Eqs. (16)- (22), can be measured from the coefficient of determination defined as where ρ orig are the input data points used to build the propagator andρ orig is the mean value of ρ orig ; because the function ρ(ω) is odd, for the evaluation of R 2 only data with ω ≥ 0 was considered. The coefficient of determination measures how the variation of the dependent variable matches the variation of the independent variable. It has a maximum value of 1, indicating a perfect fit, but is not bounded from below.
The quality of the reconstruction could have been measured using a quantity different than R 2 . For example, in a recent work [24], the authors defined the applicability of the reconstruction method as the ability to find the position of the dominant peak to within 10%. Therefore, we will provide heat maps for R 2 in ( , N res )-space onto which solid black contour lines have been drawn to indicate the accuracy in finding the dominant peak position.
All numerical analysis was performed in the Python language, using the symfit optimization package [42].

Determination of ω 0
In general, a physical cutoff ω 0 should ideally not depend on the choice of regularization, in particular it should not depend on our choice of the Tikhonov parameter α. This means that the variation of ω 0 w.r.t. α should be as small as possible in practice. Since numerically we have easier access to the variation of α w.r.t. ω 0 , the optimum value for α is more easily identified as the regions where the variation of α w.r.t. ω 0 is maximal.
This means that after generating the ω 0 v.s. α curve from our bootstrap, the point with the largest standard deviation in α w.r.t. ω 0 corresponds to the most likely physical cutoff ω 0 . This concept will guide our choices of the optimal ω 0 and hence, the corresponding α(ω 0 ).

Toy model-results and discussion
For all of the toy models described in the previous section, the p 2 -and ip-formalism derived in Equations (16)- (20) will be applied to data sets generated as described in Section 2.4, while setting Σ ij to either σ 2 i δ ij or δ ij during the inversion, where σ 2 i is the variance of G(p i ).

The Breit-Wigner Spectral Function
We start our analysis by looking at the inverse problem for the Breit-Wigner type model given in Eq. (23). This model has no IR-cutoff and the inversions performed here therefore consider ω 0 = 0. This implies that the p = 0 data point has to be excluded from the inversion for both methods. The effect of a non-zero ω 0 will be studied later on.
In Figure 1, the effect of N res at fixed noise = 0.1% is shown for the various methods. A comparison of Figure 1a with Figure 1c, and Figure 1b with Figure 1d, shows that the effects of taking into account the error on the data, i.e. setting Σ = 1, does not visibly change the quality of the spectral function computed from the inversion. Furthermore, Figures 1a and 1b show that the ip-method performs better for this toy model at low momentum scales, as the spectral function evaluated with the p 2 -method has larger oscillations at small ω. For the p 2 -method, the IR oscillations are reduced by increasing N res . The different IR behaviour of the two methods will be discussed quantitatively below.
Another important feature that can be observed from Figure 1 is that for the level of noise considered, = 0.1%, the inversions do not have a strong dependence on the number of data points, N res , taken into account in the inversion. Although the computed spectral functions coming from the inversions are not perfect, it seems that both methods capture the main features of ρ(ω) for any N res . As discussed below, it is the value of that seems to play the most important role in the inversion, with the inverted spectral function getting closer to the exact spectral function as is reduced, as expected. Note also that despite the dependence on N res of the computation of ρ is mild, the match between the computed and input spectral function slightly improves as N res increases. The effects of the noise level is shown in Figure 2 for the two methods and for N res = 128. Other values of N res show similar results. Again, taking into account the statistical errors results in a reconstructed ρ with smaller errors. In general, reducing the noise level results in a computed spectral function that is closer to the exact ρ(ω); see Figure 2a against 2c, and compare Figure 2b with 2d. Additionally, as shown in Figures 2a and 2b, the ip-method provides a spectral function that is less oscillatory in the IR. This is also evidenced by the corresponding increase in the standard deviation for ρ, suggesting that the p 2 -reconstructions are less IR stable. However, concerning the location and height of the maximum of the computed ρ(ω), the ip-method captures the location best, while the p 2 -method seems to better captures the height of the maximum.
The behaviour of the coefficient of determination R 2 as a function of the noise level and number of momenta N res is summarized in methods perform better when the statistical errors on the propagator data are taken into account in the inversion. In general, the value of R 2 is closer to unity for the ip-method. In particular for the smaller values of N res , the p 2 -method oftentimes has R 2 0.5. These results for R 2 can be viewed as an indication that overall the solution provided by the ip-method is closer to the original spectral function.
As can be seen in Figure 3 it is the value of that has a major impact on the reconstruction of the spectral function. If, for example, R 2 > 0.9 is demanded, a noise level of about 0.1% is needed for the ip-method to fulfil this condition, while the p 2 -method requires an 0.05% to achieve the same values of the coefficient of determination R 2 .
The spectral function under analysis was also investigated in [24] using various inversion techniques namely Maximum Entropy (MEM), Backus-Gilbert (BG) and the Schlessinger point method (SP). According to the authors of [24], the MEM, BG and SP are able to locate the maximum of ρ quite well, both its position and height. The SP provides the best reconstructed spectral function but not necessarily for the smallest errors (see their Figure 6). A fair comparison is difficult to perform. Their spectral functions that were reconstructed using BG are too broad and clearly quite far away from the input ρ(ω) (see their Figure 5). Their implementation of the MEM returns a spectral function with large oscillations at small momentum scales, i.e. for ω 220 MeV, and although it provides a good description of the position of the maximum of ρ, it clearly  underestimates its strength. Even though [24] does not calculate R 2 , their graphs suggest that the corresponding R 2 values would be smaller due to lesser overlap between the original and the reconstructions as N res decreases or increases. The best choice of algorithm therefore seems to depend on the desired feature of the data we are trying to capture.
Let us now discuss the effects due to a non-vanishing IR-cutoff ω 0 on the inversion. Recall that, as discussed in Section 2.5, the physical cutoff ω 0 should correspond to the point where the variation of α w.r.t. ω 0 is maximal.
As a first example, the current Breit-Wigner model has no cutoff. As can be observed in Figure 4 for N res = 128, in the low noise limit with both methods the maximal standard deviation is achieved around the maximum of the ω 0 v.s. α curve, suggesting that ω 0 = 0 is indeed the right choice.

The Bessel spectral function
In this section we discuss the results for the inversion when the input function used to generate the propagator data is the Bessel model given in Eq. (24). The interest in this type of function arises from the property ∞ 0 ρ(ω)ω dω = 0 , i.e. ρ(ω) necessarily has regions where it takes positive and negative values and, therefore, mimics a spectral function that can be associated with unphysical particles. Firstly, we consider the inversion with ω 0 = 0, which requires excluding the p = 0 data point to avoid the singularity in the matrix M . Later on we also look at the case where a finite IR cutoff is present.
In Figure 5, the effect of N res at fixed noise of = 0.1% is shown for the various methods. Similar to the Breit-Wigner model, it seems that including the Σ-matrix has no visible advantage at this noise level; compare Figure 5a   with Figure 5b shows, again, much greater oscillatory behaviour in the IR associated with the p 2method. As for the Breit-Wigner model, the inversion shows only a mild dependence on N res , with the main effect associated with an increase in the number of momenta included in the inversion being a reduction of the statistical errors. Both methods are able to locate the maximum of ρ(ω) rather well, with the ip-method performing slightly better. Moreover, both methods struggle to reproduce the oscillatory tail of the model for ω 4.
In Figure 6, the effect of the noise level on the inversion is shown for N res = 128. Again, as for the Breit-Wigner model discussed in Section 3.1, the reconstructed spectral function becomes closer to the input ρ(ω) when is reduced. By comparing Figures 6a and 6c, and Figures 6b and  6d, it can be observed that including the Σ-matrix gives a minor advantage in the IR as the noise level increases. Additionally, as Figures 6a and 6b show, the IR behaviour of the ip-method is less oscillatory than that of the p 2 -method. Figure 7 shows how R 2 , as defined in Eq. (27), changes as the resolution N res and the noise are varied. As can be observed, both methods perform better when the Σ-matrix is included. In general, in order to reproduce a given R 2 , the ip-method allows for a larger noise value than the p 2 -method. For example, if one requires an R 2 > 0.9, a noise level of about 0.1% is needed for the ip-method, while 0.05% is required for the p 2 -method. The dark areas in Figure 7 appearing for the p 2 -method with N res = 64 and where R 2 ≈ 0 are due to the large IR mismatch. Indeed, as can be seen in Figure 5, for this particular inversion the p 2 -method returns a spectral function highly oscillatory in the IR that is far from the original function. Looking at the effect of the parameters and N res , it is clear that has a much larger effect on the quality of reconstruction than N res .
As for the Breit-Wigner spectral function, we also investigated how the inversion performs when an IR-cutoff ω 0 is introduced. In Figure 8 we report the ω 0 versus α curve for N res = 128 at various noise levels. A pattern analogous to that of the Breit-Wigner case emerges, with ω 0 = 0 being suggested as a good candidate.
The spectral function (24) was designed to obey the sum rule (25). However, in practice the numerical integral over a large momentum range tends to diverge as the reconstruction at higher momenta does not go to zero faster than 1/p 2 as expected, but stays small and finite. We found that this can remedied by imposing the sum rule as a constraint on the minimizing functional. The constraint does not change the IR, but forces the UV tail to zero as p 2 → ∞ in such a way that the sum rule is satisfied. Such a constraint fit is numerically more expensive and, therefore, we choose here to focus on the bootstrap results. The analysis of the constrained fit will be the subject of a future publication.

The spectral function for a model with a cutoff
We now consider the toy model (26) which has a physical IR cutoff ω * 0 = √ 2. We start the discussion by looking at the dependence of the inversion on ω 0 . Ideally, the ω 0 we determine during the inversion will coincide with ω * 0 .
The optimal ω 0 as a function of α for N res = 128 can be estimated from Figure 9. Indeed, at low noise, α(ω 0 ) is varying rapidly around the sharpest maximum which is located close to the physical cutoff ω * 0 . This is also evident from the curve having the biggest standard deviation there, though it might not be immediately apparent from this plot. Following our earlier discussion on minimal dependence of ω 0 on α, we will use the location of this sharpest maximum as an estimator of ω * 0 . Such a maximum appears for both methods and is more clearly present for the less noisy samples.
Looking back at the curves α(ω 0 ) for the Breit-Wigner and Bessel spectral functions in Figures 4  and 8, we see that the last local maximum in the α(ω 0 ) curve is the one with the largest variance of α w.r.t. ω 0 for all toy models studied, and therefore we infer that the last maximum in the α(ω 0 ) curve provides best guess for the values of ω 0 and α.
A comparison of the reconstructed spectral function taking α at the sharpest extremum of the curve α(ω 0 ), corresponding to an ω 0 √ 2, with setting ω 0 = 0 can be seen in Figure 10  both the ip-and p 2 -method. In all these cases Σ was included. In general, for both methods the reconstructed spectral functions are quite similar outside of the IR and, clearly, the introduction of a finite cutoff gives a ρ(ω) that is closer to the original input function. Moreover, both methods are sensitive to the maximum of the spectral function even at relatively large noise levels.

Comparing the ip vs. p 2 method: numerical and analytical insights
The analysis of the toy models suggests that, in general, the ip-method outperforms the p 2 -method. In particular, the reconstructed spectral function shows a highly oscillatory behaviour in the IR for the p 2 -method which is not observed with the ip-method; see e.g. Figures 1 and 2. Moreover, the ip-method has an overall larger R 2 as shown in e.g. Figure 3.
A first hint at this difference is given by the object M for both methods. Comparing Eq. (20) with (17) and setting ω 0 = 0, we find for the ip-method. The latter is identical to that of a Laplace transform apart from the fact that it is only valid when p i p j ≤ 0.
In fact, it can be shown that the p 2 -formalism can be directly obtained by performing G = LLρ, whereas the ip-formalism is obtained by performing G = iLFρ: This helps to understand the difference in performance of the two methods. The Laplace transform is a common example of an ill-conditioned inversion problem, so doing it twice is not likely to improve the situation. On the other hand, the Fourier transform has a well-defined inversion and therefore this operation will not negatively affect the quality of the inversion.
The analysis of the condition number, defined via the ratio of the maximal and minimal singular value of the matrix 1 + 1 α 2 M Σ −1 , also helps understanding the difference between the two methods. Recall that this is the matrix which has to be inverted in order to find the residual c.
We have computed the condition number associated with the matrix for both methods at median α for the whole ( , N res )-range for the Breit-Wigner model. The results of this analysis are shown in Figure 11 and it is apparent that the condition number of the p 2 -algorithm is consistently 2-3 orders of magnitude larger, which could explain the different behaviour in the IR region. Also, the ip-method consistently reaches higher values of R 2 at lower statistical noise levels when compared to the p 2 -method. Of the two methods considered herein, one can claim that overall the ip-method performs better and, therefore, for the analysis of the lattice data for the gluon and ghost propagators we will report only the results from this method. For the record, a p 2 -method analysis of similar gluon and/or ghost data can be found in earlier work [5,43].

The Landau gauge spectral functions from lattice data
We now proceed to compute the spectral function from the Landau gauge lattice gluon and ghost propagators at T = 0. Due to rotational invariance at T = 0, it is permissible to switch between the use of p 2 4 to p 2 as the fundamental variable in Eqs. (1, 3). Indeed, as is well-known from e.g. [20,44], the standard variable in the Källén-Lehmann representation at T = 0 is p 2 . However, due to lattice effects this rotational invariance is violated, and a significant difference between p 2 4 and p 2 appears. In order to correct for this we have followed the standard technique of applying momentum cuts, as developed in the seminal papers [45,46] to deal with the breaking of rotational invariance. After these cuts the corrected p 2 provides the better measure, and is therefore used instead of the on-axis p 2 4 . In Appendix B we provide a comparison of these momentum sets to illustrate this point in detail.
The lattice data for the gluon propagator was taken from [47]. For the ghost propagator we use the data published in [48]. The propagators were obtained from simulations on an 80 4 lattice with β = 6.0, with lattice spacing a = 0.1016 (25) fm, corresponding to a physical volume of (8.1 fm) 4 . The lattice data shown below refers to renormalized data within the MOM scheme at the scale µ = 4 GeV, i.e. the scalar form factors associated with the gluon and ghost propagators are such that Details on the sampling, gauge fixing and definitions can be found in [47,48].  Figure 12: The regularization parameter α as a function of the IR cutoff ω 0 for the reconstruction of ghost and gluon propagators. Note that the αaxis is in arbitrary units.
The lattice gluon propagator was computed with a large ensemble of 550 gauge configurations and its noise level is of the order ∼ 0.5% for the all N res = 219 momentum values. On the other hand, the ghost propagator was computed using a much smaller ensemble that included only 100 gauge configurations. However, the use of several sources considerably improves the quality of the ghost lattice data and the corresponding noise level is ∼ 1% or less for all the N res = 219 momentum values. In order to estimate the variance in the reconstructions we rely on the bootstrap method, where each bootstrap sample was inverted individually, giving an ensemble of spectral density functions which was then used to calculate the mean spectral density and its variance. In total 5500 bootstrap samples were considered for the gluons, and 700 for the ghosts. However, the reconstruction and spectral function of a bootstrap sample were only included in calculating the mean and standard deviation if the Morozov criterion was met to a precision of at least 10 −10 during the inversion. This was true for about 25% of the samples.
From the results of the toy models for the noise levels and number of data points used in the inversion, see Figures 3a and 7a, one can expect an R 2 ≈ 0.9 for the gluon inversion and an R 2 ≈ 0.8 for the ghost inversion. Moreover, one also expects a good determination of the location and height of the absolute maximum of the gluon and ghost spectral functions. We call the readers attention to the fact that our analysis does not take into account possible systematics, nor do we account for correlations between the different momenta. From the technical point of view, this last remark means that we only consider variances, not covariances. This refinement, amongst other things, will be discussed in future work.
In Figure 12 we report the curves α(ω 0 ) for the inversion of the gluon data when Σ ij = σ 2 i δ ij (no sum) and Σ ij = δ ij , i.e. with or without taking into account the statistical errors during the inversion, and for the ghost inversion with Σ ij = σ 2 i δ ij (no sum). Taking Σ ij = δ ij is not displayed because for the ghosts this does not yield any interesting new information. The pattern of the α(ω 0 ) curves for the gluon and ghost inversions look rather different. For the gluon inversion, α(ω 0 ) has several extrema that can be associated with several values of the cutoff ω 0 where ∂ω 0 /∂α 0. On the other hand, for the ghost inversion the α(ω 0 ) curve has a single maximum at ω 0 = 0 and a steep decrease towards small values of α as ω 0 departs from zero towards larger values. We take this behaviour as an indication that the right cutoff value for the ghost data is ω 0 = 0 and hence we will only display the results for the ghost inversion at this particular cutoff. In the inversions of the gluons and ghosts the diagonal covariance matrix Σ ij = σ 2 i δ ij (no sum) was always included since our toy model studies showed that this gives better reconstructions than without including a covariance matrix.
We do not attempt to check the sum rule (25) as in our formulation, by construction, the correct UV asymptotic logarithmic tails of neither propagator nor spectral function are reproduced. This is best seen from Eqs. (12)- (16) in [5], and the discussion thereafter. Roughly speaking the current Tikhonov implementation gives, for µ large, ρ(µ) ∼ 1/µ and G(p) ∼ (ln p 2 )/p 2 . As mentioned earlier, the situation could be improved by including the sum rule as a constraint. We will come back to this issue in future work.

The gluon spectral function
As can be seen in Figure 12, for the gluon data inversion the curve α(ω 0 ) shows several regions where α changes quickly as ω 0 varies slightly, which can be associated with values of the cutoff ω 0 that are stable against variation of the regularization parameter α, i.e. where ∂ω 0 /∂α 0. The precise values where the derivative vanishes are more difficult to identify. From a practical point of view, and based on the observations for the toy models, we take the location of the corresponding nearby maxima as the estimated value for the IR cutoff. However, inversions are also performed at the locations of the minima in this curve to allow for a direct comparison. Figure 13 shows the reconstructions of the gluon propagator for the ω 0 identified with the extrema of α(ω 0 ) as shown in Figure 12. As for the toy models, the p = 0 data point was not included in the inversion procedure. In general, the reconstructed propagators are in very good agreement with the lattice data for all ω 0 . However, when extrapolating towards p → 0 + , it is clear that the introduction of a cutoff ω 0 greatly improves the prediction of the point G(p 2 = 0), with the minimum ω 0 = 400 MeV and the maximum ω 0 = 425 MeV performing most reliably. The minimum at ω 0 = 220 MeV and the maximum ω 0 = 250 MeV also still perform reasonably well, but the reconstructions with ω 0 = 0 MeV, ω 0 = 59 MeV and ω 0 = 111 MeV are most certainly unreliable.
With this in mind we turn our attention to the spectral densities associated with the various extrema of α(ω 0 ), as shown in Figure 14. One striking observation is the stability of the positions of the zeroes and the extrema of ρ(ω) for all reconstructions, even if the corresponding values are not exactly identical. The effect of reconstructing at a maximum or minimum of α(ω 0 ) seems minimal on these features. It follows that one can claim a global maximum for the spectral function at ω = 0.65 GeV with "Full width at half maximum"(FWHM) = 0.27 GeV, a negative minimum at ω = 1.19 GeV with FWHM = 0.49 GeV, a positive maximum at ω = 2.11 GeV with FWHM = 0.71 GeV, etc. with zeros of ρ(ω) between the quoted ω values.
The computed spectral functions reproduce the pattern observed in the preliminary study [43], where ρ reached a maximum at momentum ∼ 0.5 − 0.6 GeV and then oscillated, approaching zero at higher ω. The herein computed absolute maximum of the spectral function is roughly consistent with the predictions of [32] that used the numerical outcome from functional renormalization group (FRG) equations compatible with the scaling scenario, i.e. a gluon propagator going to zero at zero momentum following a simple power law. The FRG spectral function was calculated using     a Bayesian inspired approach that includes a dedicated guess for a basis of functions and takes into account a priori knowledge about its asymptotic behaviour. Once more, the general pattern of the spectral function computed here is in qualitative agreement with that computed in [32], an absolute maximum followed by an oscillatory behaviour towards zero, although the quantitative details differ. It should be noted, however, that contrary to [32], our method makes no explicit assumptions about the IR of the spectral function, i.e. we make no assumptions on the structure of the propagator in the IR 4 . The price paid is that we observe the oscillations ("ringing") in the IR.
It is known that different regularization recipes for ill-posed inversion problems, of which the Källén-Lehmann spectral integral calculation is an example, can yield different results. Therefore, it is important to have several toolkits to test the soundness of the computed spectral functions. This is a common observation, even applicable when gauge invariant lattice data are inverted for e.g. meson spectral functions [33]. But because the same dominant peak is found for all cutoffs ω 0 , as well as by [32], it is fair to say that this peak is meaningful.
Given a functional form for the spectral function as in Eq. (13), one can measure its derivatives. Indeed, it can also be shown [32] that in the limit of p 4 → 0 + , ∂ p4 G(p 4 ) = −∂ p4 ρ(p 4 )/2. The inset of Figure 13b shows the derivative of the various reconstructions for small momenta. For all the reconstructions where a cutoff ω 0 200 MeV has been included, the derivatives go to zero within error, which is consistent with having ρ(ω) = 0 below the cutoff, and hence ∂ ω ρ(ω) = 0, confirming the above result. As ω increases, this simple relationship between ∂ p4 G(p 4 ) and ∂ p4 ρ(p 4 ) breaks down, so no conclusions can be drawn other than that the ω → 0 behavior is correct for cutoff values ω 0 200 MeV. For the inversions at cutoffs ω 0 < 200 MeV, the relation between the derivatives of the propagator and of the spectral function is not satisfied, suggesting that a cutoff has to be included.

Ghost propagator
As discussed previously, for the ghost propagator there is no ambiguity regarding the choice of ω 0 = 0. The ghost propagator is expected to be massless and, therefore, its spectral function  should have a Dirac delta function for zero momentum. But if this is the case, the inversion of the ghost propagator data becomes rather difficult, if not impossible, to perform. Alternatively, one can rely on the ghost dressing function given by where σ(ω) is the corresponding spectral density function. Introducing the function it follows, after integratingρ, that andĜ(p) equals G(p) up to the term −g(0)/p 2 . In order to cancel this additional term, ρ(ω) has to be given by which can be checked by plugging Eq. (37) into Eq. (3) and performing the integral. A more detailed derivation of Eq. (37) can be found in Appendix C. The significance of Eq. (37) is the following: by inverting g(p) instead,ρ can be built. This is identical to ρ, apart from the fact that an additional δ (ω) has to be present at the origin. Keeping this δ (ω) in mind, we can therefore considerρ as the spectral function after subtracting the massless free ghost state. Figure 15a presents both ρ as obtained from a direct inversion of G(p), andρ as defined previously. Both reconstructions display the same minimum at p ≈ 200 MeV, butρ starts to be dominated by the 1/ω 2 behaviour for smaller momenta whereas ρ shows a maximum at p ≈ 70 MeV before going down to zero. The vertical black line in 15a indicates the smallest momentum value in the dataset, p min = 0.15 GeV. The maximum of ρ is positioned below p min , indicating that this peak could possibly be due to the attempt to reconstruct the Dirac-δ peak at p = 0. The error inρ explodes below p min , despite the fact that the error in σ(ω) stays reasonable, as can be seen from Fig. 16. This is a direct consequence of the 1/ω 2 behaviour ofρ, however the typical solution does look like the average shown in Fig. 15a. To test our approach, it has been verified thatĜ(p) + g(0)/p 2 is also a reconstruction of the original lattice data for the propagator, which is indistinguishable from a direct reconstruction of G(p). This can be seen in Figure 15b, which shows the reconstructed propagator, and clearly demonstrates thatĜ(p) + g(0)/p 2 overlays the lattice data over the full p-range where data points have been provided. To reconstruct the ghost propagator in this manner,ρ was integrated according to Eq. (36) to yieldĜ(p), and g(0) was calculated using Eq. (34). The good agreement found is an indication that the infrared ghost propagator is given essentially by its tree level value.

Conclusion
In the current paper, we improved the Tikhonov reconstruction using the Morozov discrepancy principle when applied to Källén-Lehmann inverse lattice spectroscopy, as set out previously in [5]. Based on dedicated toy models, the current research defined statistical measures for the quality of the inversion, providing a region of validity for this method. While doing so, we also supplemented the Tikhonov functional with a relative error weighting of the data.
We considered two analytically equivalent versions of the Källén-Lehmann spectral representation. However, despite this analytical equivalence, the numerical performance of the two formalisms is quite different. This difference is most notable in the IR, where the ip-formalism yields significant improvement over the previous p 2 -formalism [5]. The improvement was demonstrated most notably by a reduction in the condition number of the to-be-inverted matrix 1 + 1 α 2 M Σ −1 by 2-3 orders of magnitude. Both methods are sensitive to an IR-cutoff below which the spectral density vanishes, if such a cutoff is present.
By applying the ip-methodology to lattice SU(3) gluon data, it was found that the gluon spectral function seems to have an IR cutoff of a few hundred MeV. We also characterized a dominant peak, the location of which is fairly consistent with the findings reported in other studies like [5,32,13]. Because IR oscillations are always present in the reconstructions, we have to be careful with their interpretation. However, the dominant peak appears to be a stable prediction. For ghost data it was found that no IR-cutoff is present, consistent with the massless pole present in the ghost propagator. However, the quality of the reconstruction was greatly improved by first removing the δ-peak in the spectral density. We stress once more that our results should always be considered under the assumption that the gluon and ghost degrees of freedom have a Källén-Lehmann spectral representation to begin with. The existence of such a representation is not self-evident, since these particles are confined at zero temperature, and therefore they do not belong to the complete set of (positive norm) physical states that are usually employed to derive the spectral representation [19,20,21]. Other analytical continuations of the data are thus in principle possible. We recall that the Källén-Lehmann spectral integral allows only for branch cuts along the negative real axis, while certain analytical approaches, or fits to lattice data based thereon, entail the presence of e.g. complex conjugate poles, in se incompatible with the Källén-Lehmann spectral structure [32,49,50,51,47,52,53]. As a future improvement of our method, such complex conjugate poles could therefore be included in the inversion method. The inclusion of such poles could also lead to a reduction in the IR ringing, which would be an indication of their presence. This being said, it should also be noted that the ringing effect is also present in the presented toy model inversions, which definitely contain no complex conjugate poles' contributions. Moreover, the IR ringing also plagues other inversion strategies, like MEM [33,24], again without such poles. In recent work, [54], yet another inversion method was proposed, based on rational function interpolation. More evidence was presented for a single set of complex conjugate poles, so it would be interesting to test whether this feature prevails also within our methodology. This is currently under investigation. Additionally, imposing adherence to the sum rule as a constraint seems compatible with the ip method and will be investigated further. Such a constraint could be extended to the generalized sum rule that includes the set of complex conjugate poles, see also the recent paper [55].
Moreover, in this next phase of research we will further put our inversion strategy to the test by applying it to finite temperature lattice data for gluons, ghosts and quarks. An important question to be further addressed there is whether, and to what extent if so, the spectral functions are sensitive to the deconfinement (chiral) transition or which kind of quasi-particle behaviour can be identified [29,31,14]. This will be discussed in forthcoming work.

B Gluon propagator rotational invariance breaking
In Figure 17 we illustrate the effects associated with the breaking of rotational invariance in lattice simulations by comparing the gluon propagator lattice data for different types of momenta. In order to do so Figure 17 displays the gluon propagator G(p 2 ) at all momenta, the on-axis momenta, and the momentum cuts which are close to p µ = (1, 1, 1, 1). These momentum cuts [45,46] were devised to suppress rotational breaking effects on the propagator lattice data, and as an additional advantage they also provide access to a larger range of momenta when compared to the on-axis momenta. This is clear from Figure 17, where the momentum cuts can be seen to represent the full range of available momenta. It can also be seen from the figure that although the cuts and on-axis momenta agree in the IR region, they start to diverge slightly in the UV. Since the UV behaviour of the spectral function can already be accessed perturbatively, the most interesting physics of the spectral function for our current research is given by the IR region. Using the momentum cuts therefore gives identical information to the on-axis momenta on the IR region while providing a better connection with the UV.

C From dressing function to propagator
This Appendix will detail the relationship between the spectral density functions of the dressing function and the propagator. The dressing function is defined as g(p) = p 2 G(p).
Since ghost particles are massless, the spectral density function is expected to contain a δ function at zero momentum. Because such a δ peak complicates the numerical inversion, it is beneficial to remove it first by inverting the dressing function instead. Let where both ρ and σ are odd functions. We then definê Note that the appearance of ω 2 is motivated by the fact that both ρ and σ have to be odd functions.
Integrating overρ gives where the first term on the second line could be dropped due to the oddness of the integrand. We see thatρ nicely produces G(p) but also an extra term. This extra term can be cancelled by identifying ρ as ρ(ω) = −g(0)δ (ω) +ρ(ω).
We are forced to use the partial derivative operator δ (ω) instead of δ(ω) due to the demand that ρ(ω) should be an odd function. Integrating over Eq. (48), we find that his indeed gives G(p). Figure 18 shows this process numerically. As can be seen from the figure,Ĝ(p) can hardly be called a reconstruction of the data. However, upon adding g(0)/p 2 the result describes the data equally well as the direct reconstruction G re (p).