On the extraction of spectral densities from lattice correlators

Hadronic spectral densities are important quantities whose non-perturbative knowledge allows for calculating phenomenologically relevant observables, such as inclusive hadronic cross-sections and non-leptonic decay-rates. The extraction of spectral densities from lattice correlators is a notoriously difficult problem because lattice simulations are performed in Euclidean time and lattice data are unavoidably affected by statistical and systematic uncertainties. In this paper we present a new method for extracting hadronic spectral densities from lattice correlators. The method allows for choosing a smearing function at the beginning of the procedure and it provides results for the spectral densities smeared with this function together with reliable estimates of the associated uncertainties. The same smearing function can be used in the analysis of correlators obtained on different volumes, such that the infinite volume limit can be studied in a consistent way. While the method is described by using the language of lattice simulations, in reality it is completely general and can profitably be used to cope with inverse problems arising in different fields of research.


I. INTRODUCTION
Hadronic spectral densities are crucial ingredients in the calculation of physical observables associated with the continuum spectrum of the QCD Hamiltonian. A notable classical example is provided by the differential cross section for the process e + e − → hadrons that, at leading order in the electromagnetic coupling, is proportional to the QCD spectral density evaluated between hadronic electromagnetic currents, where E is the energy of the electron-positron pair in the center of mass frame, H and P are respectively the QCD Hamiltonian and total momentum operators and J µ em (x) is the hadronic electromagnetic current. Other important examples of observables, in which spectral densities play a crucial rôle, are the flavour-changing non-leptonic decay-rates of kaons and heavy flavoured mesons, the deep inelastic scattering cross-section, and thermodynamic observables arising in the study of QCD at finitetemperature and of the quark-gluon plasma.
It is notoriously difficult to obtain model-independent non-perturbative theoretical predictions for hadronic spectral densities. In principle this is a problem that can be addressed from first-principles within the solid framework of lattice QCD. However, in practice, one has to face highly non-trivial numerical and theoretical problems in order to extract spectral densities from lattice simulations.
The origin of these problems can be traced back to the fact that lattice results unavoidably are affected by statistical and systematic errors. More precisely, the primary observables computed in a lattice simulation are Euclidean time-ordered correlators at discrete values of the space-time coordinates and on a finite volume, e.g.
where L is the linear extent of the spatial volume V = L 3 while O andŌ are generic hadronic operators. In the following we shall not discuss cutoff effects and, therefore, we shall not indicate the dependence of the different quantities upon the lattice spacing. We shall however always assume that the correlators are known only for discrete values of the space-time coordinates. At positive Euclidean times t ≥ 0 the previous correlator can be rewritten as where, for simplicity (see below for a generalization in the case of periodic boundary conditions in time), we have assumed that the time extent of the lattice is infinite and where we have defined the associated spectral density The main problems faced during the extraction of spectral densities from lattice simulations can now be explained by starting from the previous two expressions.
The first problem is associated with the fact that the extraction of ρ L (E) from the measured lattice correlator C(t) requires an inverse Laplace-transform to be performed numerically, an ill-posed problem when the measured data are affected by uncertainties. This is the case for C(t) that unavoidably will be affected by statistical errors, particularly at large time separations where (a part from some countable exceptions) the signal-to-noise ratio degrades exponentially for Euclidean hadronic correlators.
The second problem comes from the fact that the finite volume Hamiltonian H L has a discrete spectrum. The finite volume spectral density is a distribution having support in correspondence with the eigenvalues E n (L) of H L , ρ L (E) = n w n (L) δ(E − E n (L)) .
Here it is important to notice that the finite volume spectral densities cannot be directly associated, even in the ideal case in which these can be computed exactly, to physical observables.
A big step forward in the extraction of spectral densities from lattice simulations has recently been done in Ref. [1] where it has been suggested to rely on a method originally devised to analyse geophysical observations in Ref. [2]. The central idea of the so-called Backus-Gilbert method (see also Ref. [3]) is to optimize the amount of information that can be extracted from noisy measurements, in our case C(t), by focusing on the calculation of smeared spectral densities, Notice that smeared spectral densities are smooth functions of the energy (as opposed to the distributions ρ L (E)) and that the study of their infinite-volume limit at fixed smearing function is a well-posed problem. Ideally one would like to choose the smearing functions ∆ σ (E , E) with support in a region around E of width proportional to σ and such that they become Dirac δfunctions in the limit in which the smearing radius parameter σ is sent to zero. If one can choose the same smearing function on different volumes, the infinitevolume spectral density can then be extracted by taking the double limit in the specified order.
As already done by the authors of Ref.
[1], we stress that the limit of vanishing smearing radius might not be necessary in order to compare theoretical predictions with experimental data. Indeed, it might be possible to smear experimental observations with the same function used in the theoretical calculation. In the case of e + e − → hadrons one should for example compare the theoretical predictions with ∞ 0 dΣ(E) ∆ σ (E , E). In fact one should also notice that, in order to derive results such as Eq. (1), a smearing function has to be introduced at intermediate steps of the calculations and the limit of vanishing smearing radius has then to be taken. This point is extensively discussed in Ref. [1] by making contact with the so-called Fermi's golden rule.
On the one hand, the Backus-Gilbert method is an extremely efficient algorithm 1 for controlling the statistical errors on the estimated smeared spectral densities. On the other hand, the shape of the smearing function cannot be chosen arbitrarily in this method, it is an output of the procedure. Indeed, the algorithm is designed in such a way that the width of the smearing function (having the properties of being peaked around E and of having unit area) is optimized on the basis of the number of observations and of their statistical uncertainties (in our case the number of discrete times t at which C(t) is known and the associated statistical errors). This feature of the algorithm may not represent a problem in experimental applications of the Backus-Gilbert method because, at the end of the procedure, the resulting smearing function is known and no infinite-volume limit has to be taken. It is instead a strong limitation in the context of lattice simulations where simulations at different volumes produce results with different statistical uncertainties and different number of points. In this case one gets different smearing functions at different volumes and the extraction of the infinite-volume physical observable becomes, if not impossible, extremely difficult.
In this paper we present a new method in which the shape of the smearing function is an input of the procedure and not an output. The method uses the same mechanism of the original Backus-Gilbert proposal to keep the statistical errors under control. This happens at the price of a distortion of the target smearing function induced by the presence of statistical errors and by the finite number of observations. At the end of the numerical procedure the systematic error associated with this distortion can be reliably quantified and added to the statistical uncertainties in order to provide a reliable estimate of the smeared spectral densities. Our method gives an exact reconstruction of the spectral densities, smeared with the chosen functions, in the limit of vanishing statistical uncertainties and of an infinite number of discrete lattice points along the time direction. The method is quite general and we are pretty confident that it will be useful in addressing "inverse problems" arising in different fields of research, for example the problems where the application of the Backus-Gilbert method has already proven to be useful.
The rest of the paper is organized as follows. In section II we review of the original Backus-Gilbert method and in section III we present our method. In section IV we apply our method to a benchmark system where we know the exact spectral density, while in section V we apply the method to real lattice correlators. We draw our conclusions in section VI. The paper ends with the appendix A containing the explicit formulae needed to implement our method in computer programs.

II. REVIEW OF THE BACKUS-GILBERT METHOD
In this section we review the original Backus-Gilbert method. This will help us to set the notation and to discuss, in the following sections, the similarities and the differences between our new proposal and the Backus-Gilbert approach. Although the method is general and can be applied to many different problems, in our discussion we shall use the language of lattice correlators to explain the approach. The generalization to other contexts is straightforward, see Ref. [3].
Let us consider a generic lattice correlator that, by following the same steps of the introduction, can be written as where ρ L (E) is the distribution corresponding to the finite volume spectral density. Here b T (t, E) are the socalled basis functions. These are simply given by exponentials in the limit of an infinite temporal extent of the lattice, while, in the case of periodic boundary conditions in time and by assuming that the correlator is symmetric under time-reversal, the basis functions are 2 In all what follows the time variable t is assumed to be discrete, non negative and smaller than the time extent of the lattice (0 ≤ t < T ).
The central idea of the Backus-Gilbert method is to search for a smearing function that lives in the space spanned by the basis functions, more precisely 2 The use of Eq. (10) in the case of periodic boundary conditions in time is an approximation. The approximation is much better than the naive use of Eq. (9) at finite values of T but the general spectral decomposition of a periodic hadronic correlator would require the inclusion of other contributions that vanish exponentially fast when T is sent to infinity. See Ref. [7] for a discussion of this point.
with t max < T /2. Once the coefficients g t (E ) that define the smearing function are known, the smeared spectral density can then easily be computed by starting from the correlator, The Backus-Gilbert procedure "optimizes" the choice of the smearing function, i.e. of the coefficients g t (E ), on the basis of the measured data for the correlator. In absence of statistical errors the coefficients are fixed by minimizing a deterministic functional that can be interpreted as a measure of the width of the smearing function. The functional is and the minimization is performed under the unit area constraint It is a simple exercise (see Ref. [3]) to show that the solution of this problem is given by where we have used a vector notation for the coefficients, g T (E ) = (g 0 (E ), · · · , g tmax (E )), the vector R has the following entries and the elements of the matrix A(E ) are given by It is important to notice that the matrices A(E ) tend to be nearly singular for the basis functions discussed in this paper. From a numerical point of view this might be an issue, but in fact the problem can easily be circumvented on currently available computers by using extended-precision arithmetic. In order to avoid coping with algorithmic instabilities induced by numerical rounding errors, this is what we have done in our computer programs [8].
In Figure 1 we show some examples of the smearing functions obtained by using Eq. (15). As it can be seen by comparing the plots in the different panels, the function ∆ BG (E , E) becomes more similar to a Dirac δ-function  when increasing the number of basis functions used in its definition.
In Figure 2 we show the coefficients g t (E ) corresponding to the smearing function ∆ BG (E , E) shown in the right-panel of Figure 1. The plot has been shown in order to highlight a typical feature exhibited by these coefficients. As a consequence of the nearly singular nature of the matrix A(E ) the coefficients become gigantic for some values of t and, moreover, oscillate in sign. Having noticed this feature we can now discuss the Backus-Gilbert procedure in presence of uncertainties.
The exact correlator is an idealization that is not acces-sible in the real world and, in presence of (experimental or) statistical errors, we have to consider where the index i runs over the N different statistical samples (for a lattice correlator we can think to the different bootstrap or jackknife bins),C(t) is the statistical average and δC i (t) is the deviation from the average of the i-th bin, Given the fact that the coefficients g t (E ) are huge numbers, even a tiny deviation (for example an apparently harmless rounding error) from the average gives an unacceptably large contribution to the smeared spectral function. Indeed, by applying Eq. (12) to C i (t), one gets that the sums t g t (E ) δC i (t) are also huge numbers in general and the final error on the estimated smeared spectral functions turns out to be unacceptably large. This can be viewed as a manifestation of the fact that we are dealing with a numerically ill-posed problem.
In order to keep statistical errors under control Backus and Gilbert considered another functional of the coefficients, a measure of the statistical error on the smeared spectral function, namely where Cov is the covariance matrix of the correlator, In the presence of statistical errors, the coefficients are fixed by minimizing the following functional again under the unit area constraint of Eq. (14). In this case the solution is given by where the matrix W(λ, E ) has elements with A tr (E ) already defined in Eq. (17). The real number λ is a free parameter of the algorithm, chosen in the range [0, 1].
The functional W [λ, g] is in fact a convex linear combination of the deterministic functional A BG [g] and of the error functional B[g]. The presence of the error functional in the minimization procedure forbids solutions corresponding to gigantic values of the coefficients. Statistical errors are thus kept under control at the prize of accepting that the shape of the smearing function is determined (somehow optimized) also by the statistical errors.
The tuning of the parameter λ is a subtle issue in the Backus-Gilbert procedure. Choosing λ too small 3 may result in too large statistical errors while values of λ close to one may generate smearing functions that are useless for physical applications. This point will be discussed in more details in section IV where, by using a synthetic correlator generated by starting from an exactly known spectral function and by adding random statistical noise, we shall compare the results obtained with the Backus-Gilbert method with the ones obtained by using our method.

III. THE NEW METHOD
In our method the target smearing function is an input of the algorithm. For example, the target smearing function can be chosen as a Gaussian of width σ, centred at E and normalized to have unit area in the interval The method then searches for an optimal approximation of the target smearing function in the space spanned by the basis functions, where t max < T /2. The previous formula is identical to the definition of the smearing function in the original Backus-Gilbert proposal, see Eq. (11) above. The difference is in the way the coefficients g t (λ, E ) are determined.
This is done by minimizing again a convex linear combination of a deterministic functional and of the error functional, under the unit area constraint However, in our case, the deterministic functional is chosen to be a measure of the difference between the target and the approximated smearing functions, namely 3 Notice that our parameter λ corresponds to 1 − λ in Ref. [1].
while the error functional is conveniently normalized with C(0) 2 , i.e. the square of the correlator at t = 0.
The parameter E 0 has to be chosen in such a way that the finite volume spectral function ρ L (E 0 ) = 0. This is always possible in the case of connected correlators in QCD and in the charged sectors of QCD+QED because of the presence of a mass gap. According to our experience, having E 0 > 0 is particularly convenient in the case of b T (t, E) as basis functions.
It is easy to show that the solution of the minimization procedure is given by where the vector R has already been defined in Eq. (16) and the components of the vector f (E ) are given by The matrix W has the elements Explicit expressions for R, W and f , derived in the case of the smearing function of Eq. (24), are given in appendix A.
In absence of statistical errors our procedure is a method to obtain the best approximation of the target smearing function in the space spanned by the basis functions under the norm defined by the functional A[g]. Since the target function is assumed to be analytic in the interval [E 0 , ∞] and to decay faster than any power for E that goes to infinity, the error of the approximation can be made arbitrarily small by enlarging the space spanned by the basis functions. This can be understood by looking at the deterministic functional in the case of the choice of b ∞ (t, E) as basis functions after performing the change of integration variable to x = e −E , In fact in our procedure we are just searching for the best polynomial approximation of a well-behaved function. The argument holds also in the case of b T (t, E) as basis functions because these simply reduce to b ∞ (t, E) in the limit of infinitely many times 4 . The comparison of the smearing functions∆ σ (E , E) obtained with our method in absence of statistical errors with the target function ∆ σ (E , E) is shown in Figure 3.
The different philosophy of our method with respect to the original Backus-Gilbert proposal can already be appreciated by comparing Figure 3 and Figure 1. In the Backus-Gilbert method, by changing t max one gets a different (sharper) function. In our method, by increasing t max one gets a better approximation of the target smearing function. Moreover, in our method the error of the approximation of the target smearing function is known and this information can be used to estimate the final error on the smeared spectral density as we are now going to explain.
tmax to infinity one has also to send T to infinity. On the one hand, in the presence of statistical uncertainties the difference between the target and the approximated smearing functions is due to both a finite value of t max and to the presence of the error functional B[g] in the minimization procedure, whose importance is regulated by the choice of the λ parameter. On the other hand, the quantity can always be calculated at the end of the procedure. In principle, by knowing this quantity one can write an exact expression for the bias on the smeared spectral density associated with our method, In practice we cannot use the previous formula for the obvious reasons that we don't know the true spectral density and that we cannot explore the full energy range [0, ∞]. The quality of the results obtained with our method can nevertheless be assessed by monitoring the relative deviation δ σ (E , E). This will be illustrated in the next section where, by using a benchmark system where the exact spectral density is known, we will show that a trustable estimate of the systematic error associated with our method can be obtained by using the formula The relative size of the statistical and systematic errors can be regulated by changing the parameter λ and, if the estimate of the systematic uncertainty is reliable, results corresponding to different values of λ have to be compatible within the total uncertainties. In our method the choice of the trade-off parameter can be optimized by using i.e. the function of λ obtained by evaluating the functional W [λ, g] at the solution g t (λ, E ) of the minimization procedure. This function has a characteristic shape that we show in Figure 4. At very small values λ the contribution to W coming from the error functional λB/C(0) 2 is very small for generic values of the coefficients and the minimization procedure acts on the deterministic functional (1 − λ)A in order to obtain the best approximation of the target smearing function. Conversely, at very small values of (1 − λ) the contribution of the error functional λB/C(0) 2 is dominant and the minimization procedure acts to reduce the statistical errors at the price of distorting the smearing function. The interplay between these two regimes generates a maximum in W (λ, E ), at the value λ where the deterministic and error functional are balanced. Therefore our method automatically suggests the optimal 5 choice for the trade-off parameter, i.e. λ = λ .
In all our numerical experiments we have checked that the results corresponding to values of λ smaller than λ are compatible within the corresponding total uncertainties. Indeed, statistical errors increase by decreasing λ while the relative deviation δ σ (E , E) gets smaller and smaller and, in this region, Eq. (35) can safely be used to get a reliable estimate of the systematic error. For values of λ much larger than λ , the results are instead affected by unacceptably large systematic uncertainties.

IV. THE BENCHMARK SYSTEM
In order to test our method and to compare it with the original Backus-Gilbert proposal we consider in this section a benchmark system where we know the exact spectral density. This information is used to build an exact synthetic correlator which we can then manually distort by adding random noise. We decided to consider the same benchmark system used in Ref. [1] where additional details on the model can be be found.
The benchmark system is a toy model of three scalar particles, the pion π, the kaon K and the φ meson with 5 We cannot exclude the presence of more than one maximum in W (λ, E ). We never encountered this situation in our numerical experiments but, if this happens, we suggest to define λ as the smallest value of λ where W (λ, E ) has a maximum. physical masses such that The particles are subject to the interaction Lagrangian density and the interactions are assumed to be perturbative. The authors of Ref. [1] have considered a correlator in this theory having as finite volume spectral density the following expression, where the momenta are the ones allowed by periodic boundary conditions in space, i.e. p = 2πN 3 /L, and where the energies are The infinite-volume spectral density is given by where The previous result agrees 6 with the one originally given in Ref. [1].
In our numerical experiments we have set the parameters of the model to the same values used in Ref. [1], i.e. we have set m π = 0.066, m K /m π = 3.55, m φ /m π = 7.30, g K = 1 and g π = 10 √ 8. Since we are working in lattice units the previous numbers have to be read under the formal assumption that a = 1. In evaluating the finite volume spectral density we have used a cutoff in the energy by replacing ρ L (E) with ρ L (E)θ(Λ − E) using the value Λ ≈ 19m π . In both plots the solid black curve correspond to the exact smeared spectral density that our method (not the Backus-Gilbert one) is expected to reproduce in the infinite tmax limit.

A. Exact data
In Figure 5 we compare the results obtained with our method (blue points) in absence of statistical errors with the results obtained by using the Backus-Gilbert method (orange points). The two plots in the figure have been obtained by starting from the correlator on the volumes L = 24 (first row) and L = 32 (second row) with t max = 30. In the case of our method the results have been obtained by setting E 0 = 0 and by using the target smearing function in Eq. (24) with σ = 0.1. In both plots the solid black curve corresponds to the exact smeared spectral density, which our method is expected to reproduce in the infinite t max limit. As it can be seen, in both plots the agreement between the numerical results obtained with our method and the exact result is excellent. Notice that the results obtained with the Backus-Gilbert method are not expected to reproduce the black line. In fact, because the smearing function is an output of the procedure, it can only be controlled by changing t max and, moreover, it is different at different values of E . We can only notice that our choice of setting σ = 0.1 is similar to the choice made by the Backus-Gilbert method on the volume L = 32 at high energies where the smeared spectral density is more smooth and starts to have an infinitevolume like behaviour.

B. Noisy data
In order to test our method in the case of a noisy correlator, we add uncorrelated random noise to the synthetic correlator in such a way that the signal-to-noise ratio is constant. For the results presented here, on average, the relative standard deviation of the correlator is chosen to be around 2%. For all the numerical examples shown in the rest of the paper, we only use the diagonal part of the covariance matrix in the minimization procedure.
In the reconstruction of the spectral density we estimate the statistical and systematic uncertainty independently and combine them in quadrature. For estimating the statistical uncertainty we use a bootstrapping procedure, i.e. we apply our method to a set of bootstrap samples, from which we derive the mean and standard deviation of the reconstructed spectral density. For the systematic uncertainty we use Eq. (35) such that the total uncertainty is given by where the factor 0.68 is introduced to give a consistent 1σ uncertainty on the final result.
In Figure 6 we show three plots in order to discuss the impact of the choice of the trade-off parameter λ and the significance of our estimate of the systematic error. In all the plots the black curve is the exact smeared spectral density. The data in the top-panel have been obtained at the value of λ determined from the maximum of W (λ, E ) at E /m π = 7 and, after having checked that in this case λ is not strongly dependent upon E , we have used the same value of the trade-off parameter at all energies. On the first plot, the orange band shows the statistical error, while in all plots, the blue band corresponds to the total uncertainty. As it can be seen, the The results shown in Figure 6 correspond to a relatively challenging situation because on small volumes and/or at small values of σ the smeared spectral function exhibits an oscillating behaviour induced by the fact that the energy levels are largely spaced. In cases where the smeared spectral density is very smooth, either because the smearing parameter σ is larger, or because the volume is larger, the quality of the reconstruction is much better. This can be seen in Figure 7 where we show results on the volume L = 24 corresponding to a larger value of σ with respect to the one used in Figure 6 and results on the volume L = 48 for two different values of t max . In this cases the use of Eq. (35) to quantify the systematic uncertainty results in an over-estimate of the error. This feature makes us pretty confident about the reliability of the results obtained with our method.
We close this section by illustrating the approach to the infinite-volume limit of the reconstructed smeared spectral densities in the case of this benchmark model. The plots in Figure 8 show the results obtained on different volumes by setting the smearing radius parameter to σ = 0.1 and by using b T (t, E) as basis functions with T = 2(t max + 1) and t max = 31. More precisely, the plot in the first panel (starting from the top) corresponds to L = 24, the one in the second panel to L = 30, the one shows the calculation of the pion mass, the lightest state contributing to the spectral density in this case, extracted from a standard effective-mass analysis. The bottom-panel shows the reconstructed smeared spectral density obtained by applying our method with σ = 0.1, by using bT (t, E) as basis functions with T = 48 = (2tmax + 1), by setting E0 = 0.37mπ and by using the value of λ determined at E = 3.7mπ for all the energies explored. As expected, the smeared spectral density shows a peak in correspondence of E /mπ 1 and another structure around E /mπ 3. exact smeared infinite-volume spectral density, namelŷ As it can be seen, ρ(E) is a continuous function of the energy for E ≥ 3m π but it has a cusp at E = 2m K , i.e. in correspondence of the two-kaons threshold. In the infinite-volume limit the numerical data are expected to reproduce the smeared spectral density that is instead a smooth function. This already happens, within the errors, at L = 36 for medium-high values of the energy and for all the explored energies at L = 48. Remarkably, at L = 48 the numerical data agree with the infinite volume curve up to energies of order E /m π 11 at the level of the statistical uncertainties (orange band) thus confirming that on large volumes Eq. (35) gives a conservative estimate of the systematic uncertainties.

V. LATTICE CORRELATORS
In this section, in order to show the quality of the results that can be obtained by applying our method to true data, we discuss two analyses of simulated lattice correlators from which we extract the associated smeared spectral densities.
In the first example we consider a meson pseudoscalarpseudoscalar correlator obtained by performing a lattice simulation of QCD with three degenerate flavours on a lattice volume L 3 × T = 24 3 × 48 with periodic boundary conditions in time and C * boundary conditions [9] along the spatial directions. The bare parameters of the simulation correspond to the CLS ensemble H101 and can be found in Table 1 of Ref. [10]. More precisely, the correlator is given by where u and d are the up and down quark fields that, in this unphysical simulation, have the same mass. The lightest states contributing to the finite volume spectral density associated with the correlator C QCD (t) are expected to be the pion and the three-pions states with vanishing total momentum allowed by the boundary conditions. This means that we expect the leading contributions to ρ L (E) to be proportional to δ(E − m π ) and to δ(E − E 3π ), where m π is the mass of the pion and E 3π 3m π . In the top-panel of Figure 9 we show the numerical determination of m π from a standard effectivemass analysis. In the bottom panel of the same figure we show the smeared spectral density obtained by applying our method with σ = 0.1 and by using the value of λ determined at E = 0. 5 3.7m π (see Figure 4) for all the energies explored. As it can be seen, the reconstructed smeared spectral density clearly shows a peak centred around E /m π 1 and another structure around E /m π 3.
In the second example we consider again a meson pseudoscalar-pseudoscalar correlator, but in this case it has been obtained from a QCD+QED simulation performed at the unphysical value α em = 0.05 of the electromagnetic coupling constant with dynamical up, down and strange quarks. The masses of the down and the strange quarks, having the same negative electric charge, have been set equal in this simulation and different from the mass of the positively charged up quark. The simulation has been performed on a lattice volume L 3 × T = 24 3 × 48 with periodic boundary conditions in time and C * boundary conditions [9] along the spatial directions. More details about the simulation can be found in Ref. [11]. The plot in the top-panel of Figure 10 shows the effective mass extracted from the correlator where S(x) and U (x) are the gauge-invariant interpolating operators for the strange and up quarks described in The bottompanel shows the reconstructed smeared spectral density obtained by applying our method with σ = 0.1, by using bT (t, E) as basis functions with T = 48 = (2tmax + 1), by setting E0 = 0.15m K + and by using the value of λ determined at E = 1.5m K + for all the energies explored. As expected, the smeared spectral density shows an isolated peak in correspondence of E /m K + 1 and another structure that starts in proximity of E /m K + 2.4. details in Refs. [9,11]. In this case the effective mass analysis measures m K + , the mass of the charged kaon, and this is the lightest state contributing to the finite volume spectral density. Increasing the energy, we expect contributions to the spectral density coming from states corresponding to the charged kaon plus photons and from states with three kaons. Since the volume is rather small in this case and the boundary conditions do not allow for the propagation of photons with zero momenta, after the charged kaon peak we expect a contribution to ρ L (E) proportional to δ(E − E 3K ) with E 3K m K + + 2m K 0 . By using the value of m K 0 measured in Ref. [11] we have E 3K /m K + 2.6. This expectation is confirmed by the plot in the bottom-panel of Figure 10 where the smeared spectral density shows an isolated peak in correspondence of E /m K + 1 and a structure that starts in proximity of E /m K + 2.4.

VI. CONCLUSIONS
We have presented a new method for addressing inverse problems in the presence of noisy observations. The method can be used to extract smeared spectral densities from measured correlation functions and it provides results associated with a reliable estimate of both the statistical and the systematic uncertainties.
The function used for smearing the spectral density is an input of our method, and for this reason, it can be held fixed in the analysis of data corresponding to different correlators. This feature is particularly convenient in lattice applications because it allows to study the infinite volume limit of the reconstructed smeared spectral densities in a systematic way.
The mechanism used in our method to keep statistical errors under control has been inherited from the classical method of Backus and Gilbert. However, contrary to what happens in the original Backus-Gilbert proposal, our method has a natural built-in mechanism to optimize the choice of the so-called trade-off parameter and, moreover, the significance of the estimate of the errors can be assessed by checking the compatibility of the results obtained at sub-optimal values of this parameter.
In order to illustrate the quality of the results that can be obtained with our method, we have applied it to a benchmark system where we know the exact spectral density, both on finite volumes and in the infinite volume limit. We have shown that the results obtained with our approach reproduce within the errors the expected finite volume smeared spectral densities and also that, by increasing the volume, they approach the expected infinite volume limit.
We have also applied the method to true data in the case of correlators obtained from QCD and QCD+QED lattice simulations. Using these examples we have shown that smeared spectral densities can be extracted with satisfactory accuracy and that the numerical results are compatible with the expectations coming from the knowledge of the spectrum of the two theories.
We have discussed the method by using the language of lattice correlators but, given its generality, we are pretty confident that it will be useful to address inverse problems arising in other fields of research, particularly those where the classical Backus-Gilbert method has already proven useful.