Collision energy dependence of the critical end point from baryon number fluctuations in the Linear Sigma Model with quarks

The phase structure of strongly interacting matter in the temperature and baryon density plane, has become a subject with growing attention over the last years. The interest is driven by the possibility to experimentally probe this phase structure by means of relativistic heavy-ion collisions using either current experimental facilities, in particular the RHIC-STAR detector with its Beam Energy Scan program and HADES, which already have shown results on net-proton number fluctuations Adam, J. and others (2021); Adamczewski-Musch, J. and others (2020), or dedicated facilities soon to enter into operation, such as NICA Kekelidze et al. (2017) and FAIR Senger, P. (2017). The subject has received a boost in interest after the successful detection of gravitational waves from neutron star mergers. In fact, simulations show that these mergers probe similar regions in the phase diagram as those that can be studied in relativistic heavy-ion collisions Blacker et al. (2020); Most et al. (2020).

Calculations using Lattice QCD (LQCD) have found that for finite temperature $T$ and vanishing baryon chemical potential $\mu_{B}$, the transition between the confined/chiral symmetry broken and the deconfined/partially restored chirally symmetric phases, is a crossover that happens at a pseudocritical temperature $T_{c}(\mu_{B}=0)\simeq 158$ MeV Borsanyi et al. (2020); Bazavov et al. (2019); Aoki et al. (2006). However, using effective model calculations Roessner, Simon and Ratti, Claudia and Weise, W. (2007); Ayala et al. (2020); Asakawa and Yazaki (1989); Ayala et al. (2018); Gao and Liu (2016); Gao and Pawlowski (2020) it can be shown that this transition becomes first order for low $T$ and high $\mu_{B}$. Thus, when the temperature is increased, the end of this first order phase transition line in the $T$ vs. $\mu_{B}$ plane should happen at a Critical End Point (CEP). Unfortunately, LQCD calculations cannot be used to directly determine the position of this CEP due to the severe sign problem Ding et al. (2015), although early attempts where made some time ago in Ref. Fodor, Z. and Katz, S. D. (2004). More recently, results employing the Taylor series expansion around $\mu_{B}=0$ or the extrapolation from imaginary to real $\mu_{B}$ values, suggest that the CEP has not yet been found for $\mu_{B}/T\leq 2$ and $145\leq T\leq 155$ MeV Sharma (2017); Borsanyi et al. (2020).

Microscopic model calculations indicate that an energy density of 0.5 GeV/fm³ is able to produce a transition from the hadronic phase to the phase composed of deconfined quarks and gluons and that this density can be achieved in the center of the fireball created in head-on collisions of heavy-ions at energies above $\sqrt{s_{NN}}=3-5$ GeV Mendenhall and Lin (2021). Since for low collision energies, nuclear stopping also produces a high baryon density in the interaction region, it is expected that collisions at such energies could find signals related to the appearance of the CEP.

The strategy to locate the CEP is based on finding deviations from the statistical behavior of a gas made out of resonances, which is known to also describe the crossover transition line found by LQCD at low values of $\mu_{B}$ Karsch (2017). A simple description of the transition line based on resonance degrees of freedom is provided by the Hadron Resonance Gas Model (HRGM), when the occupation numbers are given in terms of Boltzmann statistics Braun-Munzinger et al. (2003). Recall that the statistical properties of a thermal system can be described in terms of the cumulants of its conserved charges, which are extensive quantities Asakawa and Kitazawa (2016). To avoid uncertainties introduced by volume effects, the analyses consider ratios of these cumulants. The strategy is therefore linked to finding deviations in the ratios of these cumulants from those obtained in the HRGM, which are described by the Skellam distribution where, in particular, the ratios of cumulants of even order are found to be equal to 1. Of the conserved quantities that can be experimentally monitored, baryon number is accessible through the measurement of proton multiplicities Adam et al. (2020); Abdallah et al. (2021a, b); Adamczewski-Musch, J. and others (2020). The search for the location of the CEP is thus guided by the expected emergence of critical behaviour Hatta and Ikeda (2003); Stephanov et al. (1999); Bzdak et al. (2017, 2020); Athanasiou et al. (2010); Mroczek et al. (2021) on this and other conserved charges such as electric charge and strangeness, when a collision energy scan is performed. Recent analyses of the possible location of the CEP, that also resort to computing high order cumulants, include techniques such as Dyson-Schwinger equations Isserstedt, Philipp and Buballa, Michael and Fischer, Christian S. and Gunkel, Pascal J. (2019), a QCD assisted effective model Fu, Wei-jie and Luo, Xiaofeng and Pawlowski, Jan M. and Rennecke, Fabian and Wen, Rui and Yin, Shi (2021) and holographic methods Grefa, Joaquin and Noronha, Jorge and Noronha-Hostler, Jacquelyn and Portillo, Israel and Ratti, Claudia and Rougemont, Romulo (2021).

In this work we study the evolution of cumulant ratios describing baryon number fluctuations as a function of the collision energy, in the context of relativistic heavy-ion collisions, using the Linear Sigma Model with quarks (LSMq). The model incorporates a crucial feature of plasmas near transition lines, namely, the plasma screening effects. These are encoded in the contribution of the ring diagrams into the effective potential which becomes a function of the order parameter after spontaneous chiral symmetry breaking. Therefore, the statistical properties of the system can be formulated in terms of fluctuations of this order parameter Stephanov (2009, 2011) when the collision energy $\sqrt{s_{NN}}$ and thus $T$ and $\mu_{B}$, are varied. We expand on the findings of Refs. Ayala et al. (2020, 2018, 2016) where the analysis was based only on the properties of the effective potential, also up to ring diagrams order, without studying the higher order fluctuations predicted by the model and where no attempt to locate the position of the CEP as a function of the collision energy was made. In order to gain analytical insight, we resort to the high temperature approximation and work in the chiral limit. Although these approximations have limitations concerning the model accuracy for the location of the CEP, they constitute a useful guide to carry out future, more precise studies. The work is organized as follows: In Sec. II we describe the LSMq and compute the effective potential up to the ring diagrams contribution, which requires computation of the self-energies corresponding to the meson degrees of freedom. We fix the model parameters requiring that the phase transition line near $\mu_{B}=0$ corresponds to the one found by the latest LQCD calculations Borsanyi et al. (2020). In Sec. III we formulate the way the baryon number fluctuations can be described in terms of the probability distribution associated to the order parameter near the transition line and present the results of the analysis for the ratios of cumulants describing baryon number fluctuations as a function of $\sqrt{s_{NN}}$, showing that these ratios deviate from the expectations of the HRGM for energies around $\sqrt{s_{NN}}\sim 4-6$ GeV and that the CEP can be found at energies $\sqrt{s_{NN}}\sim 2$ GeV. The model thus predicts that the CEP can be found either in the lowest NICA or within the HADES energy domain. We finally summarize and conclude in Sec. IV.

Effective models are useful to identify the main characteristics of the QCD phase diagram. Although no single model can be used to describe the whole phase diagram, they can be employed to explore different regions when including different degrees of freedom. For instance, the Polyakov loop extended Nambu–Jona-Lasinio or quark-meson models have been employed in Refs. Kashiwa, Kouji and Kouno, Hiroaki and Matsuzaki, Masayuki and Yahiro, Masanobu (2008); Mao, Hong and Jin, Jinshuang and Huang, Mei (2010); Skokov, V. and Stokic, B. and Friman, B. and Redlich, K. (2010); Skokov, V. and Friman, B. and Redlich, K. (2011) to study simultaneously the deconfinement and chiral symmetry restoration aspects of QCD. Given that LQCD calculations, extended to small values of $\mu_{B}$, find coincident transition lines for the deconfinement and chiral symmetry restoration transitions, it should be possible to explore the phase diagram emphasizing independently either the deconfinement or the chiral aspects of the transition.

An effective model based on chiral symmetry is provided by the LSMq. The Lagrangian is given by

	$\displaystyle\!\!\!\!\!\!\mathcal{L}$	$\displaystyle=$	$\displaystyle\frac{1}{2}(\partial_{\mu}\sigma)^{2}+\frac{1}{2}(\partial_{\mu}\vec{\pi})^{2}+\frac{a^{2}}{2}(\sigma^{2}+\vec{\pi}^{2})$		(1)
		$\displaystyle-$	$\displaystyle\frac{\lambda}{4}(\sigma^{2}+\vec{\pi}^{2})^{2}+i\bar{\psi}\gamma^{\mu}\partial_{\mu}\psi-g\bar{\psi}(\sigma+i\gamma_{5}\vec{\tau}\cdot\vec{\pi})\psi,$

where $\psi$ is an $SU(2)$ isospin doublet of quarks,

$\displaystyle\vec{\pi}=(\pi_{1},\pi_{2},\pi_{3}),$

(2)

and $\sigma$ are an isospin triplet and singlet, corresponding to the three pions and the sigma, respectively. The squared mass parameter $a^{2}$ and the self-coupling $\lambda$ and $g$ are taken to be positive and, for the purpose of describing the chiral phase transition at finite $T$ and $\mu_{B}$, they need to be determined from conditions close to the phase boundary, and not from vacuum conditions. The original LSM, was introduced to describe the interaction between pions nucleons in the seminal work of Ref. Gell-Mann, Murray and Levy, M (1960). A review of the LSM at finite temperature is provided in Ref. Petropoulos (1999). An early study of the phase diagram using LSMq is provided in Ref. Bowman, E. Scott and Kapusta, Joseph I. (2009).

Working in the strict chiral limit, to allow for a spontaneous symmetry breaking, we let the $\sigma$ field to develop a vacuum expectation value $v$, namely

$\displaystyle\sigma\rightarrow\sigma+v.$

(3)

This vacuum expectation value can be identified with the order parameter of the theory. After this shift, the Lagrangian can be rewritten as

	$\displaystyle{\mathcal{L}}$	$\displaystyle=$	$\displaystyle\frac{1}{2}(\partial_{\mu}\sigma)^{2}+\frac{1}{2}(\partial_{\mu}\vec{\pi})^{2}-\frac{1}{2}\left(3\lambda v^{2}-a^{2}\right)\sigma^{2}$		(4)
		$\displaystyle-$	$\displaystyle\frac{1}{2}\left(\lambda v^{2}-a^{2}\right)\vec{\pi}^{2}$
		$\displaystyle+$	$\displaystyle\frac{a^{2}}{2}v^{2}-\frac{\lambda}{4}v^{4}+i\bar{\psi}\gamma^{\mu}\partial_{\mu}\psi-gv\bar{\psi}\psi+{\mathcal{L}}_{I}^{b}+{\mathcal{L}}_{I}^{f},$

where ${\mathcal{L}}_{I}^{b}$ and ${\mathcal{L}}_{I}^{f}$ are given by

	$\displaystyle{\mathcal{L}}_{I}^{b}$	$\displaystyle=$	$\displaystyle-\frac{\lambda}{4}\Big{[}(\sigma^{2}+(\pi^{0})^{2})^{2}+4\pi^{+}\pi^{-}(\sigma^{2}+(\pi^{0})^{2}+\pi^{+}\pi^{-})\Big{]},$
	$\displaystyle{\mathcal{L}}_{I}^{f}$	$\displaystyle=$	$\displaystyle-g\bar{\psi}(\sigma+i\gamma_{5}\vec{\tau}\cdot\vec{\pi})\psi.$		(5)

The terms given in Eq. (5) describe the interactions among the fields $\sigma$, $\vec{\pi}$ and $\psi$, after symmetry breaking. Since no external agent that could couple to the electric charge, such as a magnetic field, is considered, there is no distinction between neutral and charged pions. From Eq. (4) one can see that the $\sigma$, the three pions and the quarks have masses given, respectively, by

	$\displaystyle m^{2}_{\sigma}$	$\displaystyle=$	$\displaystyle 3\lambda v^{2}-a^{2},$
	$\displaystyle m^{2}_{\pi}$	$\displaystyle=$	$\displaystyle\lambda v^{2}-a^{2},$
	$\displaystyle m_{f}$	$\displaystyle=$	$\displaystyle gv.$		(6)

Notice that the square of the tree-level boson masses in Eq. (6) can vanish or even become negative as the continuous order parameter $v$ spans its domain. The physical masses, including their thermal part, are computed when $v$ takes on the value found from minimizing the effective potential. In this work we have not focused explicitly on describing the behavior of these thermal masses. Nevertheless, our results are perfectly consistent with results of previous works such as the one in Ref. Scavenius et al. (2001) when including the thermal contribution to these masses, which are provided by their self-energies (see Eq. (16) below), albeit, in our case, in the high temperature limit.

In order to determine the chiral symmetry restoration conditions as function of $T$ and $\mu_{B}$, we study the behavior of the effective potential which includes the classical potential or tree-level contribution, the one-loop correction both for bosons and fermions and the ring diagrams contribution, which accounts for the plasma screening effects.

The tree-level potential is given by

$$V^{\text{tree}}(v)=-\frac{a^{2}}{2}v^{2}+\frac{\lambda}{4}v^{4},$$

(7)

whose minimum is found at

$$v_{0}=\sqrt{\frac{a^{2}}{\lambda}}.$$

(8)

Since $v_{0}\neq 0$, we notice that chiral symmetry is spontaneously broken. However, in order to make sure that the quantum corrections at finite temperature and density maintain the general properties of the effective potential, we need to add counter-terms $\delta a^{2}$ and $\delta\lambda$ to the bare constants $a^{2}$ and $\lambda$, respectively, and write

	$\displaystyle V^{\text{tree}}$	$\displaystyle=-\frac{a^{2}}{2}v^{2}+\frac{\lambda}{4}v^{4}\rightarrow$
	$\displaystyle V^{\text{vac}}$	$\displaystyle=-\frac{(a^{2}+\delta a^{2})}{2}v^{2}+\frac{(\lambda+\delta\lambda)}{4}v^{4}+V^{\text{1;vac}},$		(9)

where $V^{\text{1;vac}}$ is the boson and fermion one-loop vacuum contribution. These counterterms need to be fixed from the conditions to keep the vacuum and its curvature fixed at their tree-value levels, namely,

	$\displaystyle\frac{1}{v}\frac{dV^{\text{vac}}}{dv}\Big{|}_{v=v_{0}}=0,$
	$\displaystyle\frac{d^{2}V^{\text{vac}}}{dv^{2}}\Big{|}_{v=v_{0}}=2a^{2}.$		(10)

The procedure is dubbed vacuum stabilization. It can be shown Ayala et al. (2021) that after this procedure is implemented, the vacuum potential $V^{\text{vac}}$ coincides with tree level potential $V^{\text{tree}}$ near the minimum and that $V^{\text{vac}}$ is independent of the chosen ultraviolet scale. Therefore, when working near the minimum, as is the case of this work, in order to avoid writing long analytical expressions, we keep working with the vacuum potential at tree level which is the practical result after implementing vacuum stabilization.

To include quantum corrections at finite temperature and baryon density, we work within the imaginary-time formalism of thermal field theory. The general expression for the one-loop boson contribution can be written as

$$V^{\text{b}}(v,T)=T\sum_{n}\int\frac{d^{3}k}{(2\pi)^{3}}\ln D_{\text{b}}(\omega_{n},\vec{k})^{1/2},$$

(11)

where

$$D_{\text{b}}(\omega_{n},\vec{k})=\frac{1}{\omega_{n}^{2}+k^{2}+m_{b}^{2}},$$

(12)

is the free boson propagator with $m_{b}^{2}$ being the square of the boson mass and $\omega_{n}=2n\pi T$ the Matsubara frequencies for boson fields.

For a fermion field with mass $m_{f}$, the general expression for the one-loop correction at finite temperature and quark chemical potential $\mu=\mu_{B}/3$ is

$$V^{\text{f}}(v,T,\mu)=-T\sum_{n}\int\frac{d^{3}k}{(2\pi)^{3}}\text{Tr}[\ln S_{\text{f}}(\tilde{\omega}_{n},\vec{k})^{-1}],$$

(13)

where

$$S_{\text{f}}(\tilde{\omega}_{n},\vec{k})=\frac{1}{\gamma_{0}(\tilde{\omega}_{n}+i\mu)+\not{k}+m_{f}},$$

(14)

is the free fermion propagator and $\tilde{\omega}_{n}=(2n+1)\pi T$ are the Matsubara frequencies for fermion fields.

The ring diagrams term is given by

$\displaystyle\!\!\!\!\!\!\!\!\!\!\!V^{\text{Ring}}(v,T,\mu)=\frac{T}{2}\sum_{n}\int\frac{d^{3}k}{(2\pi)^{3}}\ln[1+\Pi_{\text{b}}D_{b}(\omega_{n},\vec{k})],$

(15)

where $\Pi_{\text{b}}$ is the boson self-energy. The self-energies for the sigma and pion fields are in general different. However, in this work we resort to the high-temperature approximation, keeping only the leading matter effects. In this approximation the boson self-energies become independent of the boson species and we can write Ayala et al. (2021)

	$\displaystyle\!\!\!\!\!\Pi_{\text{b}}$	$\displaystyle\equiv$	$\displaystyle\Pi_{\sigma}=\Pi_{\pi^{\pm}}=\Pi_{\pi^{0}}$		(16)
		$\displaystyle=$	$\displaystyle\lambda\frac{T^{2}}{2}-N_{\text{f}}N_{\text{c}}g^{2}\frac{T^{2}}{\pi^{2}}\left[\text{Li}_{2}\left(-e^{-\frac{\mu}{T}}\right)+\text{Li}_{2}\left(-e^{\frac{\mu}{T}}\right)\right],$

where $N_{\text{c}}$ and $N_{\text{f}}$ are the number of colors and flavors, respectively, whereas ${\rm Li}_{2}(x)$ stands for the polylogarithm function of order 2 (dilogarithm).

Writing all the ingredients together, the finite temperature and baryon density effective potential, up to the ring diagrams contribution, after vacuum stabilization and boson and fermion mass renormalizaton in the $\overline{\text{MS}}$ scheme at the ultraviolet scale $\tilde{\mu}$, can be written as

	$\displaystyle V^{\text{eff}}(v)$	$\displaystyle=$	$\displaystyle-\frac{a^{2}}{2}v^{2}+\frac{\lambda}{4}v^{4}+\sum_{\text{b}=\pi^{\pm},\pi^{0},\sigma}\left\{-\frac{T^{4}\pi^{2}}{90}+\frac{T^{2}m_{\text{b}}^{2}}{24}-\frac{T(m_{\text{b}}^{2}+\Pi_{\text{b}})^{3/2}}{12\pi}-\frac{m_{\text{b}}^{4}}{64\pi^{2}}\left[\ln\left(\frac{\tilde{\mu}^{2}}{(4\pi T)^{2}}\right)+2\gamma_{E}\right]\right\}$		(17)
		$\displaystyle+$	$\displaystyle N_{\text{c}}N_{\text{f}}\left\{\frac{m_{f}^{4}}{16\pi^{2}}\Big{[}\ln\left(\frac{\tilde{\mu}^{2}}{T^{2}}\right)-\psi^{0}\left(\frac{1}{2}+\frac{i\mu}{2\pi T}\right)-\psi^{0}\left(\frac{1}{2}-\frac{i\mu}{2\pi T}\right)\right.+\psi^{0}\left(\frac{3}{2}\right)-2\left(1+\ln(2\pi)\right)+\gamma_{E}\Big{]}$
		$\displaystyle-$	$\displaystyle\left.\frac{m_{\text{f}}^{2}T^{2}}{2\pi^{2}}\left[\text{Li}_{2}\left(-e^{-\frac{\mu}{T}}\right)+\text{Li}_{2}\left(-e^{\frac{\mu}{T}}\right)\right]+\frac{T^{4}}{\pi^{2}}\left[\text{Li}_{4}\left(-e^{-\frac{\mu}{T}}\right)+\text{Li}_{4}\left(-e^{\frac{\mu}{T}}\right)\right]\right\},$

$\gamma_{E}\simeq 0.57721$ denoting the Euler-Mascheroni constant. Equation (17), with the boson self-energies given in Eq. (16), constitute the main tools to study the effective QCD phase diagram from the point of view of chiral symmetry restoration/breaking. Notice that inclusion of the plasma screening effects, encoded in the contribution of the ring diagrams, produces the possibility that the effective potential develops a cubic term in the order parameter $v$. This happens for values of $T$ and $\mu_{B}$ such that $\Pi_{\text{b}}-a^{2}\simeq 0$. This signals the possibility of a first order phase transition and of fluctuations that deviate from those of a free gas of mesons and quarks for those values of $T$ and $\mu_{B}$. Recall that the boson self-energy represents the matter correction to the boson mass. For a second order (our proxy for a crossover), namely, a continuous phase transition, these corrections should produce that the thermal boson masses vanish when the symmetry is restored. This means that at the phase transition, the effective potential not only develops a minimum but it is also flat (the second derivative vanishes) at $v=0$. This property can be exploited to find the suitable values of the model parameters $a$, $\lambda$ and $g$ at the critical temperature $T_{c}^{0}$ for $\mu_{B}^{c}=0$. Since the thermal boson masses are degenerate at $v=0$, the condition to produce a flat effective potential at $T_{c}$ for $\mu_{B}=0$ can be written as Ayala et al. (2021)

	$\displaystyle 6$	$\displaystyle\lambda$	$\displaystyle\left(\frac{T_{c}^{2}}{12}-\frac{T_{c}}{4\pi}\left(\Pi_{\text{b}}-a^{2}\right)^{1/2}\right.$		(18)
		$\displaystyle+$	$\displaystyle\left.\frac{a^{2}}{16\pi^{2}}\left[\ln\left(\frac{\tilde{\mu}^{2}}{(4\pi T_{c})^{2}}\right)+2\gamma_{E}\right]\right)$
		$\displaystyle+$	$\displaystyle g^{2}T_{c}^{2}-a^{2}=0.$

where, from Eq. (16),

$\displaystyle\Pi_{\text{b}}(T_{c}^{0},\mu_{B}^{c}=0)=\left[\frac{\lambda}{2}+g^{2}\right](T_{c}^{0})^{2}.$

(19)

Given that we have a relative freedom to chose the value of $\tilde{\mu}$, we take $\tilde{\mu}=500$ MeV, which is chosen large enough so as to consider it to be the largest energy scale in the problem. The dependence of the result is only logarithmic in $\tilde{\mu}$, therefore small variations in this parameter do not affect significantly the result. We look for solutions of Eq. (18) for the set of parameters $\lambda$, $g$ and $a$, using input form the LQCD transition curve at $T_{c}^{0}\simeq 158$ MeV from where we take the values of the curvature parameters $\kappa_{2}$ and $\kappa_{4}$. Since Eq. (18) is non-linear, the set of parameters is not unique. We work with the parametrization of the LQCD transition curve given by

$$\frac{T_{c}(\mu_{B})}{T_{c}^{0}}=1-\kappa_{2}\left(\frac{\mu_{B}}{T_{c}^{0}}\right)^{2}+\kappa_{4}\left(\frac{\mu_{B}}{T_{c}^{0}}\right)^{4},$$

(20)

with $\kappa_{2}=0.0153$ and $\kappa_{4}=0.00032$ Borsanyi et al. (2020); Guenther (2021). Notice that although Eq. (20) is provided in Ref. Borsanyi et al. (2020) as an implicit equation, for simplicity and in order to avoid a further layer of complexity into the analysis, we have normalized the variable $\mu_{B}$ to $T_{c}^{0}$. To find the values of the parameters, we first fix a value of $\lambda$, and find the solution of Eq. (18) for $g$ and $a$ that, from the effective potential in Eq, (17), produce a phase transition at the values of $T_{c}(\mu_{B})$, hereafter simply referred to as $T_{c}$, and $\mu_{B}$ given by Eq. (20). We then repeat the procedure for other values of $\lambda$. In this manner, the solutions can be expressed as a relation between the couplings $g$ and $\lambda$

$$g(\lambda)=0.31+1.94\lambda-2.06\lambda^{2}+0.97\lambda^{3}-0.20\lambda^{4},$$

(21)

and a relation between $a$ and $\lambda$

$\displaystyle\!\!\!\!\!\frac{a(\lambda)}{T_{c}}$

$\displaystyle=$

$\displaystyle 0.35+2.08\lambda-2.21\lambda^{2}+1.03\lambda^{3}-0.21\lambda^{4}.$

(22)

Figure 1 shows the values of the parameters $g$ and $a/T_{c}$ for the range of allowed values of $\lambda$ where Eq. (18) produces solutions. With all the parameters fixed, we can study the properties of the effective potential to find the values of $T_{c}$ and $\mu_{B}^{c}$ where the chiral phase transition takes place. Figure 2 shows the effective potential as a function of the order parameter. We take as examples three different sets of values of $T_{c}$ and $\mu_{B}^{c}$ along the transition curve using $a=148.73$ MeV, $\lambda=1.4$ and $g=0.88$. Notice that for $\mu_{B}=0$ and $T_{c}=158$ MeV the phase transition is second order. As $\mu_{B}$ increases, the phase transition is signaled by a flatter effective potential until the chemical potential and temperature reach values $\mu_{B}^{\text{CEP}}$ and $T^{\text{CEP}}$ where the effective potential develops a barrier between degenerate minima. For $\mu_{B}^{c}>\mu_{B}^{\text{CEP}}$ and $T_{c}<T^{\text{CEP}}$, the phase transitions are always first order. These features are illustrated in Fig. 3 where we show two examples of the effective phase diagram thus obtained. The upper panel is computed using $\lambda=1.4$, $g=0.88$ and $a=148.73$ MeV.

The lower panel is computed with $\lambda=0.4$, $g=0.82$ and $a=141.38$ MeV. The solid (red) lines represent second order phase transitions, our proxy for crossover transitions, whereas the dotted (blue) lines correspond to first order phase transitions. The star symbol represents the location of the CEP. Notice that the CEP appears only for large values of $\mu_{B}$. Also, variations of the model parameters do not change the CEP position appreciably. In fact, for the allowed range of values of $a$, $\lambda$ and $g$, the CEP location ranges between 786 MeV $<\mu_{B}^{\text{CEP}}<849$ MeV and 69 MeV $<T^{\text{CEP}}<70.3$ MeV.

Notice that the high temperature approximation is justified, in spite of the large values found for the $\mu_{B}^{\text{CEP}}$. This happens because the quantity that directly enters the calculation is the quark chemical potential $\mu$, related to $\mu_{B}$ by $\mu_{B}=3\mu$. Thus, although $\mu_{B}^{\text{CEP}}$ is rather large, $\mu^{\text{CEP}}$ is only $\sim 1.7\ T_{c}$. For this value of $\mu^{\text{CEP}}$, the combination of polylogarithmic functions and factors of $T$ in the second term of the right-hand side of Eq. (16) makes that the boson masses become $\sim\sqrt{\lambda+2g^{2}}\ T^{\text{CEP}}$. At the same time, $T^{\text{CEP}}\sim\ 0.5T_{c}$. Thus, the large temperature approximation is still valid around $(\mu_{B}^{\text{CEP}},T^{\text{CEP}})$. For even larger values of $\mu$, the absolute value of the combination of polylogarithmic functions in Eq. (16) grows faster than the decrease of the temperature and thus the high-temperature approximation is no longer valid.

Also, notice that the quark contribution, and thus the $\mu$ dependence of the thermal boson masses, enters the calculation, according to Eq. (16), as the argument of an exponential function and is controlled by $g^{2}$. This dependence translates into the effective potential and thus into the determination of the value of $\mu_{B}^{\text{CEP}}$. Thus, variations of $\mu$ into the exponent are amplified when varying $g$ and are, accordingly, more prominent than those of $T^{\text{CEP}}$, which, also according to Eq. (16), are mainly controlled by $\lambda$.

In order to study the fluctuations in the baryon number that the analysis predicts, we first look at the probability distribution. This is given by the expression

$\displaystyle{\mathcal{P}(v)}=\exp\left\{-\Omega V^{\text{eff}}(v)/T\right\},$

(23)

as a function of the order parameter around the equilibrium value in the symmetry restored phase $\langle v\rangle=0$. In Eq. (23), the factor $\Omega$ represents the system volume. For the purposes of this work, we explore large volumes compared to the typical size of the fireball crated in heavy-ion reactions, so as to mimic the thermodynamic limit. The normalized probability distribution is illustrated in Fig. 4 for different pairs of $\mu_{B}^{c}$, $T_{c}$ along the transition curve using $\lambda=1.4$, $g=0.88$ and $a=148.73$ MeV. We have extended the domain of the order parameter making $v\to|v|$ to include negative values so as to ensure that its average satisfies $\langle v\rangle=0$. Notice that for $\mu_{B}=0$ and $T_{c}^{0}=158$ MeV, for which the phase transition is second order, the probability distribution is Gaussian-like, albeit wider.

This happens because for these values of $\mu_{B}^{c}$ and $T_{c}$ the quartic terms in $v$ are important, producing a wider distribution around $\langle v\rangle$. As $\mu_{B}^{c}$ increases and the phase transition becomes first order for $\mu_{B}^{\text{CEP}}$ and $T^{\text{CEP}}$ the probability distribution becomes even wider. For $\mu_{B}^{c}>\mu_{B}^{\text{CEP}}$ and $T_{c}<T^{\text{CEP}}$, the phase transitions are always first order and the probability distribution develops secondary peaks that reflect the fact that for first order phase transitions the effective potential develops degenerate minima. We emphasize that the features of the probability distributions for first order phase transitions are due to the inclusion of the ring diagrams and thus of the plasma screening. Were these effects not to be included, the development of secondary peaks in the probability distribution would not happen and thus deviations from the Skellam statistics would not be possible.

To study the properties of the probability distribution we compute the behavior of its cumulants of third and fourth order. This is done in Fig. 5. The upper panel shows the skewness $s$ whereas the lower panel shows the kurtosis $\kappa$, as functions of $\mu_{B}$, for fixed $T$ across the corresponding critical value of $\mu_{B}^{c}$, for $\lambda=0.4$, $g=0.82$ and $a=141.38$ MeV. Notice that for second order phase transitions, both $s$ and $\kappa$ are represented by smooth curves. However, when $T$ approaches $T^{\text{CEP}}$, these functions show a peaked structure that becomes more pronounced when the transitions become first order. Notice also that when $T$ goes below $T^{\text{CEP}}$, both cumulants develop a maximum for $\mu_{B}>\mu_{B}^{c}$.

We now turn to describing the behavior of these cumulants as functions of the collision energy in heavy-ion reactions. For this purpose, we resort to the relation between the chemical freeze-out value of $\mu_{B}$ and the collision energy $\sqrt{s_{NN}}$, given by Cleymans et al. (2006)

$\displaystyle\mu_{B}(\sqrt{s_{NN}})=\frac{d}{1+e\sqrt{s_{NN}}},$

(24)

where $d=1.308$ GeV and $e=0.273$ GeV^-1. Figure 6 shows the ratio of the cumulants $C_{4}/C_{2}=\kappa\sigma^{2}$, normalized to the same ratio computed for $\mu_{B}=0$ and $T=T_{c}$, where $\sigma^{2}$ is the variance, for three different values of the volume $\Omega$. The upper panel is computed with the set of parameters $a=148.73$ MeV, $\lambda=1.4$ and $g=0.88$, whereas the lower panel corresponds to $a=141.38$ MeV, $\lambda=0.4$ and $g=0.82$. The value of $\sqrt{s_{NN}}$ for each set of parameters that corresponds to the CEP location is represented by the vertical line. Thus, we see that the CEP position is heralded not by the dip of $C_{4}/C_{2}$ but for its strong rise as the energy that corresponds to the CEP is approached.A similar result has been found in Ref. Mroczek et al. (2021). Also, as pointed out in Ref. Luo (2017), this non-monotonic behavior cannot be described by many other model calculations Xu et al. (2016); He et al. (2016). We emphasize that the reason for this failure is that other models do not incorporate long-range correlations and thus the possibility to describe critical phenomena, as we do in this work with the LSMq. Notice that in each case the ratio $C_{4}/C_{2}$ is independent of $\Omega$, as is expected for a ratio of cumulants, except near the collision energy where we find the CEP, where the high temperature approximation is less accurate. Nevertheless, the variations are of order of a few percent. Notice also that the product $\kappa\sigma^{2}$ significantly drops down for energies lower than the collision energy where the CEP is located. This result goes along recent HADES measurements of net-proton fluctuations at low energies Adamczewski-Musch, J. and others (2020).

Figure 7 shows the ratio $C_{4}/C_{2}=\kappa\sigma^{2}$, normalized to the same ratio computed for $\mu_{B}=0$ and $T=T_{c}^{0}$, for $\Omega=(200\ \text{fm})^{3}$ and for three different allowed sets of parameters $a$, $\lambda$ and $g$. Notice that although the dips have different depths, the strong rise of the ratio happens for almost the same value of the collision energy $\sqrt{s_{NN}}\sim 2$ GeV.

In conclusion, we have shown that the LSMq, with the effective potential computed up to the ring diagrams contribution and in the high temperature and chiral limits, produces an effective description of the system’s equilibrium distribution that deviates from that of the HRGM, where the ratios of cumulants of even order are always equal to 1. The plasma screening properties, encoded in the contribution of the ring diagrams, make this deviation possible, since they describe a key feature of plasmas near phase transitions, namely, long-range correlations. We have fixed the LSMq parameters using conditions at the phase transition for $\mu_{B}=0$ provided by LQCD calculations, namely the crossover transition temperature $T_{c}^{0}$ and the curvature parameters $\kappa_{2}$ and $\kappa_{4}$. The phase diagram can be directly studied by finding the kind of phase transitions that the effective potential allows when varying $T$ and $\mu_{B}$. We found that the CEP can be located in the range 786 MeV $<\mu_{B}^{\text{CEP}}<849$ MeV and 69 MeV $<T^{\text{CEP}}<70.3$ MeV. From the probability distribution obtained using the effective potential, we have computed the behaviour of the kurtosis and of the skewness and found that these cumulants show strong peaks as the CEP is crossed. We have also computed the ratio of the cumulants $C_{4}/C_{2}=\kappa\sigma^{2}$ as a function of the collision energy in a heavy-ion collision. We found that the location of the CEP coincides with the strong rise observed in this ratio and that this happens at $\sqrt{s_{NN}}\sim 2$ GeV.

Our calculations have been carried out in the high temperature limit and without including an explicit symmetry breaking term that gives rise to a finite vacuum pion mass. The high temperature approximation is valid near the (pseudo)transition lines, where what matters are the particles thermal masses. Since symmetry is being restored, the boson masses become small, of order $(\sqrt{\lambda/2+g^{2}})T$, at small values of $\mu_{B}$, whereas for fermions they vanish given that at symmetry restoration $v=0$. Therefore, a high temperature expansion is a reasonable, although of limited accuracy, description scheme for temperatures $T\sim(\sqrt{\lambda/2+g^{2}})T_{c}^{0}\sim T_{c}^{0}$ and only down to the region where the contribution from $\mu_{B}$ produces that the thermal boson masses are not that small compared to $T$.

The findings of this work are encouraging and can be extended to provide a more accurate description including an explicit chiral symmetry breaking introduced by a finite pion mass, as well as by going beyond the high temperature approximation. This is work for the near future that will be reported elsewhere.

Support for this work was received in part by UNAM-DGAPA-PAPIIT grant number IG100322 and by Consejo Nacional de Ciencia y Tecnología grant numbers A1-S-7655 and A1-S-16215. S. H.-O. acknowledges support from the U.S. DOE under Grant No. DE-FG02-00ER41132 and the Simons Foundation under the Multifarious Minds Program Grant No. 557037.