Net-baryon number fluctuations in the Hybrid Quark-Meson-Nucleon model at finite density

We study the mean-field thermodynamics and the characteristics of the net-baryon number fluctuations at the phase boundaries for the chiral and deconfinement transitions in the Hybrid Quark-Meson-Nucleon model. The chiral dynamics is described in the linear sigma model, whereas the quark confinement is manipulated by a medium-dependent modification of the particle distribution functions, where an additional scalar field is introduced. At low temperature and finite baryon density, the model predicts a first-, second-order chiral phase transition, or a crossover, depending on the expectation value of the scalar field, and a first-order deconfinement phase transition. We focus on the influence of the confinement over higher-order cumulants of the net-baryon number density. We find that whereas the cumulants show a substantial enhancement around the chiral phase transition, they are not as sensitive to the deconfinement transition.


I. INTRODUCTION
Exploration of the phase diagram of quantum chromodynamics (QCD), the theory of the strong interaction, at finite temperature and density is one of the most challenging topics in high-energy particle and nuclear physics. At vanishing density, the first-principles calculations, i.e., lattice QCD (LQCD), provide a reliable description of equation of state and a profound insight into the phase structure of QCD [1][2][3][4]. The recent LQCD findings [5,6] exhibit a clear manifestation of the parity doubling structure for the low-lying baryons around the chiral crossover. The masses of the positive-parity groundstates are found to be rather insensitive to temperature, while the masses of negative-parity states drop substantially with increasing temperature, and the parity doublers become almost degenerate with a finite mass in the vicinity of the chiral crossover T c . Despite unphysically heavy pion mass used in the study, this is likely an imprint of the chiral symmetry restoration in the baryonic sector of QCD, and is expected to occur also in cold dense matter. This can be described in a schematic framework with chiral symmetry, the parity doublet model [7][8][9]. The model has been applied to hot and dense hadronic matter, neutron stars, as well as the vacuum phenomenology of QCD [10][11][12][13][14][15][16][17][18][19][20][21].
The mechanism of quark confinement and its relation to the chiral symmetry breaking are of major importance in probing the QCD phase transition, although it is nontrivial to embed their interplay into a single effective theory. A major approach is to introduce the temporal gauge field as a static external field to a chiral Lagrangian, so that the Polyakov loop naturally appears in the thermodynamics, e.g., the Polyakov loop-extended Nambu-Jona-Lasinio (PNJL) [22][23][24][25][26] or Polyakov loop-extended * michal.marczenko@ift.uni.wroc.pl quark-meson (PQM) [27][28][29][30][31] models. While at finite temperature and low density the Polyakov loop expectation value serves as an approximated order parameter for the deconfinement, it is highly questionable that it remains so at high density.
In this work we put special emphasis on lowtemperature and finite-density region of QCD phase diagram. Hence, instead of using the Polyakov loop, we shall employ the hybrid quark-meson-nucleon (Hybrid QMN) model [21]. The model includes quark degrees of freedom on top of hadrons, but prevents the quarks from their unphysical onset at low density. This can be achieved by a new auxiliary scalar field to which the fermions are coupled non-trivially. This field serves as a momentum cutoff in the Fermi-Dirac distribution functions, and plays a role of suppressing the unphysical thermal fluctuations of fermions, depending on density.
We will study the behavior of the second-and higherorder cumulants of the net-baryon number density up to the fourth order, as well as the bulk equation of state, in the Hybrid QMN model. It is systematically examined that to what extent the thermal behaviors are dominated by the chiral criticality and the onset of deconfinement. This paper is organized as follows. In Sec. II, we introduce the parity doublet model at finite density and finite temperature, as well as its extension -the Hybrid Quark-Meson-Nucleon model -that includes the mechanism to mimic quark confinement. We discuss the obtained numerical results on the equation of state in Sec. III, and the results for the second and higher-order cumulants in Sec. IV. Finally, Sec. V is devoted to the summary and conclusions.

II. FORMALISM
In this section, we give a brief description of the Hybrid Quark-Meson-Nucleon (Hybrid QMN) model for the QCD transitions at finite temperature and density. Following Ref. [21], we first introduce the pure hadronic model of parity doublers, and then extend it incorporating the quark degrees of freedom. Throughout this paper we consider a system with N f = 2.
A. Parity Doublet model with dilaton In the conventional Gell-Mann-Levy model of mesons and nucleons [32], the nucleon mass is entirely generated by the non-vanishing expectation value of the sigma field. Thus, the nucleon inevitably becomes massless when the chiral symmetry gets restored. This is led by the particular chirality-assignment to the nucleon parity doublers, where the nucleons are assumed to be transformed in the same way as the quarks are under chiral rotations.
More general allocation of the left-and right-handed chiralities to the nucleons, the mirror assignment, was proposed in [7]. This allows an explicit mass-term for the nucleons, and consequently the nucleons stay massive at the chiral restoration point. For more details, see Refs. [7][8][9].
In the mirror assignment, under SU (2) L × SU (2) R rotation, two chiral fields ψ 1 and ψ 2 are transformed as follows: where ψ i = ψ iL + ψ iR , L ∈ SU (2) L and R ∈ SU (2) R . The nucleon part of the Lagrangian in the mirror model reads where g 1 , g 2 , and g ω are the baryon-to-meson coupling constants and m 0 is a mass parameter. Since the presence of the mass spoils scale invariance, we further impose that the theory at classical level is invariant under scale transformation. To this end, we introduce a dilaton χ [33] and replace the above mass with with a new coupling g χ . The nucleon mass emerges when the dilaton field gets condensed. The mesonic part of the Lagrangian reads where ω µν = ∂ µ ω ν − ∂ ν ω µ is the field-strength tensor of the vector field, and the potentials read The full Lagrangian of the parity doublet model is then L = L N + L M . Note that the chiral symmetry and the scale invariance are explicitly broken by the linear term in σ (Eq. (5a)). The parameters λ, λ χ and can be connected to the vacuum meson masses and the pion decay constant as where the pion mass m π = 138 MeV, the pion decay constant f π = 93 MeV and χ 0 is the vacuum expectation value of the dilaton field. We shall treat the mass of the sigma meson, m σ , as a parameter. The value of the constant B in Eq. (5c) is fixed by identifying the χ 4 with the gluon condensate G µν G µν . Since the gluon condensate is directly proportional to the non-vanishing trace anomaly, it can be transmuted into B, assuming that the QCD vacuum energy is dominated by the dilaton potential. The constant B is then estimated to be B = (273−546 MeV) 4 , and the vacuum value of the dilaton field is obtained through the following relation [21], where the lowest glueball mass m χ = 1700 MeV is taken [34,35]. The mass eigenstates N ± are related to the ψ 1 and ψ 2 fields as follows: with In the diagonal basis, the masses of the chiral partners are given by From Eq. (10), it is clear that, in contrast to the naive assignment under chiral symmetry, the chiral symmetry breaking generates only the splitting between the two masses. When the symmetry is restored, the masses become degenerate according to m ± (σ = 0) = g χ χ. Note that the chirally invariant mass needs not to be a constant in the current model as the dilaton field χ is identified with the gluon condensate G µν G µν which varies with temperature and density.
To investigate the properties of strongly-interacting dense matter, we adopt a mean-field approximation [36]. Rotational invariance requires that the spatial component of the ω µ field vanishes, namely ω = 0 #1 . Parity conservation on the other hand dictates π = 0. The mean-field thermodynamic potential of the parity doublet model reads with where γ ± = 2 × 2 denotes the spin-isospin degeneracy factor for both parity partners, and f x (f x ) is the particle (antiparticle) Fermi-Dirac distribution function, with β being the inverse temperature, the dispersion relation E x = p 2 + m 2 x and the effective chemical potential This model is composed solely of hadronic degrees of freedom. In the next subsection we give its extension including further quark degrees of freedom.

B. Hybrid Quark-Meson-Nucleon model
The hybrid approach proposed in [21] is a natural extension of the parity doublet model in twofold meaning: First, it includes the quark degrees of freedom, and the fermionic Lagrangian in Eq. (2) is extended with a quarkmeson coupling as in the standard quark-meson model, The full model Lagrangian is then Note that the nucleons and quarks are coupled to the same mean fields generated from the potential V σ . On the other hand, the nature of the repulsive interaction among quarks is still far from consensus. In general, a #1 Since ω 0 is the only non-zero component in the mean-field approximation, we simply denote it by ω 0 ≡ ω.
coupling of quarks to the repulsive ω mean field can be treated as an additional parameter. In order to avoid unnecessary complexity and to reduce the number of free parameters, we do not take it into account in the current work. Second, the model realizes the concept of statistical confinement #2 through a medium-dependent modification of the Fermi-Dirac distribution functions, where an auxiliary scalar field b (bag field) is introduced. The distribution functions for the nucleons are replaced with and for the quarks, accordingly where α is a dimensionless model parameter. The functions f ± are given in Eq. (13), and for quarks they are defined as with the quark chemical potential µ q = 1 3 µ B , and the dispersion relation E q = p 2 + m 2 q . The real scalar field b is generated from a vacuum bag potential V b . Following [21], we take the potential of the form The potential (19) develops a non-trivial vacuum expectation value at b = κ 2 b /λ b . From Eqs. (16) and (17), one finds that the nucleons favor large b, whereas the quarks do small b. The potential (19) is chosen such that, at a certain T and µ B , a transition sets in, causing the bag-field expectation value to abruptly drop. As a consequence, at low T and µ B , the quark degrees of freedom are suppressed, while the nucleons get suppressed at high T and µ B . This characteristic behavior is associated with the deconfinement transition, which is a crucial feature of the model [21]. #2 In a class of chiral models implementing the Polyakov loop [37][38][39][40][41], the suppression of unphysical quarks in hadronic sector is achieved by the expectation value of the Polyakov loop Φ, which substantially modifies the statistical distribution function. However, the quarks do have a small contribution to the pressure even at a very low temperature. The quarks are thus unconfined at any temperature in this approach. We refer to the thermodynamic suppression of the quarks as statistical confinement [42]. We emphasize that the underlying symmetry of the potential V b must be discrete. This is because otherwise additional Nambu-Goldstone modes would emerge in the low-energy sector of QCD, other than pions, and they will spoil the known phenomenology of the chiral dynamics.
The mean-field thermodynamic potential of the Hybrid QMN model reads with where γ q = 2 × N c × N f denotes the spin-color-flavor degeneracy factor for quarks. In this study, N f = 2 and N c = 3, hence γ q = 12.
The thermal values of the mean fields are determined by minimizing the thermodynamic potential: where the scalar and baryon densities read and respectively. The boundary terms in the gap equation (22d) read for the nucleons and quarks, respectively. Note that the terms (25) and (26) come into the gap equation (22d) with opposite signs. This leads to that the nucleons and quarks favor different values of the bag field. In the grand canonical ensemble, the net-baryon number density for nucleons and quarks can be calculated as follows respectively. The total net-baryon number density is then the sum of all the terms in Eq. (27), The particle-density fractions are defined as As will be seen in the next section, traces of both chiral and deconfinement transitions will be nicely pronounced in the above fractional net-baryon number densities.

C. Determination of model parameters
The original Hybrid QMN model can deal with the dilaton dynamics through a dilaton-to-fermion coupling. However, as shown in Ref. [21], due to its heavy mass, the presence of the dilaton mean-field does not affect the nuclear groundstate properties. In fact, one finds that the expectation value of the dilaton remains nearly constant in the region of interest, and does not influence the model predictions. Hence, for the sake of simplicity of the calculations and the clarity of the discussion, we assume a constant value of the dilaton field, equal to its vacuum expectation value, χ 0 . As a result, the model becomes independent of the dilaton dynamics and the corresponding gap equation (22c) is irrelevant in the following discussion. The chirally invariant mass becomes a constant, i.e., m 0 = g χ χ 0 , and is treated as a parameter. The dilaton potential becomes irrelevant as well, thus can be omitted. Consequently, the rest of the potentials (5) are rewritten as follows: where the parameters λ 2 , λ 4 , and are modified to be Note that constant terms in V σ were dropped. The gap equations can be derived similarly to Eqs. (22). The positive-parity state corresponds to the nucleon N (938), with the vacuum mass m + = 938 MeV. Its negative-parity partner is identified as N (1535), with the real part of the mass pole in the range 1490 − 1530 MeV [43]. Here, following the previous studies in Refs. [10,21], we fix its mass to be m − = 1500 MeV. The chirally invariant mass is determined to be m 0 = 790 MeV, accordingly. Such high value of m 0 is also supported by recent findings in full lattice QCD simulations [5,6], at finite temperature and vanishing baryon chemical potential, where one sees clearly that the mass of the lowest positive-parity state is largely unaffected up to the pseudo-critical temperature, while a rather strong temperature-dependence is observed in the negative-parity channel. Also, the masses of the parity doublers are found to be nearly degenerate in the vicinity of the pseudo-critical temperature. Hence, the paritydoubling structure with m 0 ∼ m + (T = 0) is approximately realized because of the chiral symmetry restoration.
The parameters g 1 and g 2 are determined by the aforementioned vacuum nucleon masses and the chirally invariant mass m 0 through Eq. (10). The nuclear matter saturation properties at T = 0 set the following two constraints The above are used to fix the value of the vector coupling g ω and the mass of sigma meson m σ . The value of the quark coupling g q is fixed by assuming that the mass of the nucleon is m + = 3m q in the vacuum. For the parametrization of the bag potential V b , we adopt the values suggested in Ref. [21]. The αdependence is to be studied in the next sections. Its value is chosen in such a way that the UV cutoff αb 0 for the nucleon distribution function and the IR cutoff b 0 for the quark distribution function do not spoil the properties of the nuclear groundstate. This roughly sets the lower bound for αb 0 to be 300 MeV. On the other hand, even though the quark degrees of freedom are suppressed through Eq. (17), they can still be thermally excited at finite temperature, even before the point where the bag field effectuates the transition. Therefore, we impose that the quark density cannot exceed 1% of the total density of the system (Y q ≤ 0.01) before the deconfinement transition point. This prescription sets an upper limit for the α parameter. In general, the limit, set by such a constraint, should decrease with increasing temperature. This is due to the fact that at higher temperature quarks are more readily excited at lower values of baryon chemical potential. Hence, the point at which the constraint is met shifts towards lower values.
The model parameters to be used in the following discussion are summarized in Table I. In Section III, we discuss the role of the α parameter and its impact on the mean-field dynamics of the Hybrid QMN model, i.e, the chiral and deconfinement transitions, as well as various thermodynamic observables, including higher-order cumulants of the net-baryon number density.

III. EQUATION OF STATE
In this section, we study the influence of the external parameter α on the thermodynamic quantities in the Hybrid QMN model. In Ref. [21], the thermodynamics of the model was studied at T = 0 MeV, already indicating a sensitivity of the order of the chiral phase transition to the choice of α parameter. In the current work, we study this dependence at finite temperature.
In Fig. 1, we show the calculated mean-field expectation values at a fixed temperature, T = 10 MeV, for two values of the α parameter, namely αb 0 = 300 MeV (left panel) and αb 0 = 390 MeV (right panel). In both cases, the results exhibit similar behavior in the vicinity of the liquid-gas phase transition since the momentum cutoff was introduced not to spoil any properties of the nuclear ground state. In the former case, the onset of chiral restoration occurs at low baryon chemical potential µ B = 1137 MeV. After this point, the σ expectation value suddenly drops to nearly zero, causing the parity-doublet nucleons to become almost equally populated. This is clearly seen in the top panel of Fig. 2, where shown are the particle abundances defined in Eq. (29). The chiral phase transition is, in this case, a first-order. At this point, quarks are still confined up to the point where the bag field expectation value develops a jump, which triggers the nucleon suppression, allowing, in turn, for the quarks to be populated. This happens at µ B = 1573 MeV.
In the case of the latter scenario (right panel of Fig. 1), the chiral transition is a smooth crossover. The transition point, defined as a peak in ∂σ/∂µ B , is located at µ B = 1305 MeV. This behavior is also resembled in the density fractions, shown in the bottom panel of Fig. 2. The deconfinement transition, as in the former case, is induced by the jump in the bag field expectation value and happens at a higher chemical potential, namely at µ B = 1868 MeV.
The chiral and deconfinement transitions are reflected in the thermodynamic pressure plotted against the netbaryon number density. This is shown in Fig. 3. While the deconfinement in both cases is pronounced in the density jump of the order of 3 − 4 ρ 0 , the chiral transition in the case of αb 0 = 300 MeV is seen as a slightly smaller density jump, roughly of the order of ρ 0 . In case of αb 0 = 390 MeV there is no true chiral transition, but a crossover, which only softens the thermodynamic pressure.
In the left panel of Fig. 4 we show different temperature profiles of the phase diagram in the (α, µ B )-plane. In the right panel of Fig. 4, we show the phase diagram in the (α, ρ B )-plane, at T = 10 MeV. The liquidgas transition, as argued before, is not affected by the α parameter. Clearly, this is so on the figure. One of the major observations is that the chiral and deconfinement transitions are always separated by about 3 − 5 ρ 0 or 500 − 600 MeV baryon chemical potential (see the left panel of Fig. 4). The deconfinement transition, driven by the non-dynamical bag field, is always of the first-order transition. This is due to the fact that the potential (20) develops a sixth-order term b 6 at low temperatures. We note that the above phase structure is modified when the chiral invariant mass m 0 changes. As m 0 is decreased, the chiral and deconfinement transitions tend to take place closer to each other and eventually become almost simultaneous with m 0 ≈ 600 MeV for higher values of α. Such transition is then necessarily of the first-order, driven by the jump in the bag-field expectation value, and triggered at lower densities.

IV. NET-BARYON NUMBER DENSITY CUMULANTS
The fluctuations of conserved charges reveal more information about the matter composition than the equation of state, and can be used as probes of a phase boundary. The critical properties of chiral models are governed by the same universality as in QCD, i.e., the chiral transition belongs to the O(4) universality class, which, at large values of the baryon chemical potential, may develop a Z(2) critical point, followed by the firstorder transition [44][45][46]. This criticality is naturally encoded in the hadronic parity doublet model, as well as in quark-based models [28,29,[47][48][49]. Recall that the Hybrid QMN model, not only includes hadron and quark degrees of freedom, but also implements the statistical confinement, through the introduction of the auxiliary bag field. It is then constructive to explore the impact of both, the chiral symmetry and the statistical confinement on the higher-order cumulants in the vicinity the two transition lines, as well as to study their asymptotic behavior.
In the grand canonical ensemble, the cumulants are defined as higher-order derivatives of the thermodynamic potential with respect to different chemical potentials. In the current work, we are interested solely in the net-baryon number cumulants. They, and their rations, are defined as follows: They are often normalized by temperature, as done in lattice QCD studies [50,51]. Alternatively, one normalizes them by the chemical potential as which are more useful in our study focusing on a cold and dense medium. In Appendix A, we summarize the differences between those two normalizations, and discuss their restrictions and applicability in probing different regimes of the QCD phase diagram.

A. Critical behavior
It is expected that, in the vicinity of phase boundaries, and especially the critical point (CP), the generalized susceptibilities of conserved charges exhibit non-monotonic  behavior. This is governed by the singular part of the free energy and should be universal in all models that embody the QCD symmetry.
In the following we focus on a particular example of the second-order cumulant χ 2 , and discuss its properties in the vicinity of phase boundaries at T = 10 MeV and various µ B in the mean-field approximation. At fixed temperature, the order of the chiral phase transition in the Hybrid QMN model is tuned by the model parameter α. As seen in Fig. 4, for a given set of temperatures, small values of the α yield a first-order chiral phase transition, which eventually goes through a critical point and turns into a crossover for larger α. In general, the position of the CP in the (α, µ B )-plane depends on the temperature. At the CP, irrespectively of the temperature, the susceptibility of the net-baryon number density, χ 2 , should diverge with the critical exponent of the 3-dimensional Ising model [52,53].
To see this in the Hybrid QMN model, we calculate the critical exponent for χ 2 around the CP driven by the chiral dynamics. We take a path with fixed temperature of T = 10 MeV and parallel to the baryon chemical po-tential axis, approaching the CP from smaller values. In the mean-field approximation one expects the critical exponent to be 2/3, as the model-independent prediction. The exponent is obtained through a logarithmic fit, where r is a constant. The value obtained numerically, 0.67, stays in good agreement with the mean-field value = 2/3. Fig. 5 shows the numerically obtained net-baryon number susceptibility near the critical point, together with the fitted function given in Eq. (36).

B. Baryon number susceptibility
Here, we compare the results obtained in the Hybrid QMN model with those in the parity doublet model, introduced in Sec. II A. For the sake of completeness of the discussion, we show the results from the Nambu-Jona-Lasinio (NJL) model in which the relevant degrees of freedom are only constituent quarks. The mean-field thermodynamic potential of the NJL model in the three-momentum cutoff scheme is given by [54,55] where β is the inverse temperature, the bare quark mass is denoted by m q , the constituent mass by M , to the order parameter M . Finally, the higher-order cumulants are calculated according to Eq. (33).
The parameters used in the parity doublet model and in the NJL model are tabularized in Table II. The parametrization in the NJL model is chosen so that it yields the same order of the chiral transition, with its strength comparable to the one obtained in the Hybrid QMN model on the level of χ 2 observable.
The characteristics of the second-order cumulants, normalized by the baryon chemical potential, χ µ B 2 , in all three models are shown in Fig. 6. In the right panel, we show the results for αb 0 = 390 MeV, compared to the results obtained in other two models. The NJL result, shown as the blue solid line, undergoes a chiral crossover transition, seen as the peak in χ µ B 2 . After this point, it directly goes to the non-interacting quark limit from above. As discussed earlier, both the parity doublet (red dashed line) and Hybrid QMN (black dash-dotted line) models feature the first-order liquid-gas phase transition, which, at T = 10 MeV, is evident around µ B = 917 MeV. The parity doublet model develops a peak-like structure, indicating the chiral crossover, and then gradually goes to zero. The asymptotic behavior is due to the monotonically increasing repulsive interactions, through increasing ω mean-field expectation value. Hence, the parity doublet model does not saturate the non-interacting quark limit at high density. The Hybrid QMN model, on the other hand, shows the chiral peak at much lower value of the baryon chemical potential, resulting in a stronger transition. This is traced back to the momentum cutoff introduced for the nucleons through Eq. (16). When the cutoff yields a non-negligible suppression of the Fermi-Dirac distribution, the Hybrid QMN result starts to visibly deviate from the parity doublet counterpart. Eventually, at higher values of the baryon chemical potential, the Hybrid QMN model result features a discontinuity connected with the change of the degrees of freedom, from nucleons to quarks. Clearly, the discontinuity is due to the jump in the expectation value of the bag field, which abruptly suppresses the nucleon density and enhances the corresponding density of the quarks, allowing for the proper asymptotic behavior to be achieved.
The left panel of Fig. 6 shows the corresponding results on the χ 2 observable for the case of αb 0 = 300 MeV. In this case, all three models yield the first-order chiral phase transition. Similarly to the previous example, the parity doublet result does not saturate the proper asymptotic behavior, while the NJL result does. The second-order cumulant in the Hybrid QMN model deviates from the result obtained in the parity doublet model much earlier, immediately after the liquid-gas phase transition. Note that such behavior is expected, since the value of the α parameter in this case is the lowest that does not spoil the nuclear groundstate properties at zero temperature. In general, lower value of α yields stronger suppression of the thermal fluctuations of the nucleon degrees of freedom, which becomes relevant already at much lower values of the baryon chemical potential. As a result, it triggers the chiral phase transition at much lower µ B . At the same time, it allows for an earlier onset of the quark degrees of freedom, and eventually the deconfinement transition occurs at lower µ B as well.
By contrast, the second-order cumulant is not sensitive to the first-order deconfinement transition in the Hybrid QMN model, as seen in Fig 6. This is connected with the fact that the deconfinement transition is driven solely by the expectation value of the b field, which is generated through the bag potential defined in Eq. (19). Since the underlying symmetry of the potential V b is discrete, the deconfinement is dominated by the massive scalar field b. To see this, let us consider T = 0 limit for simplicity, and assume that the transition from nucleons to quarks happens at some high value of the baryon chemical potential. In this limit, we expect that the mean fields σ, ω → 0. Hence, the energies of the nucleons and quarks can be ignored. In such limit, one obtains the approximated expression for the gap equation of the b field from Eq. (22d), namely The above sets an upper limit for the α parameter, α max = 2 −1/3 . This implies αb 0 452 MeV at T = 0, which is consistent with the constraint introduced in Section II C.
Moreover, Eq. (38) can be solved for b. Taking into account only the positive solution, one gets the approximate expression, and further, ∂b/∂µ B , namely (40) It is evident that Eqs. (39) and (40) are smooth functions that do not exhibit any non-monotonic behavior and asymptotically, in the limit of high µ B , go to zero.
Under the above approximation, the derivative of the gap equation (38) with respect to the bag field is From the above it is clear that the bag field is massive, and that there is indeed no additional soft mode to couple to the baryon number density other than the σ mode. Hence, the higher-order cumulants do not develop any non-monotonic behaviors at the deconfinement transition.
More details about the characteristics of the second-order cumulant can be illustrated by considering the individual contributions from the apparent µ Bdependence in the Fermi gas part and intrinsic one in the mean-field sector. We rewrite the cumulants in Eq. (33) as follows where x, y = µ B , σ, ω, b. To see the explicit dependence and sensitivity of the second-order cumulant to the mean fields, we compare the diagonal terms from Eq. (42). The contribution of the σ and b mean fields for the case of αb 0 = 390 MeV are shown in Fig. 7 normalized by µ 2 B . The diagonal term χ σσ 2 yields a positive contribution around the chiral phase transition, while it remains insensitive to the deconfinement transition (the inset plot of Fig. 7). This is because the sigma field does not play any role in activating the quarks. On the other hand, the term χ bb 2 , shows a negative contribution and much weaker sensitivity to the chiral phase transition, and develops a jump due to the first-order deconfinement transition.
We note that when the model is naively extended to higher-temperature domain, a critical point for the deconfinement transitions builds up, and the corresponding second-order cumulant exhibits an additional peak structure as well. This is because the b field becomes a soft-mode at the second-order transition. Such extension is, however, beyond the current scope of the HQMN model, as it lacks, e.g., the implementation of the thermal fluctuations of mesons and gluons, which are relevant degrees of freedom at high temperature.

C. Higher-order cumulants
In Fig. 8, we show the first four cumulants of the netbaryon number density in the Hybrid QMN model, normalized by the baryon chemical potential in the vicinity of the chiral phase transition. Shown are the cumulants for different values of the α parameter, namely for the first-order transition (blue solid line), the second-order transition at the CP (red dashed line), and the crossover (black dash-dotted line). The x-axis is normalized by the corresponding values of the chiral transition point, µ c B , in each scenario.
In general, the first two cumulants, χ 1 and χ 2 , stay positive at all values of the baryon chemical potential. This changes for higher orders, and the third-and fourthorder cumulants can both turn negative.
In the case of the first-order transition, the third-order cumulant, seen in the bottom left panel of Fig. 8, is discontinues and changes its sign from positive to small negative values at the chiral transition point. For the second-order transition, the cumulant diverges from both sides of the transition point and stays positive in the chirally broken phase, while it turns negative in the restored phase. In the case of the crossover transition, χ µ B 3 is continuous and rapidly increases in the vicinity of the transition, then abruptly drops, exhibiting a dip at negative values.
As seen in the bottom right panel of Fig. 8, the fourthorder cumulant exhibits similar characteristic structure to the third-order cumulant, while it stays positive for the first-and second-order case. For the crossover case, similarly to the third-order cumulant, in the chirally broken phase χ µ B 4 develops a peak, then suddenly drops and exhibit a dip at negative values. It then turns positive again, in contrast to χ µ B 3 , and forms the second peak in the chirally restored phase. Note that, due to the deconfinement transition at higher values of baryon chemical potential, the cumulants reproduce the appropriate Stefan-Boltzmann limit. The signal of the transition is, however, less pronounced in the higher-order cumulants. This insignificance was already evident in the second-order cumulant.
We note that the above characteristics, as expected, are qualitatively similar to the results obtained in different chiral models that incorporate the concept of statistical confinement in a different scheme, e.g., the Polyakov loop-extended quark-meson model [29]. This is because the transitions remain well separated at low temperatures, and so are the critical behaviors in the vicinity of each critical point. When the model is appropriately extended to higher temperatures, maximally three critical points could appear on the phase diagram. If they still remain separated, the cumulants, as well as their ratios, would exhibit a structure similar to the one discussed in the current work. If, on the other hand, the critical points would appear in a more complex way, e.g., closer to each other, a non-trivial structure could appear. We leave this is issue for a future work.

V. SUMMARY AND CONCLUSIONS
We have explored, under the mean-fieldapproximation, thermodynamics of the Hybrid Quark-Meson-Nucleon (Hybrid QMN) model at low temperature and various baryon chemical potential. The model describes nuclear matter in terms of the nucleon paritydoublers and mesons, and quark matter in the standard chiral-quark framework. The quark and nucleon sectors are connected via an auxiliary non-dynamical scalar field b, generated through a scalar potential V b . The role of this field is to suppress unphysical thermal fluctuations of quarks at low T and µ B and, simultaneously, those of nucleons at high T and µ B . This is achieved by the modified Fermi-Dirac distribution functions, where the α parameter is introduced to manipulate a size of the momentum cutoff. We emphasize that the inclusion of the cutoff does not spoil the universality and the critical behaviors.
We have studied the impact of the α parameter on the phase structure in the Hybrid QMN model at finite temperature. We find that the chiral and deconfinement transitions are always separated for a large m 0 . The α parameter plays a crucial role in tuning the order of the chiral phase transition. Namely, the system undergoes a first-order chiral phase transition for low values of α and goes through a critical point, to eventually turn into a smooth crossover at higher α. The deconfinement transition is, on the other hand, always of first order. This is due to the specific choice of the potential V b , which by construction exhibits a first-order transition at low temperature. The separation of the two transitions might indicate the quarkyonic phase, where the quarks are confined, but the chiral symmetry is restored [56][57][58].
We have also studied the critical behavior of the higher-order cumulants up to the fourth order, for different choice of the model parameter α. It is found that the impact of the suppression of the nucleon degrees of freedom on the chiral phase transition is twofold. First, it triggers the transition at much earlier baryon chemical potential, and second the transition is strengthened in comparison to the pure hadronic-model result. Since the changeover from nucleons to quarks is done self-consistently in the Hybrid QMN model, the observables properly saturate the asymptotic Stefan-Boltzmann limit. However, the deconfinement transition does not yield any non-monotonic behaviors.
In this study, we considered the two-flavored system, and it is a natural extension of the work to a theory including heavier flavors. The parity-doubling structure of baryons with three flavors was recently formulated in effective chiral approaches [16,59]. It was shown that the masses of the baryon parity-doublers measured on lattice [5,6] are modified to a large extent with the physical pion mass, in particular in the hyperon channels [59]. It is of great interest to establish the equations of state of nuclear/hyperon matter in a physical setup.
In view of the recent success of effective quark-hadron models in astrophysics [60][61][62], in particular in order to be consistent with current constraints, e.g., the observation of a massive two-solar-mass pulsar [63], repulsive interactions among quarks are essential. Moreover, in order to employ the Hybrid QMN model in astrophysical studies of compact stellar objects, such as neutron (or hybrid) stars, proto-neutron stars and supernovae [64][65][66], it is essential to extend the model to arbitrary isospin chemical potential. Further studies of these issues will be reported elsewhere. The potential of the b field is so far not anchored to any QCD symmetry. It is indispensable to establish its dynamics in terms of a reliable symmetry. This will lead to a better understanding of its role as an order parameter for the deconfinement phase transition. In fact, it has been suggested that such a non-trivial order parameter of the center-flavor symmetry exists in a QCD-like theory compactified on a circle [67,68]. This symmetry may be used to constrain further our model. In the following, we discuss the asymptotic behavior of the observables constructed from the net-baryon number cumulants and their ratios, depending on the normalization.
At low temperature and low baryon chemical potential the main degrees of freedom, due to confinement, are baryons, with baryon number B = ±1, while at high temperature quarks become dominant, with B = ± 1 3 . When modeling the thermodynamics of QCD, especially when including the quark degrees of freedom, it is required that, at high temperature and high baryon chemical potential, the thermodynamic quantities reach the appropriate non-interacting quark gas limit. The general expression for the thermodynamic potential of the non-interacting gas of massless particles, carrying baryon quantum number B, reads with γ being the spin degeneracy factor. The first four cumulants read Note that the fourth-order cumulant is already dimensionless and independent of the temperature and the baryon chemical potential. Moreover, it is directly proportional to the baryon quantum number B. Hence, χ 4 itself can be a useful probe of the state of matter, in both hot and dense limit.
Using the temperature-normalization scheme, defined in Eq. (34), the asymptotic behavior in the high-temperature limit is the following All of the four cumulants have well defined limits, but only the even ones depend on the baryon quantum number B. The odd cumulants vanish in the high temperature limit, owing to the baryon-antibaryon symmetry.
On the other hand, in the limit of high baryon chemical potential, the first three cumulants go to infinity, and only the fourth one is proportional to B, namely This suggests that, in order to probe the content of the QCD matter at high density, only the fourth cumulant can be utilized. The high-temperature limits of χ µ B n , defined in Eq. (35), read lim T →∞ lim T →∞ lim T →∞ and consequently, the limits of high baryon chemical potential, Here, the situation is quite opposite. Namely, the highdensity limit is well defined and are strictly proportional to the baryon quantum number B for all four cumulants. The hot limit, on the other hand, is well defined only for the third and the fourth cumulants. Very useful are also the observables formed from ratios of the net-baryon density cumulants. Here, we discuss an example of the ratio of the fourth-order to the second-order cumulant, the so-called kurtosis [50]. At low temperature, where the Boltzmann approximation can be applied, the grand canonical ensemble can be approximated by Ω −F (T ) cosh (Bµ B /T ). From this, the kurtosis reads R 4,2 = B 2 /T 2 . At high temperature, one can use Eqs. (A2) to obtain R 4,2 = 6B 2 /(3µ 2 B B 2 + π 2 T 2 ). Using the normalization by temperature, the Boltzmann approximation yields and the high-temperature limit, where the factor 6/π 2 is due to the quantum statistics. The ratio R T 4,2 is therefore roughly proportional to the square of the baryon number of the main degrees of freedom. Hence, it may serve as a very good probe for the proper determination of the state of QCD matter in the regime of low baryon chemical potential.
The same quantity, however, turns out to be less useful in the low-temperature and dense regime. In the highbaron-chemical-potential limit, the ratio is independent of the baryon quantum number B, and goes strictly to zero, regardless of which degrees of freedom are dominant. Similar result is obtained with normalization by the baryon chemical potential, where one finds that both hot and dense limits yield results independent of the content of matter, namely In fact, it can be easily checked that all the ratios, defined through Eq. (35), yield limits independent of the baryon quantum number B. Nevertheless, still useful are the cumulants themselves, which at high density should approach the non-interacting gas counterpart results. Note that these are directly proportional to the baryon quantum number B and therefore may be used to probe the state of QCD matter at high densities.
In summary, the higher order cumulants of the netbaryon number density and their ratios may be utilized as very useful probes of the state of QCD matter. However, the appropriate observables and their normalization scheme have to chosen depending on the (T , µ B )-regime of interest.