Light baryons below and above the deconfinement transition: medium effects and parity doubling

TL;DR

This study uses lattice QCD simulations to analyze how baryon properties, especially parity doubling, change across the deconfinement transition, revealing temperature-dependent mass behaviors and chiral symmetry restoration effects.

Contribution

It provides the first detailed lattice QCD analysis of parity doubling and medium effects on baryons across the deconfinement transition with multiple temperature points.

Findings

Positive-parity masses are temperature independent below Tc.

Negative-parity masses decrease as temperature approaches Tc.

Parity doubling occurs above Tc, indicating chiral symmetry restoration.

Abstract

We study what happens to the Nucleon, Delta and Omega baryons in the hadronic gas and the quark-gluon plasma, with particular interest in parity doubling and its emergence as the plasma is heated. This is done using simulations of lattice QCD, employing the FASTSUM anisotropic Nf=2+1 ensembles, with four temperatures below and four above the deconfinement transition temperature. Below Tc we find that the positive-parity groundstate masses are largely temperature independent, whereas the negative-parity ones are reduced considerably as the temperature increases. This may be of interest for heavy-ion phenomenology. Close to the transition, the masses are nearly degenerate, in line with the expectation from chiral symmetry restoration. Above Tc we find a clear signal of parity doubling in all three channels, with the effect of the heavier s quark visible.

Tables

Table 1. Table 1: Ensembles used in this work. The lattice size is N s 3 × N τ superscript subscript 𝑁 𝑠 3 subscript 𝑁 𝜏 N_{s}^{3}\times N_{\tau} , with the temperature T = 1 / ( a τ ​ N τ ) 𝑇 1 subscript 𝑎 𝜏 subscript 𝑁 𝜏 T=1/(a_{\tau}N_{\tau}) . The available statistics for each ensemble is N cfg × N src subscript 𝑁 cfg subscript 𝑁 src N_{\rm cfg}\times N_{\rm src} . The sources were chosen randomly in the four-dimensional lattice. The spatial lattice spacing a s = 0.1227 ​ ( 8 ) subscript 𝑎 𝑠 0.1227 8 a_{s}=0.1227(8) fm, the inverse temporal lattice spacing a τ − 1 = 5.63 ​ ( 4 ) subscript superscript 𝑎 1 𝜏 5.63 4 a^{-1}_{\tau}=5.63(4) GeV, and the renormalised anisotropy ξ = a s / a τ = 3.5 𝜉 subscript 𝑎 𝑠 subscript 𝑎 𝜏 3.5 \xi=a_{s}/a_{\tau}=3.5 .

$N_s$	$N_{\tau }$	$T$ [MeV]	$T/T_c$	$N_{\mathrm{src}}$	$N_{\mathrm{cfg}}$
24	128	44	0.24	16	139
24	40	141	0.76	4	501
24	36	156	0.84	4	501
24	32	176	0.95	2	1000
24	28	201	1.09	2	1001
24	24	235	1.27	2	1001
24	20	281	1.52	2	1000
24	16	352	1.90	2	1001

Table 2. Table 2: Groundstate masses m ± subscript 𝑚 plus-or-minus m_{\pm} in both parity sectors in the N 𝑁 N , Δ Δ \Delta and Ω Ω \Omega channels below T c subscript 𝑇 𝑐 T_{c} . Estimates for statistical and systematic uncertainties are included. The final column shows the T = 0 𝑇 0 T=0 values in nature Agashe:2014kda . Note that there are several candidates for m − Ω subscript superscript 𝑚 Ω m^{\Omega}_{-} . The difference δ ​ m N , Δ , Ω 𝛿 subscript 𝑚 𝑁 Δ Ω \delta m_{N,\Delta,\Omega} is defined as δ ​ m = m − − m + 𝛿 𝑚 subscript 𝑚 subscript 𝑚 \delta m=m_{-}-m_{+} and the dimensionless ratio δ N , Δ , Ω subscript 𝛿 𝑁 Δ Ω \delta_{N,\Delta,\Omega} as δ = ( m − − m + ) / ( m − + m + ) 𝛿 subscript 𝑚 subscript 𝑚 subscript 𝑚 subscript 𝑚 \delta=(m_{-}-m_{+})/(m_{-}+m_{+}) .

$T/T_c$	0.24	0.76	0.84	0.95	PDG ($T=0$)
$m_{+}^N$ [MeV]	1159(13)	1192(39)	1169(53)	1104(40)	939
$m_{-}^N$ [MeV]	1778(52)	1628(104)	1425(94)	1348(83)	1535(10)
$m_{+}^{\Delta }$ [MeV]	1459(58)	1521(43)	1449(42)	1377(37)	1232(2)
$m_{-}^{\Delta }$ [MeV]	2138(117)	1898(106)	1734(97)	1526(74)	1710(40)
$m_{+}^{\Omega }$ [MeV]	1661(21)	1723(32)	1685(37)	1606(43)	1672.4(0.3)
$m_{-}^{\Omega }$ [MeV]	2193(30)	2092(91)	1863(76)	1576(66)	2250–2380–2470
$\delta m_N$ [MeV]	619(54)	436(111)	256(108)	244(92)	596(10)
$\delta m_{\Delta }$ [MeV]	679(131)	377(114)	285(106)	149(83)	478(40)
$\delta m_{\Omega }$ [MeV]	532(37)	369(96)	178(85)	-30(79)	578–708–798
${\delta }_N$	0.211(19)	0.155(35)	0.099(40)	0.100(35)	0.241(1)
${\delta }_{\Delta }$	0.189(37)	0.110(31)	0.089(31)	0.051(28)	0.162(14)
${\delta }_{\Omega }$	0.138(10)	0.097(23)	0.050(23)	-0.009(25)	0.147–0.175–0.192

Table 3. Table 3: Operator dependence of m ± N superscript subscript 𝑚 plus-or-minus 𝑁 m_{\pm}^{N} at T / T c = 0.24 , 0.76 𝑇 subscript 𝑇 𝑐 0.24 0.76 T/T_{c}=0.24,0.76 .

$T/T_c$	operator	4	5	6
0.24	$m_{+}^N$ [MeV]	1157(13)	1156(13)	1156(13)
	$m_{-}^N$ [MeV]	1779(52)	1824(48)	1934(101)
0.76	$m_{+}^N$ [MeV]	1192(39)	1190(45)	1212(47)
	$m_{-}^N$ [MeV]	1628(104)	1698(106)	1548(201)

Equations

$$G^{\alpha {\alpha }^{\prime }}(x)=⟨O^{\alpha }(x){\overline{O}}^{{\alpha }^{\prime }}(0)⟩,$$

$$O_N^{\alpha }(x)$$

$$O_{\Delta ,i}^{\alpha }(x)$$

$$O_{\Omega ,i}^{\alpha }(x)$$

$$C^{†}C=1,{\gamma }_{\mu }^{{}^T}=-C{\gamma }_{\mu }{C}^{-1},C^{{}^T}=-C^{-1},$$

$$\mathcal{P}\psi (x){\mathcal{P}}^{-1}={\gamma }_4\psi (Px),P=\text{mbox}diag(-1,-1,-1,1).$$

$$\mathcal{P}O(x){\mathcal{P}}^{-1}={\gamma }_{4}O(Px).$$

$$P_{\pm }=\frac{1}{2}\left(1\pm {\gamma }_4\right),O_{\pm }(x)=P_{\pm }O(x),$$

$$\mathcal{P}{O}_{\pm }(x){\mathcal{P}}^{-1}=\pm O_{\pm }(Px).$$

$$\mathcal{C}\psi {\mathcal{C}}^{-1}\equiv {\psi }^{(c)}=C^{-1}{\bar{\psi }}^{{}^T},\mathcal{C}\bar{\psi }{\mathcal{C}}^{-1}\equiv {\bar{\psi }}^{(c)}=-{\psi }^{{}^T}C.$$

$$O^{(c)}=C^{-1}{\overline{O}}^{{}^T},{\overline{O}}^{(c)}=-O^{{}^T}C.$$

$$G^{\alpha {\alpha }^{\prime }}(\tau ,\mathbf{p})=\int d^{3}x{e}^{-i\mathbf{p}\cdot \mathbf{x}}{G}^{\alpha {\alpha }^{\prime }}(\tau ,\mathbf{x}),$$

$$G^{\alpha {\alpha }^{\prime }}(i{\omega }_n,\mathbf{p})={\int }_{-\infty }^{\infty }\frac{d\omega }{2\pi }\frac{{\rho }^{\alpha {\alpha }^{\prime }}(\omega ,\mathbf{p})}{\omega -i{\omega }_n}.$$

$${\rho }^{\alpha {\alpha }^{\prime }}(\omega ,\mathbf{p})=2\mathrm{Im}{G}^{\alpha {\alpha }^{\prime }}(i{\omega }_n\to \omega +iϵ,\mathbf{p}),$$

$${\rho }^{\alpha {\alpha }^{\prime }}(x)=⟨\{O^{\alpha }(x),{\overline{O}}^{{\alpha }^{\prime }}(0)\}⟩,$$

$$G^{\alpha {\alpha }^{\prime }}(\tau ,\mathbf{p})={\int }_{-\infty }^{\infty }\frac{d\omega }{2\pi }K(\tau ,\omega ){\rho }^{\alpha {\alpha }^{\prime }}(\omega ,\mathbf{p}),$$

$$K(\tau ,\omega )=T\sum _n\frac{e^{-i{\omega }_n\tau }}{\omega -i{\omega }_n}=\frac{e^{-\omega \tau }}{1+e^{-\omega /T}}=e^{-\omega \tau }\left[1-n_F(\omega )\right],$$

$$K(1/T-\tau ,\omega )=K(\tau ,-\omega ).$$

$$K(\tau ,\omega )=\frac{1}{2}\left[K_{\mathrm{e}}(\tau ,\omega )+K_{\mathrm{o}}(\tau ,\omega )\right],$$

$$\begin{array}{rl}{K}_{\mathrm{e}}(\tau ,\omega )& =\frac{\cosh (\omega \tilde{\tau })}{\cosh (\omega /2T)}=\left[1-n_F(\omega )\right]e^{-\omega \tau }+n_F(\omega )e^{\omega \tau },\\ K_{\mathrm{o}}(\tau ,\omega )& =-\frac{\sinh (\omega \tilde{\tau })}{\cosh (\omega /2T)}=\left[1-n_F(\omega )\right]e^{-\omega \tau }-n_F(\omega )e^{\omega \tau },\end{array}$$

$$K_{\mathrm{boson}}(\tau ,\omega )$$

$$G^{\alpha {\alpha }^{\prime }}(x)$$

$${\rho }^{\alpha {\alpha }^{\prime }}(x)$$

$$G_{\mu }(x)=\frac{1}{4}\text{mbox}tr{\gamma }_{\mu }G(x),G_m(x)=\frac{1}{4}\text{mbox}trG(x),$$

$$G_{\pm }(x)$$

$${\rho }_{\pm }(x)$$

$$G_{\pm }(\tau ,\mathbf{p})={\int }_{-\infty }^{\infty }\frac{d\omega }{2\pi }K(\tau ,\omega ){\rho }_{\pm }(\omega ,\mathbf{p}).$$

$$\rho (p)$$

$$G^{>}(x-x^{\prime })=⟨O(x)\overline{O}(x^{\prime })⟩,G^{<}(x-x^{\prime })=-⟨\overline{O}(x^{\prime })O(x)⟩.$$

$${\rho }_4(p)=\left(1+e^{-p^0/T}\right)\int d^{4}x{e}^{-ip\cdot x}\frac{1}{4}\text{mbox}tr⟨O(x)O^{†}(0)⟩.$$

$${\rho }_4(p)$$

$$\text{mbox}tr{P}_{\pm }O(x)\overline{O}(0)=\pm \text{mbox}tr{O}_{\pm }(x)O_{\pm }^{†}(0),$$

Full text

††arxiv: 1703.09246

1 Introduction

In the past decade the study of the quark-gluon plasma at the Large Hadron Collider at CERN and the Relativistic Heavy Ion Collider at BNL has matured into a quantitative area of research, in which more detailed questions can be asked and answered – see e.g. Refs. Braun-Munzinger:2015hba ; Armesto:2015ioy and references therein. One topic of interest concerns the changes to the spectrum of QCD, which are expected as hadrons are immersed in a hadronic gas at temperatures below the deconfinement transition, and in the quark-gluon plasma (QGP) at higher temperatures. This has been especially important for quarkonium, bound states of a heavy quark and anti-quark, as their melting/survival pattern can act as a thermometer for the temperatures reached in these collisions. Indeed, both the LHC Chatrchyan:2012lxa ; Abelev:2014nua and RHIC Adare:2014hje have reported clear suppression patterns for bottomonium states at high temperature. Ref. Andronic:2015wma contains a recent comprehensive review and Ref.  Aarts:2016hap a discussion of open questions.

For light hadrons on the other hand, the emphasis has been on the statistical properties of the hadrons emerging from the system and on the dilepton spectrum Braun-Munzinger:2015hba . Dileptons are predominantly produced by the decay of vector mesons and hence properties of their spectrum provide a connection with chiral symmetry and its restoration at high temperature. This observation has led to substantial activity on the role of chiral symmetry at finite temperature in the mesonic sector Rapp:1999ej .111Note that throughout this paper chiral symmetry will refer to SU(2)A chiral symmetry, which is spontaneously broken in the vacuum (and explicitly by nonzero quark masses).

Due to the nature of the thermal transition in QCD, studies using lattice QCD can provide important nonperturbative insight. Probably the cleanest signal with respect to chiral symmetry comes from the analysis of mesonic screening masses, which are relatively easy to compute in lattice simulations, see e.g. Ref. Cheng:2010fe , even though their relation to phenomenologically relevant quantities is not immediately clear (see however Ref. Brandt:2014uda ). Direct computation of spectral quantities in a medium, such as thermal masses, is considerably harder, due to the need to consider analytical continuation on lattices with a finite temporal extent. Recent interesting work on the pion in the hadronic gas can be found in Refs. Brandt:2014qqa ; Brandt:2015sxa . The vector meson correlator has been analysed extensively, not only due its role in the dilepton rate but also in the context of the electrical conductivity and the charge diffusion coefficient Aarts:2007wj ; Ding:2010ga ; Amato:2013naa ; Aarts:2014nba ; Brandt:2015aqk ; Ding:2016hua . Concerning quarkonia, both charmonium Umeda:2002vr ; Asakawa:2003re ; Datta:2003ww ; Aarts:2007pk ; Ohno:2011zc ; Borsanyi:2014vka and, more recently, bottomonium Aarts:2011sm ; Aarts:2013kaa ; Aarts:2014cda ; Kim:2014iga have been studied on the lattice.

Surprisingly, even though light baryons are sensitive to chiral symmetry and play an important role in the analysis of heavy-ion data, corresponding studies in the baryonic sector are very limited. In the context of lattice QCD, baryon screening masses in a gluonic medium were studied a long time ago in Refs. DeTar:1987xb ; DeTar:1987ar and, at small baryon chemical potential, in Ref. Pushkina:2004wa . More recently, screening and temporal correlators were analysed in Ref. Datta:2012fz . All these studies were carried out in the quenched approximation.

In this work we aim to improve this situation substantially. We study the $N$ (nucleon), $\Delta$ and $\Omega$ baryons, employing simulations with $N_{f}=2+1$ light flavours, with four temperatures below and four above the transition. This allows us to study the properties of these baryons and in particular in-medium modification in the hadronic gas. Chiral symmetry is closely linked to parity doubling and we analyse the emergence of parity doubling as the transition is approached. We find a qualitative difference in the response to the increasing temperature between positive- and negative-parity baryons, which may be of interest for heavy-ion phenomenology. We also contrast the behaviour in the $\Omega$ channel with the $N$ and $\Delta$ channel, to see the effect of the heavier $s$ quark.

This paper is organised as follows. In Sec. 2 we summarise the relations between baryon correlators and spectral functions, emphasising the differences with the mesonic case. We discuss positivity of the spectral functions, the role of charge conjugation, and the connection between chiral symmetry and parity doubling, both at $\mu=0$ and $\mu\neq 0$. Sec. 3 contains details of our lattice computation. The main results of our study are given in Sec. 4: we analyse the euclidean correlators and draw conclusions for both the hadronic gas and the quark-gluon plasma. These results are supported by the spectral function analysis of Sec. 5. The final section summarises and contains an outlook. We note here that our previous work in the nucleon sector, with limited statistics, can be found in Ref. Aarts:2015mma and preliminary results have appeared in Refs. Aarts:2015xua ; Aarts:2016ioq ; Aarts:2016kvt .

2 Baryonic correlators and spectral functions

We start with a brief discussion of baryonic operators and spectral relations for fermionic two-point functions. While for mesonic (bosonic) correlators the type of relations discussed below are very well known Bellac:2011kqa , for fermionic ones this is slightly less so. Moreover, it allows us to discuss how parity doubling manifests itself in correlators and spectral functions.

2.1 Baryonic operators

We consider two-point functions of fermionic operators, of the form

[TABLE]

where $\alpha,\alpha^{\prime}$ are Dirac indices and $\overline{O}=O^{\dagger}\gamma_{4}$.222We follow the conventions in Ref. Gattringer:2010zz and use euclidean gamma-matrices, $\gamma_{\mu}^{\dagger}=\gamma_{\mu}=\gamma_{\mu}^{-1}$, with $\mu=1,\dots,4$, and $\gamma_{5}^{\dagger}=\gamma_{5}=\gamma_{1}\gamma_{2}\gamma_{3}\gamma_{4}$. The simplest annihilation operators for the nucleon, $\Delta$ and $\Omega$ baryons are respectively Leinweber:2004it ; Edwards:2004sx

[TABLE]

where $C$ corresponds to the charge conjugation matrix, satisfying

[TABLE]

and hence $\gamma_{5}^{{}^{T}}=C\gamma_{5}C^{-1}$. We note here that as written Eq. (3) describes the charged $\Delta^{+}(uud)$ channel. However, since QED interactions are not incorporated and the two light quarks are taken to be degenerate (isospin limit), the operator is also relevant for the neutral $\Delta^{0}(ddu)$ channel. The $\Delta^{++}(uuu)$ and $\Delta^{-}(ddd)$ states are in principle described by an operator of the form (4), with $s\to u,d$, but again in the degenerate limit one can show that Wick contractions coming from the latter are identical to the ones derived from Eq. (3).

Under parity, elementary quark fields transform as

[TABLE]

It is straightforward to verify that this property is inherited by the baryonic operators,

[TABLE]

Hence one may introduce parity projectors and operators via

[TABLE]

such that

[TABLE]

We refer to $O_{\pm}$ as positive- and negative-parity operators.

Similarly, under charge conjugation quark fields transform as

[TABLE]

Again, this is inherited by the baryonic operators, and

[TABLE]

2.2 Spectral relations

We now derive some general spectral relations and properties of the two-point functions $G^{\alpha\alpha^{\prime}}(x)$. We work in spatial momentum space,

[TABLE]

and $\tau$ denotes the euclidean time, $0\leq\tau<1/T$, with $T$ the temperature. Fermionic fields and operators satisfy anti-periodic boundary conditions in euclidean time. A Fourier transform yields the correlator as a function of the fermionic Matsubara frequencies $\omega_{n}=(2n+1)\pi T$, $n\in\mathbb{Z}$, which can be written as a spectral integral

[TABLE]

The spectral function $\rho^{\alpha\alpha^{\prime}}(\omega,{\mathbf{p}})$ is then given by twice the imaginary part of the retarded Green function,

[TABLE]

or, in terms of the operators, by

[TABLE]

as always Bellac:2011kqa . Transforming back to euclidean time yields the integral relation

[TABLE]

with the kernel, for $0<\tau<1/T$,

[TABLE]

which can e.g. be shown by contour integration Bellac:2011kqa . Here $n_{F}(\omega)=1/(e^{\omega/T}+1)$ is the Fermi-Dirac distribution. We note that $K(\tau,\omega)$ is neither even nor odd, but satisfies

[TABLE]

A decomposition of the kernel in terms of its even and odd parts yields

[TABLE]

with

[TABLE]

where $\tilde{\tau}=\tau-1/(2T)$. Note that the normalisation is such that all kernels reduce to $e^{-\omega\tau}$ in the zero-temperature limit (for positive $\omega$). These kernels should be contrasted with the kernel appearing in bosonic spectral relations,

[TABLE]

where $n_{B}(\omega)=1/(e^{\omega/T}-1)$ is the Bose-Einstein distribution. The different denominators, $\cosh(\omega/2T)$ versus $\sinh(\omega/2T)$, reflect the quantum statistics. Note that as a consequence the problems associated with the singular behaviour of the bosonic kernel, $K_{\rm boson}(\tau,\omega)\to 2T/\omega$ as $\omega\to 0$, relevant for transport Aarts:2002cc , are absent in the fermionic case.

In order to resolve the Dirac indices, we use the decomposition (other tensor structures will not appear in our application)

[TABLE]

such that

[TABLE]

where the trace is over the Dirac indices, and similarly for $\rho_{\mu,m}$.

Below we will specialise to zero spatial momentum, for which $G_{i}$ and $\rho_{i}$ vanish. It is convenient to combine the two remaining components with the help of positive- and negative-parity projectors (8) as

[TABLE]

related via

[TABLE]

We will now prove a number of properties of $\rho_{\pm}(x)$ and $\rho_{4,m}(x)$.

2.3 Positivity

We start with positivity: we will show333Note that we use the notation $p=(p^{0},{\mathbf{p}})$ with $p^{0}=\omega$. that $\pm\rho_{\pm}(p),\rho_{4}(p)\geq 0$ for all $\omega$, while $\rho_{m}(p)$ does not have a definite sign, even when restricting to $\omega\gtrless 0$.

Suppressing Dirac indices and using the KMS condition Bellac:2011kqa , valid in thermal equilibrium, we can write

[TABLE]

where $G^{\lessgtr}$ are the usual Wightman functions Bellac:2011kqa ,

[TABLE]

We first consider $\rho_{4}(p)$ and take the trace with $\gamma_{4}$. This yields

[TABLE]

To proceed, we use the Heisenberg representation $O(x)=e^{-ix\cdot K}O(0)e^{ix\cdot K}$ and insert complete sets of eigenstates $|n\rangle$ of the translation operator $K^{\mu}$ (here $K^{0}$ is the Hamiltonian $H$ with eigenvalues $k^{0}_{n}$). Recalling that the expectation value denotes the thermal average with Boltzmann weight $e^{-H/T}/Z$, we find, after some rearrangement,

[TABLE]

where we have written the Dirac index $\alpha$ explicitly again. It is easy to see that the terms added within the summation are nonnegative and hence we arrive at positivity: $\rho_{4}(p)\geq 0$ for all $p$.

Next we consider $\rho_{\pm}(p)$ and take the trace with $P_{\pm}$. We now encounter

[TABLE]

where we used $\gamma_{4}=P_{+}-P_{-},P_{\pm}^{2}=P_{\pm},P_{+}P_{-}=0$, cyclicity of the trace, and Eq. (8). Proceeding as above then yields

[TABLE]

i.e. we find positivity of the spectral functions $\pm\rho_{\pm}(p)\geq 0$ for all $p$.

Positivity of $\rho_{4}(p)$ also follows from $\rho_{4}(p)=[\rho_{+}(p)-\rho_{-}(p)]/4\geq 0$; on the other hand, $\rho_{m}(p)=[\rho_{+}(p)+\rho_{-}(p)]/4$ does not have a definite sign and is indeed not sign-definite for $\omega\gtrless 0$, already at leading order in perturbation theory Praki:2015yua ; Praki:inprep .

To contrast, we note that for bosonic operators the spectral decomposition takes the form as above, but with a minus sign between the two thermal factors Bellac:2011kqa . In addition, if the operator satisfies $J^{\dagger}=\pm J$,444This is e.g. the case for mesonic operators of the form $J=\bar{\psi}\Gamma\psi$, where $\Gamma$ is a Dirac matrix selecting the channel, since $J^{\dagger}=\bar{\psi}\gamma_{4}\Gamma^{\dagger}\gamma_{4}\psi=\pm J$ Karsch:2003wy ; Aarts:2005hg . it follows that the corresponding spectral function $\rho_{B}(p)$ is odd under $\omega\to-\omega$, and $\omega\rho_{B}(\omega,{\mathbf{p}})\geq 0$. This can be seen by swapping $n\leftrightarrow m$ in the summation. For the fermionic operators we consider here, this argument does not apply, since $O^{\dagger}\neq\pm O$. Hence in general fermionic spectral functions are neither even nor odd.

2.4 Charge conjugation

Next we relate, in the case of vanishing baryon chemical potential (or baryon density), positive- and negative-parity correlators and spectral functions, i.e. we show that

[TABLE]

We follow Ref. Karsch:2009tp , where this is demonstrated at the level of the single-quark propagator. Here we consider baryonic (or fermionic in general) operators, transforming under charge conjugation as in Eq. (11). We assume isotropy, i.e. invariance under ${\mathbf{p}}\to-{\mathbf{p}}$, throughout.

The time-ordered correlation function is given by

[TABLE]

with the imaginary-time-ordered product

[TABLE]

Here the minus (plus) sign applies to fermionic (bosonic) operators.

At zero chemical potential, thermal expectation values are invariant under charge conjugation. We hence find, suppressing Dirac indices,

[TABLE]

From now on we take $x=(\tau,{\mathbf{x}})$, $x^{\prime}=(0,{\mathbf{x}}^{\prime})$, with $0<\tau<1/T$. Using the cyclicity of thermal expectation values Bellac:2011kqa then gives

[TABLE]

Applying this to Eq. (37), we find, in momentum space,

[TABLE]

We can now take the trace with $P_{\pm}$, which yields

[TABLE]

where we used that

[TABLE]

We have now demonstrated the first relation in Eq. (34). The second relation immediately follows, when using the integral relations (27) as well as the property (18). Physically it reflects that positive-parity states propagate forward in euclidean time, when using $G_{+}$, and backward in time when using $G_{-}$, and vice versa for negative-parity states. In terms of spectral functions, this relates the positive part of the spectrum of $\rho_{+}$ with the negative part of $\rho_{-}$, and again vice versa. Explicitly, if the spectrum is dominated by single groundstates with masses $m_{\pm}$, this implies

[TABLE]

Using the relation $\rho_{\pm}(-p)=-\rho_{\mp}(p)$, we can subsequently note that

[TABLE]

and hence these are even, respectively odd under $\omega\to-\omega$. Their spectral relations hence involve the even and odd kernels $K_{\rm e,o}(\tau,\omega)$ respectively, see Eq. (20). We remark that this only holds when there is no net density, i.e. when the density matrix is invariant under charge conjugation.

2.5 Chiral symmetry and parity doubling

The final relations we derive are for the case of unbroken chiral symmetry. Here we work in the harmonious world of thermal field theory in which chiral symmetry is simply expressed as $\{\gamma_{5},G\}=0$, sidestepping momentarily the issues related to chiral symmetry in realistic lattice QCD computations, to be discussed below.

From the anti-commutation relation of the correlator with $\gamma_{5}$, it immediately follows that

[TABLE]

and hence

[TABLE]

These relations imply that the lattice correlators are symmetric around the centre of lattice ($\tau=1/2T$), the spectral functions are even functions in $\omega$ and that identical spectral information is contained in $\rho_{\pm}(p)$. We refer to this as parity doubling. We emphasise that any of these signatures are equivalent statements of parity doubling.

An alternative proof for two massless flavours goes as follows Gattringer:2010zz . When chiral symmetry is unbroken, i.e. the quarks are massless and chiral symmetry is not broken spontaneously, the theory is unchanged when the following chiral rotation is performed on the quark fields,

[TABLE]

where $\psi=(u,d)^{T}$ and $T_{3}=\sigma_{3}/2$ acts in flavour-space. Choosing $\alpha=\pi$, we get distinct chiral transformations on the spinor fields:

[TABLE]

It is then easy to check that the $N$ and $\Delta$ operators, see Eqs. (2, 3), transform as

[TABLE]

and hence

[TABLE]

In both channels the correlator then transforms as

[TABLE]

which was to be shown.

2.6 Nonzero chemical potential

For completeness, we indicate here how the properties derived above are modified in presence of a nonzero baryon chemical potential $\mu$, such that the Hamiltonian in the Boltzmann weight is changed from $H\to H-\mu Q$, with $Q$ the baryon number.

First we consider positivity. Following the same steps as in Sec. 2.3, in which the KMS condition (31) is modified as $p^{0}\to p^{0}-\mu$, and using that the states $|n\rangle$ are simultaneous eigenstates of $H$ and $Q$ (with eigenvalues $q_{n}$), we arrive at

[TABLE]

and similar for $\rho_{\pm}(p)$. Hence positivity holds, as before.

At nonzero chemical potential, the density matrix is not invariant under charge conjugation, since baryon number changes sign. Therefore invariance is obtained by simultaneously changing $\mu\to-\mu$, which yields the relations

[TABLE]

$G_{4,m}$ are then no longer (anti)symmetric around $\tau=1/2T$, but satisfy instead

[TABLE]

Again explicitly, if the spectrum is dominated by single groundstates, Eq. (42) is modified as

[TABLE]

Finally, in the case of unbroken chiral symmetry, $G_{m}(x)=\rho_{m}(x)=0$ still holds and

[TABLE]

3 Lattice setup

We have computed baryon correlators using the thermal ensembles of the FASTSUM collaboration Amato:2013naa ; Aarts:2014nba ; Aarts:2014cda . These ensembles are generated with $2+1$ flavours of Wilson fermions on an anisotropic lattice, with a smaller temporal lattice spacing, $a_{\tau}<a_{s}$; the renormalised anisotropy is $\xi\equiv a_{s}/a_{\tau}=3.5$. The lattice action used is the Symanzik-improved anisotropic gauge action with tree-level mean-field coefficients and a mean-field-improved Wilson-clover fermion action with stout-smeared links and follows the Hadron Spectrum Collaboration Edwards:2008ja . Details of the action and parameter values can be found in Refs. Aarts:2014cda ; Aarts:2014nba . The choice of masses for the degenerate $u$ and $d$ quarks yields a pion with a mass of $M_{\pi}=384(4)$ MeV Lin:2008pr , which is heavier than in nature, while the strange quark has been tuned to its physical value. Configurations and correlation functions have been generated using the CHROMA software package Edwards:2004sx , via the SSE optimizations when possible McClendon .

We use a fixed-scale approach, in which the temperature is varied by changing $N_{\tau}$, according to $T=1/(a_{\tau}N_{\tau})$. Table 1 gives an overview of the ensembles. Access to the “zero-temperature” configurations ($N_{\tau}=128$) has been kindly provided to us by the Hadron Spectrum Collaboration. An estimate for the pseudo-critical temperature, $T_{c}=185(4)$ MeV, follows from an analysis of the renormalized Polyakov loop Aarts:2014nba and is higher than in nature, due to the large pion mass. Note that there are four ensembles in the hadronic phase and four in the quark-gluon plasma.

Concerning the baryonic correlators, Gaussian smearing Gusken:1989ad has been employed to increase the overlap with the groundstate. In order to have a positive spectral weight, we apply the smearing on both source and sink, i.e.,

[TABLE]

where $A$ is a normalisation factor and $H$ is the spatial hopping part of the Dirac operator. The hopping term contains APE smeared links Albanese:1987ds using $\alpha=1.33$ and one iteration. We tuned the parameters to the values $n=60$ and $\kappa=4.2$, by maximising the length of the plateau for the effective mass of the groundstate at the lowest temperature. Smearing is applied only in the spatial directions, equally to all temperatures and ensembles.

4 Thermal baryon correlators

In this section we present the results for the baryon correlators at all temperatures. Based on the determination of the pseudo-critical temperature $T_{c}$ via the renormalised Polyakov loop, the discussion is organised in terms of the hadronic gas ($T<T_{c}$) and the quark-gluon plasma ($T>T_{c}$). Since the transition is a crossover, it is not immediately obvious at which temperatures light and strange baryons cease to exist.555We remind the reader that the temperatures closest to $T_{c}$ are $T/T_{c}=0.95$ and $1.09$. However, below we will find clear indications that the baryonic bound states are absent at $T/T_{c}=1.09$, in the three channels we consider.

4.1 Hadronic gas

We have computed the baryon two-point functions in the $N$, $\Delta$ and $\Omega$ channels on the lattice ensembles discussed above, at zero spatial momentum ${\mathbf{p}}={\mathbf{0}}$ (we drop the momentum labels from now on). The results are shown in Fig. 1, at all the eight temperatures available. The positive- and negative-parity channels are shown separately, i.e. the negative-parity channel is obtained using Eq. (34), and the correlators are normalised to the first Euclidean time point, $\tau/a_{\tau}=1,N_{\tau}-1$ respectively, such that

[TABLE]

At low temperatures (open symbols), the correlators show exponential decay, indicating the presence of a well-defined groundstate. As the temperature is increased, some temperature dependence on the positive-parity side is observed, but considerably more temperature dependence is visible on the negative-parity side. The correlators naturally bend upwards around the minima, which are, however, not in the centre of the lattice ($\tau/a_{\tau}=N_{\tau}/2$), due to the absence of parity doubling, see e.g. Eq. (42). Above $T_{c}$ (filled symbols), the correlators appear to drop slower than exponential, indicating the absence of a well-separated groundstate.

To analyse this quantitatively, we have fitted the correlators to a combination of simple exponentials, see Eq. (42),

[TABLE]

with $m_{\pm}$ the groundstate masses in both parity channels. While Fig. 1 shows the positive- and negative-parity channels separately, the fit is carried out to the correlator $G_{+}(\tau)$ in one go. Around the minimum of the correlator, one might become susceptible to signal-to-noise problems, but we found this to be relevant at the lowest temperature only. Here we excluded points around the minimum of the correlators from the analysis, based on the quality of the fit and error analysis.666For example, in the nucleon channel we included the intervals $\tau/a_{\tau}=5-40,105-124$. In order to estimate the systematic uncertainties of the four fit parameters, we have considered various Euclidean time intervals and, to suppress contributions from excited states, we have excluded very small times. We used the so-called Extended Frequentist Method Yao:2006px ; Durr:2008zz to carry out the statistical analysis: this method considers all possible variations and weighs the final results according to the obtained $p$-value, which measures how extreme an outcome is, see Refs. Yao:2006px ; Durr:2008zz for more details. In the confined phase we found that it is possible to extract the mass parameters $m_{\pm}$, whereas above $T_{c}$ the exponential fits are no longer adequate, as can be expected in the deconfined phase (see below).

Table 2 lists the results for the masses $m_{\pm}$ in all three channels, at the four temperatures below $T_{c}$. The results are shown in units of MeV, using the estimate for the temporal lattice of $a_{\tau}^{-1}=5.63(4)$ GeV Edwards:2008ja . Also shown are the PDG Agashe:2014kda values at $T=0$. Since our light quarks are heavier than in nature, the groundstate masses in the $N$ and $\Delta$ channels at the lowest temperature are larger as well. The splitting between the positive- and negative-parity groundstate masses, denoted with $\delta m$, is of the right order, however. The strange quark mass is tuned to the physical value Lin:2008pr and the result for the $\Omega_{+}$ mass is consistent with the PDG value (within errors). Surprisingly, the $\Omega_{-}$ particle has not been unambiguously identified in the PDG and there are three candidates. The value we obtain at $T=0.24T_{c}$ seems to favour the candidate with the lowest mass, but a systematic analysis (continuum extrapolation and physical $u$ and $d$ quarks) is necessary to make a more stringent prediction. Our results for the spectrum at the lowest temperature are in agreement with those of the HadSpec collaboration for the positive-parity states Lin:2008pr ; for the negative-parity baryons the masses obtained in Ref. Edwards:2012fx on a smaller spatial lattice ($16^{3}$ instead of $24^{3}$) are somewhat lower, at the $2\sigma$ level.

As the temperature is increased, we find that the groundstate mass in the positive-parity channels is largely unaffected by temperature; the deviation from the results at the lowest temperature is always less than 5%. Very close to $T_{c}$, the values drops slightly below the ones at $T/T_{c}=0.24$. This is further illustrated in Fig. 2, where the data are plotted normalised by $m_{+}$ at the lowest temperature, in the channel under consideration. In the negative-parity channel we observe a stronger temperature dependence, which is remarkably similar in all three channels. Already at $0.75T_{c}$, the masses have dropped noticeably (see again Fig. 2) and this trend continues towards $T_{c}$. Very close to $T_{c}$ the parity channels are nearly degenerate. This is further quantified by the dimensionless ratio

[TABLE]

also included in Table 2. The smaller value of $\delta_{\Omega}$ at all four temperatures is due to both $m_{+}^{\Omega}$ being larger and $\delta m_{\Omega}$ being smaller. Both of these effects are presumably due to the $s$ quark being heavier than the $u$ and $d$ quarks, which makes the contribution to the groundstate mass due to chiral symmetry breaking less important in the $\Omega$ channel.

4.2 Quark-gluon plasma

We now turn to the temperatures above the deconfinement transition. To start, we have considered the same analysis as above, using exponential fits, assuming that the hypothesis of separated well-defined groundstates still holds. The results are shown in Fig. 3, in the $N$ and $\Omega$ channels. We observe a clear qualitative change when going from $T/T_{c}=0.95$ to 1.09 (or reducing $N_{\tau}$ from 32 to 28). The error on the would-be groundstate masses, obtained by combining systematic and statistical uncertainties, is substantially larger, which cannot be simply explained by the reduction in the number of time slices used in the fits. This, and other results presented below, lead us to conclude that bound states are absent at $T/T_{c}=1.09$, both for the light baryons and in the $\Omega$ channel. Hence even though the transition is a crossover, we find that the spectrum changes rather drastically between 0.95 and 1.09$T_{c}$.

We hence focus on the signal for parity doubling, i.e. the emergent degeneracy in the positive- and negative parity channels. Following Ref. Datta:2012fz , we study the ratio

[TABLE]

which approaches 1 in the case that separated groundstates dominate, with $m_{-}\gg m_{+}$, but vanishes in the case of parity doubling. We have previously shown $R(\tau)$ for all temperatures in the nucleon sector Aarts:2016ioq . Here we present the outcome at two selected temperatures in Fig. 4 in the $N$, $\Delta$ and $\Omega$ channels. We note the clear qualitative and quantitative difference: below $T_{c}$ the ratio is significantly different from zero,777We note that by construction $R(\tau)$ approaches zero at the centre of the lattice, $\tau/a_{\tau}=N_{\tau}/2$. while at the highest temperature it is much smaller. It should be emphasised that if chiral symmetry is exactly restored, complete degeneracy in the positive- and negative-parity channels is expected and $R(\tau)=0$. In our lattice simulations such a clear signal cannot be expected for a number of reasons. First of all we use Wilson fermions, which break chiral symmetry at short distances. We have found that smearing suppresses these contributions, yielding a better signal for parity doubling Aarts:2016ioq . Moreover, the quarks are not massless, with the two light flavours heavier than in nature. Hence this explicit symmetry breaking also affects the signal. However, this is expected to become less important at higher temperature, being suppressed as $m_{q}/T$. The effect of the finite quark mass can be seen in the splitting of $R(\tau)$ in Fig. 4 (right) between the $N,\Delta$ channels and the $\Omega$ channel at the highest temperature; this is most likely due to the larger $s$ quark mass.

In order to summarise the results for all temperatures, we show in Fig. 5 the summed ratio

[TABLE]

where $\tau_{n}=a_{\tau}n$ and $\sigma(\tau_{n})$ are the statistical uncertainties, used as weights.

Assuming that $m_{+}<m_{-}$, $R$ lies between 0 and 1, with $R=0$ corresponding to a symmetric correlator and parity doubling. We observe clear crossover behaviour in all three channels. The location of this transition is consistent with $T_{c}$, which has been determined by an analysis of the renormalised Polyakov loop. Hence it is natural to associate the transition with the approximate restoration of chiral symmetry in the quark-gluon plasma and to interpret $R$ as a quasi-order parameter. We also note that the effect is less pronounced in the $\Omega$ channel, due to the larger $s$ quark mass. It will therefore be interesting to study the effect of strangeness on parity doubling. At the highest temperature available, $R>0$ in the $\Omega$ channel; it is expected that the effect of the quark mass will eventually disappear as $m_{s}/T\to 0$.

The $N$ and the $\Delta$ baryon have the same quark content but different spin structure. In the confined phase this results in the mass splittings listed in Table 2. In the positive-parity channel the mass splitting is of the order of 300 MeV at all four temperatures, consistent with the PDG; in the negative-parity channel the mass difference is larger than in the PDG, but so is the uncertainty. In the deconfined phase, however, the quarks are quasi-free and the spin structure may become less important.888We thank Thomas Cohen for raising this question. To investigate this, we show in Fig. 6 the logarithm of the ratio of the $\Delta$ and $N$ correlators. All ratios are normalised by a single constant factor, $G_{\Delta}(0)/G_{N}(0)$ at $T/T_{c}=1.90$. As expected, this ratio is falling exponentially below $T_{c}$, due to the (approximately constant) mass difference between the $N$ and the $\Delta$ baryons (in both parity channels). Above $T_{c}$, however, we observe a flattening of the ratio, approaching 1 at the highest temperature. We interpret this as an approximate degeneracy in the $N$ and $\Delta$ channels at very high temperature, which would be of interest to study further analytically. We also note the qualitative change in the ratio immediately at $T/T_{c}=1.09$, consistent with the observed changes in the spectrum.

5 Thermal baryon spectral functions

5.1 Results

The information in the thermal correlators discussed above is also present in the corresponding spectral functions, via relation (27)

[TABLE]

As is well-known Asakawa:2000tr , a simple inversion of this type of relation, using numerically determined correlators, is not possible. Hence we use the Maximum Entropy Method MEM Asakawa:2000tr ; Bryan , which extremises a combination of the standard likelihood ($\chi^{2}$) function, determined by the data, and an entropy function,

[TABLE]

encoding prior knowledge, via the default model $m(\omega)$. The conditional probability to be extremised is of the form $\exp(-\frac{1}{2}\chi^{2}+\alpha S)$, with $\alpha$ a parameter balancing the relative importance of the data and the prior knowledge. Both $m(\omega)$ and $\alpha$ are further discussed below. In the past 15 years, this method, and related ones, have been used by a number of groups, mostly for mesonic correlators, i.e. charmonium, the dilepton rate and the electrical conductivity, see e.g. Refs. Umeda:2002vr ; Asakawa:2003re ; Datta:2003ww ; Aarts:2007pk ; Aarts:2007wj ; Ding:2010ga ; Ohno:2011zc ; Amato:2013naa ; Borsanyi:2014vka ; Aarts:2014nba ; Brandt:2015aqk ; Ding:2016hua . Applications to bottomonium, in which some simplifications occur, can be found in Refs.  Aarts:2011sm ; Aarts:2013kaa ; Aarts:2014cda ; Kim:2014iga . Here we give the first application to baryons.

Generic details of our implementation can be found in previous work Aarts:2007wj ; Aarts:2011sm ; Aarts:2014cda ; Aarts:2014nba . Here we briefly mention some differences with the bosonic (mesonic) case. We are interested in the spectrum for both positive and negative $\omega$, since $\rho_{-}(\omega)=-\rho_{+}(-\omega)$. Hence the negative part of the spectrum of $\rho_{+}$ informs us of $\rho_{-}$, and vice versa. To bring the spectral relation (68) to a numerically tractable form, we employ a cutoff $-\omega_{\rm max}<\omega<\omega_{\rm max}$, with $a_{\tau}\omega_{\rm max}=3.0$ ($\omega_{\rm max}=16.9$ GeV). The remaining finite interval is discretised using $N_{\omega}=2000$ bins. We have varied both $\omega_{\rm max}$ and $N_{\omega}$ to verify robustness. In the MEM analysis we used all the euclidean-time points, except for the time slices closest to the source and sink. At the lowest temperature, we have left out the points around the minimum of the correlators; this will be further discussed below. As default model, we use a featureless constant, $m(\omega)=m_{0}$, where the value of $m_{0}$ is determined by a fit to the correlation function using $\rho(\omega)=m_{0}$ in Eq. (68). Above $T_{c}$ we have fixed the default model to ensure a similar normalisation for all temperatures. We come back to the choice of default model below as well.

We now discuss the results. We have performed MEM on the normalised correlators $G_{+}(\tau)/(a_{\tau}G_{+}(\tau=0))$ and denote the associated dimensionless spectral functions with $\bar{\rho}(\omega)$. We note that the normalisation only affects the vertical scale but not the $\omega$ dependence. Fig. 7 contains the spectral functions in the three channels below $T_{c}$ (left) and above $T_{c}$ (right). Spectral information for the positive-parity channel can be found at $\omega>0$, whereas $\omega<0$ refers to the negative-parity channel. Below $T_{c}$, the groundstate peaks on the positive-parity side are clearly visible and their positions agree with $m_{+}$, discussed in the previous section. Excited states are suppressed, due to the choice of smearing parameters. Some broadening is observed as the temperature is increased, but given the data and resolution, it is not clear whether this is a physical effect or due to the limitations of MEM. The negative-parity groundstates are visible as well, but are considerably less pronounced. The asymmetry between the positive- and negative-parity sides below $T_{c}$ is, however, clearly visible.

Above $T_{c}$, sharp groundstate peaks are no longer discernible. The broad peaks present above $T_{c}$ are most likely a combination of physical spectral features for deconfined quarks, as seen at very high temperature in perturbation theory Praki:2015yua ; Praki:inprep , and lattice artefacts due to the finite Brillouin zone, similar to in the mesonic case. To make this statement more quantitative would require a repetition of the calculation on finer lattices, which is one of our future aims. Nevertheless, parity doubling manifests itself as $\rho_{+}(\omega)=\rho_{+}(-\omega)$, see Eq. (47). Hence the most important feature here is the emerging symmetry between the positive- and negative-parity sides as the temperature is increased. This is clearly visible for the $N$ and $\Delta$ channels, in which the position and height of the main features become comparable at positive and negative $\omega$. On the other hand, parity doubling is not yet complete in the $\Omega$ channel, as the positive-parity side is still enhanced. Nevertheless, the difference with spectral functions in the confined phase is manifest. This is consistent with the analysis of the correlators above. We note that in these plots we have not shown error bands for clarity; these will be discussed below.

The combined results in all three channels are shown in Fig. 8, at the lowest (left) and highest (right) temperature. The difference between the spectral functions in the confined and deconfined phase is clear. We also note that below $T_{c}$ the negative-parity state is best visible in the nucleon sector.

5.2 Default model and operator dependence

We now discuss some systematic effects in the construction of the spectral functions. We start with the default model dependence. The results above were obtained with a flat default model, $m(\omega)=m_{0}$. We have also used $m(\omega)=m_{0}|\omega|$ and $m(\omega)=m_{0}|\omega|^{3}$, where in each case $m_{0}$ is determined by a fit to the correlation function. The absolute value ensures positivity. In the continuum theory at leading order in weak coupling Praki:2015yua ; Praki:inprep , the spectral functions increase as $|\omega|^{5}$ for large $|\omega|\gg T,m_{q}$, but this behaviour is modified on a finite lattice Praki:2015yua ; Praki:inprep . Results are shown in Fig. 9. The error band indicates the variation with the $\alpha$ parameter using Bryan’s method Bryan and is shown for one default model only, for clarity. We observe that even though the default models are widely different, the resulting spectral functions are consistent within the uncertainty. The second peaks in the confined phase at both $\omega>0$ and $\omega<0$ are presumably a combination of excited states and lattice artefacts. Whether a structure is due to a finite lattice cutoff or represents a physical feature can ultimately be tested by repeating the computation at smaller lattice spacings. One may also test the robustness with regard to the operators used, to which we turn now.

The dependence on the operator and the amount of smearing requires some discussion. In previous studies in the mesonic sector, it has been common to use a fixed local operator of the form $\bar{\psi}\Gamma\psi$, without smearing.999For charmonium, smearing has been employed in Ref. Umeda:2002vr . Locality is well motivated when the problem under investigation is related to a symmetry, such as electromagnetism (electrical conductivity, charge diffusion, dilepton production) and in Refs.  Amato:2013naa ; Aarts:2014nba the conductivity and charge diffusion coefficient were determined using the exactly conserved lattice vector current. For spectral questions at zero temperature, smearing and optimised operators aim to increase the overlap with the ground (or other) state, in such a way that the spectrum remains invariant, but spectral weight is redistributed. On the other hand, at finite temperature, where spectral functions are broadened and bound states eventually dissolve, spectral weight will potentially be nonzero at all energies. It is then less clear which features of the spectral function are invariant (and reflect the underlying physics) and which are e.g. operator dependent.

Smearing was already discussed to some extent in Ref. Aarts:2015mma . Here we study the role of different operators. We focus on the nucleon, with the interpolator chosen to be

[TABLE]

where $n=4,5,6$ and operators $Y_{4}=\gamma_{5}$, $Y_{5}=\gamma_{4}\gamma_{5}$ and $Y_{6}=\frac{1}{2}(Y_{4}+Y_{5})$ (this nomenclature follows Chroma Edwards:2004sx ). Note that in the main part of the paper we have used operator $Y_{4}$. The operator dependence is shown in Fig. 10 for two temperatures. We observe that the correlators depend on the operator, as expected, since the overlap with ground- and excited states will differ. This manifests itself e.g. in the skewness of the correlator below $T_{c}$, while at high temperature approximate parity doubling is visible for all three operators. Below $T_{c}$, we can quantify the spectral properties more precisely by comparing the masses $m_{\pm}^{N}$ from exponential fits, see Table 3. We observe that the positive-parity mass $m_{+}^{N}$ is stable and consistent within the error. The negative-parity mass $m_{-}^{N}$ is consistent for operator 4 and 5, while for operator 6 the error is twice as large. This can be explained by noting that in Fig. 10 (top, left) the correlator is most skewed for operator 6, which leads to the smallest temporal range available on the negative-parity side, which is then reflected in the larger uncertainty.

Fig. 10 (bottom) shows the corresponding spectral functions, where at $T/T_{c}=0.76$ we observe groundstate peaks on the positive-parity side for all three operators. The position of the second peak at $\omega>0$ depends on the operator used; hence no physical relevance can be assigned to it. On the negative-parity side the overlap with the groundstate is less pronounced. In particular operator 6 seems to have especially poor overlap with low-energy features on the negative-parity side. Just as above, this finding can be understood from the asymmetric shape of the correlator: the number of data points available for MEM is very limited. At $T/T_{c}=1.52$, the approximate symmetry between the two parity sides is emerging, with the positive side still slightly enhanced, for all three operators. The fact that the overall area under the spectral curves appears different is related to the choice of normalisation. Yet the emerging symmetry, i.e. parity doubling, is present in all three cases, independent of the operator.

At the lowest temperature, we left out the points around the minimum of the correlators, both in the mass fits and the spectral function analysis, to handle a (mild) signal-to-noise problem. The effect of choosing various time ranges in the MEM analysis is shown in Fig. 11 (left). For both ranges the groundstates are clearly distinguishable and in agreement, while differences appear for the possible excited states, which is as expected.

In the results presented above, smearing was used to single out the groundstate at low temperature and suppress contributions from highly excited states at all temperatures. As a final result we show in Fig. 11 (right) the spectral function obtained at the lowest temperature in the nucleon channel, using local sources and sinks, i.e. without smearing. We observe a large contribution at higher energy, which is however not related to the low-energy states discussed above. The groundstate in the positive-parity channel is in fact still visible, as indicated in the inset, albeit much suppressed. When taken at face value, the mass is larger than found above, which is presumably due to the difficulty of extracting a signal from the local correlator. This figure therefore indicates the importance of smearing in this analysis, from a spectral function point of view.

In conclusion, we find that smearing and the choice of operator affects the correlators and hence the associated spectral functions at all temperatures. This is expected. At zero temperature, the masses of the groundstates are stable against these variations, as long as the groundstates are clearly identifiable. At nonzero temperature, the information gleaned from spectral functions is at a more qualitative level. Nevertheless, the conclusions drawn from the correlators and spectral functions are in agreement.

6 Conclusion

We studied the fate of the $N$, $\Delta$ and $\Omega$ baryons as the temperature is increased, using simulations with $N_{f}=2+1$ flavours of light quarks on anisotropic lattices. In the hadronic phase, we observed a strong temperature dependence of the groundstate masses for the negative-parity baryons, while the masses of the positive-parity baryons are stable up to the deconfinement transition. The temperature dependence is such that the positive- and negative-parity groundstates become approximately degenerate close to this transition. Degeneracy, i.e. parity doubling, is expected to coincide with chiral symmetry restoration and hence the transition from the hadronic to the quark-gluon plasma, but the precise manner in which this occurs is not known a priori. It would therefore be interesting to compare and contrast our nonperturbative predictions with model approaches, such as those discussed in Refs. Detar:1988kn ; Nemoto:1998um ; Jido:1998av ; Zschiesche:2006zj ; Steinheimer:2011ea ; Benic:2015pia ; Motohiro:2015taa ; Nishihara:2015fka ; Hohler:2015iba , to reach further insight and understanding.

In the deconfined phase, we found strong indications that the light baryons no longer exist. Here we study parity doubling directly from an analysis of the correlators, using the $R$ parameter (67), relating the positive- and negative-parity channels. We find a clear signal for the emergence of parity doubling, with the $R$ parameter acting as a quasi-order parameter. In the case of the $\Omega$ baryon, with the heavier $s$ quark, we find that parity doubling is not yet fully realised for the temperatures we considered. The effect of the quark mass is expected to vanish at higher temperatures, as $m_{q}/T\to 0$.

The conclusions from the correlator analysis are supported by the results obtained from the associated spectral functions. In the baryonic sector in vacuo, it is well understood that smearing and the use of optimised operators are essential to find clear signals for the ground- and other states. At finite temperature, with nonzero spectral weight at all energies, it is not immediately clear how to proceed with smearing and operator choice. In this paper we choose to optimise the smearing parameters and operators at zero temperature and keep them fixed as the temperature increases. With this prescription we found it is possible to obtain quantitative results from the correlator analysis and qualitative insight from the spectral functions, which are mutually consistent. It would be interesting to consider this question further and e.g. employ variational bases, widely used in vacuum, also at finite temperature, as suggested in Ref. Harris:2016usb .

As an outlook, there are various directions in which this study can be taken further, in addition to those mentioned above. From the viewpoint of lattice QCD, an important role is played by chiral symmetry. Since the Wilson-clover quarks employed here break chiral symmetry at short distances (and the two light flavours are still somewhat heavy), it would be interesting to repeat this calculation with manifestly chiral (domain wall/overlap) fermions. The signal for parity doubling should then be easily visible in the correlators, without the need to suppress short-distance contributions. A physical question is related to the role of strangeness, since a finite $s$ quark mass breaks chiral symmetry explicitly. For the $\Omega$ baryon, we indeed observed the effect of the strange quark mass in the signal for parity doubling, but a more comprehensive study of strange baryons would enlighten this further. Finally, we observed strong in-medium effects for the negative-parity baryons in the hadronic phase. It would hence be interesting to investigate whether and how this affects heavy-ion phenomenology, e.g. in the context of the hadron resonance gas or the statistical hadronisation model Stachel:2013zma .

Acknowledgments

We thank Sinéad Ryan, Tim Burns, Thomas Cohen, Claudia Ratti and Rene Bellwied for discussion. The work has been supported by STFC grant ST/L000369/1, SNF grant 200020-162515, ICHEC, the Royal Society, the Wolfson Foundation and the Leverhulme Trust, and has been performed in the framework of COST Action CA15213 THOR. We are grateful for the computing resources made available by HPC Wales. This work used the DiRAC Blue Gene Q Shared Petaflop system at the University of Edinburgh, operated by the Edinburgh Parallel Computing Centre on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This equipment was funded by BIS National E-infrastructure capital grant ST/K000411/1, STFC capital grant ST/H008845/1, and STFC DiRAC Operations grants ST/K005804/1 and ST/K005790/1. DiRAC is part of the National E-Infrastructure.