Bayesian analysis of quark spectral properties from the Dyson-Schwinger equation

TL;DR

This paper investigates the quark spectral function at finite temperature using Dyson-Schwinger equations, incorporating unquenched fermion flavors and a Bayesian reconstruction method to analyze spectral features and temperature effects.

Contribution

It introduces a Bayesian approach with a new prior for spectral reconstruction in unquenched QCD, revealing detailed quark spectral structures and their temperature dependence.

Findings

Identification of a three-peak structure in the quark spectral function.

Observation of modifications in spectral features around the pseudo-critical temperature.

Differences in temperature dependence compared to previous studies with gluon propagator models.

Abstract

We report results on the quark spectral function in the Landau gauge at finite temperature determined from its Dyson-Schwinger equation. Compared to earlier quenched results this study encompasses unquenched $N_f=2+1$ fermion flavors in the medium. For the reconstruction of real-time spectra we deploy a recent Bayesian approach (BR method) and develop a new prior in order to better assess the inherent systematic uncertainties. We identify the quark quasi-particle spectrum and analyze the (non-)appearance of zero modes at or around the pseudo-critical temperature. In both, the fully unquenched system and a simpler truncation using a model for the gluon propagator we observe a characteristic three-peak structure at zero three-momentum. The temperature dependence of these structures in case of the gluon propagator model is different than observed in previous studies. For the back-coupled and unquenched case we find interesting modifications at and around the pseudo-critical transition temperature.

Figures

Equations

$$S^{-1}(i{\omega }_p,\text{boldmathpunboldmath})$$

$$S_0^{-1}(i{\omega }_p,\text{boldmathpunboldmath})$$

$$S^{-1}(i{\omega }_p,\text{boldmathpunboldmath})$$

$$+B(i{\omega }_p,|\text{boldmathpunboldmath}|).$$

$$\Sigma (i{\omega }_p,\text{boldmathpunboldmath})$$

$$\displaystyle\times\Gamma_{\nu}(i\omega_{q},\text{\boldmath$q$\unboldmath},i\omega_{p},\text{\boldmath$p$\unboldmath})D_{\mu\nu}(i\omega_{p}-i\omega_{q},\text{\boldmath$p$\unboldmath}-\text{\boldmath$q$\unboldmath})\Big{]}\,,$$

$$g^2{D}_{\mu \nu }(i{\omega }_k,\text{boldmathkunboldmath}){\Gamma }_{\nu }(i{\omega }_q,\text{boldmathqunboldmath},i{\omega }_p,\text{boldmathpunboldmath})$$

$$=[P_{\mu \nu }^T{D}_T({\omega }_k^2,{\text{boldmathkunboldmath}}^2)+P_{\mu \nu }^L{D}_L({\omega }_k^2,{\text{boldmathkunboldmath}}^2)]{\gamma }_{\nu },$$

$$D_T$$

$$\mathcal{D}(s)$$

$$S(i{\omega }_p,\text{boldmathpunboldmath})$$

$$\rho (\omega ,\text{boldmathpunboldmath})$$

$$\displaystyle\phantom{=2\pi\big{(}}-\rho_{\rm s}(\omega,|\text{\boldmath$p$\unboldmath}|)\Big{)}\,.$$

$${\rho }_4(\omega ,|\text{boldmathpunboldmath}|)$$

$$Z_2{\int }_{-\infty }^{\infty }d\omega {\rho }_4(\omega ,|\text{boldmathpunboldmath}|)=1,$$

$${\int }_{-\infty }^{\infty }d\omega {\rho }_{\mathrm{v}}(\omega ,|\text{boldmathpunboldmath}|)=0,{\int }_{-\infty }^{\infty }d\omega {\rho }_{\mathrm{s}}(\omega ,|\text{boldmathpunboldmath}|)=0.$$

$$S_4(i{\omega }_p,|\text{boldmathpunboldmath}|)$$

$$=$$

$$S_4(i{\omega }_p,\text{boldmathpunboldmath})$$

$$-i{S}_4(i{\omega }_p)=\sum _{l=1}^{N_{\omega }}\frac{2{\omega }_p}{{\omega }_p^2+{\omega }_l^{\prime 2}}{\rho }_l.$$

$$P[\rho |S_4,I]=\frac{P[S_4|\rho ,I]P[\rho |I]}{P[\rho |I]},$$

$$P[S_4|\rho ,I]=$$

$$L=$$

$$P[\rho |I]=\exp [-\alpha S[\rho ,m]],$$

$${\frac{\delta P[\rho |S_4,I]}{\delta {\rho }_l}|}_{\rho ={\rho }^{\mathrm{Bayes}}}=0,$$

$$\displaystyle S_{\rm SJ}=\int d\omega\big{(}\rho-m-\rho{\rm log}\big{[}\frac{\rho}{m}\big{]}\big{)}\,.$$

$$\displaystyle S_{\rm BR}=\int d\omega\big{(}1-\frac{\rho(\omega)}{m(\omega)}+{\rm log}\big{[}\frac{\rho(\omega)}{m(\omega)}\big{]}\big{)}.$$

$$P[\rho |S_4,I(m)]=P[S_4|I]{\int }_0^{\infty }d\alpha P[\rho |I(\alpha ,m)].$$

$$S_{\mathrm{BR}}^{\mathrm{smooth}}$$

$$=$$

$$⟨\delta {\rho }^2{⟩}_I\approx -{\int }_{I}d\omega d{\omega }^{\prime }({\delta }^{2}Q/\delta {\rho }^2{)}^{-1}/{\int }_{I}d\omega d{\omega }^{\prime },$$

Full text

I Introduction

The wealth of data produced in heavy ion collision experiments at the Relativistic Heavy-Ion Collider (RHIC) and the Large Hadron Collider (LHC) has lead to interesting insights about the nature of the quark-gluon plasma (QGP) in various temperature regimes (see e.g. Refs. Muller and Nagle (2006); Shuryak (2009); Braun-Munzinger and Wambach (2009); Andronic (2014); Foka and Janik (2016) and references therein). Thermal and transport properties of the QGP are encoded in the correlation functions of QCD. In particular they can be assessed from real time properties of QCD’s most basic correlation functions, the quark Nickel (2007); Harada et al. (2008); Harada and Yoshimoto (2009); Karsch and Kitazawa (2007, 2009); Mueller et al. (2010); Qin et al. (2011); Qin and Rischke (2013); Gao et al. (2014) and gluon propagators Strauss et al. (2012); Haas et al. (2014); Christiansen et al. (2015); Dudal et al. (2014); Ilgenfritz et al. (2017). A prominent example is the dilepton production in a heavy ion collision. It can been related to the spectral properties of thermalized quasi-particles and specifically to the dispersion relation of quarks Braaten et al. (1990); Peshier and Thoma (2000); Arnold et al. (2002); Kim et al. (2015a). Another important example are QCD transport coefficients, that have been expressed in terms of single particle spectral functions of the fundamental fields, the quarks and gluons Haas et al. (2014); Christiansen et al. (2015). In summary, a detailed understanding of a potential quasiparticle spectrum in the QGP, in particular close to the chiral phase transition, is highly desirable.

In this work we focus on the quark spectral function encoding the quark dispersion relation and decay width in the medium. At large temperatures reliable results have been obtained in the hard-thermal loop (HTL) expansion Braaten and Pisarski (1990); Baym et al. (1992); Blaizot and Ollitrault (1993). In this regime, the quark spectral function shows two excitations in the dispersion relation, the ordinary quark with a positive ratio of chirality to helicity, and a collective ’plasmino’ mode with a corresponding negative ratio. Both have thermal masses of the order $gT$ and decay widths of the order $g^{2}T$, where $g$ is the coupling constant and $T$ the temperature. The two excitations are accompanied by a continuum contribution from a branch cut in the quark propagator due to Landau damping, i.e. the absorption of a space-like quark by a hard gluon or hard antiquark.

Beyond systematic weak-coupling expansions, there is no straightforward approach for the determination of the spectral function. Model calculations offer qualitative insights Schaefer and Thoma (1999); Kitazawa et al. (2006, 2007); Harada et al. (2008) which, however, need to be corroborated in more fundamental approaches. Models for quark spectral functions constructed along the lines of the HTL results have been fitted to data from quenched lattice QCD Karsch and Kitazawa (2007, 2009) and quenched Dyson-Schwinger calculations Mueller et al. (2010). Again, such an approach offers qualitative insights but suffers from potential biases involved in the model building. This holds in particular for temperatures around the (pseudo-)critical one, where HTL is not expected to be reliable.

In principle, functional approaches like Dyson-Schwinger equations (DSE) and the functional renormalization group (FRG) offer the possibility to determine the two-point correlators in the complex momentum plane thus allowing for a direct extraction of the corresponding spectral function. This has been performed successfully for the gluon propagator at zero temperature in Strauss et al. (2012). At finite temperatures direct computations in the complex freqeuncy plane have been carried out in matter systems, see e.g. Kamikado et al. (2014); Tripolt et al. (2014); Pawlowski and Strodthoff (2015); Strodthoff (2017). However, the additional conceptual and numerical challenges have delayed similar direct analyses in QCD so far.

There are, however, approaches that allow us to extract spectral functions from numerical data in the spatial Euclidean momentum region. These approaches utilise the fact that the spectrum is related to the Euclidean correlator via an integral transform, which needs to be inverted. This is a classic ill-posed inverse problem and Bayesian inference can be used to give meaning to it, by systematically incorporating additional prior information available. Among the different implementations of the Bayesian reconstruction strategy is the popular Maximum Entropy Method (MEM) Jarrell and Gubernatis (1996); Press et al. (1992); Asakawa et al. (2001), which originates in two-dimensional image reconstruction. It has been deployed for the study of quarks in cold and dense matter Nickel (2007) and quarks at and around the (pseudo-)critical temperature Qin et al. (2011); Qin and Rischke (2013); Gao et al. (2014). A method similar in spirit as the MEM but instead using a quadratic regulator has been applied to study gluon spectra in Dudal et al. (2014).

As has been discussed e.g. in Burnier and Rothkopf (2013) the MEM based approaches have to deal with several issues. The main difficulty is that of flat directions in the regulator functional in case of positive definite spectra. In practice this leads to very slow convergence in case that a large number of data points is supplied. The second point is related to the weighting of data and prior information, which conventionally is implemented via computing the so called evidence probability distribution. This step relies on a Gaussian approximation, which in practice is difficult to justify, see also App.C.

In order to overcome these and further difficulties, a novel implementation of the Bayesian strategy has been recently developed Burnier and Rothkopf (2013). It is specifically designed for the solution of one-dimensional inverse problems. Its generalization to arbitrary spectra Rothkopf (2017) has been applied for extracting spectral properties of gluons at finite temperature Ilgenfritz et al. (2017) in lattice QCD. In this study we both deploy the original BR method and develop in addition a new ”low-ringing” BR-type prior functional, which allows us to unambiguously distinguish between peaked structures present in the underlying correlator data and numerical ringing artifacts (see e.g. discussion in Kim et al. (2015b)) common to inverse problems (c.f. Gibbs phenomenon).

We apply our new method to temperature dependent quark propagators obtained from two different truncation schemes for the quark and gluon DSEs. On the one hand we re-analyse a model truncation for the gluon propagator and compare to previous results Qin et al. (2011); Qin and Rischke (2013); Gao et al. (2014). In these works a zero mode in the spectral function at zero momentum has been identified in addition to the two conventional symmetric peaks at finite frequency. Since such a zero mode does not appear in the HTL-studies, it has been attributed to the strong interaction physics governing the transition around the (pseudo-)critical temperature and signalling the formation of the quark-gluon plasma. While the appearance of this structure has been found to be robust under variations within a class of truncation schemes using models for the gluon Gao et al. (2014), it remained to be seen whether this is also true for the fully unquenched system. In this article we therefore study in addition a truncation based on Mueller et al. (2010); Fischer and Luecker (2013) to include results for the unquenched quark propagator with $N_{f}=2+1$. This truncation offers direct control over the Yang-Mills sector by explicitly taking quark loop effects in the gluon propagator into account. The resulting prediction for the unquenched gluon propagator at finite temperature Fischer and Luecker (2013) has been shown to agree with corresponding lattice results of Aouane et al. (2013). Furthermore, the temperature dependence of the chiral condensate evaluated on the lattice Borsanyi et al. (2010) has been reproduced. We therefore may expect realistic and quantitative results for the spectral functions as well.

The article is organised as follows. In the next section we summarize the framework to determine the quark propagator at finite temperature and chemical potential. Since all technical details have been given elsewhere, see Refs. Qin et al. (2011); Qin and Rischke (2013); Gao et al. (2014) for the model gluon and Refs. Fischer and Luecker (2013); Fischer et al. (2014a, b); Eichmann et al. (2016) for the unquenched system, we remain very brief. The section III is devoted to the Bayesian BR method reconstruction and the specific improvements we use in this work. Our results for the quark spectral functions are presented and discussed in section IV. We conclude in section V.

II The quark propagator at finite temperature and chemical potential

II.1 Quark Dyson-Schwinger equation

In order to determine the quark propagator at finite temperature and chemical potential we work with the Euclidean metric version of QCD and use the Matsubara formalism. The renormalized quark Dyson-Schwinger equation is then given by

[TABLE]

Here the inverse full quark propagator is denoted by S^{-1}(i\omega_{p},\text{\boldmathp\unboldmath}) and its inverse bare counterpart by S^{-1}_{0}(i\omega_{p},\text{\boldmathp\unboldmath}). We follow the conventions of Mueller et al. (2010) and explicitly use imaginary arguments for the energy in all functions. The dependence of all functions on the renormalisation point is left implicit and $Z_{2}$ denotes the wave function renormalization constant of the quark. The quark Matsubara frequencies are given by $\omega_{p}=(2n_{p}+1)\pi T$ with temperature $T$. The Dirac structure of the inverse propagators at finite T and $\mu=0$ can be decomposed via

[TABLE]

The dressing functions A(i\omega_{p},|\text{\boldmathp\unboldmath}|), B(i\omega_{p},|\text{\boldmathp\unboldmath}|) and C(i\omega_{p},|\text{\boldmathp\unboldmath}|) carry all non-trivial information on the quarks. The quark mass renormalization constant $Z_{m}(\mu^{2})$, the renormalized mass $m(\mu^{2})$ and $Z_{2}(\mu^{2})$ at the renormalization scale $\mu$ are determined within a MOM-renormalization scheme. The explicit form of the quark self-energy can be written as

[TABLE]

with (Landau gauge) gluon propagator $D_{\mu\nu}$, quark-gluon vertex $\Gamma_{\nu}$, gauge coupling $\alpha$ and ghost renormalization constant $\tilde{Z}_{3}$.

In order to determine the quark propagator self-consistently from its DSE, we also need to specify the dressed gluon propagator and the dressed quark-gluon vertex. The truncation used in this work has been developed over the past years Fischer (2009); Fischer et al. (2010); Fischer and Luecker (2013); Fischer et al. (2014a, b) and is guided by two main ideas. First, lattice results for the temperature dependent quenched gluon propagator are used as external input Fischer et al. (2010); Maas et al. (2012) and unquenched via adding an explicit quark-loop for each of the $N_{f}=2+1$ quark flavors used in this work. The resulting DSE for the gluon propagator is depicted in Fig. 2.

This approach ensures that the unquenched gluon propagator inherits the leading order temperature and chemical potential effects via the quark loops, giving a contribution to the thermal mass. Second order unquenching effects in the Yang-Mills diagrams are neglected. The second element of our truncation is an approximation for the full quark-gluon vertex, which combines information from the well-known perturbative behavior at large momenta and an approximate form of the Slavnov-Taylor identity at small momenta studied long ago by Ball and Chiu Ball and Chiu (1980). The explicit form of this approximate expression for the vertex is discussed in Refs.Fischer et al. (2010); Fischer and Luecker (2013); Fischer et al. (2014b) and shall not be repeated here for brevity.

A much simpler system in terms of numerical effort is obtained by substituting the dressed gluon propagator in the quark DSE by a model together with a bare quark-gluon vertex. Taking into account a Debey-like mass in the longitudinal gluon this reads $(k=q-p)$

[TABLE]

with $P_{\mu\nu}^{T,L}$ the transverse and longitudinal projection operators with respect to the heat bath. The dressing functions are given by

[TABLE]

where the functions

[TABLE]

are defined with ${\cal F}(s)=[1-\exp(-s/4m_{t}^{2})]/s$, $\tau=e^{2}-1$, $m_{t}=0.5\,$GeV, $\gamma_{m}=12/25$, and $\Lambda_{\text{QCD}}=0.234$GeV. With $\sigma D=(0.8\mbox{GeV})^{3}$ we choose $\sigma=0.5$ as a representative value for the model parameters. This is one (the simplest) example of a class of truncation schemes that have been studied in Ref. Gao et al. (2014) and found to agree qualitatively with each other in the resulting spectral functions for the quarks. As explained in the introduction we use this model as a numerically easily accessible reference, which already displays the interesting three-peak structure in the resulting spectral functions.

II.2 Quark spectral functions and representation

The spectral representation of the quark propagator is given by

[TABLE]

with spectral function \rho(\omega,\text{\boldmathp\unboldmath}) parameterized as

[TABLE]

Similar to Mueller et al. (2010) we choose conventions such that the scalar dressing functions themselves agree with those introduced by corresponding lattice studies Karsch and Kitazawa (2009). With a positive definite metric, which is not the case for gauge-fixed QCD, the components of the spectral function would furthermore obey the inequality

[TABLE]

as well as the standard sum rule

[TABLE]

Using specific projection operators for the chirally symmetric and/or static (|\text{\boldmathp\unboldmath}|=0) case one can define positive semi-definite combinations of $\rho_{4},\rho_{v}$ and $\rho_{s}$ using Eq. (9), see e.g. Mueller et al. (2010) for details. In this first study here we restrict ourselves to the reconstruction of the component $\rho_{4}$ only. The corresponding left hand side of Eq. (7) is then also given by the $\gamma_{4}$ component of the quark propagator, i.e.

[TABLE]

This component has the advantage of being antisymmetric in $\omega_{p}$ and the corresponding part of the spectral function being symmetric in $\omega$, which means that we may restrict ourselves in Eq. (7) to positive frequencies and arrive, using Eq. (8), at a purely imaginary kernel

[TABLE]

which is easily implemented numerically for the spectral reconstruction. I.e. we will directly formulate the inverse problem in imaginary frequency space, since the rational kernel in Eq. (12) suppresses spectral information much less than its Euclidean counterpart and thus is ideally suited for use in spectral reconstructions.

III Bayesian Spectral Reconstruction

III.1 Bayesian inference

In this study we use the concept of Bayesian inference to invert the relation of Eq. (12) numerically. The quark two-point function $S_{4}$ is evaluated along discretized imaginary frequencies $\omega_{p}$ with $p\in[1\ldots N_{\rm data}]$, while the spectral function is resolved along $N_{\omega}$ bins in real-time frequencies $\omega^{\prime}_{l}$

[TABLE]

Since the data are obtained from a functional QCD computation it can be provided at an arbitrary number of points $N_{\rm data}$ but contains a finite numerical error $\Delta S_{4}$ due to e.g. the evaluation of intermediate integrals in e.g. Eq. (3). Therefore for any finite $N_{\rm data}$ and $\Delta S_{4}$ a naive $\chi^{2}$ fit of the $N_{\omega}$ parameters $\rho_{l}$ would lead to an infinite number of degenerate solutions all reproducing the data within their uncertainty. Bayesian inference in the form of Bayes theorem

[TABLE]

provides a systematic prescription of how to regularize the otherwise under-determined $\chi^{2}$ fit. It does so by incorporating additional prior information $I$ on the spectrum, in addition to the correlator data $S_{4}$. The posterior $P[\rho|S_{4},I]$ for a test function $\rho$ denotes its probability to be the correct spectrum, given the computed data and prior information. It is proportional to the likelihood probability $P[S_{4}|\rho,I]$ and the prior probability $P[\rho|I]$.

The former encodes how the correlator has been obtained and in case of stochastically independent sampled data is related to the standard $\chi^{2}$ fitting functional

[TABLE]

with likelihood

[TABLE]

Here $S^{\rho}_{4}(i\omega_{p})$ denotes the correlator according to the current choice of test function via Eq. (13) and $C_{pq}$ the covariance matrix of the actual correlator data $S_{4}$. Note that we have formulated the inverse problem in the imaginary frequency domain and not based on Euclidean data. The reason is that the exponentially damped integral kernel and the information loss associated with the latter leads to much worse reconstruction results than the rational Källen-Lehmann kernel used here (for an explicit demonstration see Rothkopf (2017); Pawlowski and Rothkopf (2016)). Note also that in the case of the Euclidean kernel it was found that in order to compute the propagator from a test spectral function $\rho$ with high accuracy, one needs to deploy a logarithmic frequency grid Welzbacher (2016).

Since furthermore in our case the correlator is computed and not sampled we will assign an estimated diagonal covariance matrix to the discrete $S_{4}$ with constant relative errors $\Delta S_{4}/S_{4}={\rm const.}$. $L[\rho]$ possesses many flat directions, which translate into the non-uniqueness of the maximum likelihood $\chi^{2}$ approach.

The second and decisive term in the numerator of Eq. (14) is the prior probability, which acts as a regulator to the likelihood probability, lifting the flat directions of $L$. It tells us how compatible the test function $\rho$ is with available prior information. Conventionally it is expressed as

[TABLE]

where $\alpha$ represents a so called hyper-parameter, which weighs the influence of the data versus the prior. Prior information is encoded in $P[\rho|I]$ in two ways, on the one hand via the functional form of S itself and via the so called default model $m(\omega)$, which by definition is the most probable spectrum in the absence of data, i.e. it represents the unique extremum of $S[\rho,m]$.

Different implementations of the Bayesian strategy differ not only in the regulator $S$ which they employ, but also in how the hyperparameter $\alpha$ is treated, as well as how the most probable spectrum

[TABLE]

is determined numerically. Note that different Bayesian prescriptions can yield different results, as long as $N_{\rm data}$ and $\Delta S_{4}$ are finite. Only in the ”Bayesian continuum limit” all methods, if implemented correctly, must converge to the same result. It is therefore important to test how reconstructions change towards this limit using e.g. mock data tests.

In this study we will work with quarks at high temperature, where the spectral functions are assumed to be positive definite, a property which in turn will enter as prior information.

The popular Maximum Entropy Method for positive definite spectra proposes to use the Shannon Jaynes Entropy as regulator, a choice justified by arguments from two-dimensional image reconstruction

[TABLE]

One carries out the reconstruction multiple times with different values of $\alpha$ and then averages these intermediate results $\rho_{\alpha}$ self consistently over a probability distribution for $\alpha$, $P[\alpha|S_{4},I]$. In the standard implementation by Bryan the functional space from which the test function $\rho$ is chosen is manually limited to a dimensionality of $N_{\rm data}$. It has been shown that the extremum of $P[\rho|S_{4},I]$ in general is not contained in this choice of search space, which may result in slow convergence. In additions the determination of $P[\alpha|S_{4},I]$ relies on a Gaussian approximation.

In this study we instead use a more recent Bayesian approach to spectral function reconstruction that has been developed with the one-dimensional inverse problem of Eq. (13) in mind and its regulator functional reads

[TABLE]

The hyperparamter $\alpha$ is treated in a Bayesian fashion, in that we assume full ignorance of its values $P[\alpha]=1$ and integrate it out a priori

[TABLE]

In addition we also require that $L=N_{\rm data}$ as the correct spectrum, on average, would lead to such a value. The most probable spectrum is then obtained from carrying out a numerical optimization on $P[\rho|S_{4},I(m)]$ according to Eq. (18). In contrast to the MEM we do not restrict the search space and use a pseudo-Newton method (limited Broyden-Fletcher-Goldfarb-Shanno) to find the global extremum in the full $N_{\omega}$ dimensional search space.

The regulator $S_{\rm BR}$ was derived with the goal to limit its influence on the outcome of the reconstruction to a minimum. I.e. it shall influence the reconstruction as weakly as possible and ”let the data speak”. Comparing $S_{\rm BR}$ to $S_{\rm SJ}$ or the quadratic prior one finds that far away from the extremum $\rho=m$ $S_{\rm BR}$ shows the weakest curvature. While this makes it easy for structures encoded in the data to manifest themselves in the reconstruction it also means that common artefacts associated with inverse problems, such as Gibbs ringing, are not efficiently suppressed.

In turn, if the structures of physical interest are spectral peaks, as will be the case in the following, we have to make sure that we unambiguously distinguish ringing from genuine features. Due to its restricted search space the MEM usually produces smooth features and it is considered safe from ringing. This impression is unfortunately incorrect on the level of the regulator as illustrated in Fig. 3. For a flat default model $m=1$ we compare the penalty assigned to a function $\rho_{\rm mock}$ with a single broad feature as well as after introducing a wiggle close to its peak $\rho_{\rm test}$. The area under the spectra has been kept the same. What one finds is that both $S_{\rm BR}$ and $S_{\rm SJ}$ assign a lower penalty to the distorted curve even though the former was explicitly derived with a smoothness axiom. What distinguishes the two curves is of course their arc-length, which increases for every additional wiggly structure.

Hence we introduce a new BR-type regulator, which penalizes arc-length $\ell=\int d\omega\sqrt{1+(d\rho/d\omega)^{2}}$ explicitly. In order to leave as much of the original form of $S_{\rm BR}$ intact we add a $\ell^{2}$ term and subtract unity. The result is

[TABLE]

In the presence of the derivative term, we are not any more able to analytically compute the $\alpha$ dependent normalization of the prior probability and thus cannot marginalize $\alpha$ as in done in Eq. (21). Thus we revert to handling $\alpha$ in the same way as in the ”historic MEM” prescription, in which one adjusts $\alpha$ in order that $(L-N_{\rm data})<10^{-1}$.

As we will see also explicitly in the mock data test in the following section, this new regulator succeeds in efficiently suppressing ringing artefacts. On the other hand peaks encoded in the correlator become more washed out and it requires more data points and smaller errors to resolve these peaks to the same accuracy as with the original BR method. Hence our strategy is the following. Using the smoothening prior we will identify in which part of the spectrum actual peak structures reside and then extract their features using the standard BR approach.

The reliability of the reconstruction can be estimated in three different ways. Since in Bayes theorem in Eq. (14) data and prior information enters, we need to understand the dependence of the reconstruction result on that input. For the former we can vary the number of provided data points and carry out a Bootstrap Jackknife resampling analysis of the correlator. Since here we do not use sampled data we will instead successively lower the assigned error on the data and observe convergence toward the Bayesian continuum limit. The dependence on the prior can be estimated by repeating the reconstructions with different choices of the default model, for which we choose a flat function $m(\omega)=m_{0}$ and vary $m_{0}=\{0.1,0.5,1.0,5.0,10\}$.

Bayesian methods also provide another measure for the robustness of the reconstruction

[TABLE]

based on the curvature of the optimization functional $Q[\rho,S_{4},I]={\rm Log}[P[\rho|S_{4},I]]$. In previous works it has been found that the actual dependence on the variation of data and prior information is often larger than indicated by this quantity. One possible reason for this discrepancy is the need for two assumptions in the derivation of Eq. (23), namely that the posterior is both highly peaked and may be approximated by a Gaussian. In the following we will show error bands for the BR method that are obtained from the variation of the default model and the curvature measure, where the former dominates.

III.2 Mock data tests

Previous studies Qin and Rischke (2013) deploying a MEM-type approach hinted at the presence of a three peak structure of the quark spectral function in the Landau gauge at finite temperature. On the other hand it is was not excluded that further peaked features may be present. Therefore we must ascertain how well our reconstruction method will be able to resolve different structures given a certain quality of input data, for which we turn to a mock data analysis.

In the following we compute the Euclidean frequency correlator $S_{4}$ according to Eq. (7) based on mock spectra with different peak contents. This continuous function is then evaluated at $N^{\rm mock}_{\rm data}=128$ discretely spaced imaginary frequencies between $\omega_{p}\in[0,8]$GeV. This ideal data is fed to the reconstruction algorithm undistorted but is assigned a constant relative error $\Delta D/D$. Increasing the number of data points beyond 128 has not shown significant improvements down to $\Delta D/D=10^{-8}$.

By successively lowering $\Delta D/D$ we have investigated the approach to the Bayesian continuum limit for fixed $N_{\rm data}$ as explicitly shown in App.A. In the following we will showcase only two of these reconstructions for each mock spectrum, which are relevant for the discussion. One is the best possible outcome using a particular Bayesian implementation, i.e. using a very small $\Delta D/D=10^{-8}$. The other is the realistic case $\Delta D/D=10^{-3}$, which corresponds to the precision easily obtainable in realistic Dyson-Schwinger computations Welzbacher (2016).

As a first step let us have a look at a two-peak scenario with one structure located at the origin and a second at $\omega=2$GeV, as shown by the grey curve in Fig. 4. The two solid curves of different brightness show the reconstructions either based on the smoothed BR method (top) or its original implementation (bottom). The darker curve corresponds to the error $\Delta D/D=10^{-3}$ and the lightest color to $\Delta D/D=10^{-8}$. The light colored bands around the curves refers to the corresponding dependence on the choice of the prior.

Lowering the errors consistently improves the reconstruction of both the peak positions and widths. For $\Delta D/D\leq 10^{-3}$ the peak position of the second peak is accurate within $15\%$ for the original BR method, while the lowest structure is already captured satisfactorily. Note that for the smoothed regulator it is more difficult to reproduce the correct width of the higher lying peak at the same $\Delta D/D$ than for the original BR method. On the other hand the reconstruction based on $S^{\rm smooth}_{\rm BR}$ is free of any of the numerical ringing that appears at intermediate $\Delta D/D$ in the original BR method. Our strategy is proven valid, the smoothed BR reconstruction unambiguously tells us about the presence of two peaked feature and we can use the standard BR method to more accurately extract their properties.

Another possible scenario is the split of one of the peaks into two structures, in particular the emergence of an additional structure close to the origin. If the peak at large frequencies is well separated from these, then again neither the conventional nor the smooth BR method are challenged in identifying the three structures as seen by the results of Fig. 5.

A more difficult scenario for any Bayesian reconstruction arises in the presence of an additional peak positioned closely to the one off the origin. Here close is understood as a distance which is comparable to the width of the peaks, as shown in Fig. 6 (grey dashed).

Both methods struggle to identify the split between the two peaks even at optimal conditions $\Delta D/D=10^{-7}$ or $\Delta D/D=10^{-8}$. On the other hand it is important to note that the reconstructions of the two different methods show qualitatively different behavior in contrast to all previous scenarios. While before the shape of the peaks in both the conventional and smooth BR method agreed, here we see that at $\Delta D/D=10^{-8}$ the original BR method shows a distorted peak.

We conclude that the combination of the original and smooth BR method is promising for the investigation of quark spectral functions. The latter promises, given small enough errors, to determine the number of peaked features present in the data. The former on the other hand, while susceptible to ringing at realistic $\Delta D/D\sim 10^{-3}$ is capable of reproducing the actual peak properties more accurately based on the same quality of data.

IV Results

In this section we apply the previously tested Bayesian approaches to actual correlator data $S_{4}$ obtained from Dyson-Schwinger computations. We first re-analyse the rainbow-ladder model approach described at the end of section II.1. In this simple truncation the determination of the correlator is computationally cheap and numerical errors are easily reduced to $\Delta D/D<10^{-3}$. Subsequently we discuss our main results based on the truncation scheme with $N_{f}=2+1$ unquenched light flavors, back-coupled to the Yang-Mills sector.

IV.1 Rainbow ladder truncation

IV.1.1 Chiral limit

Let us first discuss the case of the chiral limit, which generates a second order phase transition at a critical temperature of $T_{c}=0.142\,$GeV. We have evaluated the corresponding $S_{4}$ along Matsubara frequencies at twelve temperatures above $T_{c}$ in the range of $T\in[0.15,0.60]\,$GeV. Due to the used cutoff $\Lambda=100$GeV and the discrete nature of the Matsubara frequencies this corresponds to $N_{\rm data}\in[106,27]$. The correlator computations have been checked to carry a numerical error of less than $\Delta D/D\leq 10^{-3}$, so that we can assign a corresponding diagonal correlation matrix to it. We deploy the conventional and smooth BR method with a frequency discretization of $N_{\omega}=1000$ in the interval $\omega\in[0,20]$. Note that even at low temperatures the reconstructions converged swiftly and do not show any signs of spiky defects indicating the presence of negative spectral contributions. The default model is set to to a constant $m(\omega)=m_{0}$ and we carry out the reconstruction with the different choices $m_{0}=\{0.1,0.5,1.0,5.0,10\}$. The variance in the outcome is taken as the basis for the error bands depicted in the subsequently shown plots.

In Fig. 7 we present the zero momentum reconstruction of the quark spectral function for different temperatures. Note that since we explicitly use the symmetry of the spectral function in Eq. (12) the position of the peak at positive frequency $\omega_{+}$ is mirrored exactly in the negative frequency domain. In the following it is therefore sufficient to always discuss the positive frequency band only. In the top panel the outcome of the smooth BR method is shown, while the bottom panel corresponds to the conventional one. We learn that at low temperatures at least up $T\lesssim 0.25\,$GeV two well defined peak structures exist. One is located at the origin and one above $1\,$GeV. With increasing temperature the height of the low lying peak decreases continuously and seems to asymptote to a finite value. The higher lying peak shows a clear tendency to move to higher frequencies, while broadening at the same time.

The reconstructions at different temperatures are based on dataset with different $N_{\rm data}$. Thus before we continue to a more quantitative inspection of the spectra we need to make sure that we can disentangle the actual in-medium modification from the effects of a reduction of data points. To this end we perform the following test: we take the lowest temperature dataset with $N_{\rm data}=106$ in imaginary frequencies and sparsen it by hand. Due to the discrete nature of the Matsubara frequencies this corresponds to a situation, where the $T=0.15\,$GeV spectrum would be encoded in a correlator evaluated at $T=0.3\,$GeV or $T=0.6\,$GeV respectively. In Euclidean time this corresponds to constructing the reconstructed correlator Ding et al. (2012). From similar tests performed in previous Bayesian studies we expect that with decreasing number of data points the resolution of peaks diminishes, eventually inducing changes in the position and width of the reconstructed features.

And indeed, as shown in Fig. 8 we find that going from $N_{\rm data}=106$ to $N_{\rm data}=53$ (dark grey dashed) leads to a visible weakening of the peak structures, while their position remains unaffected. Comparing to the reconstruction from actual $T=0.3\,$GeV with the same lower number of data points (green solid) however shows clear differences. Both the diminishing of the lowest lying peak height, as well as the shifting of the second peak to higher frequencies can thus be attributed to genuine in-medium effects. Interestingly the direction of change in the peak position is opposite to that sketched in Qin and Rischke (2013); Gao et al. (2014). At $T=0.6\,$GeV, i.e. $N_{\rm data}=26$ the sparsened data do not allow the reconstruction of the two peaks at all and we must assume that the reconstruction is not trustworthy at this point.

The investigation of the effects of finite momentum on the quarks does not suffer from a similar ambiguity, since for fixed $T$ neither the number of data points nor the relative errors change. In Fig. 9 we show reconstructed spectra at $T=150\,$MeV respectively over a range of momenta |\text{\boldmathp\unboldmath}|/T=\{0,\frac{1}{2},1,\frac{3}{2},2,\frac{5}{2},3\}.

As seen before at low temperature both reconstruction approaches unambiguously show the presence of two peaks. One is located around the origin, another one positioned close to $\omega\approx 1\,$GeV. Increasing the momentum to |\text{\boldmathp\unboldmath}|/T=3 induces changes in the spectrum that are of the same qualitative nature as those from increasing temperature. The peak at the origin decreases significantly in area, while not extending further towards higher $\omega$. The second peak diminishes much more weakly and is seen to shift to higher frequencies, as expected from the naive momentum dependence of the dispersion relation.

For visualization purposes we present in Fig. 10 the reconstructions at a fixed intermediate momentum |\text{\boldmathp\unboldmath}|/T=1 (top) and |\text{\boldmathp\unboldmath}|/T=3 (bottom) for all different temperatures investigated. All the effects on the peak around the origin, as well as the second peak at finite frequencies as discussed above are clearly visible here.

We continue with a quantitative analysis of the in-medium modification of the quark spectral features. Two quantities of interest here are the position of the higher lying peak denotes with $\omega_{+}$ and the height of the peak around the origin referred to as $Z_{0}$. Both are shown in Fig.11. For $\omega_{+}$ the expectation from resummed hard-thermal loop perturbation theory at small |\text{\boldmathp\unboldmath}|/T\ll 1 is a linear dependence on the temperature \omega_{\pm}^{\rm HTL}=m_{T}\pm\frac{|\text{\boldmathp\unboldmath}|}{3} with thermal quark mass $m_{T}\propto T$. While at low temperatures and small momenta we see a rise stronger than linear, at intermediate $T$ our $\omega_{+}$ indeed shows a behavior compatible with a linear increase. Consistent with our conclusions from the sparsening test, for temperatures much higher than $T=0.3\,$GeV, the reconstruction becomes unreliable and at the same time we see that the linear rise abates and goes over to a constant.

The behavior found here is again different from that presented in Gao et al. (2014). Firstly, we find that the peak at non-zero frequency $\omega_{+}$ monotonously moves to larger frequencies with increasing temperature in contrast with the previously reported behavior, where the function $\omega_{+}(T)$ shows a minimum shortly above $T_{c}$. The second difference is related to the height of the zero frequency peak, which in the temperature range investigated does not go to zero but stays finite. Since the high temperature regime is not reliably captured due to the sparseness of the Matsubara frequencies we cannot yet conclude whether it eventually vanishes after all.

In all quantitative statements we need to keep in mind that our current numerical precision for the correlator data is limited to $\Delta D/D=10^{-3}$. While we are confident, judging from the mock data tests, that the number of peaks and the direction of changes with temperature are correctly captured, it is fathomable that the peak position may still change by up to $20\%$ if the errors are further reduced. The width of the peaks carries at least the same uncertainty at this point along the ”Bayesian continuum limit”.

IV.1.2 Finite mass case

We proceed with a second set of correlators $S_{4}$, for which the current quark mass has been set to $m_{q}=3.7\,$MeV at a renormalization point of $19\,$ GeV. Here we restrict ourselves to the $p=0$ case. The question to answer is how the spectral structures differ in contrast to the chiral case. We use the same temperature regime and discretization of the correlator data and leave the errors unchanged.

The results for the zero momentum spectral reconstructions at different temperature are given in Fig. 12, with the smooth BR method in the top panel and the conventional one in the bottom one. Qualitatively the figures are very similar to the results in the chiral case. There exist two peaks, one at the origin and one above $1\,$GeV. The position of the second peak moves to higher frequencies as temperature increases, while the height of the lowest lying peak decreases continuously. The only visible difference is that at low temperature the height of the peak around the origin is discernibly smaller than in the chiral case.

We again checked that the changes between the outcome at different temperature is indeed attributable to in-medium effects by manually coarsening the $T=0.15\,$GeV correlator data and repeating the reconstruction with it. The results of this procedure are similar to the ones in the quenched case: the deletion of every second data point, i.e. $N_{\rm data}=56$, weakens the peak strength while leaving the peak position unaffected. The reconstruction based on only every fourth data point corresponding to the situation at $T=0.6\,$GeV fails to identify the encoded two peak structure and is therefore deemed not trustworthy.

Just as in the chiral case, we determine the position of the second peak $\omega_{+}$ and the height of the peak around the origin $Z_{0}$, shown in top and bottom panel of Fig. 13 respectively111 Note that the quantity $\omega_{+}$ shown here does not fully coincide with the conventional thermal mass at |\text{\boldmathp\unboldmath}|=0, as it is reconstructed from $S_{4}$, instead of $S_{4}\pm S_{s}$. Since the latter makes the spectral function non-symmetric we postpone its analysis to future work.. As is to be expected from the close resemblance of the reconstructed spectra, the values obtained do not differ markedly between the finite mass (triangle) and the chiral (circle) case. The only differences are found at low temperatures. The values of $\omega_{+}$ for finite mass actually show a linear behavior down to $T=0.15\,$GeV, whereas in the chiral case they deviate from a straight line at that point. $Z_{0}$ is smaller at low temperatures in the presence of a finite quark mass but already at $T=0.2\,$GeV no significant difference remains.

IV.2 Unquenched truncation with back-coupled $N_{f}=2+1$ quark flavors

Our main result concerns the quark spectrum computed from the quark and gluon Dyson-Schwinger equations in a modern truncation incorporating $N_{f}=2+1$ light quark flavors with physical masses. In this case we encounter a chiral crossover at a pseudo-critical temperature of $T=0.155\,$GeV via the inflection point of the quark condensate and $T=0.160\,$GeV via the chiral susceptibility in agreement with corresponding lattice results. We have computed $S_{4}$ along Matsubara frequencies at vanishing spatial momentum for eleven temperatures in a larger temperature range of $T\in[0.125,1.0]\,$GeV, i.e. also for temperatures below the pseudo-critical one. At the same cutoff of $\Lambda=100$GeV as before this now corresponds to $N_{\rm data}\in[127,16]$. The correlator computations have been checked to carry a numerical error of less than $\Delta D/D\leq 10^{-3}$, so that we can assign a corresponding diagonal correlation matrix to it. Both the smooth and original BR method are carried out with a frequency discretization of $N_{\omega}=1000$ in the interval $\omega\in[0,20]$. The default model is set to to a constant $m(\omega)=m_{0}$ and we carry out the reconstruction with the different choices $m_{0}=\{0.1,0.5,1.0,5.0,10\}$. The variance in the outcome is taken as the basis for the error bands depicted in the subsequently shown plots. In order to keep the presentation of the reconstructed spectra clear, we show in Fig.14 only a subset of the reconstruction in a temperature range pertinent to the discussion below. For completeness the full results are plotted in the appendix B in Fig. 20.

The first interesting result is that the reconstructions at low temperatures show unstable behavior that hints at a failure of the reconstruction limited to positive definite spectra. We see in Fig. 14 that at $T=0.125\,$GeV two sharp peaks appear at positions very different from those at higher $T$. The results at and above the transition region, $T=0.150\,$GeV and $T=0.175\,$GeV appear to be in better shape at first glance. However, truncating these datasets (e.g. from $N_{\rm data}=106$ to $N_{\rm data}=94$ for the lowest temperature) actually changes the behavior of the spectral functions around $\omega=0$. The same tests for the two analyses in the previous subsections had shown virtually no effect on the reconstruction, which is what is expected for a well converged result. We thus conclude, that in the truncation scheme with back-coupled quarks positivity violations characteristic for the low temperature spectral functions of the quarks persist for much larger temperatures than in the simple model case and prohibit convergence. Only the spectra at and above $T=0.2\,$GeV do not show such artificial behavior and are therefore deemed trustworthy.

This is also indicated in Fig. 15, where we show the reconstruction outcome of taking the correlator at $T=0.2\,$GeV and $N_{\rm data}=80$ and sparsening it by factor two (dark grey dashed) or factor four (light grey dashed). For $N_{\rm data}=40$ the reconstruction is only very weakly affected, while for $N_{\rm data}=20$ a sole peak at the origin remains. We conclude as before that the reconstruction eventually becomes unreliable at high temperature but that at intermediate $T$ we are able to observe genuine in-medium modification.

We now analyze the position of the peaks. Although we find the same number of peaks present as in the model calculation, their behavior under variations of temperature seems to be different. Whereas the amplitude of the lowest lying peak still decreases with increasing temperature, the position of the second peak does not show a clear pattern. In particular it appears to not move monotonously to higher values of frequency with increasing temperature.

This is obvious from Fig. 16, where we plot again the location $\omega_{+}$ of the second peak (upper panel) and the amplitude $Z_{0}$ of the zero mode (lower panel). For purpose of comparison we also depict the results from the finite mass rainbow ladder truncation as circles. Clear qualitative and quantitative differences are visible. Instead of monotonously rising in value, $\omega_{+}$ appears to decrease first up to around $T=0.4\,$GeV before then increasing again in an almost linear fashion towards higher temperatures. Interestingly the initial downward trend starts in the low temperature regime, where we did not deem the reconstruction reliable due to the possible presence of positivity violation. Then we must further clarify whether the behavior of $\omega_{+}$ up to $T=0.4\,$GeV might still suffer from the influence of residual positivity violation. This will require the application of a reconstruction algorithm for general spectral functions, which is foreseen as next step in this line of study.

$Z_{0}$ on the other hand behaves at least qualitatively similar in the region where we trust the reconstruction. Below $T=0.15\,$GeV it also shows a clear dip, which is related to the appearance of the artificial spiky structures at intermediate frequencies there. Above $T=0.2\,$GeV it decreases monotonously.

Compared to the values reported in Gao et al. (2014) the behavior of $\omega_{+}$ here is quite similar. If the reconstructions between $T=0.2\,$GeV and $T=0.4\,$GeV are reliable, in particular as they do not show any obvious pathologies, then we also observe a dip in $\omega_{+}$ at intermediate temperatures. The height of the central peak on the other hand never fully vanishes in our case.

V Conclusions

We have investigated the spectral properties of quarks in the Landau gauge, based on Dyson-Schwinger equations according to two different truncation schemes. In the rainbow-ladder approximation model both the chiral and finite current quark mass case have been considered, while our main result concerns quark spectra in a modern truncation with $N_{f}=2+1$ unquenched flavors of light medium quarks.

The reconstruction of the spectral functions is based on a recently developed Bayesian approach, the so called BR method, formulated in imaginary frequencies. We further developed in this study a low-gain variant of the BR method, which successfully suppresses numerical ringing, which can affect the original BR method and in turn helps us to unambiguously determine the number of physical peaks in the spectrum. The accuracy of the reconstruction further benefits from the use of the Källen-Lehmann kernel instead of the Euclidean one.

In mock data tests we have shown the capabilities and limitations of our Bayesian reconstruction approach for either a best-case scenario with correlator precision $\Delta D/D=10^{-8}$ and a real-world setting with $\Delta D/D=10^{-3}$. In cases with two or three peak features, which are expected to be relevant for our study the combination of the conventional and smooth BR method allowed us to unambiguously identify the number of encoded peaks and estimate their properties. In the most challenging but least likely case that two rather broad peaks at high frequency are located close to each other it required the best case scenario to infer the presence of all features.

The reconstructed spectra for the rainbow ladder truncation model with vanishing and finite current quark mass showed very similar behavior. At low temperatures two peaks are present, one at the frequency origin, another one above $\omega=1$GeV. Changing the temperature or changing the spatial momenta induced qualitatively similar changes. The lowest lying peak height diminishes but does not vanish up to the highest parameter values investigated. The second peak both broadens and moves to higher frequencies.

A quantitative analysis of the height of the low lying peak $Z_{0}$ and position of the second peak $\omega_{+}$ revealed a different behavior as reported in previous studies. We did not find any indication of a non-monotonicity in $\omega_{+}$ with respect to temperature and our value for $Z_{0}$ always took on finite values in contrast to a vanishing $\omega=0$ peak in Gao et al. (2014).

We have made sure that the observed changes in $Z_{0}$ and $\omega_{+}$ with temperature can be attributed to the thermal medium. To disentangle the effects from a degrading of the reconstruction due to less available data points at high temperature, we repeated the reconstructions with manually sparsened correlator data sets and identified the regime, where the Bayesian method is reliable. And indeed, in the region where the reconstruction can be trusted we find that $\omega_{+}$ shows a linear rise with $T$ qualitatively compatible with hard-thermal loop predictions. At the same point where the reconstruction becomes unreliable we also see that the linear rise begins to artificially flatten off.

In the unquenched truncation scheme with $N_{f}=2+1$ flavors of light medium quarks the positive definite Bayesian approach is challenged at low and high temperature. For $T<0.2\,$GeV we find indications that non-positive spectral contributions are present, which lead to artificially spiky structures, while a sparsening test shows that for $N_{\rm data}<20$ the reconstruction also becomes unreliable. In the intermediate temperature window we observe again two peak structures with the lower one decreasing monotonously in height. The second peak however behaves very differently than before as it now appears to exhibit a dip in $\omega_{+}$, similar to the behavior reported in Gao et al. (2014).

In order to unambiguously determine, whether the non-monotonous behavior of $\omega_{+}$ can be attributed to physics encoded in the correlator, we will have to extend the analysis of this study in the future to non-positive spectral functions. In the context of gluon spectral functions in lattice QCD, a generalization of the BR method has been proven a useful tool Ilgenfritz et al. (2017). Implementing a smooth version of this generalized BR method will constitute an important step towards a robust and quantitative picture of the low temperature regime of quark spectra, which is work in progress.

Acknowledgement

C.F. and C.W. were supported by the German Federal Ministry of Education and Research (BMBF) under Contract No. 05P15RGFCA and the Helmholtz International Center for FAIR within the LOEWE program of the State of Hesse. J.P. and A.R. acknowledge support by EMMI, the grants ERC-AdG-290623, BMBF 05P12VHCTG as well as that this work is part of and supported by the DFG Collaborative Research Centre ”SFB 1225 (ISOQUANT)”.

Appendix A Mock test Bayesian continuum limit

In this appendix we present figures for the explicit approach of the mock spectral reconstructions of Sec.III.2 towards the Bayesian continuum limit at fixed $N_{\rm data}=128$. Each of the figures contains seven solid curves, denoting the reconstruction according to an assigned relative error $\Delta D/D=[10^{-8},10^{-2}]$. In the upper panel the smooth BR is deployed, while in the lower one it is the original BR method.The darkest curve corresponds to the largest error. In Fig. 17 the two-peak scenario is shown, while Fig. 18 contains the results for three peak scenario with the second peak lying close to that at the origin. The last Fig. 19 shows the outcome for a three peak scenario with two closely places structure at finite $\omega$.

Appendix B Complete spectral reconstructions for the unquenched case

Here we plot in Fig.20 for completeness the spectral reconstructions for the unquenched truncation with back-coupled $N_{f}=2+1$ quark flavors. One observes the appearance of artificial peaked structures at the lowest teperature $T=0.125$GeV, while at $T=0.2$GeV and above we obtain the same number of peaks as in the model computations. As was discussed in the context of Fig.16 the peak position of the peak located at finite frequencies however displays a qualitatively different behavior here than in the model truncation.

Appendix C Test of the evidence probability computation

In contrast to the BR method, where the hyperparameter $\alpha$ is self-consistently integrated out apriori, MEM-like approaches marginalize $\alpha$ at the end of the reconstruction procedure Jarrell and Gubernatis (1996); Asakawa et al. (2001). I.e in MEM one conventionally computes the corresponding spectrum $\rho_{\alpha}$ for many different values of $\alpha$ and then determines the probability distribution $P[\alpha|\rho,D,I]$. The individual $\rho_{\alpha}$ are subsequently averaged, weighted by $P[\alpha|\rho,D,I]$. In order to compute $P[\alpha|\rho,D,I]$ one however relies both on the assumption that the posterior probability is highly peaked and that it allows for a Gaussian approximation. Both are not tested in practice.

Common lore states that $P[\alpha|\rho,D,I]$ will have a peak at a finite $\alpha$ for which one particular spectrum contributes most strongly. If the maximum were at $\alpha=0$ the method reverts to an under-determined $\chi^{2}$ fit and no unique extremum exists. Here we give numerical evidence that the existence of a peak in the approximated $P[\alpha|\rho,D,I]$ depends on the choice of search space used. Furthermore if the search space is extended to the full size of the problem (in which there still exists a unique Bayesian answer) we find that only a maximum at $\alpha=0$ remains.

We use the same mock data as in the two-peak scenario in sec.III.2 and compare four different scenarios. We deploy the MEM with limited search space $S_{\rm BR}$ and $N_{\rm base}=N_{\rm data}$ according to Bryan and compare with (solid line) an implementation without restriction, where $N_{\rm base}=N_{\omega}$ (dashed line). The Bayesian result is selected with a step tolerance in the minimizer of $\Delta=5\times 10^{-8}$. In addition we replace the Shannon-Jaynes Entropy by the BR prior and repeat the reconstruction with $S_{\rm BR}$ and restricted search space $N_{\rm base}=N_{\rm data}$ (solid line) or without, i.e. using $S_{\rm BR}$ and $N_{\rm base}=N_{\omega}$ (dashed line). We have of course adapted the computation of $P[\alpha|\rho,D,I]$ to this new prior. The results for the probabilities are shown in Fig. 21.

We find indeed that only for the restricted search space a peak at finite values of $\alpha$ remains. This issue is independent of the actual regulator used, both $S_{\rm SJ}$ and $S_{\rm BR}$ show the same trend. We believe that the underlying reason is that in the presence of a restricted search space the minimizer is at some point not able to lower the value of $L$, while in the full search space it can be brought very close to zero. This finite minimal value of $L$ then prohibits the probability to rise further. Since the Bayesian answer is unique if it exists Asakawa et al. (2001), we conclude that it is the approximations made to determine $P[\alpha|\rho,D,I]$, which prevent us from obtaining that unique result in the full search space.

As a consequence we revert to the historic MEM choice of setting $\alpha$ such that $L=N_{\rm data}$ in case of the smooth BR method, where an apriori marginalization of the hyperparameter is not analytically feasible.