Pseudo-scalar meson spectral properties in the chiral crossover region of QCD

Determining the type of excitations that can exist in a thermal medium is key to understanding how hadronic matter behaves at extreme temperatures. In this work we study this question for pseudo-scalar mesons comprised of light-strange and strange-strange quarks, analysing how their low-energy spectral properties are modified as one passes through the high-temperature chiral crossover region between $T=145.6 \, \text{MeV}$ and $172.3 \, \text{MeV}$. We utilise the non-perturbative constraints satisfied by correlation functions at finite temperature in order to extract the low-energy meson spectral function contributions from spatial correlator lattice data in $N_{f}=2+1$ flavour QCD. The robustness of these contributions are tested by comparing their predictions with data for the corresponding temporal correlator at different momentum values. We find that around the pseudo-critical temperature $T_{\text{pc}}$ the data in both the light-strange and strange-strange channels is consistent with the presence of a distinct stable particle-like ground state component, a so-called thermoparticle excitation. As the temperature increases this excitation undergoes collisional broadening, and this is qualitatively the same in both channels. These findings suggest that pseudo-scalar mesons in QCD have a bound-state-like structure at low energies within the chiral crossover region which is still strongly influenced by the vacuum states of the theory.


Introduction
The nature of excitations in a thermal medium is an open question that underpins many physical phenomena.In order to address this question it is essential to determine the non-perturbative structure of correlation functions in quantum field theories (QFTs) at finite temperature, since these encode the in-medium dynamics.Significant advancements in the understanding of stronglyinteracting phenomena at finite temperatures have been achieved by using lattice quantum chromodynamics (QCD) simulations to compute the form of correlation functions, and to infer dynamical features based on the properties of these quantities [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15].Lattice simulations are performed in imaginary time τ , and the correlation functions of particular interest for the study of spectral properties at finite temperature are the two-point functions of gauge-invariant operators O Γ (τ, ⃗ x) where Γ denotes a specific set of quantum numbers, and the expectation value is with respect to the thermal background state at temperature T = 1/β.Assuming that the background state is in thermal equilibrium, it follows from taking the spatial Fourier transform of Eq. (1.1) that [16,17] where ρ Γ (ω, ⃗ p) is the Fourier transform of the thermal commutator of O Γ (τ, ⃗ x).The spectral function ρ Γ (ω, ⃗ p) is a quantity of particular importance as it contains information about all possible spectral excitations with quantum number Γ that can occur in the thermal ground state.Equation (1.2) therefore implies that imaginary-time temporal correlator C Γ (τ, ⃗ p) data can in principle be used to extract spectral information.This pursuit has received significant focus in the literature but requires one to overcome an ill-posed inverse problem, namely that discrete temporal correlator data and the structure of Eq. (1.2) are not sufficient to uniquely reconstruct the form of ρ Γ (ω, ⃗ p).For this reason, all of the strategies for obtaining information about ρ Γ (ω, ⃗ p) require additional input, either from perturbative calculations or phenomenological modelling [18,19].
In the specific case ⃗ p = 0, the temporal correlator C Γ (τ, ⃗ p) coincides with the spatial integral over C Γ (τ, ⃗ x).If one instead fixes a spatial direction z, and integrates C Γ (τ, ⃗ x) over τ and the remaining directions x and y, one obtains the so-called spatial correlator C Γ (z).For hadronic systems, particular progress has been made in establishing their properties in recent years [20][21][22][23][24][25][26][27][28].At large distances z, it is expected that these correlators have an exponential behaviour where the coefficient m scr Γ is referred to as the screening mass.In the zero-temperature limit m scr Γ approaches the mass of the lowest-energy vacuum state created by O Γ , and so in this sense the screening mass provides a measure of the in-medium modifications experienced by this state.A significant advantage of spatial as opposed to temporal correlators is that they can be calculated on the lattice at arbitrarily large distances, whereas the temporal correlator extent is limited by the inverse temperature β due to its periodicity.It is therefore computationally easier to extract information about the thermal properties of a system from C Γ (z), particularly at large temperatures.As with the temporal correlator, C Γ (z) is also directly related to ρ Γ (ω, ⃗ p).In this case, one finds e ipzz ∞ 0 dω πω ρ Γ (ω, p x = p y = 0, p z ). (1.4) The fact that C Γ (z) and C Γ (τ, ⃗ p) are both governed by ρ Γ (ω, ⃗ p), but the dependence is significantly different, means that lattice simulations of these correlators provide distinct ways in which to probe the spectral structure of the theory.Therefore, any extraction of ρ Γ (ω, ⃗ p) from lattice correlator data must be consistent with both sets of data.
In Ref. [29] it was demonstrated that spectral function components can also be extracted from spatial correlator data.This approach was applied to N f = 2 flavour lattice QCD data for the light-quark pseudo-scalar meson correlator at different temperatures above the pseudo-critical temperature T pc , and used to extract the low-energy contributions to ρ PS (ω, ⃗ p).By computing the temporal correlator predictions from these components via Eq.(1.2), it was found that the corresponding temporal correlator data does indeed impose non-trivial constraints on these components.The goal of the present work is to generalise the approach developed in Ref. [29] to N f = 2 + 1 flavour lattice data for pseudo-scalar meson correlators involving both one and two strange-quark components.The remainder of this paper is structured as follows: in Sec. 2 we outline the theoretical foundations of the approach developed in Ref. [29], in Sec. 3 we discuss the lattice data analysis and physical implications of the results, and finally in Sec. 4 we summarise our findings.

Spatial correlators at finite temperature
As outlined in Sec. 1, spatial correlators C(z) contain important information about the properties of QFTs at finite temperature.In Ref. [29] it was shown that the spatial correlator actually possesses a general non-perturbative representation, the structure of which encodes the connection between the spectral degrees of freedom and the behaviour of C(z).In this section we will briefly outline the origin of this representation, and how it can be used to extract spectral information from spatial correlator data.

Causality constraints
In Refs.[30][31][32][33][34][35] the authors developed a non-perturbative framework for describing the structure of scalar correlation functions at finite temperature, generalising the rigorous vacuum-state QFT formulations that have led to numerous foundational insights over the last few decades [36][37][38].In particular, it was demonstrated that the constraint of causality, namely that the fields ϕ(x) satisfy the condition: [ϕ(x), ϕ(y)] = 0 for (x − y) 2 < 0, imposes significant constraints on the structure of thermal correlation functions.In the case of the spectral function ρ(ω, ⃗ p), this constraint implies that ρ(ω, ⃗ p) satisfies the general representation [30] ρ where the quantum number label Γ is dropped here to indicate that this holds for any type of scalar field ϕ(x).Equation (2.1) corresponds to the finite-temperature generalisation of the Källén-Lehmann spectral representation [39,40], a foundational result of vacuum-state QFT.A significant characteristic of Eq. (2.1) is that all of the dynamical and temperature-dependent effects are encoded in D β (⃗ u, s), the so-called thermal spectral density, and this controls the correlation between the energy and momentum of the system.Establishing the basic properties of D β (⃗ u, s) is therefore key to understanding in-medium phenomena at finite temperature.For the purposes of this analysis we are interested in the general structure of the spatial correlator C(z), as defined in Eq. (1.4).It turns out that Eq. (2.1) imposes various constraints on the structure of the corresponding Euclidean two-point function [41], and this implies that C(z) can be written [29] where

Spectral properties from spatial correlators
It is clear that establishing the properties of the thermal spectral density is key for understanding how in-medium effects manifest themselves in correlation functions.In the limit of vanishing temperature full Lorentz invariance is restored, which implies: or equivalently: D β (⃗ x, s) → ρ(s), where ρ(s) is the spectral density of the zero-temperature theory.
As is well known, the analytic properties of the two-point correlation functions are then fixed by ρ(s), the singularities of which reflect the existence of particle states in the spectrum of the theory.
In particular, the appearance of a singular δ(s − m 2 ) component implies the presence of an onshell state with mass m.At vanishing temperature, the spectral characteristics of the theory are directly related to the singular structure of ρ(ω, ⃗ p) in the variable p 2 = ω 2 − |⃗ p| 2 .However, as soon as the temperature becomes non-vanishing D β (⃗ u, s) will possess a non-trivial ⃗ u dependence, and so the connection between the structure of ρ(ω, ⃗ p) and the spectral properties of the theory is less clear cut.Establishing this connection requires one to understand the dispersive properties of the excitations of the thermal ground state.

Thermoparticles
In Ref. [30] it was suggested that the thermal medium may contain particle-like excitations.A characteristic feature of these excitations is that they are formed from the same discrete δ(s − m 2 ) component as mass m vacuum particle states, and their contribution to the position-space thermal spectral density is given by [34] where the second term D c,β (⃗ x, s) parametrises all other spectral contributions.These excitations were subsequently given the name thermoparticles [31] in order to draw a sharp distinction with other particle-like contributions considered in the literature such as collective quasi-particle modes, which vanish when T → 0. As outlined in Refs.[30][31][32][33][34] and summarised recently [42], there are several well-motivated reasons for why thermoparticle components D m,β (⃗ x) δ(s − m 2 ) provide a natural description for particle-like states at finite temperature.Perhaps the most compelling is that due to the representation in Eq. (2.1), the appearance of a non-trivial damping factor D m,β (⃗ x) causes the zero-temperature peak in ρ(ω, ⃗ p) at p 2 = m 2 to become broadened, since the momentum space factor D m,β (⃗ u) is no longer restricted to the point ⃗ u = 0. Physically, this broadening describes the effect of collisions with the medium, and has been shown to depend on the underlying dynamics of the theory [34], as one would expect.
Another significant feature of thermoparticles is that they describe intrinsically stable excitations which dominate the large-time evolution of the thermal state [34].The probability for their detection is therefore reduced only by there collisions with the medium.In the zero-temperature limit these collisions become increasingly suppressed, and one recovers the original vacuum-state contribution to the spectral function.The factorised structure of the thermoparticle representation also means that it can be directly generalised to describe intrinsically unstable excitations, by replacing δ(s − m 2 ) in Eq. (2.4) with a suitable resonance-type function, such as a relativistic Breit-Wigner.This description then enables one to distinguish between reductions in particle-detection probability that arise from two completely different physical origins: 1. Statistical effects -caused by collisions with the medium due to the non-vanishing temperature.These effect are parametrised by the damping factor D m,β (⃗ x).

2.
Mixing effects -brought about by decay processes that occur due to the intrinsic instability of underlying vacuum states.
We will not discuss further here the specific structure of thermoparticle contributions and their corresponding damping factors, only to mention that their properties have been explored in different models [34], and that more recently in Refs.[43] and [44] it was shown that these quantities can in fact be used to perform non-perturbative analytic calculations of different in-medium observables, including the shear viscosity.

Thermoparticle signatures
Now that we have a physically well-motivated description of how stable particle-like excitations could potentially appear at finite temperature, one can use this to look for their signatures in thermal correlation functions.In the case of the spectral function, it follows from Eq. (2.1) that the corresponding thermoparticle contribution ρ TP (ω, ⃗ p) always has a discrete energy threshold at the zero-temperature particle mass m, and thus one can write where the structure of ϱ depends on the specific properties of the damping factor.From a physical perspective, Eq. (2.5) implies that the thermal background must be supplied with an energy of at least m in order to have a non-vanishing probability of creating a thermoparticle state.The threshold in Eq. (2.5) is in fact a special case of a more general situation: Of course, given the decomposition in Eq. (2.4) the full spectral function would also contain additional components from D c,β (⃗ x, s), and this could give spectral contributions both above and below the thermoparticle threshold |ω| = m.The interplay between these different spectral contributions and their relative size depends on the dynamics of the specific theory.
For the spatial correlator C(z) it was shown in Ref. [29] that if thermoparticle components dictate the low-energy behaviour of the spectral function, just as stable particle states do at zero temperature, then these contributions will dominate the behaviour of the spatial correlator in the following manner and this domination will be especially pronounced at large distances.Taking the derivative of Eq. (2.6), one finds that the damping factor can then be extracted from the leading large-z behaviour of the spatial correlator derivative Given the form of the damping factor, one can then use the spectral representation in Eq. (2.1) to calculate the thermoparticle contribution to ρ(ω, ⃗ p).This procedure can also be generalised to situations in which there exists more than one thermoparticle-type state [29].
The discussions in this section demonstrate that because of the constraints imposed by causality, if thermoparticles exist and provide a dominant low-energy spectral function contribution, their complete analytic structure can be computed directly from the corresponding spatial correlator.In Ref. [29] this procedure was applied to QCD in order to explore the thermal properties of light meson states.This study focussed in particular on two-flavour lattice QCD data for the spatial pseudoscalar meson correlator, and found evidence for the existence of a pion ground state consistent with a thermoparticle-type behaviour at a temperature 220 MeV above the pseudo-critical temperature T pc .A natural question is whether this type of spectral characteristic is also displayed in meson states comprised of heavier quarks.In Sec. 3 we will explore and confirm this possibility using lattice QCD data for N f = 2 + 1 flavour pseudo-scalar meson correlators.

Lattice setup and analysis strategy
For the analysis in this work we studied the Euclidean correlators of pseudo-scalar meson operators O PS in QCD, focussing in particular on the meson channels comprised of one light and one strange quark: O PS = lγ 5 s, and two strange quarks: O PS = sγ 5 s.For this purpose we calculated both spatial and temporal correlators in N f = 2 + 1 flavour QCD with two mass-degenerate flavours of light quarks and a strange quark, in two different configurations: {β = 7.010, T = 36.4,145.6 MeV}, and {β = 7.188, T = 43.1,172.3 MeV}.We used gauge field configurations generated with the highly improved staggered quark (HISQ) action with physical values for the light and strange quark masses.Details of these configurations, including the choice of scale setting, can be found in Ref. [27].The logic for choosing these configurations is that it enabled a self-consistent computation of the correlators at an effectively vanishing temperature and at temperatures both below and above the pseudo-critical temperature T pc = 156.5 MeV [45], hence allowing one to investigate meson spectral properties within the crossover region.The calculation of the correlators were performed using Möbius domain wall fermions (MDWF) [46] and a mixed action approach, with domain wall valence quarks and staggered sea quarks.More details regarding the numerical setup, including the choice of lattice action and parameter tunings, can be found in Appendix A. The analysis strategy was based on the proceedure adopted in Ref. [29], and used the following steps: 1. Fit the form of the spatial correlator C PS (z) at each temperature.In both channels and in each configuration these data were consistent with a sum of pure exponentials, and so the following functional form was fitted with m scr i the screening mass of the different components.The results of these fits are listed in Tables 1 and 2. The second term is needed in order to take into account the periodic boundary conditions of the lattice meson fields along the z direction, with L the symmetric spatial lattice extent.Since m scr i → m i for vanishing temperature, where m i is the corresponding vacuumstate mass, analysis of the low-temperature configuration samples (36.4 MeV and 43.1 MeV) enabled us to determine an estimate for m i .
2. Use the best-fit parameters {c i , m scr i } at the higher temperatures (145.6 MeV and 172.3 MeV) and the extracted zero-temperature masses m i , together with the analytic relations in Sec.2.2, to compute the damping factors of the hypothesised thermoparticle components.As demonstrated in Ref. [29], if these components are present, and dominate at low energies, it follows that the summed pure exponential behaviour in Eq. (3.1) corresponds to N thermoparticle states, each with damping factors of the form The low-energy contributions from other spectral components in Eq. (2.4) are assumed to be negligible in this hypothesis.The validity of this hypothesis is tested in the final step.
3. Apply the spectral representation in Eq. (2.1) to calculate the thermoparticle contributions to the full spectral function ρ PS (ω, ⃗ p).As with the light-light pseudo-scalar meson case studied in Ref. [29], these contributions ρ 4. Use these spectral function expressions to compute the thermoparticle contribution to the corresponding temporal correlator C PS (τ, ⃗ p) via Eq.(1.2), and compare these predictions with the lattice data.This comparison can be directly made at each temperature by using the nonrenormalised lattice data.To compare the spectral functions across different temperatures one must then make a choice of renormalisation prescription.
In step 1 the functional form in Eq. (3.1) was fitted to the lattice data points C PS (z = n z a; T ), with a the lattice spacing and n z the integer spatial site number, where 0 ≤ n z ≤ N s − 1 with L = N s a.These fits were performed in the data range [n min , N s /2] for various values of n min , and for N = 1, 2 and 3 states in order to assess their stability.The best-fit parameter values {c i , m scr i } were obtained by minimising the specific Chi-squared statistic χ 2 (n min , N )/d.o.f., and the fits were deemed stable if these values remained insensitive to variations of n min in some non-trivial range.The final parameter values were determined by averaging the best-fit values over this plateau-like region, and their corresponding errors were obtained using a bootstrap analysis in which fits were performed using random samples of the data in order to estimate the standard deviation.The plotted results of this stability analysis for the best-fit parameters {c i , m scr i } are displayed in Appendix B in Figs.11-14.The horizontal lines in the plots indicate the estimated parameter range from the bootstrap analysis.The plateau-averaged value of the best-fit parameters, and their bootstrap-calculated errors, are listed in Tables 1 and 2.
Having determined the best estimates for the parameters {c i , m scr i } and their respective errors in step 1, in the remaining steps these values were used to compute the total thermoparticle contributions to ρ PS (ω, ⃗ p), and compare the corresponding predictions for C PS (τ, ⃗ p) with the lattice data.These comparisons are discussed in Secs.3.2 and 3.3.Associated with the computations of ρ PS (ω, ⃗ p) and C PS (τ, ⃗ p) are uncertainties due to the error in the fit parameters.These uncertainties were computed using a bootstrap-type analysis and are displayed as error bands in the figures of Secs.3.2-3.4.The fit parameter errors also enter in the analysis of the prediction quality in Appendix C. By following the procedure set out in steps 1-4 we were able to assess whether the QCD lattice data in both the lγ 5 s and sγ 5 s channels were consistent with the presence of thermoparticle-type components in the spectral function ρ PS (ω, ⃗ p) both below and above T pc .

Light-strange channel
Based on the analysis procedure set out in Sec.3.1, the results for the fit parameters and their respective errors are summarised in Table 1.For the pseudo-scalar meson operator O PS = lγ 5 s the lowest-lying particle state at zero temperature is the kaon K, which is stable in QCD and has a neutral charge state K 0 of mass 497.6 MeV [47].At higher energies the situation is significantly less certain, but a recent analysis from the LHCb experiment suggests that the K(1460) is the next heaviest state, with a mass of 1482.4MeV, and corresponds to the first radial excitation of the K [48].From Table 1 one can see that the low temperature fit for β = 7.188 agrees quite well with the experimentally measured m K and m K(1460) values, but the β = 7.010 case gives an overestimation of the masses, particularly for the excited state.In both configurations there is a significant increase in the screening mass of the kaon between the low and high temperatures, and the values at 145.6 MeV and 172.3 MeV are consistent with the continuum extrapolated results obtained in Ref. [27].For the first excited state the change in m scr runs in opposite directions as the temperature is increased in the different configurations.For β = 7.188 the m scr value increases  In Fig. 1 is displayed the low-energy contribution to the spectral function using the T = 145.6MeV parameters in Table 1, together with the resulting prediction for the temporal correlator C PS (τ, ⃗ p) and its comparison with the corresponding lattice data for different values1 of |⃗ p|.In the left plot one can see that the kaon state undergoes collisional broadening but the first excited state remains on its vacuum mass shell.In the right plot one finds that the temporal correlator prediction agrees well 2 with the data at large τ , and then starts to underestimate the data as τ becomes smaller.The point at which this underestimation occurs coincides with the approximate boundary where the two-state spatial correlator fit starts to break down, which is indicated by the vertical dashed line in the right plot.This is consistent with the fact that higher excited states, which should increasingly dominate at smaller τ , are not included in the prediction.It should be noted that for the largest values of |⃗ p| there start to arise sizeable lattice cutoff effects, and so the comparison at these values becomes significantly less certain.are significantly broadened, and that the temporal correlator predictions agree with the data, but only for the largest values of τ and smaller values of momentum |⃗ p|, as quantified in Appendix C. It is difficult to make firm conclusions here since the deviations at larger |⃗ p| may be as a result of cutoff effects, and the slight overestimation in prediction between the smallest τ points and the dashed line, where the two-state fit starts to break down, could be due to systematic effects.As τ is decreased below the dashed line the prediction immediately begins to undershoot the data, which is consistent with the absence of higher-energy spectral function contributions in the prediction.So far we have considered the form of the non-renormalised spectral functions and corresponding temporal correlator predictions separately for each lattice setup.In order to compare the spectral functions at the temperatures below and above T pc one must choose a renormalisation prescription.By making the choice: C PS (τ = 0.5 fm, ⃗ p = 0) = 1 MeV 3 at both T = 36.4MeV and T = 43.1 MeV, the renormalised spectral function at ⃗ p = 0 has the form in Fig. 3.One can see that as the temperature increases through the crossover region both the ground and first excited states become collisionally broadened, with the amplitude of their peaks significantly reduced.

Strange-strange channel
Table 2 summarises the results obtained for the fit parameters and their respective errors after applying the procedure set out in Sec.3.1 to the strange-strange correlator data.Unlike the lγ 5 s channel, for the operator O PS = sγ 5 s the lattice correlator can contain a disconnected component, and this could in principle be non-negligible.However, for the purposes of this study we follow the approach of Refs.[26,27] in treating the ground state to be a hypothetical unmixed projection η ss of the η meson which is generated by the purely connected correlator component.The physical characteristics of higher excited states in this channel are less certain, but it has been suggested that η(1295) could be the first radial excitation of the η, and that it is mainly comprised of an ss component [47].In both lattice configurations there is a significant increase in the screening mass of the η ss between the lower and higher temperature data fits, and the values at all tempera- T    tures are consistent with the continuum-extrapolated results obtained in Ref. [27].For β = 7.010 the screening mass of the first excited state does not change significantly as the temperature is increased, but for β = 7.188 the value of m scr decreases between 43.1 MeV and 172.3 MeV.In Fig. 4 is displayed the low-energy contribution to the spectral function using the T = 145.6MeV parameters in Table 2, together with the resulting prediction for the temporal correlator and its comparison with the corresponding lattice data for different values of |⃗ p|.Just like in the lγ 5 s channel, in the left plot one can see that the η ss ground state broadens, but the first excited state remains on its vacuum mass shell.In the right plot one finds that the temporal correlator prediction also has the same characteristic features, matching the data well3 at large τ , and then starting to increasingly underestimate it as τ becomes smaller.This coincides with the region where the two-state fits start to break down, and is therefore consistent with the absence of higher-energy spectral function contributions.Although there are sizeable lattice cutoff effects for larger values of momentum, the prediction still provides a good description of the large τ data for all values of |⃗ p|.
Due to the decrease in the excited state screening mass between 43.1 MeV and 172.3 MeV there are two possible interpretations: either this decrease is not significant once systematic effects are taken into account, or this decrease has a physical origin.Given the first interpretation, one should therefore simply set the vacuum mass of the excited state equal to its screening mass at 172.3 MeV.The corresponding spectral function and temporal correlator predictions are plotted in Fig. 5.Here one can see very good agreement between the prediction and lattice data at large τ , as is the case at 145.6 MeV.In the second interpretation, the decrease in the excited state screening mass could be explained by the breakdown in the assumption that the thermoparticle-like components are entirely dominant at low energies.If one fixes the excited state mass to the value obtained at 43.1 MeV, and assumes that this state does not broaden, one obtains the results in Fig. 6.In this case one can see that the predictions underestimate the data at all momentum values.If true, this would give support to the hypothesis that the spectral function contains a significant non-thermoparticle-like low-energy contribution.By definition, this contribution would appear in the thermal spectral density like the second term in Eq. (2.4) and ultimately behave like a collective thermal excitation.

Model discrimination
A well-known problem with standard reconstruction approaches is that different reconstructed spectral functions may display large differences, and yet still provide a good description of the correlator data from which they were extracted.However, by comparing both spatial and temporal correlator predictions one can break this degeneracy, since the correlators have a different dependence on the spectral function, as discussed in Sec. 1.To demonstrate the utility of this approach, in the remainder of this section we will consider a different causal spectral function which also gives rise to a purely exponential spatial correlator contribution, namely the relativistic Breit-Wigner where m is the zero-temperature mass, and Γ a temperature-dependent width.Although this contribution is causal [49] it does not describe a thermoparticle-type state, since its corresponding thermal spectral density D β (⃗ x, s) cannot be written: D m,β (⃗ x) δ(s − m 2 ).This implies that a Breit-Wigner excitation of the form in Eq. (3.4) describes an intrinsically unstable state, in contrast to the thermoparticle case.Combining Eqs.(3.4) and (1.4), the spatial correlator takes the form  broadening at this temperature, these predictions correspond to spectral functions with the same on-shell excited state contribution as the thermoparticle case, but with a Breit-Wigner ground state component.In Figs. 7 and 8 there are clear differences with the thermoparticle predictions, which for ease of comparison are plotted together in Fig. 9 for the largest τ data points.The Breit-Wigner contribution gives rise to a flatter temporal correlator behaviour at large τ , and this leads to a poorer description of the data in both meson channels.This is most pronounced in the sγ 5 s channel, where the prediction is significantly above the largest τ data points for all values of |⃗ p|.
Overall, these findings suggest that the lowest-lying spectral function components in each channel are not consistent with a ground state contribution with a causal Breit-Wigner structure.It is interesting to see here that even at this relatively low temperature of 145.6 MeV, the comparison of both spatial and temporal correlators results in a non-trivial discriminating power between different spectral function models.

Spectral implications for QCD
The analyses in Secs.3.2 and 3.3 are consistent with the hypothesis that both the lγ 5 s and sγ 5 s meson spectral functions contain a stable bound-state-like thermoparticle component which dominates the low-energy behaviour of ρ PS (ω, ⃗ p) at T = 145.6MeV, and has the form of Eq. (3.3).In Appendix D it is demonstrated that these components have a purely non-perturbative character.Although 145.6 MeV is a relatively low temperature, the contribution of these components to the temporal correlator is distinct enough to be differentiated from other casual spectral function models, as outlined in Sec.3.4.There are two particularly significant aspects to these findings: 1.The thermoparticle picture appears to provide a good description of how particle states are modified in the presence of a thermal medium.Instead of behaving like collective thermal excitations, these states are intrinsically stable, and their amplitudes dissipate solely due to their interactions with the medium.
2. The qualitative features of ρ PS (ω, ⃗ p) are the same in both meson channels, and in fact also agree with those found in Ref. [29] in the lγ 5 l channel.This suggests that the low-energy spectral functions of pseudo-scalar mesons in QCD have a common structure which is fixed by the vacuum dynamics of the theory.
Both of these results give support to the non-perturbative finite-temperature QFT framework put forward in Refs.[30][31][32][33][34].The second observation also appears to align with the findings of Ref. [34], where it was demonstrated that the behaviour of damping factors associated with asymptotic states, and hence low-energy spectral components, can actually be fixed by the equations of motion.Due to the significantly larger errors associated with extracting the parameters of excited states, the conclusions regarding the status of these spectral contributions remain more uncertain.But within the estimated uncertainties, the T = 145.6MeV data in both lγ 5 s and sγ 5 s channels are consistent with the first excited states undergoing no thermal modifications, and hence remaining on their vacuum mass shells.Physically, this is perhaps not so surprising, since the vacuum masses of these states are still significantly larger than the temperature scale, which is below T pc .At T = 172.3MeV, which is above T pc , the excited states have more of an effect on the structure of the spectral function.In the lγ 5 s channel the results are consistent with the first excited state undergoing collisional broadening, whereas for sγ 5 s it appears that the excited state remains on shell, but that there may be additional non-thermoparticle-like contributions.This suggests that as the temperature passes through the crossover region and starts to increase above T pc , the flavour dependence of the meson states has a more significant bearing on the low-energy spectral properties of the theory.

Conclusions
By utilising the general non-perturbative constraints satisfied by correlation functions at finite temperature, in this work we analysed how the spectral properties of pseudo-scalar mesons comprised of light-strange and strange-strange quarks are modified as one passes through the high-temperature chiral crossover region.A key component of this analysis was to test the hypothesis that stable particle-like excitations, so-called thermoparticles, exist in thermal media and provide a natural description for low-energy meson states at finite temperature.For this purpose we followed the strategy developed in Ref. [29], where it was established that thermoparticle contributions to spectral functions can be directly extracted from the behaviour of the spatial correlators, and then used to predict the form of the corresponding temporal correlator C PS (τ, ⃗ p).Since the spatial and temporal correlators depend on the spectral function in different ways, this provides a highly nontrivial test of the validity of the extracted thermoparticle contributions.
Applying this procedure to N f = 2 + 1 flavour lattice QCD data we show that these data are consistent with the appearance of a stable bound-state-like thermoparticle contribution in the spectral function for both lγ 5 s and sγ 5 s flavour mesons, and that this contribution is present for temperatures around the pseudo-critical temperature T pc .This generalises the results of Ref. [29], where evidence of such a component was found in the lγ 5 l channel for N f = 2.We tested the robustness of the thermoparticle contribution by comparing the C PS (τ, ⃗ p) prediction with lattice data at multiple momentum values, and also with a different causal spectral function model.Overall, these findings suggest that the thermoparticle hypothesis provides a consistent description of how vacuum states undergo collisional broadening at finite temperature, and that non-perturbative effects continue to play an important role for hadrons in the chiral crossover region.This picture aligns with the expectations of chiral spin symmetry, where chiral symmetry is effectively restored but a dominant fraction of quarks remain bound [50][51][52].Although the focus of this study was correlators involving pseudo-scalar meson operators, the approach itself can in principle be generalised to other types of hadronic states of different spin, as well as to systems with non-vanishing baryon density.This could help provide a better understanding of in-medium effects in different regions of the QCD phase diagram.We leave these generalisations to future works.authors gratefully acknowledge the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) for funding this project by providing computing time through the John von Neumann Institute for Computing (NIC) on the GCS Supercomputer JUWELS [53] at the Jülich Supercomputing Centre (JSC).We also acknowledge the EuroHPC Joint Undertaking for awarding this project access to the EuroHPC supercomputer LUMI, hosted by CSC (Finland) and the LUMI consortium through a EuroHPC Extreme Scale Access call.For the calculations on these machines we were using Grid [54,55] and the Grid Python Toolkit (GPT) [56].Part of the analysis was performed on the GPU cluster at Bielefeld University.We thank the Bielefeld HPC.NRW team for their support.
All of the lattice data used in the analysis of this paper and presented in each of the figures can be found in Ref. [57].

A Mixed action approach and mass tuning
Ideally, any lattice QCD action should fulfil nearest neighbour locality, be free of doublers, and have a massless Dirac operator that is chirally symmetric.However, it has been proven in even dimensions that it is not possible to simultaneously satisfy all three of these properties [58][59][60].One therefore has to weigh up the advantages and disadvantages of different lattice actions and select the one that best suits the problem at hand.For the purposes of chiral symmetry preservation, common choices are the highly-improved staggered quarks (HISQ) [61,62], or a form of domain wall fermions (DWF) [63,64].The HISQ action is computationally cheap compared to other lattice actions, but there only remain traces of chiral symmetry, whereas DWF actions have improved chiral characteristics, but their realisation requires the introduction of a fifth lattice dimension which makes them computationally more expensive.As a compromise, one can mix these actions to benefit from both the cheapness of the HISQs and the good chiral properties of the DWF action.
For the purposes of this work we used Möbius domain wall fermion (MDWF) observables [46], and computed them on HISQ configurations.DWF actions are known to only fulfil chiral symmetry approximately when the lattice extent of the fifth dimension L 5 is finite.The reliability of this approximation can be measured by the so-called residual mass [46] m res = x ⟨ Q(t, x)γ 5 Q(t, x)q(0, 0)γ 5 q(0, 0)⟩ x ⟨q(t, x)γ 5 q(t, x)q(0, 0)γ 5 q(0, 0)⟩ where Q(t, x) is a quark field in the midpoint of the fifth dimension, and q (t, x) is a quark field bound to the domain wall.m res depends on the DWF parameters, which in the specific case of MDWFs are L 5 , the domain wall height M 5 , and b 5 , c 5 , which are related to the lattice spacing of the fifth dimension a 5 via: a 5 = b 5 − c 5 .In the limit L 5 → ∞ the residual mass vanishes, and hence chiral symmetry is restored.
A drawback of the mixed action approach is that an additional mass term ∆m mix [65] appears which causes the bare sea and valence quark masses to differ, thus requiring the valence quark mass to be tuned.We performed this tuning on the lowest-temperature configurations detailed in Tables 1 and 2 in two steps: 1. Adjust the MDWF parameters to reduce m res such that its value is to within 10% of the light valence bare quark mass m l,val .
2. Vary m q,val for q ∈ {l, s} until the pseudo-scalar meson ground state mass m q q is consistent with the corresponding physical neutrally charged particle mass, as listed in the PDG [47].The parameter values from this procedure are listed in Tab. 3, and in Fig. 10 is shown the comparison between the physical ground state masses and those extracted using the tuned parameters.
Figure 10: Pseudo-scalar meson ground state masses calculated with the tuned bare quark masses in Tab. 3, and the corresponding physical neutrally charged particle masses from the PDG [47].Since there is no pure pseudo-scalar ss bound state, we used the value mη = 2(m K 0 ) 2 − (m π 0 ) 2 in order to tune the strange bare quark mass.

C Prediction quality
To provide a quantitative measure of the agreement between the temporal correlator data C data where n * τ is the first τ point above the dashed line of Figs 1, 2, 4-6 where the prediction is applicable.Furthermore, we define ∆ stat to be the average of the statistical error bars on the lattice data, and ∆ theor the average of the bootstrap-computed error band on the theoretical prediction.Details of the error band computation can be found in Sec. 3

D Perturbative comparison
In Fig. 15 we compare the leading-order perturbative predictions for the temporal correlator at zero momentum with the lattice data at T = 145.6MeV and 172.3 MeV in both channels.To remove unknown renormalisation constants we perform this comparison for the unit-normalised correlator at n τ = 1.Since the quark masses are small relative to those of the meson ground states, the prediction is computed for massless quarks using the spectral function derived in Ref. [5].From Fig. 15 one can see that the perturbative predictions significantly underestimate the data for large τ .This demonstrates that these predictions only affect the large-energy behaviour of the spectral function, which dominates the small τ structure of C PS (τ, ⃗ p = 0).It should be noted that the perturbative predictions displayed in Fig. 15 are continuum results, and therefore do not take lattice discretisation effects into account.However, such effects are generally small in the low-energy regime [5], and would not significantly alter the large-τ predictions.These results demonstrate that the perturbative contributions to ρ PS (ω, ⃗ p) are highly sub-dominant at low energies, and hence the spectral components analysed in Sec. 3 are of a non-perturbative nature.One expects that this sub-dominance will change at higher temperatures due to the increasing importance of nonthermoparticle-like components, as outlined in Sec.3.3.

Figure 1 :
Figure 1: lγ5s channel spectral function extracted from spatial correlator lattice data at 145.6 MeV (left), and the temporal correlator prediction and corresponding lattice data (right).The vertical dashed line in the right plot indicates the approximate boundary below which knowledge of further excited states is necessary.The effects of the vacuum mass mi uncertainties are not included in the left plot in order to improve its clarity.

Figure 2 Figure 2 :
Figure2shows how the spectral function and corresponding temporal correlator prediction change as the temperature is increased above T pc .One can see that both the kaon and first excited state

Figure 3 :
Figure 3: Renormalised lγ5s channel spectral function for T = 145.6MeV and 172.3 MeV.The effects of the uncertainties in the vacuum masses mi on the spectral function are not included in the plot.This explains why the excited state components at T = 145.6MeV and 172.3 MeV appear not to possess the same energy threshold, although the extracted vacuum masses are consistent within uncertainties.

Figure 4 :
Figure 4: sγ5s channel spectral function extracted from spatial correlator lattice data at 145.6 MeV (left), and the temporal correlator prediction and corresponding lattice data (right).The vertical dashed line in the right plot indicates the approximate boundary below which knowledge of further excited states is necessary.The effects of the vacuum mass mi uncertainties are not included in the left plot in order to improve its clarity.

Figure 5 :
Figure 5: sγ5s channel spectral function extracted from spatial correlator lattice data at 172.3 MeV (left), and the temporal correlator prediction and corresponding lattice data (right).The vertical dashed line in the right plot indicates the approximate boundary below which knowledge of further excited states is necessary.In this plot the excited state vacuum mass is set equal to the screening mass.The effects of the vacuum mass mi uncertainties are not included in the left plot in order to improve its clarity.

Figure 6 :
Figure 6: sγ5s channel spectral function extracted from spatial correlator lattice data at 172.3 MeV (left), and the temporal correlator prediction and corresponding lattice data (right).The vertical dashed line in the right plot indicates the approximate boundary below which knowledge of further excited states is necessary.In this plot the excited state mass is set equal to its value at 43.1 MeV.The effects of the vacuum mass mi uncertainties are not included in the left plot in order to improve its clarity.

. 5 )Figure 7 :
Figure 7: lγ5s spectral function extracted from spatial correlator data at 145.6 MeV assuming a Breit-Wigner ground state form (left), and the temporal correlator prediction and corresponding lattice data (right).The vertical dashed line in the right plot indicates the approximate boundary below which knowledge of further excited states is necessary.

Figure 8 :
Figure 8: sγ5s spectral function extracted from spatial correlator data at 145.6 MeV assuming a Breit-Wigner ground state form (left), and the temporal correlator prediction and corresponding lattice data (right).The vertical dashed line in the right plot indicates the approximate boundary below which knowledge of further excited states is necessary.

Figure 9 :
Figure 9: Comparison of the temporal correlator predictions with the 145.6 MeV lattice data assuming either a thermoparticle (dark shading, solid line) or Breit-Wigner (light shading, dashed line) ground state in the lγ5s channel (left), and the sγ5s channel (right).This is plotted in the reduced τ range for which the ground and first excited states are dominant.
is the position-space thermal spectral density, which only depends on |⃗ x| due to the assumed isotropy of the thermal ground state.Since the properties of the spectral function are entirely determined by the structure of D β (⃗ x, s), it is not surprising that this is also the case for C(z).Since Eq. (2.2) is a general representation, essentially relying only on the causality of the theory, one can use it to investigate how model-specific characteristics are reflected in the properties of C(z).This will be discussed in Sec.2.2.

Table 1 :
Parameter fit values and their estimated uncertainties from the OPS = lγ5s spatial correlator data.

Table 2 :
Parameter fit values and their estimated uncertainties from the OPS = sγ5s spatial correlator data.

Table 4 :
.1.The values of MAPE, ∆ stat , and ∆ theor are displayed in Table4for each channel and spatial momenta.If the MAPE is comparable or smaller than the combined lattice and theoretical errors, namely: MAPE ≲ ∆ stat + ∆ theor , the prediction provides a good description of the data.Channel T [MeV] p [MeV] MAPE [%] ∆stat [%] ∆ theor [%]The mean absolute percentage error (MAPE) of the CPS(τ, ⃗ p) predictions, average statistical lattice error, ∆stat, and theoretical uncertainty, ∆ theor , across the prediction range.In descending order these values correspond to the data and predictions in Figs. 1, 2, 4, 5, and 6, respectively.