Bayesian study of relativistic open and hidden charm in anisotropic lattice QCD

We present the first combined study of correlators and spectral properties of charmonium and open-charm-mesons at finite temperature using a fully relativistic lattice QCD approach. The QCD medium is captured by second generation anisotropic 24^3xN_t lattices from the FASTSUM collaboration, including 2+1 flavors of clover discretized quarks with m_pi~380 MeV. Two Bayesian methods are deployed to reconstruct the spectral functions, the recent BR method as well as the Maximum Entropy Method with Fourier basis. We take particular care to disentangle genuine in-medium effects from method artifacts with the help of the reconstructed correlator. Consistent with the direct inspection of correlators, we observe no significant in-medium modification for J/Psi and eta_c around the crossover, while the chi_c states on the other hand show clear changes around the transition. At the highest temperature, T=352 MeV, J/Psi and eta_c also exhibit discernible changes compared to the vacuum. For D mesons around $T_c$, no significant modifications are observed, but we find clear indications that no bound state survives at the highest temperature of T=352 MeV. Above T_c we discover a significant difference between D and D* mesons, the latter being much more strongly affected by the medium.

We present the first combined study of correlators and spectral properties of charmonium and open-charm-mesons at finite temperature using a fully relativistic lattice QCD approach. The QCD medium is captured by second generation anisotropic 24 3 × Nτ lattices from the FASTSUM collaboration, including 2+1 flavors of clover discretized quarks with mπ ≈ 380MeV. Two Bayesian methods are deployed to reconstruct the spectral functions, the recent BR method as well as the Maximum Entropy Method with Fourier basis. We take particular care to disentangle genuine inmedium effects from method artifacts with the help of the reconstructed correlator. Consistent with the direct inspection of correlators, we observe no significant in-medium modification for J/Ψ and ηc around the crossover, while the χc states on the other hand show clear changes around the transition. At the highest temperature, T = 352MeV, J/Ψ and ηc also exhibit discernible changes compared to the vacuum. For D mesons around Tc, no significant modifications are observed, but we find clear indications that no bound state survives at the highest temperature of T = 352MeV. Above Tc we discover a significant difference between D and D * mesons, the latter being much more strongly affected by the medium.

PACS numbers:
Heavy quarks and the mesons formed by binding these to pairs with either light (open heavy-flavor) or other heavy quarks (hidden heavy-flavor, heavy quarkonium) are an essential part of our toolset to probe the physics of the quark-gluon plasma in relativistic heavy-ion collisions at RHIC and LHC. A close collaboration between experiment [1] and theory [2,3] has emerged to shed light on the various physics aspects which need to be understood in order to produce a comprehensive understanding of the heavy particle yields measured in experiment.
In the early stages of the collision, the production of heavy quark-antiquark pairs from partonic processes, such as gluon splitting, needs to be described. Due to the large momentum transfers involved, QCD perturbation theory provides important insights here (see [1] for a more detailed overview).
Subsequently the heavy quarks may propagate through the collision center, interacting strongly with the bulk matter and in turn lose a significant amount of their original kinetic energy via radiative or collisional processes. The physics of in-medium propagation and energy loss has been investigated via a multitude of models often based on a Boltzmann equation [4] or, if further simplified, on a Langevin or equivalently Fokker-Plank equation [5,6]. The transport coefficients characterizing the evolution have been estimated either via resummed perturbation theory [7,8], strong coupling holography [9,10] or lattice QCD simulations [11]. Direct perturbative analysis of energy loss in a thermal medium [12], computations based on a Dyson-Schwinger like T-matrix resummation [13], as well as modeling with microscopic kinetic equations [14] have provided further insightful results.
Just after the production of a quark-antiquark pair these two particles may actually form a quarkonium bound state. Even in the absence of a QCD medium, this binding is a genuinely non-perturbative process, and its theoretical description has benefitted greatly from the separation of scales between the heavy quark mass and other relevant scales in the system, such as the (local) temperature or the intrinsic scale of quantum fluctuations in QCD, Λ QCD . It allows the construction of non-relativistic field theories (EFT) [15], which provide a simplified language of the physics relevant to heavy quark binding and quarkonium production. In turn the in-medium modification of heavy quarkonium properties can be related to both screening effects in the thermal medium, captured in the Debye mass of the system, as well as scattering effects, such as e.g. Landau damping [16] or color singlet to octet transitions [17].
While for beauty quarks the scale separation is indeed pronounced and EFT approaches are well established, the situation for charm quarks is less straight forward. Estimates of initial temperatures in the collision center from relativistic hydrodynamics give values of T LHC init ≈ 600 MeV, the charm quark mass being barely double this. Thus for charmonium and its open heavyflavor cousins D and D * a genuinely relativistic description is usually preferred. Standard lattice QCD simulations, which routinely provide vital insight into the properties of light hadrons at zero and finite temperature, can be directly applied to the charm sector too, as the required lattice spacings, while small, remain computationally manageable.
We may ask what insight can be gained from performing a thermal lattice simulation of charmonium and open-charm mesons? As the simulated system is in genuine equilibrium, we may learn about the probability of heavy-quark binding being able to sustain bound states, either open or hidden, in the late stages of the collision.
And since in most models the production of open heavyflavor occurs via a hadronization prescription at freezeout, from the theory point of view this is just the relevant regime to look at. The success of the statistical model of hadronization to predict charmonium yields from RHIC to LHC energies [18] further supports a focus on the late stages of the collision to understand not only open, but also hidden charm production.
Considering the binding properties of both quarkonium and open-charm together is valuable, since the latter system provides the natural baseline for the measurements of suppression in the former. In model computations in a hot pion medium e.g., D meson modifications were found to strongly influence J/ψ production during the expansion of the fireball [19].
Recent measurements of a finite D meson flow [20] and even J/Ψ [21,22] flow have unambiguously shown that the charm quarks at the LHC are in at least partial kinetic equilibration with the surrounding medium. In turn the charm quarks may have already shed most of their knowledge about the earlier stages of the collision and follow the bulk in its expansion. Hence a thermal description seems to be a good starting point.
Thermal lattice simulations obviously cannot tell us about possible cold nuclear matter effects that may affect the production of heavy mesons. Recent measurements by the ALICE Collaboration of the nuclear modification factor R pP b for D mesons in proton-lead collisions at the LHC have however shown that the strong suppression of heavy quarks they observed at transverse momenta above p T > 2GeV is not due to cold nuclear matter effects, but to the strong coupling of the charm quarks to the QGP medium [23,24].
What we can measure on the lattice are meson correlation functions, which are distinct from the actual particle correlation functions determined in experiment. They provide us with insight on global in-medium modification if one considers e.g. the appropriate ratios of vacuum correlators with those at T > 0. In order to extract information about the in-medium modification of individual meson states, e.g. the J/Ψ or D * , the spectral function encoded in the correlator needs to be computed via an unfolding procedure, which, as will be discussed, presents a formidable numerical challenge. Once we have access to the spectra, we can learn about medium induced mass shifts and thermal broadening, which also for heavy flavors can be translated to changes in the production yields and thus to measurable observables (see e.g. [25,26]).
The study of quarkonium spectral functions in a thermal medium is an active field of research with several different strategies present on the market. On the one hand there are efforts that have successfully used nonrelativistic effective field theories to extract bottomonium spectral information in-medium [27][28][29][30][31][32]. Recently also charmonium at finite temperature has been considered in this fashion [33,34]. As a plus, the non-relativistic setting allows one to access more relevant datapoints on the lattice, since the corresponding correlators are not symmetric in Euclidean time. At the same time the effective field theory approach introduces additional uncertainties, related e.g. to the (non-perturbative) matching to QCD and the absence of a naive continuum limit. Nonrelativistic quarkonium spectra have also recently been computed based on new results on the complex lattice QCD heavy-quark potential [25,26,36].
On the other hand there are several groups that aim at computing quarkonium in a fully relativistic framework. Those who focus on beauty require extremely finely spaced lattices, rendering the computation feasible only in the quenched approximation [37,38], which has also been used recently to study charmonium in [39]. The charm mass however is already light enough that it is possible to consider charmonium also in a full QCD medium, as has been done previously in Refs. [42,43].
In contrast to quarkonium, so far there have been only very few lattice studies carried out on in-medium open charm mesons [44][45][46], as they require a relativistic formulation at least for the light quark partner. Previous works used spatial correlators and cumulants to gain insight into the survival and melting patterns of D mesons around the transition temperature T c . In contrast to spatial correlation functions, which are related to the spectral function of the system via a double integral transform, the temporal correlators connect via a single integral transform. Thus spectral information manifests itself more directly on the latter than the former. In the present study we will, for the first time, present results for the temporal correlators and spectral functions of open charm mesons. Preliminary results have been presented in [47,48]. Let us also note an alternative QCD based approach to thermal in-medium D-meson properties, the sum-rule approach [49][50][51].
In the following section we briefly discuss the lattice QCD setup (Sec. I A) and introduce the Bayesian approach to spectral function reconstruction in (Sec. I B). Our presentation of numerical results start at T = 0 in Sec. II A. The in-medium correlator and spectral analysis for charmonium is contained in Sec. II B and for open charm in Sec. II C. We close with a summary and conclusion in Sec. III.

A. Lattice Setup
In this study we deploy lattices generated by the FASTSUM collaboration to describe the QCD medium in which the open and hidden heavy flavor mesons are immersed. These second generation ensembles [31,52] feature 2+1 flavors of anisotropic clover fermions and a mean-field improved anisotropic Szymanzik gauge action. With an anisotropy parameter of ξ = a s /a τ = 3.5 the spatial lattice spacing is a s = 0.123 fm. While the strange quark mass is tuned to its physical value, the pion mass takes on a value of m π ≈ 380 MeV. The temper- ature is changed in the fixed scale manner, by varying the number of points in the Euclidean time direction. The available temperatures and the corresponding number of configurations are listed in table I. The configurations were sampled every 10 HMC trajectories. The integrated Polyakov loop autocorrelation time for the ensembles above T c was about 20 trajectories, while the plaquette autocorrelation time below T c was found to be about 30 trajectories. The pseudocritical temperature was determined from the inflection point of the Polyakov loop to be T c = 185(4) MeV, for more details see Ref. [52].
Our finite temperature simulations utilize the same bare action parameters as those by the Hadron Spectrum Collaboration [53], who also kindly provided the zero temperature ensemble used for calibration at N τ = 128. The relativistic charm quark action and its parameters are taken over from the vacuum study [54,55] and both the configurations and the correlators were generated using the Chroma software package [56].
The Euclidean correlators G(τ ; T ) we compute in the lattice simulation are related to the spectral function ρ(ω; T ) via the integral relation where Information about the hadronic states, including their energies, widths and thresholds in the medium, is encoded in the peak and continuum structure of ρ. Our strategy to extract spectral functions from Euclidean time correlators via Bayesian inference is detailed in the next subsection. We note that the correlators have not been renormalised, and we can therefore only determine the shape of the spectral function, and not the absolute value, which would be required to determine for example transport coefficients.
We may already learn about the overall in-medium modification by considering the changes induced by thermal effects in the correlators themselves. Comparing correlators at different temperatures in the relativistic theory requires us to disentangle two effects. From Eq. (1) we see that besides the spectral function, also the kernel carries a temperature dependence. Hence we cannot simply truncate the low temperature correlator and divide it with the one at higher temperature, as is done in nonrelativistic settings, where the kernel does not contain a temperature dependence.
Instead we should consider the so-called reconstructed correlator where T r denotes a reference temperature at which the spectral function can be reliably constructed and which is usually chosen to be the lowest available temperature. Then we can obtain the Euclidean correlator that ensues if the low temperature spectral function were present at the higher temperature. I.e. in the absence of inmedium modification we have ρ(ω; T ) = ρ(ω; T r ), and hence G(τ ; T ) = G r (τ ; T, T r ). The converse is not necessarily true, since the loss of information occurring in Eq. (1) may lead to a situation where multiple changes in the underlying spectral function compensate each other in the Euclidean correlator. We will use the reconstructed correlator in the following to form direct ratios with the actual finite temperature correlators and as a benchmark for the systematic uncertainties of the spectral reconstruction algorithm. It is important in this regard to note that the reconstructed correlator can actually be computed directly from the underlying correlator G(τ ; T r ) without the need for a spectral reconstruction [41]. Based on the following identity for hyperbolic functions one finds

B. Bayesian Spectral reconstruction
The extraction of spectral functions from (1) represents an ill-posed inverse problem. The data obtained from lattice simulations are stochastic estimates of the correlator G(τ ) with finite precision ∆G/G, evaluated at N τ discrete points G(τ j = a t j) ≡ G j . For a numerical application we also need to discretize the spectrum using N ω points, in anticipation of the many different structure we may find: narrow bound state peaks, wide resonances, continuous structures above the threshold and even a transport peak at low frequencies. This necessitates the use of N ω N τ , so that attempting a naive χ 2 fit of the parameters ρ l to G i leads to an infinite number of solutions reproducing the simulation data within their errors.
In this study we deploy Bayesian inference [57,58] to give meaning to the ill-posed inversion task. This well established branch of mathematics allows us to systematically incorporate additional, so called prior information about the spectrum into the reconstruction task. This in turn regularizes the naive χ 2 fit and leads to a unique Bayesian answer. The starting point is Bayes' theorem, which states that the probability for a test function ρ to be the correct spectrum given simulation data and prior information is proportional to two ρ-dependent terms.
The first, P [G|ρ, I] = exp[−L], is called the likelihood and encodes how the data are generated. In the case of stochastically sampled data, such as in lattice QCD, the likelihood may be expressed in terms of the quadratic distance between the simulation data G i and the Euclidean correlator G ρ i arising from the current test function ρ inserted in (7) Here C ij denotes the usual covariance matrix. Note that in Monte Carlo simulations on the lattice individual correlators may not be perfectly uncorrelated. This leads to residual autocorrelations, which reduces the effective number of available configurations N cfg . A χ 2 fit wishes to determine the maximum of the likelihood function, which however contains many degenerate extrema. The second term in the numerator, the prior probability P [ρ|I(m)] = exp[αS[ρ, m]], provides the necessary regularization of the likelihood in the Bayesian approach. It is conventionally formulated with a hyperparameter α which weights the influence of data and prior information and which has to be determined self-consistently. Prior information enters into the functional S in two distinct ways. On the one hand the form of S itself favors certain spectra. On the other hand S depends on the so called default model m, which by definition is the correct spectrum in the absence of data. I.e., m acts as the unique extremum of the regulator S.
In the end one carries out a numerical optimization on the posterior P [ρ|G, I] to find the most probable spectrum given data and prior information Via the competition between L and S the reconstructed spectrum will be partially fixed by data and partially constrained by prior information.
In this study we will use two different implementations of the Bayesian strategy: the well known Maximum Entropy Method [59][60][61][62], as well as the more recent BR method [63]. While both are Bayesian in nature they differ in the deployed prior functional, in the way how the hyperparameter α is handled and in the implementation of the optimization problem of Eq. (10). Note that two Bayesian approaches in general will yield different outcomes as long ∆G/G and N t is finite. Only in the "Bayesian continuum limit" will both methods agree.
Based on arguments from two-dimensional image reconstruction the MEM proposes to use the Shannon-Jaynes entropy, as regulator. Following the state of the art implementation by Bryan [64] we work with a restricted search space in order to reduce the computational burden of the optimization in Eq. (10). While Bryan uses the N τ dimensional SVD basis of the Kernel we here deploy a different set of basis functions, the so called Fourier basis, which has been shown to provide improved frequency independent resolution [65]. The hyperparameter α is treated in the following fashion: One repeats the reconstruction based on S SJ for many different values of α and then computes an approximation of the probability P [α|ρ, D, I], which is closely related to the evidence probability P [D|I]. The final spectrum is the average over all individual reconstructions weighted with P [α|ρ, D, I].
The second method we deploy is the BR method, which has been developed with one-dimensional reconstruction problems in mind. Combining besides positive definiteness a smoothness criterion and the requirement of independence of the units used, the BR method regulator reads The absence of the factor ρ in front of the logarithm allows one to analytically determine the α-dependent normalization of the corresponding prior probability P [ρ|I].
In turn it becomes possible to marginalize the hyperparameter α a priori by assuming full ignorance of its values, i.e., P [α] = 1: The optimization of Eq. (10) is then implemented with the integrated posterior P [ρ|G, I(m)].
The numerical implementation for this study uses the MPFR library to implement arbitrary precision datatypes, which in practice are set to 512 bits of precision. Due to the presence of the ρlog[ρ] term in S SJ the MEM converges relatively slowly compared to the BR method. Therefore in the latter, as is common practice, we accept an extremum if Eq. (10) is satisfied to an accuracy of ∆ MEM = 5 × 10 −8 , while for the BR method we use ∆ BR = 10 −60 .
The robustness of the outcome of a Bayesian reconstruction depends on both the precision of the input data G i and the reliability of the prior information I. To quantify the effect of the former we carry out a ten-bin Jackknife, where the reconstruction is repeated, each time with a different subset of simulation datasets removed. The dependence on the prior information can be assessed by repeating the reconstruction with different default models. As a proxy for both effects Bayesian methods provide an internal measure of robustness, which under certain assumptions is related to the curvature of the minimization functional Q = log[P [ρ|G, I]]. In previous studies we have seen that this internal measure, possibly due to its reliance on assumptions, often underestimates the actual uncertainty of the reconstruction, as determined from the Jackknife and default model variation.
In the following we will use an equidistant frequency discretization of N MEM ω = 1000 and N BR ω = 3000 for the spectral reconstructions in the interval ω num ∈ [0, 5]. The reason why we can use a smaller number of points for the MEM is related to the fact that for a given number of datapoints the number of Fourier basis functions is fixed. As long as all features of the basis function are resolved, the result becomes virtually independent of the actual N MEM ω . The Euclidean times are taken to be τ = {0, . . . , N τ }. Since the naive finite temperature kernel diverges at the origin we instead deploy the regularized variant K(ω, τ ) = cosh[ω(τ − β/2)]/ cosh[ωβ/2], so that raw output of the reconstruction is ρ/ arctan[ωβ/2], which we convert to ρ/ω afterward. This also means that the default model provided is the default model for ρ/ arctan[ωβ/2]. We will show reconstructions for the choice m(ω) = 1 but vary the default model both in amplitude m 0 = {0.1, 1} and form m(ω) = m 0 (1 + ω) 2 to assess the residual influence of the default model.

II. NUMERICAL RESULTS
Before embarking on an investigation of the properties of open and hidden charm mesons at finite temperature, let us inspect first the situation at zero temperature, as obtained on the FASTSUM ensembles. Not only will the obtained vacuum results allow us to judge the accuracy of our computations but also provides us with the necessary baseline for interpreting the Bayesian spectral reconstructions at finite temperature.
Charmed states @ T≃0 Charmed states @ T≃0 We show the vacuum correlation functions at N τ = 128 of all charmed states investigated in this study in the top panel of Fig.1. The anticipated differences in the ground state masses manifest themselves visibly in the different slopes at intermediate Euclidean times. To access the value of the mass quantitatively we compute the so called effective mass using a ten-bin Jackknife resampling. From the middle panel of Fig.1 we see that it asymptotes to the ground state mass at late times. The higher the ground state mass, the lower the signal to noise ratio, i.e. in particular the scalar and axial-vector quarkonium channel still suffer from significant uncertainty. Fitting the late time behavior reveals the ground states masses listed in Tab.II. As our study aims at investigating the in-medium spectral structure of charmed mesons, where information beyond masses, i.e. thermal widths, are of interest, we do not deploy distillation or smearing techniques for our correlators. Thus it is not surprising that a multiexponential fit did not result in stable results for the masses of the first excited states, which usually require a more refined variational approach.
We continue by carrying out first Bayesian spectral reconstructions on the T = 0 correlators using both the BR method and the MEM. Results are shown in Fig. 2. Errorbands in the plots arise from the combined variance obtained via a ten-bin Jackknife (statistical uncertainty), as well as from varying the underlying default model (systematic uncertainty). At the lowest temperature a relatively large number of datapoints is available, i.e. here the statistical error dominates the uncertainty budget compared to the default model dependence.
As shown in Fig. 2 both methods are able to locate the position of the lowest lying structure well from the currently available data. As expected from the use of naive interpolation operators, the overlap with the excited states is not strong enough for an accurate determination of the next higher lying structure. Even though the BR method in general does better in reproducing isolated ground state peaks, with the current data quality neither method has a clear advantage in reproducing excited state peaks. For all states except P-wave charmonium a second bump at higher frequencies is hinted at in the reconstructions, which however is not significant.
Note that except for the small signal-to-noise case of the P-waves the BR method produces a sharper ground state peak than the MEM. Furthermore the variance of the ground state structure is smaller in the BR method.
Interestingly the MEM appears to show a more pronounced second structure than the BR method, which however also carries a larger uncertainty band.
In preparation for the in-medium investigation we need to elucidate how increasing the temperature, i.e. the diminishing of available Euclidean datapoints, affects the spectral reconstruction. To this end we wish to take the lowest temperature result and encode it in a correlator which corresponds to a lattice at higher temperatures. In studies of non-relativistic mesons using the effective field theory NRQCD this is straightforwardly achieved by simply truncating the T = 0 dataset to the number of points available at T > 0. In a relativistic setting where the correlation function shows periodicity related to the inverse temperature a simple truncation is not appropriate and instead we have to turn to the reconstructed correlator, which takes into account the explicit temperature dependence of the integral kernel in Eq. (1).
We first prepare reconstructed correlators of several different temporal extents by applying Eq. (6) to the low temperature correlator at N τ = 128. Since that equation is only well defined for an Euclidean extent which divides N τ = 128 without remainder, G rec is computed here for N τ = 64, N τ = 32 and N τ = 16 to be used in the Bayesian reconstruction. For the subsequent comparison with finite temperature correlators, we have resorted to padding the N τ = 128 correlators with zeros to compute also extents such as N τ = 40.
In Fig. 3 several representative examples of the outcome of the Bayesian analysis using the reconstructed correlator are shown. In the top row the BR method has been deployed, while in the bottom row the MEM was used. While quantitatively different, similar qualitative trends emerge. One finds that a change from N τ = 128 to N τ = 64 only induces a small broadening and a minor shift of the ground state structure to higher frequencies. Going to even smaller N τ = 32 however already induces significant modifications to the low lying peak with ∆m = 0.2 − 0.5GeV. At N τ = 16 the well defined ground state peak cannot be resolved anymore with the currently available statistics. In general the MEM produces more strongly washed out features than the BR method here. We have to keep in mind that the actual spectrum encoded in all the different reconstructed correlators is exactly the same and the differences in the outcome are simply a manifestation of the degradation of the resolution of the Bayesian method with limited available data. Now that we understand the inherent method deficiency related to limiting the number of available datapoints, we can proceed to inspecting genuine finite temperature data.  B. Quarkonium at finite temperature

Correlation functions
We first take a look at what information we can glean from the thermal lattice correlators themselves. Figure 4 shows the ratio of the thermal correlator G(τ ; T ) to the reconstructed correlator G r (τ ; T, T r ) defined in Eq. (6), with T r = 0.24T c corresponding to N τ = 128. In the absence of any thermal modifications, this ratio equals 1. We show data for the pseudoscalar (top), vector (second to top), scalar (second to bottom) and axial vector (bottom) channels. As is expected, the results for the axialvector channel are virtually the same as for the scalar channel. Note that in the scalar channel, we have only accessed the correlation function up to T = 1.52GeV.
In the pseudoscalar channel, we find that the thermal modifications of the correlator are small, 2%, at all temperatures. With the current level of precision no significant in-medium modification can be observed up to T c . However, for T > T c , we find that the correlator ratio begins to deviate systematically from unity, and at the highest temperature (1.9T c ), the deviation is genuinely significant. Interestingly there is a qualitative jump between the data at T = 1.52T c and T = 1.90T c , the former following the monotonic upward bending above 1, similar to lower temperatures, while the latter lies fully below unity. A rather significant modification of the spectrum must occur in between these two temperatures.
Comparing this result to the same quantity evaluated in [41], we find an interesting difference. While our ratios show a clear intermediate upward bend between T c < T < 1.5T c , Ref. [41] depicts ratios which decrease monotonically within errors. At around 2T c , once we observe the jump in our results the Euclidean time dependence again agrees qualitatively. We note that the reconstructed correlator ratios determined from first generation FASTSUM ensembles [66] exhibit the same qualitative features as those of [41], suggesting that the difference is largely due to the limited statistics in the present study.
Let us turn to the vector channel. Similarly to the pseudoscalar case, we again see a clear difference between the data below and above T c , with G/G r remaining roughly consistent with 1 below T c and deviating systematically from 1 above T c . The latter deviation manifests itself in a monotonic upward bending, consistent with what has been observed in previous lattice studies using both relativistic [41] and non-relativistic formulations [33,34]. The maximum change in the ratio here lies at 5% at T = 1.90T c and contrary to the pseudoscalar channel, no sudden changes in the correlator occur. The strength of the modification in our case is slightly smaller than what was observed in previous studies.
The P-wave charmonium states behave rather differently compared to the pseudoscalar and vector case, A strong modification is present, already below T c , and may reach up to a factor of 3.5 at T = 1.90T c . The change again takes the form of an upward bending, which increases in strength with temperature. The qualitative  behavior is consistent with findings of previous studies, both based on relativistic and non-relativistic heavy quarks. The strength of the modification is similar to those observed in the relativistic case but is much more pronounced than in the non-relativistic study, where it remains below 2% for T ≈ 2T c . One possible explanation is related to the appearance of zero modes in the correlator [35], which are not resolved in the effective theory approach.
The observed upward bending behavior in the correlator ratio has been interpreted in the context of nonrelativistic potential computations to be consistent with the disappearance of excited state peaks and the steady reduction of the threshold down to lower energies, as temperature rises. Extracting these intricate changes in the spectrum directly via Bayesian inference is extremely challenging, since the relevant structures are located relatively close to each other compared to the maximum lattice momentum. In addition, the small overall change in the ratios of the pseudoscalar and vector channel are only well resolved for N τ /4 < τ < N τ /2 rendering the number of effective relevant datapoints comparatively small.

S-wave (J/Ψ and ηc) spectra
We begin our investigation of in-medium spectral properties with the J/Ψ and η c particle. In Fig. 5 we show their spectral functions, as obtained by the BR method and in Fig. 6 by MEM with the Fourier basis, each time at eight different temperatures. Errorbands are again taken as the combined variance arising from statistical and default model dependence. With a diminishing number of datapoints at increasing temperatures the default model dependence starts to become the dominating source of error at the T = 1.52T C and T = 1.90T C . The gray dashed lines, inserted as reference guide, denote the positions of the vacuum ground and first excited state peak, as obtained from a variational computation by the Hadron Spectrum Collaboration [54].
A first naive inspection of the BR method result by eye tells us that with increasing temperature the former ground state structures monotonically move to higher frequencies until above T ∼ 1.5T c no discernible structure can be found. The strength of the peak reduces concurrently. The MEM shows qualitatively similar results for the spectra at T T c but no well pronounced peaks appear above T c . Such differences between the methods have been observed before in the study of bottomonium in lattice NRQCD [31,32]. On the one hand the resolution of the MEM reduces with the smaller number of available datapoints due to a smaller number of admitted basis functions. On the other hand the BR method is known to be susceptible to ringing artifacts, which may lead to more pronounced peak structures than actually present in the data.
In light of these different methods artifacts, how can we disentangle them from genuine in-medium effects? It is here that the reconstructed correlators again play an important role. We already investigated what happens to the reconstruction of the T ≈ 0 correlator if it is consistently restricted to a smaller number of Euclidean time steps and found that shifts to higher frequencies and broadening ensued. Thus in Figs. 7 and 8 we compare the reconstructed in-medium spectra at three temperatures corresponding to N τ = 128, N τ = 32 and N τ = 16 with  the reconstruction obtained from the corresponding reconstructed correlator. Fig. 7 depicts the outcome based on the BR method, Fig. 8 based on the MEM. The comparison is illuminating, as it reveals that the shift and broadening observed at T = 0.95T c , contrary to what a simple inspection by eye might suggest, is fully consistent with the vacuum spectrum reconstructed from the same number of N τ = 32 datapoints. Only at T = 1.9T c or equivalently N τ = 16, does the in-medium spectrum show a more strongly washed out behavior than the reconstruction from G rec . In the vector channel these differences are already significant, going beyond the combined statistical and systematic errorbars. In the pseudovector channel the difference is hinted at but is not yet significant.

P-wave (χc0 and χc1) spectra
We now turn to the P-wave charmonium spectra. Figure 9 and Figure 10 show the scalar (χ c0 ) and axial-vector (χ c1 ) spectral functions obtained by the BR method and by the Fourier basis MEM respectively.
A first inspection by eye suggests that for both channels already at T = 0.76T c there occurs a substantial shift and weakening of the ground state structure. From T ≈ T c on, no sign of the ground state peak remains. Interestingly in the χ c1 channel, it is the MEM result that shows slightly more peaked ground state features below T c but also with larger uncertainty. To ascertain whether this first impression holds true, let us carry out the comparison to the spectra from the reconstructed correlator again, shown in Fig. 11. As was the case with the S-waves, the restricted temporal range leads to an upward shift and gradual weakening of the ground state peak. However, at T = 0.95T c (N τ = 32), the ground state peak is clearly present in the reconstructed correlator spectral function, while the inmedium spectral function obtained from the thermal correlators shows no sign of such a peak. Within the current data quality we conclude that the 1P state disappears from the spectrum below or near T c . Since the MEM comparison for the P-wave carries very large errorbands, we do not gain any further insight from it and omit it from our discussion here.
On the one hand the P-wave channels are those with the smallest signal to noise ratio, leading to sizable errorbands. On the other hand the magnitude of the inmedium modification already observed in the correlator ratio is highly significant and can be therefore be picked up by the Bayesian reconstruction. We are currently increasing the statistics on the P-wave correlators to bring them to a similar level of relative precision to those of the S-wave to make sure that the relatively early disappearance of the ground state peak is not simply due to a lower precision.

C. Open Charm at finite temperature
After considering the charmonium states, we now turn to the study of the open charm mesons. Here we consider the pseudoscalar (D, D s ) and vector (D * , D * s ) channels, as these are the most relevant for phenomenology.

Correlatior ratios
Up to this point open heavy flavor mesons at finite temperature were solely investigated by use of spatial correlation functions or bulk observables. Here we present for the first time the actual Euclidean time correlator ratios, which have a direct connection with the in-medium spectral function.
In figure 12 we show the correlators G(τ ) divided by the reconstructed correlators G r (τ ) at the same temperature, for all four channels. We observe that the correlator at T = 0.76T c is still consistent with no modifications in all channels, while some nontrivial trough structure is already hinted at at intermediate τ .
Starting from T = 0.84T c genuine departures from unity emerge. The maximum changes we find are around 10%, which is stronger than what is observed in the corresponding charmonium channels, consistent with a lower binding energy.
The shape of the changes is also distinct from that observed in the charmonium sector. Instead of an upward bending that increases with temperature, a deeper and deeper trough emerges with a minimum at around  11: Comparison of the BR reconstructed cc scalar (top) and axial-vector (bottom) channel T > 0 spectra with those based on the T ≈ 0 reconstructed correlator at the same number of datapoints. All reconstructions are obtained using τ ∈ [2, Nτ /2 − 1] and the same statistics. Already at 0.95Tc a significant in-medium modification is observed, and no discernible peak structures are found.  τ ≈ 0.35fm. One possible origin of this difference may lie in the absence of a densely populated regime of excited states, which cannot melt as the temperature increases.
Since the potential picture is not readily applicable for open-heavy flavor, a similarly intuitive explanation as for the modification of charmonium is not at hand. We note that the pseudoscalar D meson ratio behavior differs clearly from that of pseudoscalar charmonium discussed in the previous section but is reminiscent of what has been observed for the η c in [41] and in the first generation FASTSUM ensembles. Up to around T c the trough minimum seems to remain almost at the same position in all channels. In the D meson case, the qualitative behavior remains the same also above T c and the trough deepens. For the D * mesons however, a sharp upward rise sets in at τ /a > N τ /4.
Let us take over some intuition from the potential computations of quarkonium, where the upward bending arises once the lower lying peak structures in the spectrum start to be affected by in-medium effects. While the ground state masses of D and D * mesons differ by less than 100MeV, (the latter being heavier than the former) the distinct patterns in the ratio suggest that their stability against in-medium effects differs significantly. D * mesons appear much more sensitive to thermal fluctuations than D mesons, which in vivo should lead to to a stronger suppression.
Several heavy-ion experiments, i.e. STAR [67], CMS [68] and ALICE [69][70][71] have conducted measurements of the nuclear suppression factor of open charm mesons. While the former two collaborations have presented mainly D results, ALICE has determined both D and D * production yields separately. The currently available data quality from Run 1 and Run 2 at the LHC has allowed them to determine R AA at intermediate p T and a large number of participants corresponding to the 10 − 20% centrality class. It is in this regime where the charm quarks have the highest probability to become equilibrated with the medium and thus we may attempt to connect to our fully thermal results.
In Run 1 no indications have been found that the D or D * species show different suppression patterns. There is a slight tendency of D 0 to show a bit stronger suppression than the D * even though the difference is not significant. The preliminary Run 2 data on the other hand show a tendency for the D * meson to be more strongly suppressed at low p T compared to the D's, but again due to the relatively large errorbars the effect is not significant yet. A possible discrepancy between our observation of a significantly different behavior of the D and D * correlator ratios and the absence of such differences in the observed suppression needs to be further understood.
On the theory side our interpretation of the upward bending of the ratios is based on a naive intuition borrowed from potential model computations of charmonium and surelys needs refinement. A direct determination of the spectral properties of course would be most illuminating, which is what we will attempt in the following section.
On the experimental side, the Run 2 data are still being investigated and it will be very interesting to see whether the tendency of the low p T suppression of D * will eventually prevail in the final analysis.
Note also that except for the highest temperature (T = FIG. 13: Spectral functions for the lc (left) and sc (right) pseudoscalar channels, obtained using the BR method. Also shown are the ground and first excited state energies determined from a variational analysis [55]. 1.9T c ), the modifications in the strange-charm sector are smaller than those in the light-charm sector, which is consistent with the the hypothesis that D s yields may be increased relative to D yields. Experimentally this hypothesis is supported by the preliminary ALICE Run 2 analysis, which clearly showed that the R AA for D s is larger than the average R AA for D and D * .

Spectral functions
We saw that the open-charm correlator ratios behave very different to their charmonium counterparts. In order to give this theory observation a physical interpretation we need to translate it into a statement about the spectral properties of open charm mesons. As a first step we may ask whether the in-medium changes of the correlator are an indication of an equally different behavior of the ground state spectral structure? To this end we carry out the corresponding Bayesian spectral reconstructions.
In Fig. 13 and Fig. 14 the results based on the BR method and the MEM are plotted respectively. The gray dashed lines are the vacuum ground and first excited state masses obtained from a variational approach in Ref. [55].
Both methods show that below T c , the D mesons exhibit consistently more pronounced structures, compared to their D * cousins. The naive inspection by eye again finds a broadening and shifting of peaks with temperature. The BR method exhibits remnant peak structures up to T ≈ 1.5T c , while no sign of any remnant structure appears at T = 1.90T c . The MEM on the other hand shows overall more washed out structures, so that at T > T c one is hard pressed to identify a genuine peak. For both methods we also see that above T > T c the strength of the D meson structures is slightly stronger than those of D * , which is consistent with the observation of stronger deviations from unity in the D * correlator ratio.
The mandatory comparison with the spectra coming from the corresponding reconstructed correlators, shown in Figs. 15, 16, tells us that in agreement with the correlator ratios, up to T c no significant modification of the ground state occurs, while at 1.9T c neither (u, d) + c nor s + c channels show any sign of a ground state remnant. The currently available data quality however does not yet allow us to distinguish whether also in terms of the ground state spectral peak D * mesons are more susceptible to in-medium effects than the D's. Within our current accuracy, no significant difference is found between the spectral functions of D and D s mesons.

III. SUMMARY AND CONCLUSION
We have presented a combined analysis of charmonium and open-charm meson properties in a thermal medium based on fully relativistic lattice discretized heavy quarks. On the second generation FASTSUM ensembles we scrutinized the ratios of the in-medium correlators to the so called reconstructed correlators, as well as the corresponding in-medium spectral functions obtained from both the BR and the Maximum Entropy methods. Method artifacts related to the diminishing of the available number of correlator points in Euclidean time at higher temperatures were crosschecked by reconstructing at the same time the low temperature spectrum from the reconstructed correlator itself.
For the charmonium S-wave particles η c and J/Ψ we found that the correlator ratios show a qualitatively and quantitatively very similar behavior to that observed in studies deploying non-relativistic QCD. At temperatures above T c , where the in-medium modification is statistically significant, an upward bending occurs with strength increasing with temperature. One exception is the η c channel at T = 1.9T c , where compared to T = 1.5T c a significant change appears, the ratio moves from bending above unity to below unity. Interestingly, only at this high temperature is the ratio consistent with what has been observed in previous studies of relativistic charmonium, where the upward bending in the η c channel was absent. The other difference to previous relativistic charmonium studies lies in the strength of the vector channel deviations from unity, which in our case are much milder than those observed previously.
The χ c0 and χ c1 channel ratios show a much stronger in-medium modification. Even though the absolute errors in these channels are much larger than for the Swaves, the deviation from unity is easily discernible. The strong upward bend is consistent with previous findings in relativistic studies but much larger than the corresponding one observed in non-relativistic QCD.
The Bayesian spectral reconstructions show qualitatively similar results: peaks seem to shift in frequency and broaden with increasing temperature. In general, the BR method shows stronger peaked features than the MEM, which is understood to arise from the restricted search space of the latter and the susceptibility of the former to numerical ringing. To distill genuine in-medium effects, we compare the in-medium spectra to the ones obtained from the reconstructed correlator. Here also the BR method and MEM give consistent results, showing that around T = T c no in-medium modification is visible in the S-wave, while significant weakening of the peak structure occurs for the P-waves. At T = 1.9T c also in the S-wave channel a weakening of the peak structure can be observed.
For the D and D * mesons we carry out the same analysis starting with the correlators. This is the first time that instead of spatial correlation functions, the actual Euclidean time correlators are considered, which are related to the in-medium spectral function via a single integral. Already below T c significant in-medium modifications are observed, which instead of manifesting themselves in an upward bend, take the form of a deeper and deeper trough. The absence of a densely populated excited states regime is proposed as reason for the difference to the charmonium ratios. Above T c the D and D * mesons start to show strong differences. The former continue monotonically to deepen the trough, while the latter show a sudden rise of the ratios, which hints at the D * ground state peak starting to be affected earlier than for the D.
The spectral reconstructions, while being consistent with the behavior observed in the correlator ratios below The lc channels are on the left, and the sc channels on the right. All reconstructions are obtained using τ ∈ [2, Nτ /2 − 1]. At 0.95Tc no significant in-medium modification is observed, while at T = 1.9Tc there is no evidence of any surviving bound state.  T c , are not yet precise enough to distinguish any differences between the modification of the D and D * ground state. In all considered channels there is no significant ground state modification visible at T c , while at 1.9T c no more structure appears to remain. In terms of phenomenological relevance the open charm results provide two main insights: We find that even though the D and D * mesons only differ by around 100MeV in vacuum binding energy, the latter show a much stronger in-medium modification above T c in the correlator ratio than the former. Taking intuition from potential model computations of charmonium, such behavior would translate into stronger suppression for the D * . On the other hand none of the currently available open charm data are able to resolve any differences between the different D and D * nuclear modification factors. The newest preliminary results from Run 2 by ALICE may contain a (not yet significant) hint towards stronger D * suppression at very low p T though. Further experimental efforts into reducing the uncertainty of the measured R AA in P b + P b collisions are thus very welcome.
The second insight is related to the systematically smaller medium modifications observed in the D s compared to the D correlators. Consistent with intuition, the system with the heavier of the light quarks appears more strongly bound and thus more stable. In turn the suppression of the D s in a purely thermal setting is expected to be weaker than that for the D. In the new Run 2 results from ALICE this ordering of suppression among different heavy-flavor mesons with and without s-quarks is clearly observed in the corresponding R AA .
In order to investigate in more detail the different behavior of individual in-medium D and D * meson states above T c , hinted at by their correlator ratios, we will need to significantly improve the robustness of the spectral reconstruction. Before the arrival of the third generation of FASTSUM ensembles, which will feature twice the number of Euclidean datapoints, we are currently working on increasing the statistics of the meson correlator mea-surements on the second generation ensembles, which at the same time will help to further constrain validity of a disappearance of the P-wave charmonium ground state peak at around T c .
This study has been carried out at a single lattice spacing and with relatively heavy up and down quarks. The quark mass effects will be investigated in the near future when a new ensemble with m π ≈ 230 MeV becomes available. This is not expected to have a substantial impact on the charmonium sector, but may allow us to determine with more certainty whether there is a difference in the behaviour of D and D s mesons. The relatively large lattice spatial lattice spacing together with the use of a tree-level spatial clover coefficient was found in [54,55] to result in an underestimate of the S-wave hyperfine splitting in both charmonium and open-charm systems of 20-40 MeV, which is beyond the precision of the current study. No significant volume dependence was found at zero temperature in these studies, but finite volume effects will in the future be studied directly by including results also for a larger volume of ∼ (4 fm) 3 . In addition to this, it is planned to improve on the current study by halving the temporal lattice spacing, hence reducing lattice artefacts while doubling the number of available data points.