Plasma screening and the critical end point in the QCD phase diagram

In recent years, the study of hadron matter under extreme conditions of temperature and baryon density has become a subject of great interest. Of particular importance is the possibility of experimentally exploring the QCD phase structure by means of relativistic heavy-ion collisions. Currently, this exploration is carried out in experimental facilities such as RHIC, with the STAR Beam Energy Scan program, and HADES Adam, J. and others (2021); Adamczewski-Musch, J. and others (2020). Dedicated experiments soon to come on line, such as the NICA-MPD Kekelidze et al. (2017); Abgaryan et al. (2022) and FAIR-CBM Senger, P. (2017), are expected to widen the energy range for the exploration of the phase diagram.

From Lattice QCD (LQCD) calculations, its known 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 occurs at a pseudocritical temperature $T_{c}(\mu_{B}=0)\simeq 158$ MeV Borsanyi et al. (2020); Bazavov et al. (2019); Aoki et al. (2006). Also from the calculations made using effective models 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 found that this transition becomes first order at low $T$ and high $\mu_{B}$. Therefore, the point at which the first-order phase transition line in the $T$ vs $\mu_{B}$ plane ends and the crossover begins, as the temperature increases, is called the 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) but 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); Borsányi et al. (2021).

The Hadron Resonance Gas Model (HRGM) describes the crossover transition line for low values of $\mu_{B}$ found by LQCD Karsch (2017) when the occupation numbers are given in terms of Boltzman statistics Braun-Munzinger et al. (2003). Therefore, the strategy to locate the CEP consists in finding deviations from the statistical behavior of the HRGM. The statistical properties of a thermal system are characterized in terms of the cumulants of its conserved charges, that are extensive quantities Asakawa and Kitazawa (2016). To avoid uncertainties introduced by volume effects, the analyses involve ratios of these cumulants. The strategy narrows down to find deviations in the ratios of these cumulants from those obtained in HRGM, which are described by the Skellam distribution where the ratios of cumulants of even order are equal to 1. The baryon number is a conserved quantity that can be experimentally probed by means of measurements of proton multiplicities Adam et al. (2020); Abdallah et al. (2021a, b); Adamczewski-Musch, J. and others (2020). Therefore the location of the CEP can be identified by the appearance of critical behavior Hatta and Ikeda (2003); Stephanov et al. (1999); Bzdak et al. (2017, 2020); Athanasiou et al. (2010); Mroczek et al. (2021) in this and other conserved charges such as electric charge and strangeness, when a collision energy scan is performed.

In this work we study baryon number fluctuations as a function of the collision energy looking at the evolution of cumulant ratios, in the context of relativistic heavy ion collisions, using the Linear Sigma Model with quarks (LSMq) including the plasma screening effects. The latter are encoded in the ring diagram contribution to 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. To provide analytical insight, we employ the high temperature approximation and work in the chiral limit. Although these approximations have limitations in terms of the accuracy for the CEP localization, they are a useful guide for future more precise studies. The remaining of this 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. The parameters of the model are fixed by requiring that the phase transition line in the vicinity of $\mu_{B}=0$ corresponds to the one found by the most recent LQCD calculations Borsanyi et al. (2020). In Sec. III, we formulate how baryon number fluctuations are described in terms of the probability distribution associated with the order parameter near the transition line and report the results obtained from the analysis of the cumulant ratios used to describe baryon number fluctuations as a function of $\sqrt{s_{NN}}$. We show 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 in Sec. IV. The complete analysis can be found in Ref. Ayala et al. (2022).

The QCD phase diagram can be partially described by effective models; they can be used to explore different regions in parameter space depending on the degrees of freedom in the model. Given that LQCD calculations find that coincident deconfinement and chiral symmetry restoration transitions lines, it should be possible to explore the phase diagram emphasizing the chiral aspects of the transition only, like LSMq does. The Lagrangian of the LSMq 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})$
			$\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 a $SU(2)$ isospin doublet of quarks,

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

(2)

and $\sigma$ are isospin triplet and singlet, corresponding to the three pions and the sigma meson, 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. In order to allow for a spontaneous symmetry breaking, we work in the strict chiral limit and we let the $\sigma$ field to develop a vacuum expectation value $v$, namely,

$\displaystyle\sigma\rightarrow\sigma+v,$

(3)

which can later be taken as the order parameter of the transition. 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}+\frac{a^{2}}{2}v^{2}-\frac{\lambda}{4}v^{4}$
		$\displaystyle+$	$\displaystyle 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=-\frac{\lambda}{4}\sigma^{4}-\lambda v\sigma^{3}-\lambda v^{3}\sigma-\lambda\sigma^{2}\pi^{+}\pi^{-}-2\lambda v\sigma\pi^{+}\pi^{-}$
		$\displaystyle-\frac{\lambda}{2}\sigma^{2}(\pi^{0})^{2}-\lambda v\sigma(\pi^{0})^{2}-\lambda(\pi^{+})^{2}(\pi^{-})^{2}$
		$\displaystyle-\lambda\pi^{+}\pi^{-}(\pi^{0})^{2}-\frac{\lambda}{4}(\pi^{0})^{4}+a^{2}v\sigma$
	$\displaystyle{\mathcal{L}}_{I}^{f}$	$\displaystyle=-g\bar{\psi}(\sigma+i\gamma_{5}\vec{\tau}\cdot\vec{\pi})\psi.$		(5)

The expressions in Eq. (5) describe the interactions among the fields $\sigma$, $\vec{\pi}$ and $\psi$, after symmetry breaking.

From Eq. (4) we can see that the sigma, the three pions and the quarks have masses given by

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

respectively. We study the behavior of the effective potential in order to analyze the chiral symmetry restoration conditions in terms of temperature and quark chemical potential. The effective potential 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$ then chiral symmetry is spontaneously broken. To include quantum corrections at finite temperature and 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},$$

(9)

where

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

(10)

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}],$$

(11)

where

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

(12)

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})],$
			(13)

where $\Pi_{\text{b}}$ is the boson self-energy. The self-energies for the sigma and pion fields are in general different. However, we are working in the high-temperature approximation, keeping only the leading matter effects. In this approximation the boson self-energies become mass independent and therefore independent of the boson species. Additionally, because no external agent that might couple to the electric charge, such as a magnetic field, is considered, there is no distinction between neutral and charged pions and we can write Ayala et al. (2021)

	$\displaystyle\Pi_{\text{b}}$	$\displaystyle\equiv$	$\displaystyle\Pi_{\sigma}=\Pi_{\pi^{\pm}}=\Pi_{\pi^{0}}$
		$\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, while ${\rm Li}_{2}(x)$ stands for the polylogarithm function of order 2.

The matter correction to the boson mass is encoded in the boson self-energy. For a second order (our approach to a crossover due the strict chiral limit) these corrections should cause the thermal boson masses to vanish when the symmetry is restored. This means that in the 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$. So, 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.$		(16)
		$\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. (LABEL:Pileading),

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

(17)

We take $\tilde{\mu}=500$ MeV, which is large enough to be considered the largest energy scale in our problem. Notice that the dependence in $\tilde{\mu}$ is only logarithmic, therefore small variations on this parameter do not change significantly the result. To fix $\lambda$, $g$ and $a$, we look for a set of parameters such that the solutions of Eq. (16) produce a curve comparable to the LQCD transition curve at $T_{c}^{0}\simeq 158$ MeV. In order to compare the curves, we use a common parameterization 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},$$

(18)

with $\kappa_{2}=0.0153$ and $\kappa_{4}=0.00032$ Borsanyi et al. (2020); Guenther (2021). In other words, fixing the parameters is reduced to finding the set of $\lambda$, $g$ and $a$ such that the coefficients $\kappa_{2}$ and $\kappa_{4}$ of their transition curve are comparable to those given by LQCD. Explicitly, we first fix a value of $\lambda$, and find the solution of Eq. (16) for $g$ and $a$ that, from the effective potential in Eq, (15), 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. (18). 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},$$

(19)

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}.$

(20)

Since Eq. (16) is non-linear, the set of parameters is not unique.

With all the parameters fixed, we can study the properties of the effective potential and find the values of $T_{c}$ and $\mu_{B}^{c}$ where the chiral phase transition takes place. Figure 1 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. This is shown more clearly in Fig. 2, where we show the phase diagrams 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. These are the results directly obtained from our analysis. In all cases, we locate the CEP in a region of low temperatures and high quark chemical potential. 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. Now that we have the complete analysis of the phase diagram from the effective potential, in the next section we introduce the elements necessary to study the phase diagram in terms of baryon number fluctuations.

To study the fluctuations in the number of baryons predicted by the analysis, we start by looking at the probability distribution, which is a function of the order parameter around the equilibrium value in the restored symmetry phase $\langle v\rangle=0$ given by

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

(21)

where, the factor $\Omega$ represents the system volume. For this work, we explored large volumes compared to the typical fireball size created in heavy ion collisions, to simulate the thermodynamic limit. Using $\lambda=1.4$, $g=0.88$ and $a=148.73$ MeV, we illustrate in Fig. 3 the normalized probability distribution for different pairs of $\mu_{B}^{c}$, $T_{c}$ along the transition curve. We have extended the domain of the order parameter making $v\to|v|$ to include negative values and ensure that its average satisfies $\langle v\rangle=0$. Notice that for $\mu_{B}=0$ and $T_{c}^{0}$, for which the phase transition is second order, the probability distribution is Gaussian-like, albeit wider. This is due to the fact that quartic terms in $v$ are important for these values of $\mu_{B}^{c}$ and $T_{c}$ producing a wider distribution around $\langle v\rangle$. As $\mu_{B}^{c}$ increases and the phase transition becomes first order 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. If these effects were not included, the development of secondary peaks in the probability distribution would not occur and, therefore, deviations from the Skellam statistic would not be possible. To understand the properties of the probability distribution we calculate the behavior of the fourth order cumulant, also known as kurtosis. Figure 4 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, $\kappa$ is represented by smooth curves. However, when $T$ approaches $T^{\text{CEP}}$, this function shows a peaked structure that becomes more pronounced when the transitions become first order. Notice also that when $T$ goes below $T^{\text{CEP}}$, the kurtosis develops a maximum for $\mu_{B}>\mu_{B}^{c}$. Let us next proceed to describe the behavior of the kurtosis as functions of the collision energy in heavy-ion reactions. To this end, 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}}},$

(22)

where $d=1.308$ GeV and $e=0.273$ GeV^-1. Figure 5 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$. Additionally 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 include long-range correlations, as we do in this work with the LSMq, which is crucial to describe critical phenomena. Notice that, as is expected for a ratio of cumulants, in each case the ratio $C_{4}/C_{2}$ is independent of the volume $\Omega$ except near the collision energy where we find the CEP, where the high temperature approximation is less accurate.

Notice also that the product $\kappa\sigma^{2}$ significantly drops down for energies lower than the collision energy where the CEP is located, in agreement with recent HADES measurements of net-proton fluctuations at low energies Adamczewski-Musch, J. and others (2020). Finally, Fig. 6 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$. It is worth mentioning that, even though the dips have different depths, the sharp increase in the ratios occurs for almost the same value of collision energy $\sqrt{s_{NN}}\sim 2$ GeV.

In this work we have used the LSMq in the high temperature and chiral limits to explore the QCD phase diagram from the point of view of chiral symmetry restoration. We have computed the finite temperature effective potential up to the contribution of the ring diagrams to account for the plasma screening effects. This model makes an effective description of the system equilibrium distribution that deviates from that of the HRGM, where the ratios of cumulants of even order are always equal to 1. When we include the plasma screening properties, encoded in the ring diagrams contribution, we find a deviation from the HGR model, since the screening properties describe a key feature of plasma in the transitions, namely, the long-range correlations. We fix 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 obtained 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 found that these cumulants show strong peaks as the CEP is crossed. Finally, we describe 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 conclude that the CEP location coincides with a sharp rise in this ratio at $\sqrt{s_{NN}}\sim 2$ GeV. Our studies were performed at the high temperature limit and without including an explicit symmetry breaking term that gives rise to a finite vacuum pion mass. Due to this, the encouraging findings of this work 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 relaxing 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. JJCM acknowledges financial support from the University of Sonora under grant USO315007861. 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. METY is grateful for the hospitality of Perimeter Institute where part of this work was carried out and this research was also supported in part by the Simons Foundation through the Simons Foundation Emmy Noether Fellows Program at Perimeter Institute. Research at Perimeter Institute is supported in part by the Government of Canada and by the Province of Ontario.