Meson Screening Masses in (2+1)-Flavor QCD

We present lattice QCD results for mesonic screening masses in the temperature range 140 MeV $\lesssim T \lesssim$ 2500 MeV. Our calculations were carried out using (2+1)-flavors of the Highly Improved Staggered Quark (HISQ) action, with a physical value for the strange quark mass and two values of the light quark mass corresponding to pion masses of 160 MeV and 140 MeV. Continuum-extrapolated results were obtained using calculations with a variety of lattice spacings corresponding to temporal lattice extents $N_\tau = 6 - 16$. We discuss the implications of these results for the effective restoration of various symmetries in the high temperature phase of QCD, as well as the approach toward the perturbative limit.


I. INTRODUCTION
At high temperatures the properties of stronginteraction matter change from being controlled by hadronic degrees of freedom to deconfined quarks and gluons. While the thermodynamics in the low temperature phase of QCD resembles many features of a hadron resonance gas, with hadrons keeping their vacuum masses, this quickly changes at temperatures close to and above the crossover transition to the high temperature phase. In fact, the zero temperature hadronic degrees of freedom seem to provide a quite satisfactory description of thermal conditions close to the transition to the high temperature phase [1], although there is evidence of thermal modification of the spectrum [2]. At high temperature, however, quarks and gluons deconfine, which also is reflected in properties of hadron correlation functions and the thermal masses extracted from them (see e.g. [3]). Resonance peaks in spectral functions, which enter the integral representations of thermal hadron correlation functions, broaden and shift with temperature [4]. In spatial correlation functions [5] the finite temporal extents, 0 ≤ τ ≤ 1/T , of the Euclidean lattice acts on spatial quark and anti-quark propagators like a finite volume effect, which influences the long-distance behavior of these correlation functions. Their exponential decay at large distances defines screening masses, which differ substantially from the pole masses at zero temperature, and approach multiples of πT at high temperature, * Deceased. which is characteristic for the propagation of free quark quasi-particles in a thermal medium.
The chiral cross-over separating the low and high temperature regimes for non-vanishing quark masses is characterized by a smooth but rapid change of the chiral condensate around T pc . The pseudo-critical temperature T pc , for the physical value of the ratio of light and strange quark masses, has recently been determined from fluctuations of various chiral observables : T pc = (156.5 ± 1.5) MeV [6].
Despite a small explicit breaking of the chiral symmetry by the residual light quark masses, the chiral symmetry, which is spontaneously broken in the hadronic phase, gets effectively restored above T pc . The deconfinement of the light quark and gluon degrees of freedom is believed to be strongly related to the drop of the chiral condensate and the resultant effective restoration of the chiral symmetry. If chiral symmetry is restored then the excitations of the plasma are also expected to carry those informations in spatial hadron correlators. In fact, the analysis of spatial hadron correlation functions and their asymptotic large distance behavior [5] is found to be a sensitive tool for studies of different patterns of chiral symmetry restoration at high temperature. Generally it is found in calculations at physical values of the quark masses that the temperature dependence of screening masses, differs significantly in quantum number channels sensitive to the restoration of the SU L (2) × SU R (2) chiral flavor symmetry and the anomalous axial U A (1) symmetry, respectively. While the former will be restored completely at chiral transition temperature in the chiral limit, the latter remains broken also at high temperature by the Adler-Bell-Jackiw anomaly [7][8][9]. However, with the thermal suppression of non-perturbative breaking effects, which at zero temperature arise, for instance, from the presence of topologically non-trivial gauge field configurations [10], the anomalous axial symmetry may be "effectively restored". It has been argued that the question whether or not the chiral symmetry and anomalous axial symmetry get effectively restored at the same temperature may have significant qualitative consequences for the structure of the QCD phase diagram in the chiral limit [11].
Calculations with staggered fermions [12,13] show evidence for U A (1) symmetry breaking also above T pc and provide evidence for the close relation between axial symmetry breaking and the density of near-zero eigenmodes [14]. However, to what extent the flavor singlet anomalous axial U A (1) symmetry gets effectively restored at the chiral phase transition temperature, T 0 c = 132 +3 −6 MeV [15], which defines the onset of a true phase transition in the chiral limit, is still an open question [16][17][18][19].
Several recent lattice QCD calculations performed in 2 and (2+1)-flavor QCD with physical quark mass values utilizing overlap and Möbius domain wall [20][21][22][23][24][25] as well as Wilson [26] fermions observe an effective restoration of the U A (1) symmetry at temperatures above the pseudocritical temperature T pc , i.e. at about (1.2−1.3)T pc . This is in accordance with earlier findings in calculations of screening masses with staggered fermions, where effective U A (1) restoration has been observed through the degeneracy of scalar and pseudo-scalar correlation functions and screening masses at temperatures T 1.3T pc [12].
One of the motivations of this study is to also determine the extent to which U A (1) is effectively restored at the chiral crossover temperature through screening masses for which we have performed continuum extrapolation not yet performed in earlier studies. At the level of screening correlators, U A (1) restoration will lead to a degeneracy between the scalar (S) and pseudoscalar (P S) correlators, while chiral symmetry restoration yields a degeneracy between the vector (V ) and axial vector (AV ) correlators. We calculate mesonic correlation functions numerically using (2+1)-flavor lattice QCD for all the possible flavor combinations including light and strange quarks, namely, light-light (ūd), light-strange (ūs) and strange-strange (ss). Within each flavor combination, we determine scalar, pseudoscalar, vector and axial vector ground sate screening masses. The temperature dependence of this set of meson correlation functions has been analyzed before [12], including also charmonia [27], on coarse lattices using the p4 discretization scheme for staggered fermions. With this calculation we substantially improve over earlier work by using the Highly Improved Staggered Quark (HISQ) action with physical values for the light and strange quark masses and by performing calculations in a wide range of lattice spacings, 0.017 fm ≤ a ≤ 0.234 fm that allows us to perform controlled extrapolations to the continuum limit in the temperature range 140 MeV ≤ T ≤ 974 MeV. Albeit not continuum extrapolated, we extend the calculation of screening masses to temperatures as large as 2.5 GeV. Results for screening masses for charmonia, open strange-charm as well as forss channels, with the HISQ action but for only a single lattice spacing corresponding to N τ = 12, have been reported before [28]. This paper is organized as follows: In the next section, we briefly review properties of spatial meson correlation functions and their evaluation using the staggered fermion discretization scheme. We describe the staggered fermion set-up for our calculations in Sec. III. We then present our results in Sec. IV where we start with updating our scale setting in Sec. IV A and present some zero-temperature meson masses. Staggered fermion specific cut-off effects, so-called taste splittings, for T = 0 are shown in Sec. IV B. We present results for temperatures around the chiral crossover regime in Sec. IV C where we also discuss effective U A (1) restoration. In Sec. IV D, we present our results for the screening masses at high temperatures compared to chiral crossover temperature and compare these with predictions from resummed thermal perturbation theory. Finally we state our conclusions in Sec. V. For completeness we have appendices where we start with an update of the parametrization for scale setting in Appendix A and then in Appendix B and Appendix C, we summarize our statistics and tabulate the continuum-extrapolated values of the screening masses, respectively.

II. SPATIAL CORRELATORS AND SCREENING MASSES
Properties of the hadron spectrum at zero and non-zero temperature are commonly determined from an analysis of two-point correlation functions M Γ (x)M Γ (y) , where the operators M Γ project on to a specific set of quantum numbers and x, y are Euclidean space-time coordinates. At zero temperature the lowest excitation (mass) in a given quantum number channel is conveniently extracted from the asymptotic large Euclidean time behavior of the correlation function. At finite temperature, the calculation of correlators separated in Euclidean time is limited by the limited extent of this direction that determines the inverse temperature of the system, β = 1/T . In contrast there are no such restrictions for spatially separated correlators, also known as screening correlators.
In QCD, the finite temperature meson screening correlators, projected onto zero transverse momentum (p ⊥ ≡ (p x , p y ) = 0) and lowest Matsubara frequency of a bosonic state (p 0 ≡ ω 0 = 0), are defined by where M Γ ≡ψΓψ is a meson operator that projects onto a quantum number channel Γ selected by Γ = Γ D ⊗ t a with Dirac matrices Γ D and a flavor matrix t a . The an- gular brackets · · · , denote the expectation value over the gauge field ensemble. The correlators decay exponentially for large z, which defines the corresponding screening mass m Γ (T ). As already mentioned, for T → 0, the screening masses tend to the mass of the T = 0 meson with the same quantum numbers. For T → ∞, they approach the common value m Γ = 2πT irrespective of the spin and flavor [5], which indicates that the dominant excitations consist of two almost free fermionic excitations (quarks) which each have a lowest Matsubara frequency (energy) ω 0 = πT . For non-zero T , the relation between screening mass and pole mass could be highly non-trivial due to the emergence of non-analytic structures in the spectral function [29]. On the lattice, the continuum Dirac action must be replaced by a suitable discrete variant. Staggered fermions, which we use in this work, are described by onecomponent spinors rather than the usual four-component spinors. Because of this, they are relatively inexpensive to simulate. However the price to be paid is that the relation to the continuum theory is subtle. The continuum limit of the theory is the Dirac theory of four fermions rather than one. As a result, each meson too comes in sixteen degenerate copies which are known as tastes and the corresponding operators are of the form is the 16-component hypercubic spinor and Γ D and Γ T are Dirac matrices in spin and taste space respectively. Although different tastes are degenerate in the continuum, on the lattice this degeneracy is broken by gluonic interactions. The masses of the taste partners can be determined from the decay of correlation functions of staggered meson operators M(x) = n,n φ(n, n )χ(x + n)χ(x + n ), where x is the hypercube co-ordinate and n and n point to the various vertices of the unit hypercube and φ is a sitedependent phase factor whose form depends on the spin and taste quantum numbers of the meson [30][31][32].
In this work, we will only consider local operators, i.e., operators with n = n . In Table I we list the eight local staggered meson operators that were studied in this work and their mapping to the familiar mesons of QCD. We note that the operators M3, M4 and M5 (respectively M6, M7 and M8) refer to the x, y and τ components of the same axial vector (respectively vector) meson. In the spatial correlation functions the meson operators were separated along the z direction. One thus may average over the M3 and M4 (respectively M6 and M7) components in order to improve the signal. Note however, that unlike at T = 0, at finite temperature one cannot average over all three transverse directions due to absence of Lorentz invariance in the definition of the correlators [33]. In the vector and axial vector channels we thus deal with two distinct correlation functions and resulting screening masses, denoted as transverse and longitudinal.
A typical staggered meson correlator, for a fixed separation (in lattice unit) between source and sink, is an oscillating correlator that simultaneously couples to two sets of mesons with the same spin but with opposite parities: where n σ = z/a denotes the spatial separation of the source and sink operators M φ . For large enough distances the correlator of Eq. 3 may be constrained to a single term, i.e. i = 0. In Eq. 3 we also replaced the large distance exponential fall-off given in Eq. 2 by a hyperbolic cosine which arises due to the periodic nature of correlators on lattices with finite spatial extent N σ .

A. Data sets
We calculated the six distinct mesonic correlators, constructed from local staggered fermion operators introduced in the previous subsection, numerically using (2+1)-flavor gauge field ensembles generated with the HISQ action and a Symanzik improved gauge action. The HISQ action [34][35][36] is known to have the least amount of taste-splitting [37], due to which it has been used in several precision studies both at T = 0 as well as at finite temperature [35,[37][38][39][40]. The gauge ensembles for β ≤ 7.825, have been generated by HotQCD collaboration and previously had been used to study the QCD equation of state of strongly interacting matter [41,42].
For β > 7.825, we have used the gauge ensembles from TUMQCD collaboration, generated for the study of the expectation values of the Polyakov loop and its correlators [43,44]. Gauge configurations have been generated on lattices of size N 3 σ × N τ , where N τ = 6, 8, 10, 12 and 16, and N σ = 4N τ . Most of the data for these five different values of the temporal lattice size, corresponding to five different values of the lattice spacing a at fixed value of the temperature T = 1/(N τ a), have been collected in a temperature range 140 MeV ≤ T ≤ 172 MeV using physical values of the light (m l ) and strange (m s ) quark masses, i.e. a quark mass ratio 1/27. On lattices with temporal extent N τ = 8, 10 and 12 we also used data sets obtained with a slightly larger quark mass ratio, 1/20. These data sets cover a larger temperature range up to about 2.5 GeV. The Goldstone pion masses for these two quark mass ratios are 140 MeV for m l /m s = 1/27 and 160 MeV for m l /m s = 1/20.
All the above-mentioned gauge configurations used in this analysis have been generated with a strange quark mass tuned to its physical value by tuning the mass of the ηs s meson, M ηss = 686 MeV. This value is based on leading order chiral perturbation theory relation, M ηss = 2m 2 K − m 2 π , between the ηs s , π and K masses. Once the strange quark mass had been determined, the light quark mass was set to either m l = m s /27 or m l = m s /20, as already discussed. The former choice of quark mass were used for temperatures below and near the chiral cross-over temperature, T pc , while the higher quark mass was used at higher temperatures (T 172 MeV) where quark mass effects are negligible. The tuning of the strange quark masse, which leads to our line of constant physics, is also discussed in detail in Ref. [41]. All our simulation parameters and the number of gauge field configurations analyzed are summarized in Appendix B.
The conversion of hadron masses, calculated in lattice units, into physical units as well as the determination of our temperature scale requires the calculation of one physical observable that is used for the scale setting. For this purpose we use the kaon decay constant, f K = 156.1/ √ 2 MeV, also used in other thermodynamics studies with the HISQ action. We give the updated The purpose of the new calibration of the parametrization of f K a(β) in Appendix A is to improve on the scale at the larger β-values in this study. Note that when compared to the previous scale [40,41], this leads to a small ∼ 1% decrease of the lattice spacing at the largest β-values compared to the previous scale determination [40,41], while the differences are negligible for β 7.0.

B. Hadron correlation functions
A general meson correlator M(x)M(y) consists of quark line connected and disconnected parts. In this work we only focus on flavor non-singlet mesonic correlators which do not have disconnected contribution.  The analysis of chiral symmetry restoration, including the U A (1) restoration, can be performed using flavor nonsinglet correlators alone [21,45]. The (fictitious) ηs s meson, whose mass was used to fix the bare quark masses, also does not receive any contributions from disconnected diagrams [28]. We generally had to retain up to 2-3 states in Eq. 3. Such multi-state fits present a challenge as a straightforward fit is often highly unstable. For this purpose we developed a routine to guess the initial parameters directly from the data [46] for different terms of the sum in Eq. 3. We also developed [46] a fit parameter estimation routine that works directly on the oscillating correlators. This method relies on the fact that the mass of the oscillating and non-oscillating part are usually roughly of similar size and thus assumes their equality in the first step: 1. At a small fit interval [n σ,min : n σ,max = N σ /2], perform one state fits on all even points of the correlator and we call the resulting fit parameters say A even φ,0 and m even φ,0 . Repeat the same for the odd points (A odd φ,0 , m odd φ,0 ).

2.
Assuming similar size of the non-oscillating and oscillating mass, the fit parameters for the combined fit may be estimated with 3. Using the parameters from step 2 as initial guess, perform a full one state fit with oscillating and nonoscillating part.
4. Increase the fit interval. Guess the mass of the next excited state of either the even or the odd part (we used m . Adjust the corresponding amplitude (A − φ,1 or A + φ,1 ) such that the first point of the correlator in the fit interval is reproduced. 5. Perform a full fit with higher states. Use the parameters from steps 3 and 4 as initial guess.
6. Repeat steps 4 to 5 until the desired number of states is reached.
Having developed a method to perform automated multiple state fits, we still have to find which set of fit parameters is the most reasonable one for a given fit interval. For that purpose we have used the corrected Akaike information criterion (AICc) [47,48]: For each fit interval we have performed different multiple state fits (maximum 3 states for oscillating correlators and maximum 4 states for non-oscillating correlators) and selected the one which has the smallest AICc. In Fig. 1 a comparison between the different multiple state fits and the result that is selected by the AICc is shown. In contrast to the one state fit, this results in an early onset of a stable plateau. After the fits have been performed the plateaus were manually selected for each correlator. The final value for the screening mass and its uncertainty are determined by Gaussian bootstrapping. More technical details about the automated fitting procedure can be found in Ref. [46].
We calculated screening correlation functions using point as well as corner wall sources. The point source is the simplest type of source that one can use to calculate mesonic screening functions and we have used one source for each color. However it does not suppress the excited states, therefore isolating the ground state can be difficult unless the states are well-separated or the lattice extent is large. The use of extended (smeared) sources can often help to suppress excited state contributions, allowing to extract the ground state mass and amplitude even on smaller lattices. Here we have used a corner wall source which means putting an unit source at the origin of each 2 3 cube on a chosen (in our case) z-slice [49][50][51]. In Fig. 2, we show a comparison of a mass calculation using point and corner wall sources at two different temperatures. As discussed earlier, in both cases we found that it is necessary to take into account contributions from higher excited states to obtain reliable results for the ground state screening masses. In Fig. 2 we have only shown the fit results where we have taken one state for both oscillating and non-oscillating part of the correlator (denoted by self-explanatory notation '(1,1)') and the AICc selected plateaus for the corresponding fit interval. We found that the use of a corner wall source provided advantages only for the noisy correlators, which in particular are the vector and axial vector channels at low temperatures. In the bottom panel of Fig. 2  Numbers appearing in parenthesis corresponds to number of states taken in the fit for non-oscillating and oscillating part of the correlator. The method of using AICc selection criterion to find a plateau among various fits has been described in the main text.
better signal as compared to a point source and one gets a longer plateau with smaller uncertainty when the minimum distance for the fit, n σ,min , is varied. Therefore, we used the corner wall source only where it was necessary, i.e. for the vector and axial vector channels below T = 300 MeV. For all the other cases however, we found that higher state fits for the point source worked just as well and that their results agreed with the corner wall results. We also found that in the case of a corner wall source, the excited state often has a negative amplitude and therefore, the influence of the higher states is to shift the result for the screening mass downward rather than upward as can be seen from the top panel of Fig 2. IV. RESULTS

A. Scale setting and line of constant physics
As the scale setting calculations as well as the determination of the line of constant physics had been performed prior to our current screening mass analysis we tried to re-confirm the scales used in our calculation through additional zero temperature calculations performed on lattices of size 64 4 . We performed calculations at three values of the gauge coupling, β = 7.01, 7.13 and 7.188. Using the parametrization of f K a(β) given in Appendix A this corresponds to lattice spacings a = 0.085 fm, 0.076 fm and 0.072 fm, respectively. The strange quark mass has been fixed using m s a(β) from Ref. [41] and the light quark mass was taken to be m l = m s /27. The resulting zero-temperature meson spectrum is shown in Fig. 3. The solid horizontal lines in the figures correspond to the experimentally determined values of the respective masses [52]. The slight mismatch for M ηss (m K ), arising from the slight mistuning of the strange quark mass, is visible in right (middle) panel of Fig. 3. We note that results for most of the P S, V and AV mesons agree well with the physical zero temperature spectrum within errors. The scalar meson, in theūd sector however, seems to have twice the pseudoscalar mass rather than the true scalar mass forūd sector. This is a well-known staggered artifact [53][54][55] and we will also discuss its effect for non-zero temperatures in Sec. IV C.
However, such definite trend is absent in heavierūs and ss sectors. A slight mismatch can also be observed for the AV masses inūd sector with no definite trend with decreasing lattice spacing.

B. Taste splittings at T = 0
Although use of staggered quarks will lead to tastesplitting in every hadronic channels, its effects are particularly severe in the pseudoscalar sector (π, K and ηs s ), since these are the lightest states in the theory. In Fig. 4, we plot the masses of the sixteen different tastes of each of the three pseudoscalar mesons, i.e., π, K and ηs s , for three different values of the lattice spacing. The correlators for the different taste partners are constructed using non-local operators [32] with Γ D = γ 5 and various Γ T , as shown in Fig. 4. In each case, the lightest meson is the meson with the quantum numbers Γ T = Γ D = γ 5 . This meson is the only Goldstone boson that is massless in the chiral limit at finite lattice spacing and the masses of the other fifteen mesons approach its mass in the continuum limit. The masses of the other partners have been normalized to the mass of the corresponding Goldstone boson for that particular lattice spacing. Our results extend the previous HISQ results for taste-splitting to smaller lattice spacings. A more detailed discussion on the taste splitting effects, also in comparison to other staggered discretization schemes can be found in [38,40].
One can define the root mean square (RMS) pion mass m P S RM S , as a measure of the taste splitting [56]: The γ-matrix suffixes in Eq. 4 refer to the taste structure of the mesons. The RMS mass approaches the Goldstone mass in the continuum limit; hence its deviation from the Goldstone mass at a given lattice spacing is a way of quantifying the taste-breaking effects. The sixteen tastes group into different multiplets, in a way understood from staggered chiral perturbation theory [56]. This is the reason for the factors of 3 in Eq. (4). We find that the RMS taste splitting is of the order of 15-25 % for the light-light(ūd) sector but decreases to about 4-8 % for the strange-strange(ss) sector. We also see that this splitting decreases as the lattice spacing decreases, as expected.
Lastly we note that the masses plotted here are consistent with the trend observed in Fig. 2 of Ref. [40] where the taste-splitting was calculated, with the same action but for coarser lattices and a slightly heavier quark mass.

C. Screening masses around the cross-over region
We now present our results for screening masses calculated in a range of temperatures going from just below the chiral cross-over temperature, T pc = 156.5(1.5) MeV, to about 2T pc , namely 140 MeV ≤ T ≤ 300 MeV. This temperature range is important both from the phenomenological point of view as well as regarding the restoration of chiral SU A (2) and axial U A (1) symmetries. As already mentioned earlier, our screening masses were calculated at two values of the light quark mass viz. m l = m s /27 for T 172 MeV and m l = m s /20 for all higher temperatures. It is worth to mention here that we have also calculated screening masses with m l = m s /20 for T 172 MeV but we do not show them here because we have fewer statistics compared to that for m l = m s /27. For higher temperatures, the quark mass dependence is negligible and the heavier quark mass can be used without affecting any of the conclusions. Horizontal lines correspond to the physical values of the masses [52]. The scalar meson mass is 2mπ instead of ma0 (or mπ +mη) due to a staggered artifact at finite lattice spacing. This discrepancy vanishes when the continuum limit of the correlator would be taken before calculating the screening mass [53,54] (see Sec. IV C).    Using the fitting procedure described in Sec. III, we calculated screening masses for five different values of the lattice spacings corresponding to N τ = 6, 8, 10, 12 and 16, which allow for a continuum extrapolation. As the temperatures do not agree among the different lattices, the screening masses have to be interpolated between the different temperature values. In our extrapolation method, the interpolation and the extrapolation are performed in one single fit: For the interpolation we use simple splines. Then, the extrapolation is performed by replacing the spline coefficients by a function linear in 1/N 2 τ and performing a joint fit, that includes all the data. The spline's knot positions are placed according to the density of data points. The knots are positioned in such a way that the same number of data points lies between each pair of subsequent knots. This means in particular that more knots are used at the low temperature region, where the interpolation is more curvy. To stabilize the spline, we use some of its coefficients, to constrain the spline's derivative with respect to T at some points. These constraints are placed far outside of the actual region where the extrapolation is performed [46]. The error band are computed using Gaussian bootstrapping and by performing the extrapolation on each sample. Final values and errors are calculated using median and 68%-percentiles of the bootstrap distribution. In Fig. 6 we show two examples of continuum extrapolations following the above-mentioned procedure in the P S and   S sector for a limited temperature range. More technical details of the continuum extrapolations can be found in Ref. [46]. Continuum extrapolated masses of all four channels for all three flavor combinations have been tab-ulated in Appendix. C.
We plot the screening masses for 140 MeV ≤ T ≤ 300 MeV, for the different flavor sectors and for all lattice spacings, in Fig. 5. The mesons with angular momentum J = 0 (S and P S) were easier to determine, especially for lower temperatures, as compared to the J = 1 mesons (V and AV ). We find some cut-off dependence in the scalar sector, especially for smaller N τ . For the other sectors, the cut-off dependence was indistinguishable within the statistical error. We perform the continuum limit for all the sectors, using data from five different values of the cut-off corresponding to our five different values of the temporal lattice extent, mentioned earlier. The resulting continuum extrapolated bands are plotted in Fig. 7. In Fig. 5 and Fig. 7 we also show the pseudo-critical temperature region as a gray vertical band. The massless infinite temperature limit m free scr = 2πT is shown as a dashed line in each of the plots.
For T T pc the screening masses are expected to approach the mass of the lightest zero temperature meson with the same quantum numbers, e.g., theūd pseudoscalar screening mass should approach the pion mass m π . We see that this behavior is readily realized for the P S and V sectors. Already for T 0.9T pc the corresponding zero temperature masses are approached in thē ud,ūs andss sectors to better than 10%. Although the zero temperature limits are not yet reached in the AV channel at this temperature, we see clear indications for a rapid approach to the corresponding zero temperature masses for all combinations of heavy and light quarks. These values are in all cases approached from below, i.e., at the pseudo-critical temperature the AV screening masses are smaller than the corresponding zero temperature masses. In thess sector the screening mass of the f 1 -meson is about 15% lower than the f 1 -mass around T pc and reduces to about 7% already at T 0.9T pc . The situation is similar in theūs sector. However, thermal effects are substantially larger in theūd sector. Here we find that the screening mass of the a 1 mesons at T pc differs by about 35% from the corresponding zero temper-  ature mass and the two masses still differ by about 20% at T 0.9T pc . Note that also from our calculations for m l = m s /20, where we have results at even lower temperatures, we found that the screening masses go towards corresponding zero-temperature masses steeply. Similar behavior was also found in calculations with staggered fermions utilizing the p4 discretization scheme [12].
The situation is far more complicated in the S sector for finite lattice spacings. In nature, the lightest flavored scalar meson is either the a 0 (980) or the a 0 (1450). Rather than either of these values, as can be seen from the left panel of Fig. 5 and Fig. 7, the scalar screening mass approaches the value 2m π instead. The reason for this is that for staggered fermions, the scalar can decay into two pions at finite lattice spacing [53]. This decay is forbidden in nature due to parity, isospin and G-parity (I G ) conservation. The unphysical behavior in the staggered discretization comes from the contribution of the different tastes in the intermediate states of loop diagrams. If one takes the continuum limit for the correlator before calculating the screening mass, then the contribution from different tastes cancels out and the physical behavior is recovered [53][54][55]. Since we, however, calculate the screening masses first and then take the continuum limit, we obtain the unphysical ππ state rather than the true scalar ground state or the physically allowed πη decay. The unphysical decay only occurs for mesons with isospin I = 1. For theūs case (I = 1/2), the decay to Kπ actually occurs in nature. In Figs. 5 and 7, we see that the scalar screening mass indeed tends to m π + m K as T → 0.
As the cross-over temperature is approached, the vector and axial vector screening masses should become equal due to effective restoration of chiral symmetry. At T = 0, the axial vector meson a 1 is about 500 MeV heavier than the vector meson ρ. As the temperature is increased, the AV screening mass decreases while the V mass increases slightly until the two masses become degenerate right at the pseudo-critical temperature (left panel of Fig. 7). In contrast, in theūs andss sectors, AV and V masses become equal at higher temperatures compared to T pc . Moreover, the relative change of AV masses w.r.t. V masses from low temperature towards degeneracy temperature progressively decreases when one goes fromūd toss. It must be noted that the approach is nevertheless smoother compared to previous results that were obtained using the p4 discretization scheme for staggered fermions [12]. cross-over temperature, as noted from Fig. 7, is quite similar to what has been seen in the calculation of nucleon masses, where the mass of one particular parity (the one with higher zero-temperature mass) of nucleon changes a lot and comes close to its parity partner, which on the contrary, hardly changes from low temperature towards chiral cross-over temperature [2,57,58].
In Fig. 7, we also see that the scalar and pseudoscalar screening masses in theūd sector become degenerate around T ∼ 200 MeV. Unfortunately, one cannot immediately draw any conclusions about an effective U A (1) restoration from this due to the pathology of theūd scalar correlator that we have discussed above. However, as we have already mentioned, the unphysical contribution cancels out if one would take the continuum limit for the correlator first. Moreover, as the pion screening mass increases around the cross-over region while the continuum scalar screening mass is expected to decrease around T pc before rising again at higher temperatures, this unphysical decay channel might be closed around T pc due to lack of phase space. Therefore the degeneracy of the screening masses in the S and P S channel around T ∼ 200 MeV is an indication towards an effective restoration of the U A (1).
Despite the above argument, we may nevertheless try and estimate the effective U A (1) restoration temperature directly from the correlators. Although it is difficult to calculate the continuum limit of staggered correlators due to their oscillating behavior, one may instead consider the corresponding susceptibility, which is given by the integrated correlator, and calculate its continuum limit instead. The staggered π and δ susceptibilities are de-  fined as We plot our results, along with the continuum extrapolations, for the difference of the scalar and pseudoscalar susceptibilities for theūd sector in Fig. 8. In order to be able to take the continuum limit, we have renormalized the quantity with m 2 s . We have also normalized these numbers to T 4 pc . For reference, we also show the pseudocritical temperature region by a gray band in the figure. For theūd sector we see that the difference is non-zero around the pseudo-critical temperature and only goes to zero for T ∼ 200 MeV. There are some theoretical arguments in favor of effective U A (1) restoration at the chiral phase transition [59] in the chiral limit. On the other hand lattice calculations, performed away from the chiral limit, have found evidence in favor of this scenario [14,20,23,60].
Before moving on, we note that the behavior of the screening masses and susceptibilities in theūs andss sectors is qualitatively the same although the degeneracies discussed above occur at progressively higher temperatures [46]. This mass ordering of degeneracy temperatures is in complete accordance with what has been observed for even heavier sectors [28], although one has to keep in mind that the mass effects in the susceptibilities for heavier sectors are expected to be much larger than the U A (1) breaking effects due to quantum fluctuations.

D. Screening masses at high temperatures
In the previous subsection we have seen that the temperature dependence of the screening masses at T > 250 MeV qualitatively follows the free theory expectations, namely the screening masses are proportional to the temperature, with proportionality constant not very different from 2π. Furthermore, the AV and V screening masses are close to the free theory expectations, while the P S and S screening masses are 10-20% smaller. In this subsection we will study the screening masses at higher temperature with the aim to see how the degeneracy of P S(S) and AV (V ) screening masses expected in the infinite temperature limit sets it. We would like to see if contacts to the weak coupling calculations can be made at high temperatures.
Although attempts have been made [5,58,[61][62][63][64][65][66] to compare screening masses from lattice QCD to those from weak coupling calculations, it is not clear in which temperature range weak coupling results can be reliable. For this reason it is important to perform lattice calculations at as high temperatures as possible. Therefore, we extended the calculations of the meson screening masses to T = 1 GeV using four lattice spacings corresponding to N τ = 6, 8, 10 and 12, and performed the continuum extrapolations. The results are shown in Fig. 9. We find that the lattice spacing dependence is very small for T > 300 MeV, and within errors the N τ = 8 results agree with the continuum extrapolated values. Therefore, for 1 GeV < T < 2.5 GeV we calculated the screening masses using only N τ = 8 lattices. The results of these calculations are also shown in 9. We clearly see from the figure that the AV and V screening masses overshoot the free theory value around T = 400 MeV and are almost constant in temperature units. The P S and S screening masses overshoot the free theory expectation only at temperature larger than 1 GeV and remain smaller than the AV and V screening masses up to the highest temperature considered.
The behavior of the screening masses in the weak coupling picture beyond the free theory limit can be understood in terms of dimensionally reduced effective field theory, called electro-static QCD (EQCD) [67]. This approach turned out to be useful for understanding the lattice on the quark number susceptibilities [68,69], the expectation value of Polyakov loop [43] and the Polyakov loop correlators [44]. It is interesting to see if deviation of the screening masses at high temperature from 2πT can be understood within this framework.
In EQCD the correction to the free theory value for the screening masses is obtained by solving the Schrödinger equation in two spatial dimensions with appropriately defined potential [70][71][72]. At leading order the potential is proportional to the coupling constant of EQCD, g 2 E [72], which in turn can be expressed in terms of the QCD coupling constant g 2 = 4πα s . At leading order g 2 E = g 2 T , and g 2 E has been calculated to 2-loops [73]. Moreover, at leading order the potential and the correction to the free theory value are independent of the spin, i.e. the P S(S) and AV (V ) screening masses receive the same correction, that has been calculated in Ref. [72]. This correction is positive in qualitative agreement with our lattice results. In Fig. 9 we show the corresponding weak coupling result from EQCD. We used the 2-loop result for g 2 E and the optimal choice for the renormalization scale µ/T = 9.08 [73]. We varied the scale µ by factor of two around this optimal value to estimate the perturbative uncertainty, which turned out to be very small (the uncertainty corresponds to the width of the weak coupling curve in Fig. 9. We see that the weak coupling results from EQCD are slightly larger than the AV (V ) screening masses and significantly larger the the lattice results for P S(S) screening masses. This is not completely surprising because the EQCD coupling constant g 2 E is not small except for very high temperatures and thus higher order corrections may be important. Beyond O(g 2 E ) the correction will be spin dependent [70,71]. Since the coupling constant decreases logarithmically the screening masses will approach 2πT only for temperatures many orders of magnitude larger than those considered here, when the AV (V ) and P S(S) screening masses become degenerate. It would be interesting to calculate the O(g 4 E ) correction to meson screening masses and see whether EQCD predictions work quantitatively.

V. CONCLUSIONS
We have performed an in-depth analysis of mesonic screening masses in (2+1)-flavor QCD with physical (degenerate) light and strange quark masses. In the vicinity of the pseudo-critical temperature for chiral symmetry restoration, T pc and up to about 1 GeV we could perform controlled continuum extrapolations, using input from five different values of the lattice cut-off. Comparing screening masses for chiral partners, related through the chiral SU L (2) × SU R (2) and the axial U A (1) trans-formations, respectively, we find in the case of light-light mesons evidence for the degeneracy of screening masses related through the chiral SU L (2) × SU R (2) at or very close to T pc while screening masses related through an axial U A (1) transformation start becoming degenerate only at about 1.3T pc . In particular, the V and AV mesons (J = 1), which are related by chiral SU L (2) × SU R (2) transformations, become degenerate at T T pc , while the S and the P S (J = 0) mesons, which are related by axial U A (1) transformations, only become degenerate around 1.3T pc . The onset of these degeneracies also occurs in the light-strange and strange-strange meson sectors, but at higher temperatures.
At high temperatures the screening masses overshoot the free theory expectations in qualitative agreement with the weak coupling calculations at O(g 2 E ). While mesonic screening masses in given angular momentum (J) channels become degenerate, screening masses in channels with different J, e.g. J = 0 and J = 1, stay well separated even up to the highest temperature, T = 2.5 GeV, that was analyzed by us. We argued that it is necessary to go beyond O(g 2 E ) calculations in order to understand this feature within the EQCD framework. This non-degeneracy has also been observed in Ref. [74], where it was also shown that these two sets of mesons only become degenerate at asymptotically high temperatures. This conclusion is in agreement with the results that we have presented in this paper in Sec. IV D (Fig. 9). For the scale setting in this project we used the Kaon decay constant, i.e. f K a(β). Including the measurements up to β = 7.373, listed in Ref. [41], we have updated the parametrization used in Ref. [40] : In Fig. 10 we have compared the fit described with Eq. (A1) to the same from Ref. [40]. It can be seen from the plot that one overestimates f K a(β) with the old parametrization for β 6.9 by ∼ 1%. One can look in Ref. [40,41] for more details on this kind of parametrization. Here we summarize our data sets and the number of configurations on which point and wall source correlators have been calculated are given in the last two columns of the tables which are labled 'point'and 'wall', respectively. Here we have tabulated the continuum extrapolated screening masses of P S, S, V and AV channels and in each channel for all three flavor combinations i.e.ūd,ūs andss.  (7) 6.29 (7) 6.02(8) 6.31(6) 1.000 6.3(1) 6.6(1) 6.3(2) 6.60(8)  (7) 6.0(2) 6.31(6) 1.000 6.3(1) 6.5(2) 6.3(2) 6.59 (9) Table XII. Continuum-extrapolated values of the strangestrange screening masses.