Pion spectral properties above the chiral crossover of QCD

Spectral functions encode a wealth of information about the dynamics of any given system, and the determination of their non-perturbative characteristics is a long-standing problem in quantum field theory. Whilst numerical simulations of lattice QCD provide ample data for various Euclidean correlation functions, the inversion required to extract spectral functions is an ill-posed problem. In this work, we pursue previously established constraints imposed by field locality at finite temperature $T$, namely that spectral functions possess a non-perturbative representation which generalises the well-known K\"{a}ll\'{e}n-Lehmann spectral form to $T>0$. Using this representation, we analyse lattice QCD data of the spatial pseudo-scalar correlator in the temperature range $220-960 \, \text{MeV}$, and obtain an analytic expression for the corresponding spectral function, with parameters fixed by the data. From the structure of this spectral function we find evidence for the existence of a distinct pion state above the chiral pseudo-critical temperature $T_{\text{pc}}$, and contributions from its first excitation, which gradually melt as the temperature increases. As a non-trivial test, we find that the extracted spectral function reproduces the corresponding temporal lattice correlator data for $T=220 \, \text{MeV}$.


Introduction
The phases of QCD under extreme conditions and the nature of its associated effective degrees of freedom are among the most pressing problems of theoretical physics, affecting experimental programs from heavy ion collisions, to astro-particle and gravitational wave physics. Of particular interest is the question, deeply related to the confinement problem, of how ordinary hadronic matter gets modified in medium to eventually dissolve into the expected quark gluon plasma. In principle, the answer is provided by the spectral properties of Euclidean two-point functions of gauge-invariant operators O Γ (τ, x) where Γ denotes a set of quantum numbers, and the expectation value is over a thermal ensemble at temperature T = 1/β, which for the purposes of this work we specialise to zero baryon density, µ B = 0. The spatial Fourier transform of the correlation functions take the universal form where the associated spectral functions ρ Γ (ω, p) contain the desired information about the possible excitations in a given quantum number channel.
In order to fully describe a thermal system of strong-interaction particles one ultimately requires a non-perturbative approach. Lattice QCD is a powerful such approach, and numerical computations of the correlation functions in Eq. (1.1) have led to significant advancement in the understanding of finite-temperature phenomena [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17]. On the lattice it is the spatial correlators, integrated over the orthogonal directions, that are accessible over the largest distances and hence contain the most information about the thermal system. Particular progress has been made in establishing the properties of these correlators in recent years [18][19][20][21][22][23][24][25], which are defined as where the second equality follows from Eq. (1.2). Equation (1.4) demonstrates that the structure of C Γ (x 3 ) is entirely controlled by the spectral function of the corresponding hadronic operator, and thus directly probes the spectral properties of QCD for T > 0. A fundamental issue with Eq. (1.2) though is that the extraction of the corresponding spectral function is an ill-posed inverse problem. For this reason, most strategies to obtain QCD spectral functions from lattice data require intricate statistical methods, combined with input based on either perturbative calculations or phenomenological modelling. A general discussion of such strategies can be found in Refs. [26,27].
In this work we pursue a different approach, which was developed in Refs. [28][29][30][31][32] and is based on the T > 0 generalisation of axiomatic formulations of local vacuum-state QFT, whose applications over the years have led to numerous foundational insights, including the relationship between spin and statistics, the generality of CPT symmetry, and the rigorous connection of Minkowski and Euclidean QFTs [33][34][35]. The non-perturbative T > 0 framework of Refs. [28][29][30][31][32] focussed on the simplest case of Hermitian scalar fields φ(x), and established that characteristic features such as the loss of Lorentz symmetry can be incorporated by defining a thermal background state |Ω β at temperature T = 1/β, which is no longer invariant under the full Poincaré group. Together with the standard constraints brought about by the assumption of thermal equilibrium [36,37], it was demonstrated that the locality of the fields 1 alone imposes particularly significant constraints, and implies the following representation of the scalar spectral function 2 [28] ρ(ω, p) In the T → 0 limit Eq. (1.5) reduces to the well-known Källén-Lehmann spectral representation 3 , and hence represents its T > 0 generalisation. From the structure of Eq. (1.5) one can see that the effects of the background state are entirely captured by the thermal spectral density D β ( u, s).
Determining the properties of this quantity is therefore essential for describing the characteristics of scalar particles in thermal media.
Equation (1.5) is completely general and holds for any scalar field satisfying locality. In order to understand the characteristics of the thermal spectral density for specific theories, additional information is necessary. In Ref. [28] the authors proposed that the singular structure of D β ( u, s) in the variable s is preserved relative to the vacuum theory, as long as no phase transition is met. This implies that discrete and continuous contributions are decomposed, and hence if a theory contains a stable particle state of mass m at T = 0 one can write [32] where D c,β ( u, s) is continuous in s. As outlined in Refs. [28][29][30][31][32], there are several reasons for why the discrete component in Eq. (1.6) provides a natural description of a particle state in finitetemperature QFT. Firstly, for T > 0 the so-called damping factor D m,β ( u) is non-trivial, which due to the structure of Eq. (1.5) causes ρ(ω, p) to have contributions outside of the mass shell p 2 = m 2 , and hence the T = 0 mass peak of the particle becomes broadened, as one would expect.
Moreover, the precise nature of this broadening has been shown to be controlled by the underlying interactions between the particle state and the constituents of the thermal medium [32]. Although Eq. (1.6) assumes the T = 0 particle state to be stable, this can in principle be generalised to unstable states by replacing δ(s − m 2 ) with a suitable resonance-type function in s, such as a relativistic Breit-Wigner. The factorisation of the ( u, T ) and s dependence therefore ensures that this representation can distinguish between particle decays of different physical origin, namely those brought about by dissipative thermal effects, controlled by ( u, T ), and those due to any intrinsic instability of the T = 0 particle states. The structure of damping factors in specific models was first 1 By locality we mean: [φ(x), φ(y)] = 0 for (x − y) 2 < 0, which is simply the physical assumption that all measurements respect causality. 2 As is standard in the literature, the spectral function refers to the Fourier transform of the thermal two-point commutator 3 For the vacuum spectral function the Källén-Lehmann representation has the momentum-space form [38,39]: 2π (ω) ∞ 0 ds δ p 2 − s ρ(s), where ρ(s) is the spectral density whose singularities capture the presence of stable particle states. explored in Ref. [32], and more recently in Refs. [40][41][42], where it was shown that these quantities can in fact be used to perform non-perturbative analytic calculations of in-medium observables, including the shear viscosity.
The goal of the present work is to apply the observations from Refs. [41,42] to lattice data for pseudo-scalar correlators in order to compute the thermal spectral density D β ( u, s), and from it the spectral function ρ(ω, p). In Sec. 2 we establish an analytic connection between the spatial correlator and the spectral representation in Eq. (1.5). In Sec. 3 we use these results to extract the properties of the light-quark pseudo-scalar spectral function from lattice QCD data, and discuss the physical implications of these results. Finally, in Sec. 4 we summarise our main findings.

Spectral representation of spatial correlators
We begin by discussing the relation between the thermal spectral density D β ( u, s) and the imaginarytime two-point function C(τ, x). In Ref. [41] it was shown for scalar fields that this connection is defined via the following integral relation where D β ( x, s) is the inverse Fourier transform of D β ( u, s). Using the fact that D β ( x, s) depends only on | x| in an isotropic medium, it follows by combining Eqs. (1.3) and (2.1) that the spatial correlator has the general spectral representation As outlined in Sec. 1, in the simplest case where there exists a single stable particle state at T = 0, there is evidence to suggest [32] that the thermal spectral density has the decomposition in Eq. (1.6). However, hadronic correlators in QCD in general receive contributions from multiple states. In this case, the decomposition in Eq. (1.6) has the natural position-space generalisation where m i are the T = 0 masses associated with these states. Using analogous arguments to those in Ref. [41], if D c,β ( x, s) is non-vanishing for s ≥ s c , and the gap between the scale s c and some subset of the masses {m 1 , ..., m n } is sufficiently large, the damping factors associated with these states will dominate the behaviour of C(x 3 ) 4) and this domination will be especially pronounced at large x 3 . Equation (2.4) demonstrates not only that the vacuum particle states have a significant impact on the behaviour of C(x 3 ), but also how this behaviour is connected to the dynamical properties of the theory. In the situation that one of the masses m 1 is significantly smaller than the others, it follows from Eq. (2.4) that the corresponding damping factor D m 1 ,β associated with this state can be directly calculated from the leading large-x 3 behaviour of the spatial correlator derivative In principle, the damping factors associated with the higher mass states {m 2 , ..., m n } can also be extracted, in this case from the sub-leading large-x 3 behaviour of C(x 3 ). The practical feasibility of this depends on the relative size of the masses and the behaviour of D m i ,β ( x). Once explicit particle damping factors are available, one can then use Eq. (1.5) to calculate their analytic contributions to the full spectral function ρ(ω, p).

Preliminary considerations
Before applying the relations derived in the last section to lattice QCD data, let us recall a result from the literature that explicitly supports the validity of the ansatz in Eq. (1.6) in QCD, at least for low temperatures and in the chiral limit. In Refs. [43,44] QCD vector and axial-vector correlators were evaluated using the well-tested PCAC current algebra, which to leading order in = T 2 /(6f 2 π ) allows one to express the low-temperature correlators by the vacuum correlators Whilst these are not scalar correlators, they provide particular examples for the more general assumption leading to Eq. (1.6), namely that certain analytic structures of vacuum correlators are preserved at finite temperatures. In the case of Eq. (3.1), the poles of the vacuum correlators are mixed, but explicitly present in the finite-temperature correlators. Similar considerations apply, in less detail but more generality, to the spectral decomposition of lattice meson correlators, where Z is the partition function, and H has the discrete spectrum H|n = E n |n . Since the Hamiltonian is independent of T , this spectrum is identical to that in the vacuum, as long as no non-analytic phase transition is crossed. Again, it is evident that the finite-temperature correlator inherits the vacuum particle states and their analytic structures, which get mixed with weight factors determined by temperature 4 . Since there is no non-analytic phase transition in QCD at µ B = 0, the proposition in Eq. (1.6) appears well suited for lattice pion correlators over a large temperature range, and we shall now use it to extract the corresponding spectral information.  Table 1.

Pseudo-scalar lattice correlator analysis
We now focus on spatial pseudo-scalar correlators defined as ) isospin doublets, and τ a the Pauli spin matrices. For the remainder of the paper the isospin index a will be dropped since the spatial correlator is independent of a for degenerate quark masses. We analyse the lattice data generated in Ref. [45] with two mass-degenerate flavours of domain wall fermions, which have good chiral symmetry properties also at finite lattice spacing. Thermal ensembles for T = 220, 320, 380, 480, 660, and 960 MeV were generated with a fixed bare quark mass of am u = am d = 0.001, which is compatible with the physical pion mass at T = 0. The correlator in Eq. (3.3) was evaluated in the range x 3 ∈ [x * 3 , L/2], where L is the spatial lattice extent, and x * 3 = 1/(T N τ ) is the smallest non-vanishing separation by one lattice spacing, with N τ the temporal lattice extent. The calculation in Ref. [45] was also performed such that the correlator was normalised to one at x 3 = x * 3 . To avoid ambiguity we will refer to this normalised correlator as C PS (x 3 ). More details regarding the precise lattice setup can be found in Ref. [45] and a brief summary is given in Table 1.
We begin our analysis by parametrising the functional form of the correlator C PS (x 3 ). In order to do so, one must first understand the particle states that are generated by the light quark pseudo-scalar operator O a PS . For masses below 2 GeV, these states consist of the pion π, π(1300), and π(1800) [46], the latter two of which are generally interpreted as being the first two radial excitations of the pion, the π * and π * * [47]. In accordance with numerous previous investigations, we find the large x 3 behaviour of C PS (x 3 ) to be consistent with a pure exponential. Taking into account the periodic boundary conditions along the x 3 direction, we thus started by fitting the (7) 776 (17) (8) 8604(77) 1.243 Table 1: Lattice parameters for the spatial pseudo-scalar correlator data from Ref. [45], together with the best-fit parameters pertaining to the dashed lines in Fig. 1.
standard functional form where m scr π defines the so-called screening mass of the lowest energy state, in this case the pion. The lower boundary of the fit range x min 3 was chosen such that the inclusion of an additional data point would significantly reduce the quality of the fit, producing an O(1) increase in the χ 2 /d.o.f. By adding an additional exponential component to our fit ansatz, the overall quality of the fit was significantly improved, enabling excellent fits over the entire data range [x * 3 , L/2] for T = 660 MeV and 960 MeV, and excluding only the smallest separation point for T = 220, 320, 380, and 480 MeV. These data and the fits are shown in Fig. 1. The corresponding fit parameters, together with their uncertainties, are listed in Table 1.
In Eq. (3.5) the exponents of the second pair of terms have the interpretation of a screening mass for the second-lightest pseudo-scalar meson state, the π * . As opposed to the pion, the π * is unstable under strong interactions, but this is not observed in the data, presumably due to an insufficient lattice resolution. The exact correlator can also contain additional contributions from higher excited states, including the π * * , although these states are increasingly exponentially suppressed. One could in principle extract the properties of these states from the short distance behaviour of the correlator on finer lattices with higher resolution and more lattice data points, but at the current resolution these states do not play any role. As an independent check, the resulting pion screening masses (listed in Table 1) can be compared with the continuum extrapolated values from N f = 2+1 simulations [48], and good agreement is observed.
Since the QCD pion in vacuum is a stable massive scalar particle, and corresponds to the lightest T = 0 state, the considerations of Sec. 2 apply, and one can use Eq. (2.5) to extract the damping factor associated with this state. Doing so, one obtains where α π and γ π are T -dependent parameters. Note that a damping factor with this functional form was already observed in Ref. [42] for pion states in the quark-meson model, in this case with correlators generated using a functional renormalisation group (FRG) approach 5 . The parameter γ π was shown to control the thermal broadening of the pion states, and hence corresponds to the interaction energy between the pions and the thermal background. For the sub-leading π * contribution it must also be the case that the corresponding damping factor has the form D m π * ,β ( x) = α π * e −γ π * | x| . Therefore, the resolution of two independent exponential contributions to C PS (x 3 ) is consistent with the two lowest-lying particle states in the T = 0 spectrum having exponential damping factors. Applying Eq. (2.4), contributions of this form generate the general spatial correlator behaviour which implies the following relation between the screening masses, vacuum masses, and damping factor exponents for each particle contribution Due to the restoration of Lorentz invariance in the T → 0 limit, the parameters γ i must vanish, and hence as expected for the T = 0 spatial correlator of a massive particle with weight α i (T = 0). Having determined the functional form of the π and π * damping factors, one can now use the corresponding thermal spectral density together with the general representation in Eq. (1.5), to establish how the π and π * states contribute to the pseudo-scalar spectral function ρ PS (ω, p). One finds that which is the two-state generalisation of the result obtained in Ref. [42]. Equation (3.11) therefore represents the spectral function that best describes the lattice data in the range [2x * 3 , L/2]. Substituting Eq. To interpret the physics of this result, let us focus on the zero-momentum spectral function ρ PS (ω) := ρ PS (ω, p = 0), which has the form (3.12) 5 Recent analyses of pion spectral properties in the quark-meson model were also performed in Refs. [49][50][51], based on an FRG analytic continuation approach.
Particularly distinctive characteristics of Eq. (3.12) are that the π and π * contributions have sharp thresholds at their respective T = 0 masses, and each particle contribution has a single peaked behaviour, with the peak maxima located at Mathematically, the thresholds are a direct consequence of the discrete δ(s − m 2 i ) components in the ansatz of Eq. (1.6). It reflects the physical situation that the thermal background state |Ω β is itself composed of pions, and therefore one has to inject an energy of at least m π or m π * in order to create a zero-momentum excitation with π or π * quantum numbers, respectively. The shift of the peaks to energies larger than m π and m π * corresponds to the fact that in addition one has to supply the interaction energy of the particles with the thermal background in order to excite those states with maximal probability. Hence, the thresholds have a microscopic origin in the mass gaps of the constituents of the thermal system, whilst the peak positions are a collective feature of the macroscopic system.
To complete our knowledge of the pseudo-scalar spectral function, we still need values for the coefficient parameters α π and α π * . In principle, these are fixed by the absolute normalisation of the correlation functions. However, this would require renormalised correlators, which are not available at present. Fortunately, by taking the ratio of the fitted coefficients this normalisation cancels, and one obtains the following expression for the relative weight of the π and π * coefficients (3.14) Since the spectral function has the functional form in Eq. (3.12), this implies that it satisfies the integral relation ∞ 0 dω π ω ρ PS (ω) = α π + α π * .
(3.15) By choosing to normalise the spectral function to one, it follows from Eqs. (3.14) and (3.15) that the coefficients in the normalised spectral functionρ PS (ω) can be calculated directly from the fitting parameters asα (3.17) Using these relations, together with the fit parameters and Eq. (3.12), in the left plot of Fig. 2 we evaluateρ PS (ω) with its respective 1σ error band for m π = 140 MeV at different temperatures. In the right plot of Fig. 2 we also showρ PS (ω)/ω 2 , which is often studied in the literature. One sees that the rescaling by ω 2 causes a separation of the π and π * peaks, and their relative heights to change order.

Systematics
Before discussing the physics implications of our results, it is useful to assess the systematic uncertainties involved in their extraction. Starting with the lattice data we analysed here, it is worth repeating that domain wall fermions have good chiral properties at finite lattice spacing, and there are no doublers involved. This is important since the data we consider are not extrapolated to the continuum. The range of lattice spacings considered here (cf. Table 1) is comparable to those used in a recent continuum extrapolated calculation of screening masses [48], so finite lattice spacing effects are expected to be small. Similarly, for the spatial correlators, the box size in units of the largest correlation length (cf. Table 1) is m src π L 12, and based on experience, finite-size effects in the correlators are exponentially small.
A larger quantitative effect is caused by the fact that the bare quark mass is held fixed at am ud = 0.001 in lattice units. Since there are different lattice spacings involved at different temperatures (cf. Table 1), the system is not strictly staying on a so-called line of constant physics. Using the relation m 2 π ∼ m q from chiral perturbation theory this amounts to a roughly 50% uncertainty on the vacuum pion mass, and hence we set m π = 140 ± 70 MeV. By using a linear extrapolation based on the T = 0 excited pion lattice calculations of Ref. [52] we were able to estimate how this lattice uncertainty propagates to the mass of the first excited state, finding that m π * = 1412 ± 38 MeV. In order to establish the effect that these uncertainties have on the spectral function we varied the masses within their uncertainty ranges. Whilst we found that this produced some quantitative modifications of the peak heights on the order of 20%, the overall qualitative characteristics ofρ PS (ω) were not significantly altered. Since the analysed lattice data correspond to N f = 2 QCD, one would also expect the spectral function to be quantitatively modified to some extent once strange quarks are included.
Besides discrete particle components, the spectral function also contains continuous contributions, represented by D c,β ( x, s) in Eq. (2.3). However, since the spatial correlator is directly related to the momentum-dependent spectral function via Eq. (1.4), if these components did provide a significant contribution, particularly at small energies, they should be detectable in the fits. Together with the controllable effects discussed here, the largest systematic uncertainty in the present approach is the decomposition ansatz of Eq. (1.6) proposed in Ref. [28], which encodes the significant influence of the vacuum states on the analytic structure of the spectral function. This ansatz is motivated by various physical considerations, as outlined in Secs. 1 and 3.1, and a general investigation will be the subject of future work. In the context of this study there exists a highly non-trivial test of the consistency of the extracted spectral functions, which is independent of how they are obtained. This will be discussed in Sec. 3.4.

A non-perturbative test of the spectral function
It is apparent from Eqs. (1.2) and (1.4) that the spectral function governs the behaviour of both spatial and temporal correlation functions. This therefore provides a straightforward and nonperturbative test of the consistency of any extracted spectral function: the spectral function determined from spatial correlators must predict the corresponding temporal ones, and vice versa. In our case, Ref. [53] provides lattice data of the temporal correlator C PS (τ ) := C PS (τ, p = 0) at T = 220 MeV with the same lattice action and parameters as the spatial correlator data analysed here. From Eqs. (1.2) and (1.5) it follows that and hence the extracted thermal spectral density in Eq. (3.10), from which our ρ PS is calculated, can be used to predict the form of C PS (τ ), and compared with the data of Ref. [53].
The lattice temporal correlator data from Ref. [53] are normalised as C PS (τ = τ * ) = 1, where τ * = 1/(T N τ ) is lattice separation one. In order to perform a comparison we need to ensure that both the spatial and temporal lattice correlators are normalised such that the zero-momentum spectral function has unit normalisation, as in Fig. 2. To do so, consider the derivative of C PS (τ ) for τ → 0 + , which due to Eq. (3.18) gives Rescaling the temporal lattice data by −2 dC PS dτ (τ → 0 + ) therefore ensures the proper normalisation, and we denote the corresponding correlator byĈ PS (τ ). A final obstacle is the estimation of the derivative itself, since the data does not extend down to τ = 0. For this purpose we performed interpolations of the data based on Hermite polynomials and cubic splines, which were then extrapolated to τ = 0. Using the spread in results across the different interpolation methods we estimated the error on the derivative to be of the order of 10%. In Fig. 3 we plot theĈ PS (τ ) prediction from the spatial correlator data against the appropriately rescaled temporal lattice data.
We see that theĈ PS (τ ) prediction matches well with the rescaled lattice data. As τ becomes small the prediction begins to underestimate the lattice points. This is entirely expected, sinceρ PS (ω) only contains information about the two lowest-lying states π and π * , whereas for small values of τ the higher-excited states such as the π * * exert a greater influence over the behaviour ofĈ PS (τ ). Indeed, we observe quantitative agreement for τ 0.0007 MeV −1 ≈ m −1 π * , which is the smallest screening length resolved by the spatial correlator. Overall, the analysis in this section provides robust evidence that the spectral function in Eq. (3.11) is fully consistent with the available stateof-the-art lattice correlator data from Refs. [45,53] on length scales m −1 π * at T = 220 MeV.

Dynamical implications
Apart from the appearance of thresholds related to the vacuum spectrum, which we discussed already, the most prominent feature in Fig. 2 (left) is that the pion component of the spectral function gives a pronounced peak at the lowest temperature T = 220 MeV, which corresponds to approximately 1.2 T pc . As the temperature increases, this peak eventually disappears due to two combined effects. Firstly, there is a broadening of the peak, as expected from the increasing collision rates in the medium. Secondly, the peak location, defined in Eq. (3.13), begins to exceed the π * threshold at ω = m π * , causing the pion peak to effectively disappear. The π * features a distinct peak up to around T = 660 MeV, although the inclusion of higher-excited pion states, which are not captured by the lattice data, would modify this picture. Nevertheless, because the contributions of the higher-excited states have increasingly larger thresholds, their behaviour would not affect the pion contribution toρ PS (ω). The emerging physical picture is that of sequential melting of the pion and its excitations with increasing temperature, and that this process happens gradually and at temperatures significantly larger than the chiral symmetry restoration scale T pc .
Note that the general characteristics of Fig. 2 (right) are similar to those seen in some earlier reconstruction studies [54,55], where instead the spectral function was extracted from data of temporal correlators C Γ (τ ) via Eq. (1.2). In particular, these studies also see separated peaks of decreasing size for larger values of ω, even when T > T pc . A contrast with Fig. 2 though is the absence of discrete thresholds, which is most likely because the reconstructed spectral functions are constrained to be smooth. Further meaningful comparisons can also be made with non-perturbative studies in the O(4) model [56][57][58], and analyses of the pseudo-scalar spectral function in low-temperature QCD [59][60][61][62]. In the O(4) model the spectral functions also appear to possess peaks above the transition temperature, as in Fig. 2, and the peak locations are related to the screening masses [56]. In QCD analyses [59][60][61][62] the pion screening mass m scr π is similarly found to increase with temperature, and the notion of a quasi-particle mass is introduced in the low-temperature regime, motivated by the observation that the spatial correlator has a purely exponential decay at large distances. This exponential behaviour is explained by the presence of a discrete delta component in the spectral function whose argument is proportional to m scr π at p = 0. However, in our case we demonstrate that an exponential spatial correlator can also be achieved with exponential damping factors (cf. Eq. (3.7)), and a resulting spectral function (Eq. (3.11)) with a quite different analytic structure.
Finally, we note that the picture of pions melting at temperatures significantly above the effective chiral symmetry restroration scale is fully consistent with the expectations based on chiral spin symmetry, which was observed to be approximately realised in a temperature range T pc T 3T pc on the same lattice configurations [45,53]. As explained in those studies, chiral spin symmetry can only emerge dynamically in a regime where the chromoelectric quark gluon interactions dominate the quantum effective action of QCD. Refs. [45,53] therefore concluded that in this intermediate temperature regime the thermal medium consists of hadron-like states where chiral symmetry is nearly restored, and quarks are still effectively bound. The spectral function of Fig. 2 fully supports this picture. In Ref. [63] it was discussed how this chiral spin symmetric band smoothly extends to finite baryon density, where it may connect to quarkyonic matter in the cold and dense regime [64,65], which is expected to contain chirally symmetric (parity-doubled) baryons. In this context, it would be most interesting to extend the present study to lattice correlators evaluated at imaginary baryon chemical potential, which can be computed without a sign problem, and understand the influence of small baryon densities on the spectral functions.

The scalar spectral function and chiral symmetry restoration
In Ref. [45] the correlators for several other iso-vector meson operators are given for the temperatures listed in Table 1, in particular the Lorentz scalar, which contains contributions from the a 0 meson. The physical a 0 has a vacuum mass of roughly 980 MeV, but this value will be lower in N f = 2 QCD because the a 0 is expected to have a sizeable ss component [46]. Moreover, the a 0 is unstable under strong decays in vacuum, and the effect of this instability also needs to be taken into account. For these reasons, we postpone the analysis of the other iso-vector meson correlators to future work. Nevertheless, a few qualitative observations can already be made. Due to the effectively restored U(1) A symmetry for T 220 MeV, the scalar and pseudo-scalar spatial correlators are degenerate within the available accuracy of the data [45]. One might therefore expect the corresponding spectral function ρ S to be identical to ρ PS in this temperature regime. This, however, would only be the case for a theory with exact chiral symmetry, whereas for physical QCD, some care is in order.
Although scalar and pseudo-scalar screening masses are approximately degenerate for T > T pc , Eq. (3.8) demonstrates that this can be achieved even if the damping factor exponents γ i (T ), and hence the spectral functions, are different, due to the unequal vacuum masses. If ρ S also has the same qualitative properties as ρ PS , this implies that the low-energy behaviour of these spectral functions also cannot be equal since the pion and a 0 components will have different thresholds. It is well known from reconstruction studies that enormous accuracy of the correlation functions is necessary to resolve details of the spectral function, due to its folding with the thermal kernel in Eq. (1.2). In other words, slightly different spectral functions may well result in approximately degenerate correlators. This leads to the following natural picture: starting with different vacuum properties of the pion and a 0 , the spectral functions ρ PS and ρ S will begin to move towards one another for higher temperatures on account of the thermal modifications encoded in their respective damping factors. In the temperature region T T pc investigated here, one still expects some differences, in particular around the pion threshold region, but plausibly these differences will only integrate to a small deviation between the respective correlators. As the temperature is further increased, the vacuum properties of the particles will be overpowered by the thermal modifications, and the spectral functions will effectively merge in the limit T → ∞, as expected. The nonperturbative approach presented in this work therefore offers a viable path to unveil the details of chiral symmetry restoration and deconfinement in QCD.

Conclusions
Understanding the non-perturbative structure of Euclidean correlation functions at finite temperature is essential for obtaining a correct physical interpretation of how particles behave in the presence of a thermal medium. In this work, we utilise the general spectral representation satisfied by thermal two-point functions to derive a corresponding representation for the spatial correlator C(x 3 ) of real scalar fields. We find that if there is a mass gap to a stable vacuum particle state comprising the bulk of the medium, and a further gap to the continuum threshold, the large-x 3 behaviour of the correlator is dominated by discrete particle contributions, as expected from basic physical arguments. In such a regime, the spectral function can be extracted directly, avoiding the well-known inverse problem.
Using these analytic results, we analyse lattice QCD data for the light-quark pseudo-scalar meson operator in the temperature range 220 − 960 MeV, and extract the form of the corresponding spectral function ρ PS (ω, p). As a non-trivial test, we demonstrate that the extracted spectral function reproduces the corresponding temporal lattice correlator data for T = 220 MeV. We find that the pion π and its first-excited state π * dominate the behaviour of ρ PS , and that the π is clearly distinguishable in a range of temperatures above the chiral pseudo-critical temperature T pc . Similar qualitative features should also be realised for the a 0 meson and its excitations, due to effective chiral symmetry restoration. Explicit relations between the vacuum masses, screening masses, and spectral function parameters are derived, and insights into the distinction between decays due to intrinsic vacuum instability and thermal broadening are also discussed. These findings suggest that non-perturbative effects continue to play a significant role above T pc even for hadronic states composed of light quarks, and are consistent with the expectations based on the approximate realisation of chiral spin symmetry in this temperature regime. Although this study focused on the spectral properties of correlators involving pseudo-scalar meson operators, our approach can in principle be generalised to meson states with higher spin, other hadronic states, as well as to regimes with non-vanishing baryon density. This work represents a step towards the analytic characterisation of non-perturbative in-medium effects in QCD.