Hagedorn bag-like model with a crossover transition meets lattice QCD

Thermodynamic functions, the (higher-order) fluctuations and correlations of conserved charges at $\mu_B = 0$, and the Fourier coefficients of net-baryon density at imaginary $\mu_B$, are considered in the framework of a Hagedorn bag-like model with a crossover transition. The qualitative behavior of these observables is found to be compatible with lattice QCD results. Fair quantitative description of the lattice data is obtained when quasiparticle-type quarks and gluons with non-zero masses are introduced into the bag spectrum. The equation of state of the model exhibits a smooth and wide crossover transition.


I. INTRODUCTION
The empirically known spectrum of hadrons suggests a rapid, possibly exponential, increase of the density of states at large masses [1]. An exponentially rising hadron mass spectrum was first proposed by Hagedorn in the 1960s [2] in the framework of the statistical bootstrap model, long before the advent of QCD as the fundamental theory of strong interactions. Evaluations within the MIT bag model [3] similarly suggest an exponentially increasing mass spectrum [4]. The Hagedorn mass spectrum is characterized by the Hagedorn temperature T H ∼ 170 MeV, above which a transition to a new state of matter occurs according to the early ideas [5]. The asymptotic freedom property of QCD suggests a transition to the quark-gluon plasma phase at high temperatures and densities. First-principle lattice QCD simulations at zero baryochemical potential are consistent with a smooth crossover transition between hadronic and partonic matter [6], characterized by the pseudocritical temperature T pc 155 MeV [7,8] obtained from the analysis of chiral observables.
Hagedorn states are possibly created in multi-particle reactions, e.g. during heavy-ion collisions [9], most abundantly close to the Hagedorn temperature, as was discussed in [10,11]. Their appearance can explain the fast chemical equilibration of the hadronic gas, this especially concerns the abundances of (multi-)strange baryons and their anti-particles.
The presence of Hagedorn states in a hadron resonance gas model provides also a lowering of the speed of sound and of the shear viscosity over entropy density ratio [12,13]. Recently also a microcanonical bootstrap description of Hagedorn states with explicit conserved quantum numbers has been developed: The energy spectra of the resulting hadrons from the decays of such exotic states follow an exponential law akin with the Hagedorn temperature and thus look thermal by itself [14]. Incorporating those states into a microscopic hadronic transport model, again fast equilibration of strange and multistrange baryons and mesons has been shown [15], and the full dynamics of heavy-ion collisions within such an unorthodox picture has also been developed [16].
In the simplest Hagedorn model all hadrons are treated as point particles. Due to the exponentially increasing hadron mass spectrum, the Hagedorn temperature T H becomes the limiting temperature above which the partition function diverges -a behavior which cannot be reconciled with lattice QCD results. Early on, it has been suggested that hadrons in a statistical system should be treated as spatially extended objects [17][18][19][20][21], which essentially corresponds to a van der Waals type excluded volume correction in the partition function of a hadron gas. The inclusion of the spatial size of hadrons leads to a disappearance of the "limiting" temperature under certain conditions [21]. Moreover, as shown in Refs. [20,22], there is a possibility of a first-order, second-order, or a crossover transition in the gas of spatially extended quark-gluon bags, with thermodynamic properties at high temperatures being similar to the MIT bag model equation of state [17]. Both phases are described within a single partition function. Different possibilities, such as a crossover transition at zero baryon density and a first-order phase transition at finite baryon density in the gas of quark-gluon bags were explored in various works [23][24][25][26][27][28].
The temperature dependence of thermodynamic functions at zero chemical potential within the gas of extended quark-gluon bags with a crossover transition was considered in Ref. [27] in the context of lattice QCD equation of state. General qualitative features were found to be compatible with lattice QCD, although a quantitative description is lacking.
Recently, a lot of lattice data has appeared on correlations and fluctuations of conserved charges. These observables correspond to the derivatives of the partition function with respect to chemical potentials, and they have long been considered sensitive to phase transitions [29,30]. Nowadays, the susceptibilities are actively being used to formulate, test, and constrain various effective QCD models for equation of state at non-zero baryon density [31][32][33][34][35][36][37]. In the present work we explore to what extent the behavior of the conserved charges susceptibilities in the gas of extended hadrons and quark-gluon bags is compatible with lattice QCD. We also consider an extension of the conventional quark-gluon bag model by introducing the effects of a quasiparticle-like finite, constituent quark and gluon masses, which lead to a substantially improved agreement with the lattice data.

A. The partition function
The model assumes a multi-component system of color neutral objects. These objects, henceforth referred to as particles, have finite sizes -the eigenvolumes. The particles can carry three abelian charges -baryon number, electric charge, and strangeness. These three charges are characterized by the corresponding chemical potentials µ B , µ Q , µ S . For convenience, we will employ the fugacities, λ i ≡ e µ i /T . First, we consider a system with a finite number of different components f , a generalization to an infinite number of components will follow later. The particles under consideration can have arbitrary integer values of baryon charge, electric charge, and strangeness.
It is assumed that particles are non-overlapping -this constraint is modeled through the excluded-volume correction [20,38]. The grand canonical partition functions reads Here and s i , are the baryon charge, electric charge, and strangeness of particle species i, d i is its degeneracy, and The presence of the theta function in (1) causes certain technical difficulties. These technical difficulties can be overcome by considering the isobaric (pressure) ensemble [20].
The isobaric partition function,Ẑ(T, s, λ B , λ Q , λ S ), is given as the Laplace transform of In the thermodynamic limit, V → ∞, the grand canonical partition function behaves The isobaric partition function in the thermodynamic limit has the form It follows from Eq. (5) thatẐ(T, s, λ B , λ Q , λ S ) has a singularity at s = s * = p/T , and no singularities at s > s * . The integral representation (5) is unconditionally divergent at s < s * , and, therefore, does not directly provide information about the behavior of the isobaric partition function in that region. Thus, there is a possibility thatẐ has singularities at s < s * . These considerations lead to the conclusion that the system pressure p(T, λ B , λ Q , λ S ) in the thermodynamic limit is defined by the farthest-right singularity s * of the isobaric partition functionẐ, i.e.
The exact nature of the farthest-right singularity s * depends on the input particle spectrum.

B. The mass-volume density of states
In the isobaric ensemble, the input particle spectrum enters through Eq. (4) only. This permits a generalization to a system with an infinite number of different components, characterized by some density function. First, let us rewrite Eq. (4) in the following form Here the sum i goes through all particles which carry the specific baryon number b, the electric charge q, and strangeness s. One can now introduce a mass-volume density of states ρ(m, v; b, q, s), which determines the number of particle states carrying fixed quantum numbers b, q, and s, in the mass-volume interval [m, v; m + dm, v + dv]. Equation (7) is Equation (8) is reduced to (7) for ρ(m, v; b, q, s) = i∈{b,q, Finally, let us rewrite Eq. (8) as follows where we have introduced the generalized fugacity dependent mass-volume density of states: Note that ρ(m, v) ≡ ρ(m, v; λ B = 1, λ Q = 1, λ S = 1) is the mass-volume density of all states irrespective of their quantum numbers.
We follow the picture presented in Refs. [22,23]. The particle spectrum is assumed to consist of two contributions: Tables [1]; 2. an exponential Hagedorn spectrum of the heavy quark-gluon bags.

the established, low mass hadrons and resonances listed in Particle Data
Therefore, where ρ H corresponds to the established hadrons listed in PDG, and ρ Q corresponds to the quark-gluon bags. The PDG hadrons form a discrete spectrum, therefore ρ H is given as a finite sum of δ functions: 1 Each of the PDG hadrons is assumed to have a finite eigenvolume parameter v i . In the spirit of the bag model, we assume here that the eigenvolumes of the PDG hadrons are proportional to their mass: v i = m i /ε 0 , where ε 0 = 4B unless stated otherwise. Here B is the bag constant. In principle, one can consider a different parametrization of v i , e.g. based on phenomenological knowledge of some hadron-hadron interactions, the only requirement here is that all eigenvolumes are non-vanishing, i.e. v i > 0. 1 We neglect here the effects of finite resonance widths. These can have an important effect in precision thermal model applications, such as thermal fits [39], but are not very relevant for the mostly qualitative aspects of the equation of state studied here.

C. Quark-gluon bags
The mass-volume spectrum of the quark-gluon bags, ρ Q , is the crucial ingredient of the model, which determines some of its most qualitative features. The form of ρ Q depends strongly on the assumptions regarding the internal color-flavor structure of the bags (see, e.g., Refs. [24,25]). In the region where both m and v are large, the spectrum can be described in the framework of the bag model [3]. The mass-volume density of states was computed assuming bags filled with non-interacting massless quarks and gluons, at zero chemical potentials [4,40,41], and also for finite baryon chemical potential [42]. One Here V 0 is a model parameter which is given a sufficiently large value, B is the bag constant, and M 0 > 0 is a parameter which regularizes the mass-volume density close to the lower threshold [23]. As will be shown, the exact value of M 0 has no significance for applications considered in this paper. σ Q corresponds to the energy density (or three times the pressure) of the non-interacting gas of massless quarks and gluons. Here this quantity is taken as a function of all three chemical potentials in (2+1)-flavor QCD: with (13) implies that the eigenvolume of a QGP bag with a given mass is fluctuating. These fluctuations are given by the distribution ρ Q (m, v; λ B ). The presence of volume fluctuations is crucial for obtaining a transition to quark-gluon plasma: assuming the fixed mass-volume relation, e.g. m = 4Bv from the MIT bag model [3], leads to a constant energy density at high temperature, which is incompatible with lattice QCD (see Refs. [24,25,27] for more details).
The values of parameters C, γ, and δ in the pre-exponential factor in Eq. (13) depend strongly on the details of the bag model calculation. For example, they depend on whether the colourlessness constraint for the bags is implemented [41][42][43], as well as on other internal symmetry constraints considered [24,25]. Therefore, these parameters are usually treated as free model parameters. Such a philosophy is considered in the present work as well. Similar arguments apply to the possible dependence of C, γ, and δ on fugacities λ B , λ Q , λ S . In the absence of a detailed knowledge, we omit the possible dependence of these parameters on fugacities in the present study.
The Hagedorn temperature. The exponential hadron mass spectrum was first introduced by Hagedorn [2]. It has the general form with the parameter T H called the Hagedorn temperature. The mass-volume relation (13) for quark-gluon bags employed in the present work also implies this same exponential asymptotic mass spectrum. The mass spectrum for the quark-gluon bags reads The integral converges as long as M 0 > 0. The function w(v; m) has a peak for large values of m. Therefore, the integral in (13) can be approximated for large m with the Laplace's method. One has w(v; m) ≈ w(v 0 ; m) + Equation (17) is the familiar relation between the mass of a quark-gluon bag and its average volume from the MIT bag model with massless quarks [3]. Applying Laplace's method to Eq. (16) one obtains the asymptotic mass spectrum This form coincides with the Hagedorn mass spectrum (15) with .
where f is given by Eq. (4). FunctionẐ has a pole singularity at s H = f (T, s H , λ B , λ Q , λ S ).
Another possibility is a singularity s Q in the function f (T, s, λ B , λ Q , λ S ) itself.
Let us compute the function f for the particular mass-volume density given by (11)- (13).
First, we split it into two parts The quark-gluon bags in Eq. (24) are heavy, m 2 GeV, therefore one can use the non-relativistic approximation: This approximation has a relative accuracy of 10% or better for m/T > 20, a condition which is realized in the applications considered in this work. The error due to the nonrelativistic approximation for quark-gluon bags in the resulting pressure is even smaller, see The expression for f Q then simplifies to where The integration over m in Eq. (26) can be carried out approximately using Laplace's method. One has One can invert the above relation to obtain One can see that v(m 0 ) is a decreasing function of temperature. This elucidates the so-called effect of thermal compression of bags.
Let us note that Applying Laplace's method to Eq. (26) one obtains Here s B corresponds to s B = p B /T , where p B (T, λ B , λ Q , λ S ) coincides with the pressure in the MIT bag model equation of state [17]: Recalling the definition of the "upper" incomplete gamma function one can perform the integration over v in Eq. (33) explicitly Note that dependence on λ B , λ Q , λ S in the above relation enters through . Also note that the application of Laplace's method eliminates the dependence of the final result on the parameter M 0 from Eq. (13).
The function f has a singularity at s Q = s B , as follows from (36). Thus, the system pressure is defined at given temperature and chemical potentials as The model may contain a phase transition, defined as a "collision" of singularities s H and s Q at particular values of the thermodynamic parameters, i.e. s H (T c ) = s Q (T c ) at the critical temperature T c of the phase transition. This mechanism was first described in Ref. [20]. In this case s H (T ) > s Q (T ) for T < T c and s H (T ) < s Q (T ) for T > T c . A detailed analysis performed in Ref. [28] reveals that a phase transition as described above is only realized when γ + δ < −3 and δ < −7/4. For other values of these parameters a crossover-type transition is realized, i.e. s H (T ) > s Q (T ) for all T . If the crossover transition i.e. the system proceeds to the phase which has thermodynamic properties similar to that of the quark-gluon plasma described by the MIT bag model equation of state.
Whether this is the case depends on the values of γ and δ [28].
Lattice QCD simulations at physical quark masses reveal that the transition from hadronic to partonic degrees of freedom at µ B = 0 is of a crossover type [6][7][8]. Therefore, in this work we focus only on the case where the crossover transition is realized. The possibility of a real phase transition at finite µ B will be considered in a separate publication.
The farthest-right singularity of the isobaric partition functionẐ is equal to s H for all possible values of thermodynamic parameters when the crossover scenario is realized. The pressure, p = T s H , satisfies the following transcendental equation: The first term in Eq. (38) corresponds to the discrete part of the particle spectrum. It equals the sum of the partial pressures evaluated self-consistently for an ideal Boltzmann The second term corresponds to the contribution of the quark-gluon bags. The two terms are not independent -they both depend self-consistently on the total system pressure to which they both contribute. Equation (38) is solved numerically in the present work. The energy density, the entropy density, the speed of sound and the various susceptibilities are obtained from the pressure function as derivatives with respect to T or λ B,Q,S , through the standard thermodynamic relations.
One should note that the application of the Laplace's method used to perform the mass integration in the derivation of Eq. (38) is, strictly speaking, most accurate in the limit of large masses, i.e. when m 0 [Eq. (30)] is large. The method is expected to be less accurate if contributions of "small" quark-gluon bags close to the mass-volume density threshold V 0 , are significant. This situation can take place in the vicinity of the crossover transition, where the bags start to appear in addition to the PDG hadrons. In Appendix V B we show numerically that, for the parameter sets used in the present study, Laplace's method allows to evaluate the pressure with a precision of better than 2%, for all temperature values considered. We therefore adopt Laplace's method in all our subsequent calculations, as that method remedies some technical difficulties and presents a clear physical picture. In a more elaborate study one may omit the Laplace's method approximation altogether.
It should be pointed out that the mass-volume density (13) of quark-gluon bags is obtained for asymptotically large masses and volumes, whereas the lower end of the massvolume spectrum is regulated by the cut-off parameters V 0 and M 0 only. One may ask how this lower end of the mass-volume spectrum matches with the spectrum of the PDG hadrons at even lower masses. In Appendix V A we show that the spectra of PDG hadrons and quark-gluon bags can indeed be merged rather smoothly for the model parameter sets under consideration here.

Particle number density
The particle number density is the average number of particles per unit volume. It is the sum (integral) of individual densities corresponding to different particle species: where the sum goes over all species. As the considered particle spectrum includes the continuous spectrum of the quark-gluon bags, the sum in Eq. (39) in general corresponds to the integral over the mass-volume density of states, i.e.
The individual densities can be computed by introducing the fictitious fugacities λ i for all species into the partition function. n i are given by the standard expression for the multi-component excluded-volume model [44]: Here and p is the system pressure.
The total hadron density reads with For the particle spectrum (11) consisting of the PDG hadrons and quark-gluon bags the above quantities can be computed explicitly. The calculation proceeds in the same fashion as done for the pressure function in Sec. II D: the PDG part of the mass-volume density is computed explicitly, whereas for the quark-gluon part one first applies Laplace's method to perform the integration over the mass in Eqs. (43) and (44), the remaining integrals over the volume can be expressed in terms of the incomplete Gamma function. The result is:

Filling fraction
Another interesting quantity is the filling fraction (f.f.) -the ratio between the average total volume occupied by hadrons over the system volume. The definition of this quantity The average eigenvolume of a particle in the thermal system can be computed as follows: (48)

Average particle mass
The average mass of a particle in the thermal system is given by The explicit calculation, employing the Laplace integration over m and the incomplete Gamma function to express the integral over v, yields

III. CALCULATION RESULTS FOR BAGS WITH MASSLESS QUARKS AND GLUONS
Calculations here are performed for the following set of parameters: All model parameters are fixed, the only exception is the δ exponent. In the present study we fix γ = 0, for simplicity. Non-zero γ values can be considered equally well. For γ = 0, the crossover-type transition is realized if δ ≥ −3 [28]. Therefore, the δ exponent is varied here in the range −3 ≤ δ ≤ − 1 2 , in the steps of 1 2 . This variation of δ corresponds to the range 0 ≤ α ≤ 5 2 for the exponent α in the Hagedorn mass spectrum [see Eq. (21)]. The scan in δ performed here is, therefore, similar to the study presented in Ref. [27]. The temperature dependence of the scaled pressure p/T 4 and the scaled energy density ε/T 4 is depicted in Fig. 1. The model shows a crossover transition, the functions plotted approach the Stefan-Boltzmann limit of massless quarks at high temperatures. The results are compared with the lattice QCD data of the Wuppertal-Budapest collaboration 3 (blue bands) [45]. On a quantitative level, the agreement of the model with the lattice data is not very good. This especially true for the energy density: the model predicts a peak in the temperature dependence of ε/T 4 -a qualitative feature not seen in lattice simulations.
The main reason for this disagreement is that the QGP phase is described by the MIT bag model with massless quarks, which is known to provide only a rough description of QCD thermodynamics at large temperatures.
The auxiliary quantities introduced in Sec. II E are depicted in Fig. 2. The filling fraction (Fig. 2a) shows a monotonic increase with temperature, from small values (f.f. 0) at small temperatures towards f.f. 1 at high temperatures. This implies that almost the whole volume is occupied by the finite-sized particles at high temperatures.
The particle chemistry at different temperatures can be clarified by studying the temperature dependence of the mean hadron mass m (Fig. 2b). The non-trivial behavior of v with respect to δ similarly implies a non-trivial behavior of the particle number density n (Fig. 2d). Indeed, as follows from Eqs. On the other hand, at δ < −7/4 one has v → ∞ which implies n → 0. The numerical results shown in Fig. 2d are consistent with these considerations. In fact, the vanishing hadron number density for δ < −7/4 implies that an arbitrary large but finite subvolume of the system is occupied at high temperatures by a single bag filled with quark-gluon plasma.
Fluctuations and correlations of conserved charges are another observables, accessible with lattice QCD, suggested long ago to be sensitive to the parton-hadron transition [29,30]. These observables, henceforth referred to as susceptibilities of conserved charges, are defined by the derivatives of the pressure function with respect to the corresponding chemical potentials: The matrix of the second order conserved charges susceptibilities has been studied in lattice QCD simulations at the physical point in Refs. [47,48]. Lattice simulations agree well with the predictions of the ideal HRG model at temperatures below the pseudocritical one, and show a behavior consistent with an approach towards the Stefan-Boltzmann limit at high temperatures. The temperature dependencies of the second order diagonal susceptibilities, χ B 2 , χ Q 2 , and χ S 2 , are shown in Fig. 3 in comparison with the lattice QCD data. The behavior of the susceptibilities in the model is qualitatively compatible with lattice QCD results. From a quantitative point of view, one sees that the approach to the Stefan-Boltzmann limiting values is too fast in the model compared to the lattice data.

A. Modification of the model
While the simple bag model picture above appears to describe many qualitative features seen in lattice data, the quantitative description of the main thermodynamical functions, such as pressure, energy density, interaction measure, and the speed of sound, is obviously not very good. This description cannot be notably improved solely by a variation of the These values are taken here as a representative case, other combinations of the thermal quark and gluon masses are certainly possible.
Equation (14) for σ Q should be modified to reflect massive quarks and gluons in the bag model equation of state. The modified σ Q is temperature-dependent and reads With this modification of σ Q , Eq. (34) reproduces the heavy-bag model studied in Ref. [56].
Here, for simplicity, we only consider temperature-independent quark and gluon masses.
A more involved model may take into account their temperature dependence, as typically done in quasiparticle models. This can be achieved by employing, e.g., a hard-thermal-loop description [57,58] for the intrinsic thermal pressure of the bags. An explicit temperature dependence of the quark and gluon masses, however, will require careful considerations regarding the thermodynamic consistency of the model [50], and will be considered elsewhere.

B. Modification of the parameters
The Equation (21) which relates the Hagedorn temperature T H to the bag constant B should be modified for quarks and gluons with finite thermal masses. This is because the quantity σ Q is modified and now depends on the temperature, i.e. σ Q = σ Q (T ). The Hagedorn temperature can therefore be obtained as the solution of the following transcendental equation: The decreased value of the bag constant necessitates an increase of the value of the parameter V 0 , to avoid a large overlap of the spectrum of stable hadrons and quark-gluon bags. We, therefore, adopt the value V 0 = 8 fm 3 in the following.
As before, we set γ = 0 and C = 0.03 GeV −δ+2 . We do not vary the value of δ but settle here for the value δ = −2. In this case one expects a crossover transition to a gas of infinitely large quark-gluon bags in the limit of high temperatures, as elaborated in the previous section. were taken to be four times smaller, i.e. ε 0 = 16B. The lattice QCD data of the Wuppertal-Budapest collaboration [45] are shown by the blue bands.

C. Thermodynamic functions
The temperature dependence of the scaled pressure p/T 4 , the scaled energy density ε/T 4 , the interaction measure (ε − 3p)/T 4 , and the speed of sound squared c 2 s = dp/dT at µ B = 0 is depicted in Fig. 5. All quantities are in a rather good agreement with the lattice QCD data of the Wuppertal-Budapest collaboration [45]. It is particularly notable that the scaled energy density shows a monotonic behavior, consistent with lattice QCD, in contrast to the previous simple model where the bags are filled with massless quarks and gluons. The behavior of the auxiliary quantities from Sec. II E is found here to be quite similar to the one shown in Fig. 2 for the massless quarks and gluons case.
The model describes the lattice data on a semi-quantitative level. Sizeable deviations are seen close to T H only for the trace anomaly and the speed of sound squared, which is sensitive to the second temperature derivative of the pressure function. An improved description of the data can be achieved by variations of the free parameters of the model. However, in light of the general limitations of the present model this appears to be rather unnecessary. We proceed, instead, by studying, in this model, the behavior of those lattice observables, which are considered to be sensitive probes of the nature of the quark-hadron transition.

D. Second order correlations and fluctuations of conserved charges
We turn now to the behavior of the susceptibilities of conserved charges. The temperature dependence of the matrix of the second order conserved charges susceptibilities in the present model is depicted in Fig. 6. These are compared to the lattice QCD data of the Wuppertal-Budapest [47] and HotQCD collaborations [48]. The model predictions agree qualitatively with the lattice data for all observables considered. A notable underestimation of the lattice data is observed for χ Q 2 and χ BQ 11 in the vicinity and also above the crossover temperature region. We argue that these observables are sensitive to the eigenvolume values assumed for the discrete, PDG part of the hadronic spectrum. In order to illustrate this sensitivity, we additionally depict in Fig. 6 , computed in the Hagedorn model with quark-gluon bags filled with massive quarks and gluons. Lattice QCD data of the Wuppertal-Budapest [47] and HotQCD collaborations [48] are shown by the blue and green bands, respectively.
for charged pions, owing to their small masses and to the fact they carry electric charge [59]. The resulting effect on χ Q 2 is depicted in Fig. 6(b) by the dash-dotted line. The inclusion of the Bose statistics for pions improves the agreement with the lattice data for χ Q 2 in the crossover temperature region, other observables are almost unaffected.
We also consider the baryon-strangeness correlator ratio, followed by a quick saturation just above the Hagedorn temperature T H in the model. This behavior is consistent with the lattice QCD data.

E. Higher-order susceptibilities
Higher-order susceptibilities are expected to be particularly sensitive to crossing the crossover transition [61]. The higher-order susceptibilities, such as χ B 4 /χ B 2 , χ S 4 /χ S 2 , χ B 6 /χ B 2 , and χ B 8 , were recently computed in lattice QCD [62][63][64][65]. Here we study the above observables in the Hagedorn bag-like model with massive quarks and gluons. The results are depicted in Fig. 8, together with the lattice data.
The so-called kurtosis of the net baryon number fluctuations -the χ B 4 /χ B 2 ratio -shows a rapid decrease from unity towards the Stefan-Boltzmann limiting value of 2/(3π 2 ) in the temperature range T = 150 − 200 MeV, see Fig. 8(a). In the conventional ideal HRG model this quantity is equal to unity, owing to the fact that no multi-baryon hadrons are known to exist. This scenario is reasonable for low temperatures, where a dilute hadron gas is expected, and it is realized in the present model. The deviation of this observable from unity in the vicinity of the pseudocritical temperature, seen in lattice QCD, is often interpreted as a signal for a rapid hadron melting and transition to a deconfined phase [66]. However, this onset of deviations from unity is well captured also in a HRG model with repulsive excluded volume interactions [32,33,35,67]. The present Hagedorn model extends the HRG model to include the exponentially increasing Hagedorn bag spectrum with massive quarks and gluons, as well as the excluded volume corrections. The combined effect of these two extensions leads to the behavior shown in Fig. 8(a) which is consistent with the lattice data on a quantitative level. At this point it is necessary to emphasize the importance of including both the exponential Hagedorn spectrum and the excluded-volume interactions.
Only when both of effects are included simultaneously is it possible to obtain a smooth transition to the quark-gluon plasma type behavior of χ B 4 /χ B 2 at high temperatures. Taking into account the Hagedorn spectrum, but not the excluded volume corrections, would lead to a monotonic increase of χ B 4 /χ B 2 with temperature, with a subsequent divergence at the Hagedorn temperature [68] -a behavior inconsistent with lattice QCD.
The temperature dependence of the net strangeness kurtosis -the χ S 4 /χ S 2 ratio -shows a peak at T 160 MeV, both in the model and in the lattice data. It is also consistent with an approach to unity at low temperatures, and an approach to the Stefan-Boltzmann limit at high temperatures. The presence of the peak in the temperature dependence of χ S 4 /χ S 2 is an interplay of two effects: (i) the presence of the multi-strange hyperons in the list of known hadrons causes the initial increase of χ S 4 /χ S 2 to above unity [62], whereas (ii) the presence of the excluded-volume corrections suppresses this ratio at higher temperatures [33]. The inclusion of the Hagedorn quark-gluon bag states enforces the quark-gluon plasma type behavior of this observable at high temperatures. As shown by the dashed line in Fig. 8(b), this observable is rather sensitive to the magnitude of the eigenvolumes taken for the PDG hadrons. This sensitivity is less pronounced for χ B 4 /χ B 2 . The behavior of the sixth-and eight-order net baryon susceptibilities χ B 6 /χ B 2 and χ B 8 , shown in Fig. 8 (c) and (d), shows a strong non-monotonic temperature dependence. As far as the present level of accuracy in the lattice data is concerned, the Hagedorn model provides a reasonable quantitative description of these data. As within its present formulation the model exhibits a crossover transition at both zero and non-zero baryon density, this agreement suggests that strong non-monotonic behavior seen in lattice data is not unambiguously related to possible critical phenomena, at least not directly.

F. Fourier coefficients at imaginary µ B
The model can also be applied to study observables at imaginary chemical potentials. This is achieved through the analytic continuation. These observables can then be compared with lattice QCD data at imaginary chemical potentials. Such a comparison with an independent set of lattice observables provides an important cross-check of the model validity.
Here we consider the behavior of the model at the imaginary baryochemical potential µ B = iθ B T , the electric and strangeness chemical potentials are set to zero. The QCD partition function is an even function of µ B because of the CP-symmetry, and it is periodic in the imaginary µ B /T direction with the period of 2π -the Roberge-Weiss symmetry [69].
Therefore, it is sufficient to apply the model in the interval 0 < θ B < π, where it exhibits analytic behavior 4 .
The QCD pressure at imaginary baryochemical potential can be written in terms of the Fourier series with the Fourier coefficients The net baryon density at imaginary µ B reads with b k ≡ k p k .
The leading four Fourier coefficients b k of the net baryon density at imaginary µ B were recently calculated in lattice QCD and presented in Ref. [70]. The coefficients were used to constrain the parameters of various phenomenological models, such as the excluded volume HRG model [70] or the cluster expansion model [36]. an uncorrelated gas of hadrons at low temperatures. The higher-order coefficients start to visibly depart from zero as the temperature is increased to above the Hagedorn temperature T H . The coefficients show an alternating sign structure: the odd-order coefficients, b 1 and b 3 , are positive, whereas the even-order ones, b 2 and b 4 , are negative. The emergence of this structure was explained in Ref. [70] in terms of a baryonic excluded-volume, an alternating sign structure is also expected for an uncorrelated massless gas of quarks at high temperatures.
The present model agrees qualitatively with the available lattice data. It does appear to underestimate b 1 and overestimate the higher-order coefficients at certain temperatures.
The quantitative description can be improved by a variation of the model parameters. On the other hand, the chiral transition, taking place with increasing temperature in the crossover region, has not been discussed, as the present model is not suited for a straightforward calculation of the chiral order parameter -the chiral condensate qq . The phenomenological picture of the chiral transition, however, can be given. The bag-like states start to rapidly occupy nearly the whole system volume during the (crossover) transition within a small temperature interval [see the behavior of the filling fraction in Fig. 2 (a)]. As the (MIT-)bags inside are chirally restored with qq = 0, the total overall order parameter should rapidly decrease towards 0, mimicking the chiral transition.
The temperature dependence of the chiral transition in the presented picture can be characterized by the available volume fraction, 1 − f.f., which is depicted in Fig. 10 by the solid line. In a chirally "broken" phase at very low temperatures the particle densities are small and this quantity is close to unity. In a chirally "restored" phase at high temperatures, where the bags occupy almost the whole volume, the available volume fraction is close to zero.
The chiral transition takes place in a relatively narrow T ∼ 150 − 180 MeV temperature range. A related (but not identical) quantity computed on the lattice is the subtracted chiral condensate ∆ l,s [71], which is expected to show a similar qualitative behavior as one traverses the chiral transition. The corresponding lattice QCD data of the Wuppertal-Budapest collaboration [7], computed for physical quark masses, is depicted in Fig. 10  If a true second order chiral phase transition occurs at vanishing baryon chemical potential at the critical temperature, the parameters γ and δ appearing in Eq. (13) would have to be fine-tuned. A theoretical explanation for the behaviour of those parameters in the chiral limit, m q → 0, is not obvious. Also then, as m π → 0, the hadronic masses in the present thermal and extended hadronic model have to be shifted too. Mean field type studies, as e.g. in elaborated chiral models [72], cannot be straightforwardly overtaken. We leave these ideas for future studies. The Hagedorn quark-gluon bag-like model, introduced in Refs. [20,50], is historically one of the first models for a hadron-parton transition in QCD. The quantitative aspects of such a model have not been studied extensively before, the present model results and their comparison to the lattice data suggest that this approach is compatible with first-principle lattice QCD results. One remarkable feature of such a model is that the whole transition, be it a real phase transition or a crossover, is described within a single partition function.
This is quite different from many conventional phenomenological models for the equation of state, where the hadronic and partonic phases are usually described by different partition functions that are then being matched, either via the Maxwell construction [73,74] or a smooth switching function [31,32].
In the present work we have considered only the case where the crossover transition is  The picture for the "massless quarks" parameter sets [Eq. (51)] turns out to be quite similar and not shown here.

B. Accuracy of the approximations
In the present work the pressure has been determined as the solution of the transcendental equation (38). Two approximations were used in order to obtain the quark-gluon bag part Both approximations are expected to be rather accurate for the heavy quark-gluon bags, as discussed in the main text. Nevertheless, it can be useful to quantify the error introduced by these approximations. In order to do that, we consider the exact transcendental equation for the model pressure without approximations: Here is φ(T, m) is taken in the exact, relativistic form (2). At each temperature value, we first obtain the approximate solution p(Approx) by solving Eq. (38). This corresponds to  the non-relativistic approximation and Laplace's method, for calculating the pressure is well justified, at least for the parameter set considered. This can be attributed to two reasons: 1. At large temperatures the system is dominated by heavy bags. As already elaborated, the heavier the bags are, the more accurate are the approximations.
2. At low temperatures the system is dominated by the PDG hadrons, which are evaluated without approximations. Therefore, possible inaccuracies in evaluating the con-