Chiral restoration temperature at finite spin density in QCD

TL;DR

This study explores how finite spin density influences the chiral phase transition temperature in QCD, revealing a transition from crossover to first-order and identifying a critical end-point through an extended phase diagram.

Contribution

It introduces a method to extend the QCD phase diagram at finite spin density using the linear sigma model and lattice data comparison, addressing regularization ambiguities.

Findings

Chiral transition temperature decreases with increasing spin potential.

Transition changes from crossover to first-order at a critical end-point.

Chiral symmetry restoration occurs at high spin potentials regardless of temperature.

Abstract

We investigate the impact of a uniform spin density on the critical temperature of the chiral phase transition in finite-temperature QCD in the scope of the linear sigma model. We demonstrate that at a finite spin potential $\mu_\Sigma$, corresponding to a finite spin density, the predictive power of the model is challenged by an ambiguity associated with a contribution of the vacuum renormalization term to the free energy. Eliminating the regularization freedom through comparison with recent low-$\mu_{\Sigma}$ lattice data, we extend the phase diagram of QCD at finite spin density to regions inaccessible to lattice simulations. We show that, as the spin potential increases, the temperature of the chiral crossover transition diminishes and the chiral crossover turns into a first-order transition at a second-order critical end-point $(T,\mu_\Sigma)_\mathrm{CEP}\simeq (0.142,0.098)$ GeV. With increasing spin potential, the critical temperature touches the zero-temperature axis at $\mu_{\Sigma} = 0.310$ GeV, implying that the chiral symmetry is restored at higher potentials at any temperature.

Equations

$$\mathcal{L}={\mathcal{L}}_{\mathcal{M}}+{\mathcal{L}}_q,$$

$${\mathcal{L}}_{\mathcal{M}}=\frac{1}{2}\left({\partial }_{\mu }\sigma {\partial }^{\mu }\sigma +{\partial }_{\mu }\vec{\pi }\cdot {\partial }^{\mu }\vec{\pi }\right)-V_{\mathcal{M}},$$

$$V_{\mathcal{M}}=\frac{\lambda }{4}({\sigma }^2+{\vec{\pi }}^2-v^2{)}^2-h\sigma ,$$

$$V_{\mathcal{M}}=V_0+\frac{1}{2}{m}_{\sigma }^2\delta {\sigma }^2+\frac{1}{2}{m}_{\pi }^2\delta {\vec{\pi }}^2+\cdots ,$$

$$m_{\sigma }^2=\lambda (3⟨\sigma {⟩}^2-v^2),m_{\pi }^2=⟨\sigma {⟩}^2-v^2.$$

$${\mathcal{L}}_q=\overline{\psi }\left[i/\partial -g(\sigma +i{\gamma }_5\vec{\tau }\cdot \vec{\pi })\right]\psi +{\delta }_{\Sigma }{\mathcal{L}}_q.$$

$${\delta }_{\Sigma }{\mathcal{L}}_q={\mu }_{\alpha ,\mu \nu }{\mathcal{S}}^{\alpha ,\mu \nu },$$

$$\displaystyle{\cal S}^{\alpha,\mu\nu}=\frac{1}{2}\bar{\psi}\bigl{\{}\gamma^{\alpha},\Sigma^{\mu\nu}\bigr{\}}\psi\,,\qquad\Sigma^{\mu\nu}=\frac{i}{4}\bigl{[}\gamma^{\mu},\gamma^{\nu}\bigr{]}\,.$$

$${\delta }_{\Sigma }{\mathcal{L}}_q={\mu }_{\Sigma }{\psi }^{†}({\mathbf{n}}_{\Sigma }\cdot \mathbf{\Sigma })\psi ,$$

$${\mathcal{H}}_q=-i{\gamma }^0\mathbf{\gamma }\cdot \nabla +{\gamma }^{0}g\sigma -{\mu }_{\Sigma }({\mathbf{n}}_{\Sigma }\cdot \mathbf{\Sigma }).$$

$${\mathcal{H}}_q(\mathbf{p}){\psi }_{\varsigma }^{(s)}(\mathbf{p})=\varsigma E_{\mathbf{p}}^{(s)}{\psi }_{\varsigma }^{(s)}(\mathbf{p}),$$

$$(\varsigma E_{\mathbf{p}}^{(s)}-g\sigma -\frac{1}{2}{\mu }_{\Sigma }{\mathbf{n}}_{\Sigma }\cdot \mathbf{\sigma }){\psi }_u$$

$$(\varsigma E_{\mathbf{p}}^{(s)}+g\sigma -\frac{1}{2}{\mu }_{\Sigma }{\mathbf{n}}_{\Sigma }\cdot \mathbf{\sigma }){\psi }_d$$

$$(\mathbf{A}\cdot \mathbf{\sigma }-B){\psi }_u=0,$$

$$\mathbf{A}$$

$$B$$

$${\epsilon }^4={\mu }_{\Sigma }^2[(\mathbf{p}\cdot {\mathbf{n}}_{\Sigma }{)}^2+g^2{\sigma }^2].$$

$$[E_{\mathbf{p}}^{(s)}{]}^2={\mathbf{p}}^2+g^2{\sigma }^2+s^2{\mu }_{\Sigma }^2-2s{\mu }_{\Sigma }\sqrt{(\mathbf{p}\cdot {\mathbf{n}}_{\Sigma }{)}^2+g^2{\sigma }^2}.$$

$$E_{\mathbf{p}}^{(s)}$$

$$=\sqrt{(\mathbf{p}\times {\mathbf{n}}_{\Sigma }{)}^2+[\sqrt{g^2{\sigma }^2+(\mathbf{p}\cdot {\mathbf{n}}_{\Sigma }{)}^2}-s{\mu }_{\Sigma }{]}^2}.$$

$$E_{\mathbf{p}}^{(s)}$$

$$=\sqrt{E_p^2+s^2{\mu }_{\Sigma }^2-2s{\mu }_{\Sigma }{E}_{p_z}},$$

$${\mathcal{Z}}_q=\int [id{\psi }^{†}][d\psi ]e^{-S_E},S_E=\underset{0}{\overset{\beta }{\int }}d\tau \int d^{3}x{\mathcal{L}}_E,$$

$$\psi (X)=\frac{1}{\sqrt{{\mathcal{V}}_{3d}}}\sum _n\sum _{s,\varsigma }\sum _{\mathbf{p}}{e}^{i({\omega }_n\tau +\mathbf{p}\cdot \mathbf{x})}{\psi }_{n,\varsigma }^{(s)}(\mathbf{p}){\hat{b}}_{n,\varsigma }^{(s)}(\mathbf{p}),$$

$$[{\hat{\psi }}_{n,\varsigma }^{(s)}(\mathbf{p}){]}^{†}{\hat{\psi }}_{n,{\varsigma }^{\prime }}^{(s^{\prime })}(\mathbf{p})={\delta }_{s{s}^{\prime }}{\delta }_{\varsigma {\varsigma }^{\prime }}.$$

$$S_E=i\beta \sum _n\sum _{s,\varsigma }\sum _{\mathbf{p}}({\omega }_n-i\varsigma E_{\mathbf{p}}^{(s)})[{\hat{b}}_{n,\varsigma }^{(s)}{]}^{†}(\mathbf{p}){\hat{b}}_{n,\varsigma }^{(s)}(\mathbf{p}),$$

$${\mathcal{Z}}_q={\left[\prod _{n,\mathbf{p}}\prod _{s,\varsigma }i\beta ({\omega }_n+i\varsigma E_{\mathbf{p}}^{(s)})\right]}^{N_f{N}_c},$$

$$\ln {\mathcal{Z}}_q=N_f{N}_c\sum _{n,\mathbf{p}}\sum _{s,\varsigma }\ln [\beta (-i{\omega }_n-\varsigma E_{\mathbf{p}}^{(s)})].$$

$$F_q^{\mathrm{ZP}}$$

$$F_q$$

$$V_{\mathcal{M}}\to V_{\sigma }=\frac{\lambda }{4}({\sigma }^2-v^2{)}^2-h\sigma .$$

$$F_q^{\mathrm{ZP}}$$

$$I_1=-\frac{g^4{\sigma }^4}{32{\pi }^2}\left[\frac{1}{ϵ}+\frac{3}{2}-{\gamma }_E+\ln \left(\frac{4\pi M^2}{g^2{\sigma }^2}\right)\right],$$

Full text

I Introduction

Spin polarization phenomena in the quark–gluon plasma (QGP) have recently attracted significant attention Becattini (2022). Observations of global and local spin polarizations in non-central heavy-ion collisions Adamczyk et al. (2017) reveal the transfer of angular momentum from colliding nuclei to hadronic matter and provide sensitive probes of the thermodynamic, topological, and electromagnetic properties of QCD matter Huang et al. (2021); Fukushima (2019); Becattini et al. (2024).

Measurements of the global polarization of $\Lambda$ hyperons and of local spin polarization—including the spin alignment of various hadrons—suggest that, at a certain stage of its evolution, the vortical QGP develops a finite spin density of quarks. A precise analytical formulation of spin degrees of freedom, however, remains subtle due to the so-called “pseudogauge symmetry ambiguity” Becattini et al. (2019); Speranza and Weickgenannt (2021); Fukushima and Pu (2021), which reflects the non-uniqueness of a local definition of the spin operator. This ambiguity can potentially be resolved within interacting field theories, as recently demonstrated in Ref. Buzzegoli and Palermo (2024), through the canonical spin formulation of quarks in an effective infrared model of QCD.

Based on this idea, in our paper we explore the effect of finite spin quark density on the QGP thermodynamics and chiral transition properties within the Linear sigma model coupled to quarks Gell-Mann and Lévy (1960), LSMq, using the canonical formulation of spin operator. The spin density is described by the spin potential, $\mu_{\Sigma}$, which plays the role of a thermodynamic quantity conjugated to the spin density and also enters the hydrodynamical description of the system Montenegro et al. (2017); Florkowski et al. (2020); Gallegos et al. (2021); Ambrus et al. (2022); Abboud et al. (2025); Becattini and Singh (2025). A similar incorporation of spin density has been recently employed in first principle lattice simulation with two dynamical light quarks Braguta et al. (2025) where it was shown that a finite spin density leads to a decrease of both chiral and deconfinement temperatures in the limit of the physical quark masses.

It is worth mentioning that the spin polarization, being a part of the total angular momentum of the system, is not a separately conserved quantity because the angular momentum can be shared between the orbital and spin contributions. Therefore, the spin potential should be distinguished from genuine chemical potentials associated, for example, with the conserved electric, baryonic or isospin charges. Still, the lattice simulations of the phase diagram of QCD at a finite spin potential have demonstrated the excellent continuum scaling of the physical quantities Braguta et al. (2025). The absence of infrared and ultraviolet artifacts in these simulations suggests that the implementation of the spin potential is free from pathologies that could be associated, for example, with the spin non-conservation and the pseudogauge freedom. Therefore, the first-principle lattice simulations suggest that the canonical spin potential for quark degrees of freedom has a physical sense in QCD.

However, similar to baryon chemical potential ($\mu_{B}$), due to the sign problem in the implementation on lattice of real $\mu_{\Sigma}$, the lattice results reported in Ref. Braguta et al. (2025) were limited to low spin potentials, effectively restricting the exploration of regimes with substantial quark spin polarization. In our paper, we study the influence of the spin density on the thermodynamics of QCD within the effective model framework which, theoretically, has no such limitation. Thus, the aim of the current project is twofold: (i) to match the predictions of the model for the chiral transition line to the lattice data at small values of the spin potential; and (ii) to explore the region of arbitrary values of the spin density, which is inaccessible to lattice simulations.

One advantage of using the LSMq model over other effective QCD models is that it is renormalizable. In this work, we have explored the effect of vacuum renormalisation on the model thermodynamics. We have observed that the incorporation of the $\mu_{\Sigma}$-dependent finite counter-terms in the vacuum free energy is crucial to obtain a qualitative agreement with the lattice results of Ref. Braguta et al. (2025). Within that setup, we investigated the influence of spin polarization on the phase diagram of QCD and also studied the behavior of spin density as a function of the thermodynamic parameters $T$ and $\mu_{\Sigma}$.

We mention that the nonconservation of the canonical spin tensor has profound implications, as will be explicitly demonstrated through our analysis. First of all, it does not commute with the Dirac Hamiltonian, therefore giving rise to a nonlinear dispersion relation that leads to a $\mu_{\Sigma}$-dependent zero-point free energy. Second, unlike the baryon chemical potential, the spin potential has no Fermi level. Third, at large temperatures, the thermal part of the quark free energy becomes independent of $\mu_{\Sigma}$ and the spin expectation value is due solely to the renormalized meson potential. None of the above (unnatural) features are present in the case when the polarization is induced via the conserved helicity operator, considered in our previous work Chernodub and Ambrus (2021). It would be interesting to test this hypothesis with other conserved spin operators, such as the Pryce operator introduced in Ref. Pryce (1948) and revisited recently in Refs. Bauke et al. (2014); Cotaescu (2022).

The structure of this paper is as follows. Section II introduces the details of the model. Section II.1 is devoted to the mesonic sector, while Section II.2 addresses the quark sector, where the finite spin density term is incorporated into the quark Lagrangian and the corresponding dispersion relation is derived. The thermodynamics of the model is discussed in Section II.3. In Section III, we present a detailed analysis of vacuum renormalization, which plays a central role in achieving a consistent description of QCD thermodynamics with finite spin density within the LSMq framework. Section IV highlights the subtleties involved in solving the saddle-point equations with the modified dispersion relation in the presence of a spin potential. The calculation of the phase diagram and determination of the ambiguity of the renormalization procedure with the help of the lattice data are presented in Section V. Our conclusions are summarized along with a brief discussion in Section VI. Three appendices are devoted to the discussion of particularities of the dimensional regularization (Appendix A), the infrared divergences (Appendix B) and the small-mass limit (Appendix C). Throughout this work, we adopt the Minkowski metric $g={\rm diag}(1,-1,-1,-1)$.

II Model

The linear sigma model coupled to quarks Gell-Mann and Lévy (1960) describes interacting mesonic and quark degrees of freedom. For a 2-flavor system considered in this work, the quark sector is given by the isospin doublet $\psi=(u,d)^{T}$, where $u$ and $d$ denote the up and down quark spinors, respectively. The mesonic fields include the $\sigma$ and $\vec{\pi}$ mesons, where $\vec{\pi}=(\pi^{+},\pi^{-},\pi^{0})$ is the pion scalar isotriplet. The model Lagrangian can be written as follows:

[TABLE]

where $\mathcal{L}_{\mathcal{M}}$ and $\mathcal{L}_{q}$ correspond to the mesonic and the quark parts of the Lagrangian, respectively, which we will briefly discuss in the following.

II.1 Mesonic sector

The mesonic sector incorporates the chiral properties of the model. In the QCD with finite quark masses, the chiral symmetry is explicitly broken, and the chiral phase transition manifests as a smooth crossover rather than a true phase transition. Consequently, the field $\sigma$ plays the role of the approximate order parameter associated with the chiral symmetry breaking: it takes a non-vanishing value in the low-temperature, chirally-broken phase and gradually approaches zero as the temperature increases.

The mesonic part of the Lagrangian,

[TABLE]

incorporates, besides the kinetic terms for the meson fields, the phenomenological interaction potential:

[TABLE]

with $\lambda$, $v$ and $h$ being model parameters.

The breaking of the chiral symmetry in the vacuum is implemented at the level of the $\sigma$ meson, which acquires a non-vanishing vacuum expectation value. In contrast, the pion field retains a vanishing vacuum expectation value, $\langle{\vec{\pi}}\rangle=0$, as it can always be removed by a global chiral rotation of the pion and sigma fields (see, for example, Ref. Schechter and Ueda (1971)). The first term in the potential (3) leads to a strong spontaneous breaking of the chiral symmetry, supplemented by weak explicit symmetry breaking due to the second term.

The potential (3) can be expanded over the small fluctuations $\delta\sigma$ and $\delta\pi$ of the meson fields as:

[TABLE]

which allows the meson masses to be computed as:

[TABLE]

The masses of the charged and the neutral pions do not split due to the absence of electromagnetic interactions.

The parameter $h$ ensuring the explicit symmetry breaking is obtained by enforcing the PCAC relation $h=f_{\pi}m_{\pi}^{2}$, where $f_{\pi}=93$ MeV is the pion decay constant and $m_{\pi}=138$ MeV is the pion mass in the vacuum. Demanding $\langle\sigma\rangle=f_{\pi}$ Gell-Mann and Lévy (1960) and $m_{\sigma}=600$ MeV, the remaining model parameters $v$ and $\lambda$ are fixed through the relations $v^{2}=f_{\pi}^{2}-m_{\pi}^{2}/\lambda$ and $m_{\sigma}^{2}=2\lambda f_{\pi}^{2}+m_{\pi}^{2}$.

II.2 Quark sector

We extend the usual quark sector of the linear sigma model coupled to quarks (LSMq) which includes both the quark-meson interaction and the kinetic term for quark fields, to account for a non-vanishing spin density:

[TABLE]

The interaction term involving the $\sigma$ and $\vec{\pi}$ fields enables dynamical mass generation for the quarks through spontaneous chiral symmetry breaking, with $\langle\sigma\rangle\neq 0$ leading to a constituent quark mass $m_{q}=g\langle\sigma\rangle$. Taking $m_{q}=307$ MeV yields a Yukawa coupling $g\simeq 3.3$.

We express the additional term $\delta_{\Sigma}\mathcal{L}_{q}$ as Braguta et al. (2025)

[TABLE]

where the canonical spin part of the angular momentum density for spin-1/2 Dirac fermions is given by

[TABLE]

The associated background relativistic spin field $\mu_{\alpha,\mu\nu}$ introduces a finite spin density and a finite spin current in the system. For the study of polarization, we consider the case when only the rotation part of the spin potential is non-vanishing. Defining $\mu_{\alpha,\mu\nu}=-\frac{1}{2}g_{\alpha 0}\varepsilon_{0\mu\nu k}\mu_{\Sigma}^{k}$, with $\varepsilon^{0123}=-\varepsilon_{0123}=1$, we arrive at

[TABLE]

where we took $\mathbf{n}_{\Sigma}=\bm{\mu}_{\Sigma}/\mu_{\Sigma}$ to be the direction of the vector spin potential, with $\mu_{\Sigma}=\sqrt{\mu_{1}^{2}+\mu_{2}^{2}+\mu_{3}^{2}}$ being the three-norm of $\bm{\mu}_{\Sigma}$. We also used the convention $\Sigma^{i}=\frac{1}{2}\varepsilon^{0ijk}\Sigma_{jk}$, i.e. $\Sigma^{3}=\Sigma^{12}=\frac{i}{2}\gamma^{1}\gamma^{2}$, with the other components given by circular permutations. The spin term (7) in the Lagrangian (6) does not break explicitly the global chiral symmetry, given by the transformations $\psi\to e^{i\omega_{5}\gamma_{5}}\psi$ and ${\bar{\psi}}\to{\bar{\psi}}e^{i\omega_{5}\gamma_{5}}$ with a coordinate-independent, real-valued parameter $\omega_{5}$.

In the mean-field approximation, the meson fields are replaced by their (homogeneous) expectation values in the quark Lagrangian (6). Taking into account that $\langle{\vec{\pi}}\rangle=0$, the Dirac Hamiltonian operator reads as follows:

[TABLE]

To study the energy branches supported by the system, we construct the eigensystem of $\mathcal{H}_{q}$ at the level of Fourier modes $e^{i\mathbf{p}\cdot\mathbf{x}}\psi_{\varsigma}^{(s)}(\mathbf{p})$,

[TABLE]

with $\varsigma=\pm 1$ labeling particle ($+1$) and antiparticle ($-1$) modes, while $s=\pm 1/2$ denotes the polarization degree of freedom. Denoting $\psi_{\varsigma}^{(s)}=(\psi_{u},\psi_{d})^{T}$, we have

[TABLE]

Multiplying the top relation by $\bm{\sigma}\cdot\mathbf{p}$ and using $(\bm{\sigma}\cdot\mathbf{p})^{2}=\mathbf{p}^{2}$, $\psi_{d}$ can be eliminated in favor of $\psi_{u}$, leading to

[TABLE]

where we denoted

[TABLE]

with $\varepsilon^{2}=[E^{(s)}_{\mathbf{p}}]^{2}-\mathbf{p}^{2}-g^{2}\sigma^{2}-\mu_{\Sigma}^{2}/4$. We now multiply Eq. (14) by $\mathbf{A}\cdot\bm{\sigma}+B$. Using the relation $(\mathbf{A}\cdot\bm{\sigma}-B)(\mathbf{A}\cdot\bm{\sigma}+B)=\mathbf{A}^{2}-B^{2}$ for any C-numbers $\mathbf{A}$, $B$, non-trivial solutions require that $\mathbf{A}^{2}=B^{2}$. After some simplifications, we arrive at

[TABLE]

We employ the polarization parameter $s=\pm 1/2$ to keep track of the sign of the square root of the above expression, such that

[TABLE]

Taking another square root and demanding that $E^{(s)}_{\mathbf{p}}{>}0$, we arrive at

[TABLE]

Specializing the above result to the case when the spin polarization is directed along the $z$ axis, i.e. $\mathbf{n}_{\Sigma}=\mathbf{e}_{z}$, we obtain

[TABLE]

where $E_{p}=\sqrt{p^{2}+g^{2}\sigma^{2}}$ and $E_{p_{z}}=\sqrt{p_{z}^{2}+g^{2}\sigma^{2}}$ are, respectively, the total and longitudinal energies in the limit of vanishing spin potential. The very same spatially-anisotropic form of the dispersion relations also appears in condensed matter systems, such as in Weyl semimetals Burkov et al. (2011) or in Dirac semimetals subjected to rotating electric fields Ebihara et al. (2016). In addition, the double-square-root structure of the energy dispersion (19) emerges also in QCD in the presence of chiral Farias et al. (2016); Azeredo et al. (2024) and isospin Avancini et al. (2019); Lopes et al. (2021); Ayala et al. (2023) chemical potentials.

It is instructive to clarify the physical meaning of the spin potential $\mu_{\Sigma}$. This quantity enters the Lagrangian, Eqs. (6) and (9), as a Lagrange multiplier in a way formally-analogous to a conventional chemical potential. However, there is a fundamental difference between the two quantities. A standard chemical potential quantifies the energy required to add or remove a particle carrying a conserved charge. For instance, the baryonic chemical potential $\mu_{B}$ corresponds to the energy cost of inserting a unit of baryonic charge (e.g., a proton) or, equivalently, removing a unit of negative baryonic charge (e.g., an antiproton) from the ensemble.

By contrast, the mentioned analogy fails for the spin potential: first, spin is not a conserved quantum number, and second, the quasiparticle energy dispersion in Eqs. (18) or  (19) does not take the form of a simple shift of a momentum-dependent energy by a chemical potential. Thus, while $\mu_{\Sigma}$ plays a role formally reminiscent of a chemical potential, it cannot be interpreted in the same thermodynamic sense. Instead, it should be regarded as an external parameter that controls the degree of spin polarization in the system.

Nevertheless, there exists a region of parameter space in which $\mu_{\Sigma}$ may be interpreted, albeit in a rather loose sense, as playing a role analogous to that of a chemical potential. To see this point, let us consider the energy dispersion (18) for a quasiparticle with a momentum aligned with the direction of the spin, ${\bf p}=p\,{\bf e}_{\Sigma}\equiv p\,{\bf e}_{z}$. Then, in the notations of Eq. (19), $E_{\mathbf{p}}^{(s)}=|E_{p_{z}}-s\mu_{\Sigma}|$, implying that a relatively small spin potential, $|\mu_{\Sigma}|<2E_{p_{z}}$, acts like a chemical potential for spin-polarized fermions propagating along the direction of the spin polarization, $E_{\mathbf{p}}^{(s)}=E_{p_{z}}-s\mu_{\Sigma}$.

Consequently, the spin potential $\mu_{\Sigma}$ may be interpreted, similarly to a chemical potential, as the energy cost associated with adding (or removing) a unit of spin ($\hbar$ in dimensional units) provided two conditions are satisfied: (i) the single-particle energy exceeds $|\mu_{\Sigma}/2|$ and (ii) the momentum of the particle is aligned (or anti-aligned, respectively) with the direction of the spin polarization so that ${\bf p}\|{\bf n}_{\Sigma}$. Since fermions and anti-fermions carry a half of the unit of spin, $\hbar/2$, the spin potential adds or removes the energy $\mu_{\Sigma}/2$.

II.3 Thermodynamics of the spin system

Under the path integral formalism, we have the total partition function $\mathcal{Z}=\mathcal{Z}_{\mathcal{M}}\mathcal{Z}_{q}$, where $\mathcal{Z}_{\mathcal{M}}=e^{-\beta\mathcal{V}_{3d}V_{\mathcal{M}}}$ is the mesonic contribution, with $\mathcal{V}_{3d}$ being the space volume. The quark contribution $\mathcal{Z}_{q}$ is given as

[TABLE]

where the Euclidean Lagrangian density satisfies $\mathcal{L}_{E}=-\mathcal{L}_{q}=\psi^{\dagger}(\partial_{\tau}+\mathcal{H}_{q})\psi$. The functional integration can be performed using the normal mode decomposition of the fermion field,

[TABLE]

with $X_{\mu}=(\tau,\mathbf{x})$ denoting the Euclidean four-coordinates and $\omega_{n}=(2n+1)\pi T$ being the Matsubara frequencies, ensuring the anti-periodicity of the fermionic wavefunction with respect to the imaginary time $\tau=it$. The modes $\psi_{n,\varsigma}^{(s)}(\mathbf{p})$ correspond to the Hamiltonian eigenfunctions $\psi_{\varsigma}^{(s)}(\mathbf{p})$ introduced in Eq. (11), satisfying the orthogonality relation

[TABLE]

The Fourier coefficients $\hat{b}_{n,\varsigma}^{(s)}(\mathbf{p})$ are represented as $N_{f}\times N_{c}$ Grassmann variables.

After carrying out the functional integration, the Euclidean action $S_{E}$ defined in Eq. (20), can be expressed as

[TABLE]

where the energy branches $\varsigma E^{(s)}_{\mathbf{p}}$ correspond to the eigenvalues of the Hamiltonian introduced in Eq. (11).

The quark contribution to the partition function can then be obtained using the Gauss integration formula for Grassmann variables, $\int d\hat{b}^{\dagger}d\hat{b}\,e^{a\hat{b}^{\dagger}\hat{b}}=a$, yielding

[TABLE]

where the degeneracy factor $N_{f}N_{c}$ accounts for $N_{f}$ quark flavours and $N_{c}$ color degrees of freedom. Taking the logarithm of the above expression, we obtain

[TABLE]

The total thermodynamic potential of the spin system $\mathcal{F}$ subsequently has contributions from both the mesonic and quark sectors. The mesonic contribution can be obtained in a straightforward way as $-T\ln\mathcal{Z}_{\mathcal{M}}=\mathcal{V}_{3d}V_{\mathcal{M}}$. On the other hand, the quark contribution of the free energy $\mathcal{F}_{q}=-T\ln\mathcal{Z}_{q}=\mathcal{V}_{3d}(F^{\rm ZP}_{q}+F_{q})$ can be decomposed into its zero-point (ZP) and thermal parts:

[TABLE]

III Vanishing temperature: Zero-point renormalization

In the mean-field approximation of a homogeneous system, discussed in Sec. II.3, the thermodynamic potential density $F=\mathcal{F}/\mathcal{V}_{3d}$ receives a contribution from the mesonic potential (3) with vanishing pion field, $\vec{\pi}=0$:

[TABLE]

Keeping the field $\sigma$ as a background scalar field, the potential gets a divergent one-loop contribution from the quark sector, denoted by $F_{q}^{\rm ZP}$ in Eq. (26), which we treat via a self-consistent renormalization. This section outlines the procedure of renormalization in the presence of a finite spin density.

We begin with the standard LSMq model in Subsec. III.1 and subsequently incorporate the axial-vector interaction to isolate and address spin-potential-dependent divergences. In Subsections III.2 and III.3, we introduce two schemes: a naïve implementation in Subsec. III.2, leading to a spin-independent renormalized potential; and a spin-dependent deformation in Subsec. III.3, featuring a free parameter $\ell$, allowing the spin to modify the vanishing temperature limit of the theory.

III.1 Standard LSMq

In standard LSMq, the zero-point quark contribution to the free energy in Eq. (26) reads

[TABLE]

Employing dimensional regularization, as explained in Appendix A, $I_{1}$ evaluates to Skokov et al. (2010)

[TABLE]

where $\gamma_{E}$ is the Euler-Mascheroni constant and $M$ is the renormalization scale introduced in the minimal subtraction scheme.

In order to self-consistently account for this vacuum contribution coming from the fermion loop, we employ a bare mesonic potential of the form

[TABLE]

Compared to Eq. (28), the above expression displays a logarithmic term, introduced to balance the analogous term appearing in $I_{1}$ above. We now demand that, at zero temperature, the renormalized potential $V^{\rm ren}_{\sigma}=V_{\sigma}^{\rm bare}+F_{q}^{\rm ZP}$ reduces to the classical mesonic form given in Eq. (28):111Ref. Skokov et al. (2010) argues that the $\sigma^{4}\ln\sigma$ term is required to ensure that the chiral transition is a crossover in the chiral limit of the theory. With our choice of parameters, the transition is a crossover even without this term. Therefore, for simplicity, we do not keep it in the renormalized potential.

[TABLE]

where we subtract the value of the potential in the vacuum to ensure that $P=-(V_{\sigma}^{\rm ren}+F_{q})$ represents the mechanical pressure. It is straightforward to identify:

[TABLE]

where we took into account that, while $\lambda_{0}v_{0}^{2}=\lambda v^{2}$ is finite, $\lambda_{0}v_{0}^{4}=0$, because $v_{0}^{2}=O(\epsilon)$ vanishes when the limit $\epsilon\to 0$ is taken.

III.2 Spin-deformed LSMq: a naïve approach

Next, we turn our focus to the modifications due to the inclusion of the spin potential $\mu_{\Sigma}$. Contrary to the case without spin, the zero-point energy $F_{q}^{\rm ZP}$ in Eq. (26) exhibits an explicit dependence on the spin potential. As pointed out in Ref. Farias et al. (2016), the arbitrariness of the regularization methods should apply only to vacuum integrals, i.e., integrals that do not depend on the parameters of the medium, $T$ and $\mu_{\Sigma}$. We therefore write

[TABLE]

and subsequently express the zero-point quark contribution to the free energy as

[TABLE]

where $I_{1}$ is defined in Eq. (29) and we introduce two more integrals, explicitly given as

[TABLE]

while $\delta I_{\Sigma}$ denotes the regular remainder. After dimensional regularization, these integrals reduce to (see Appendix A)

[TABLE]

while the remainder takes the form (cf. Appendix A)

[TABLE]

with $p_{\Sigma}=\sqrt{\frac{1}{4}\mu_{\Sigma}^{2}-g^{2}\sigma^{2}}$ and $\theta_{\Sigma}=\theta(\tfrac{1}{2}|\mu_{\Sigma}|-g|\sigma|)$.

A naïve regularization scheme would be to again enforce that the regularized potential agrees with $V_{\sigma}$ in Eq. (28) for any value of $\mu_{\Sigma}$ when $T=0$. This requires us to employ medium-dependent bare parameters, which become functions of $\mu_{\Sigma}$. This procedure is in principle allowed in effective models that arise as limits of more complicated, fundamental models. Hence, we introduce

[TABLE]

Compared to Eq. (31), the above expression has an extra logarithmic term, to compensate the analogous term coming from $I_{2}$ in Eq. (37), as well as a term involving the step function $\theta_{\Sigma}=\theta(\tfrac{1}{2}|\mu_{\Sigma}|-g|\sigma|)$, ensuring the cancellation of the part that depends on $p_{\Sigma}$ in Eq. (38). Imposing Eq. (32) again, we arrive at

[TABLE]

while $\lambda_{\rm n}=\lambda_{0}$ from Eq. (33) and the step function $\tilde{\theta}_{\Sigma}=1-\theta_{\Sigma}$ is non-vanishing only when $g|\sigma|>\tfrac{1}{2}|\mu_{\Sigma}|$. In contrast to Eq. (33), the term $\lambda_{\rm n}v_{\rm n}^{4}$ no longer vanishes, as $v_{\rm n}^{2}=-\mu_{\Sigma}^{2}/2g^{2}+O(\epsilon)$. Hence, the offset $V_{\rm n}$ diverges as $\epsilon\to 0$. In the scheme discussed above, the bare parameters of the meson potential exhibit an explicit medium dependence, such that the potential after absorbing the zero-point terms becomes medium-independent. However, as we will show in Sec. V, this regularization scheme leads to results which are inconsistent with lattice data. We therefore consider an extension of this procedure, described below.

III.3 Renormalization with spin effects

It is conceivable that, after renormalization, the meson potential acquires a dependence on $\mu_{\Sigma}$. We retain the restriction that at $\mu_{\Sigma}=0$, Eq. (32) is satisfied. Furthermore, at vanishing temperature, we impose that the point $\sigma=f_{\pi}$ corresponds to a (local) minimum and that the $\sigma$ meson mass remains unchanged:

[TABLE]

where $M_{\sigma}=\sqrt{\lambda(3f_{\pi}^{2}-v^{2})}=600$ MeV (the vacuum mass of the $\sigma$ meson). To satisfy the above equations, we write the bare potential as

[TABLE]

In this second renormalization prescription, we consider that the parameters with the subscript “s” represent a $\mu_{\Sigma}$-dependent deformation of those in Eq. (39). For simplicity, we control this deformation by allowing the coefficient $L^{(2)}_{\rm s}=L_{\rm n}^{(2)}+\Delta L_{\rm s}^{(2)}$ to deviate from $L_{\rm n}^{(2)}$ as

[TABLE]

where $\ell$ is a free parameter. Note that, contrary to $L^{(2)}_{\rm n}$, $\Delta L^{(2)}_{\rm s}$ is non-vanishing regardless of the relation between $\mu_{\Sigma}$ and $g\sigma$. The parameters $(\lambda_{\rm s},v_{\rm s}^{2},V_{\rm s})$ in Eq. (42) receive $\ell$-dependent corrections with respect to $(\lambda_{\rm n},v_{\rm n}^{2},V_{\rm n})$ in Eq. (40), denoted $(\Delta\lambda_{\rm s},\Delta v_{\rm s}^{2},\Delta V_{\rm s})$, i.e.

[TABLE]

Imposing Eqs. (41) gives

[TABLE]

In our proposed scheme, $L^{(4)}_{\rm s}=L^{(4)}_{\rm n}$. The resulting renormalized meson potential $V_{\sigma}^{\rm ren}=V_{\sigma}^{\rm bare}+F_{q}^{\rm ZP}$ reads:

[TABLE]

An important consequence of Eq. (46) is that the coefficient of the $\sigma^{4}$ term in the meson potential becomes $\mu_{\Sigma}$-dependent. Crucially, the thermodynamic stability of the system requires that this coefficient be positive. This translates into the constraint $\ell>0$. If, on the contrary, $\ell<0$, then the sign of the coefficient of $\sigma^{4}$ becomes negative when the spin potential $\mu_{\Sigma}$ exceeds the following threshold:

[TABLE]

The naïve scheme discussed in Subsec. III.2 is recovered when $\ell=0$. In Sec. V, we will view $\ell>0$ as a free parameter of the scheme, whose value we will employ in Subsec. IV.3 to fit the lattice data reported in Ref. Braguta et al. (2025).

III.4 Phase structure

In the mean-field approximation, the expectation value of the $\sigma$ meson is obtained by demanding the minimization of the free energy. Figure 1 shows the thermodynamic potential density at vanishing temperature, $F=V^{\rm ren}_{\sigma}$. In panel (a), we fixed $\ell=1$ and varied the spin potential, $\mu_{\Sigma}$. At low values of $\mu_{\Sigma}$, the dominant minimum of $V^{\rm ren}_{\sigma}$ is at $\sigma=f_{\pi}$ and the system is in the chirally-broken state. As $\mu_{\Sigma}$ increases, a second minimum of $V^{\rm ren}_{\sigma}$ close to $\sigma\simeq 0$ becomes dominant and chiral symmetry is restored. Panel (b) shows that, at fixed $\mu_{\Sigma}=0.5$ GeV, increasing $\ell$ promotes chiral restoration, leading to the appearance of the second minimum, which eventually becomes dominant.

To pinpoint the value of $\sigma$ that minimizes the thermodynamic potential density at $T=0$, we impose the saddle point equation

[TABLE]

Figure 2 shows the variation of $\sigma/f_{\pi}$ with respect to the spin potential $\mu_{\Sigma}$. For $\mu_{\Sigma}<\mu_{\Sigma}^{c}$, the system remains in the chirally-broken phase. When $\mu_{\Sigma}>\mu_{\Sigma}^{c}$, a first-order phase transition to the chirally-restored phase can be observed. Since at $T=0$, $\partial V^{\rm ren}_{\sigma}/\partial\sigma$ depends on $\ell$ and $\mu_{\Sigma}$ only through the combination $\ell\mu_{\Sigma}^{2}$, the critical spin potential obeys the scaling $\mu_{\Sigma}^{c}\sqrt{\ell}={\rm const}$. We test this hypothesis in panel (b) of Fig. 2 and we observe that this constant equals $\mu_{\Sigma}\sqrt{\ell}\simeq 0.5$ GeV.

We can get a better understanding of the chiral restoration transition by considering a small-$\sigma$ expansion of the renormalized meson potential $V_{\sigma}^{\rm ren}$ in Eq. (46),

[TABLE]

where

[TABLE]

The saddle point equation reduces to

[TABLE]

from where we can deduce that $\sigma\ln(\sigma/f_{\pi})=-[h+(2a+b)\sigma]/(2b)$. At the minimum of the thermodynamic potential, we have

[TABLE]

Demanding that the phase transition happens at the above point is equivalent to imposing $V^{\rm ren}_{\sigma}=0$, which leads to

[TABLE]

Imposing now that $\sigma\ln(\sigma/f_{\pi})=-[h+(2a+b)\sigma]/(2b)$ yields the solution $\mu_{\Sigma}\sqrt{\ell}\simeq 0.50$ GeV, in agreement with panel (b) of Fig. 2.

III.5 Pressure and spin density

We now discuss the system pressure, $P=-F=-V^{\rm ren}_{\sigma}-F_{q}$. At vanishing temperature, the thermal quark contribution vanishes, $F_{q}=0$, and the pressure receives only the zero-point contribution, $P=P^{\rm ZP}=-V^{\rm ren}_{\sigma}$. In the chirally-broken state, $\sigma=f_{\pi}$ and $V^{\rm ren}_{\sigma}$ vanishes, see Eq. (46), therefore $P^{\rm ZP}(\sigma=f_{\pi})=0$. As discussed in Sec. III.4, the chiral restoration takes place when $\mu_{\Sigma}\sqrt{\ell}\simeq 0.5$ GeV. By the Gibbs construction, at this point the pressure must be equal to its vacuum value, $P=0$, increasing as $\mu_{\Sigma}^{2}\ell$. Since in the chirally-restored phase, $\sigma\simeq 0$, the scaling with $\mu_{\Sigma}\sqrt{\ell}$ of the chiral restoration transition transfers also to the pressure, as can be seen in Fig. 3(a).

At a finite spin potential $\mu_{\Sigma}$, the system can develop a non-vanishing spin polarization $Q_{\Sigma}=\langle\mathcal{S}^{z}\rangle$ corresponding to an average spin density $\mathcal{S}^{z}=\psi^{\dagger}\Sigma^{z}\psi$ along the $z$ direction. In the grand canonical ensemble, the spin density can be obtained via

[TABLE]

In general, this quantity can be represented as the sum of two terms, $Q_{\Sigma}=Q_{\Sigma}^{\rm ZP}+Q_{\Sigma}^{q}$, corresponding to the contributions from the renormalized meson potential at zero temperature, $Q_{\Sigma}^{\rm ZP}=-\partial V_{\sigma}^{\rm ren}/\partial\mu_{\Sigma}$, and from the quark thermal potential, $Q_{\Sigma}^{q}=-\partial F_{q}/\partial\mu_{\Sigma}$, respectively. Focusing on the contribution at vanishing temperature, when $F=V_{\sigma}^{\rm ren}$ and

[TABLE]

it can be seen that $Q^{\rm ZP}_{\Sigma}$ vanishes in the chirally-broken phase, when $\sigma=f_{\pi}$. In the chirally-restored phase, the condensate is very small, $\sigma\ll f_{\pi}$, and Eq. (55) implies that $Q_{\Sigma}^{\rm ZP}$ increases linearly with the spin potential:

[TABLE]

Contrary to the pressure and the point of the first-order phase transition, which scale as $\mu_{\Sigma}\sqrt{\ell}$, the spin density $Q_{\Sigma}$ scales as $\mu_{\Sigma}\ell$. Figure 3(b) shows the $Q_{\Sigma}$ as a function of $\mu_{\Sigma}\sqrt{\ell}$, for various values of $\ell$.

IV Finite temperature

After renormalizing the zero-point contribution to the free energy, in this section we turn our attention to its thermal part and examine certain subtle aspects of the energy dispersion relation in the presence of a finite spin potential.

IV.1 Thermal contribution to the thermodynamic potential

We start by writing the thermal part $F_{q}$ of the quark free energy, given in Eq. (27), with respect to spherical coordinates:

[TABLE]

In order to study the thermodynamic properties of the system, the mean value of the $\sigma$ meson is obtained by enforcing the saddle point equation:

[TABLE]

where $V_{\sigma}^{\mathrm{ren}}$ is introduced in Eq. (46). We now write down the expression of $\partial F_{q}/\partial\sigma$ with respect to the spherical coordinates introduced above:

[TABLE]

The presence of $E_{\mathbf{p}}^{(s)}$ in the denominator of Eq. (59) reveals an unexpected infrared divergence. As can be seen from Eq. (19), $E_{\mathbf{p}}^{(s)}$ vanishes when $s\mu_{\Sigma}=E_{p_{z}}=\sqrt{p^{2}\cos^{2}\theta+g^{2}\sigma^{2}}$ and $p_{\perp}=0$. This integrable singularity can be removed by changing the integration variable from $p$ to $E_{\mathbf{p}}^{(s)}\to E$, as discussed in Appendix B. Without repeating the details here, several mathematical manipulations allow Eqs. (57) and (59) to be put in the following compact form:

[TABLE]

where $s=\pm 1/2$, $p=\sqrt{E^{2}-g^{2}\sigma^{2}}$ and $\widetilde{E}_{s}=E-s\mu_{\Sigma}$.

The modulus $|\widetilde{E}_{s}|$ appearing in Eq. (60a) indicates that the spin potential $\mu_{\Sigma}$ does not introduce a Fermi level. It is easy to check that $\lim_{T\to\infty}F_{q}=0$, regardless of the value of $\mu_{\Sigma}$. To gain more insight into the properties of $F_{q}$, we employ integration by parts to arrive at

[TABLE]

Keeping in mind that $\widetilde{E}_{s}=E-s\mu_{\Sigma}$, the following relation,

[TABLE]

shows that $\mu_{\Sigma}$ acts as the usual baryon (vector) chemical potential, while the step function serves to subtract the $T\to 0$ contribution, thereby eliminating the Fermi level contributions from the corresponding quantities. The first two terms in Eq. (61) are easily recognizable as the pressure $P_{\rm FD}(\mu,T)$ and charge density $Q_{\rm FD}=\partial P_{\rm FD}(\mu,T)/\partial\mu$ of a Fermi gas with temperature $T$ and usual (vector) chemical potential equal to $\mu=\mu_{\Sigma}/2$, allowing $F_{q}$ to be expressed as

[TABLE]

where $\Delta P_{\rm FD}(\mu,T)=P_{\rm FD}(\mu,T)-P_{\rm FD}(0,\mu)$ refers to the Fermi-Dirac pressure with the Fermi-level contribution (corresponding to $T=0$) subtracted, and similar for $\Delta Q_{\rm FD}(\mu,T)$, while

[TABLE]

The term $\delta F_{q}$ plays a role also in the evaluation of the fermion condensate $\langle\bar{\psi}\psi\rangle$, defined via $\partial F_{q}/\partial\sigma=g\,\langle\bar{\psi}\psi\rangle$:

[TABLE]

with $\Delta\langle\bar{\psi}\psi\rangle_{\rm FD}=-g^{-1}\partial\Delta P_{\rm FD}/\partial\sigma$.

We now discuss the small-mass limit of $F_{q}$. The Fermi-Dirac pressure is known:

[TABLE]

where terms of combined sixth order in $g\sigma$ and $\mu_{\Sigma}$ were neglected. Subtracting the $T\to 0$ part of $P_{\rm FD}$ in the same limit,

[TABLE]

we arrive at

[TABLE]

where we assumed the hierarchy $g\sigma<|\mu|$, valid at large temperatures, where $g\sigma=O(T^{-2})$. Differentiating the pressure (68) with respect to $\mu$, we get

[TABLE]

Finally, for the term $\delta F_{q}$ at small $g\sigma$, we have, to leading order in $\mu_{\Sigma}/T$ (see Appendix C):

[TABLE]

Substituting the small-mass expansions from Eqs. (68), (69) and (70) into Eq. (63) gives

[TABLE]

It is worth remarking that, when the mass is negligible, the quark contribution to the thermodynamic potential becomes independent of $\mu_{\Sigma}$.

IV.2 Phase structure

We now consider the chiral restoration phase transition at finite temperature and spin potential. Since the renormalized free energy depends on the parameter $\ell$, in this subsection we investigate the phase structure of the model as a function of $\ell$. These results will then allow us to fix the parameter $\ell$ in the following section using the lattice data reported in Ref. Braguta et al. (2025).

At a vanishing spin potential, the system undergoes a crossover transition to the chirally-restored phase at the pseudocritical temperature $T_{\rm pc}=0.147$ GeV. It is worth noticing that this temperature does not depend on the renormalization parameter $\ell$, which is associated with a finite spin potential $\mu_{\Sigma}$.

In Fig. 4, we show the thermodynamic potential as a function of $\sigma/f_{\pi}$ for various values of $\mu_{\Sigma}$ with two specified combinations of $\ell$ and $T$. In panel (a), we study the effect of the spin potential in the case $\ell=0$, corresponding to the naïve renormalization approach, discussed in Sec. III.2 with the temperature set at $T=0.18$ GeV (above the crossover transition temperature). The absolute minimum of $F$ (denoted by the black blobs) shifts to higher values of $\sigma$ as $\mu_{\Sigma}$ is increased, signaling the breaking of chiral symmetry. In panel (b), we consider a larger value $\ell=1$ and fix $T=0.14$ GeV (below the pseudocritical temperature). In this case, increasing $\mu_{\Sigma}$ causes the minimum of $F$ to shift towards lower values of $\sigma$, indicating chiral restoration. For the particular case shown in panel (b), the chiral transition is of first order.

In what follows, we always set the mean-field value of the $\sigma$ meson corresponding to the minimum of the total thermodynamic potential. In Fig. 5, we show $\sigma/f_{\pi}$ as a function of $T$ for various values of $\mu_{\Sigma}$, with two fixed values of $\ell$. For the case $\ell=0$, shown in panel (a), the renormalized meson potential $V_{\sigma}^{\rm ren}$ is independent of $\mu_{\Sigma}$, such that $F$ receives $\mu_{\Sigma}$ contributions only from the thermal quarks via $F_{q}$. It can be seen that increasing $\mu_{\Sigma}$ causes the transition point to migrate towards larger temperatures, inhibiting chiral restoration. In contrast, for the case $\ell=1$, shown in panel (b), the thermodynamic potential receives $\mu_{\Sigma}$ contributions also via $V_{\sigma}^{\rm ren}$. We observe two effects: first, the transition point shifts to lower temperatures, indicating that the spin potential promotes chiral restoration; second, the transition becomes sharper, eventually turning into a first-order transition. Comparing panel (a) with panel (b), we see that the thermal quark contribution $F_{q}$ inhibits and softens chiral symmetry restoration, while the meson contribution $V_{\sigma}^{\rm ren}$ favors and sharpens it.

IV.3 Fit to lattice data

We are now in a position to make the connection between the effective model study presented in this paper and the lattice QCD analysis reported in Ref. Braguta et al. (2025). Due to the infamous sign problem, that work considered an imaginary spin potential, $\mu_{\Sigma}^{I}=-i\mu_{\Sigma}$, and found that the pseudocritical transition temperature, $T_{\rm pc}(\mu_{\Sigma}^{I})$, increases quadratically with $\mu_{\Sigma}^{I}$. This behavior is characterized by the curvature $\kappa_{\Sigma}$, as defined in Ref. Braguta et al. (2025):

[TABLE]

valid at small $\mu_{\Sigma}^{I}/T$. In QCD, the above formula can refer to the pseudocritical temperature of the deconfinement or chiral transitions, when $\kappa_{\Sigma}\to\kappa_{\Sigma}^{L}$ or $\kappa_{\Sigma}^{\psi}$, respectively. In this paper, we focus on the latter coefficient, extrapolated to physical pseudoscalar and vector meson masses, which was reported in Ref. Braguta et al. (2025) as

[TABLE]

where the number in the parentheses indicates the numerical uncertainty in the last two digits corresponding to a combined statistical and systematic error of Monte Carlo simulations.

Moving now to the case of a real spin potential, Eq. (72) can be analytically continued as follows:

[TABLE]

where $\kappa_{\Sigma}^{(2)}$ is expected to retain its value in Eq. (73). Equation (74) indicates a decreasing dependence of $T_{c}$ on $\mu_{\Sigma}$. As shown in Fig. 5, this behavior is consistent with the effect of the renormalized meson potential, $V_{\sigma}^{\rm ren}$ (visible when $\ell=1$), and opposite to the effect of the thermal quark contribution, $F_{q}$ (visible when $\ell=0$).

The above discussion indicates that the parameter $\ell$ allows us to change the regimes between the thermal quark-dominated behavior ($\ell=0$) to the renormalized meson-dominated behavior (at $\ell\gtrsim 1$). In panel (a) of Fig. 6, we represented the ratio $T_{c}(\mu_{\Sigma})/T_{c}(0)$ with respect to $\mu_{\Sigma}/T_{c}(0)$ for various values of $\ell$. The points marked “Lattice fit” correspond to the values of $T_{c}(\mu_{\Sigma})/T_{c}(0)$ computed using Eq. (74) up to the second order, using the curvature $\kappa_{\Sigma}^{(2)}=\kappa_{\Sigma}^{\psi}\simeq 0.06$, compatible with the lattice result (73) within the statistical uncertainty of the numerical simulations. As expected, when $\ell$ is small, the pseudocritical temperature increases with $\mu_{\Sigma}$, in contrast to the lattice result. The critical value of $\ell$ corresponding to vanishing curvature is $\ell_{0}\simeq 0.70$, while for $\ell>\ell_{0}$, the curvature becomes negative.

Panel (b) of Fig. 6 shows the fitted values of $\kappa_{\Sigma}^{(2)}$ and $\kappa_{\Sigma}^{(4)}$ as functions of the renormalization parameter $\ell$. The curvature $\kappa_{\Sigma}^{(2)}$ shows a monotonically-increasing trend. At $\ell=0$, one finds $\kappa_{\Sigma}^{(2)}\simeq-0.022$ and $\kappa_{\Sigma}^{(4)}\simeq 0.005$. The value $\kappa_{\Sigma}^{(2)}=0$ is reached at $\ell=0.70$, as indicated in panel (a). The range given in Eq. (73), shown with a thin, blue horizontal band around $\kappa_{\Sigma}^{\psi}=0.0595$, corresponds to the following interval of the parameter $\ell$:

[TABLE]

which is shown as a thin, vertical purple band around the value $\ell=2.61$. Our fitting procedure also allows us to predict the value for the quartic coefficient:

[TABLE]

which was not reported in Ref. Braguta et al. (2025). For definiteness, we will employ this central value of $\ell=2.61$ as the “lattice-compatible” $\ell$ in what follows, which corresponds to the rounded value $\kappa_{\Sigma}^{\psi}=0.06$. The fit was performed in the range $\mu_{\Sigma}\in[0,0.1]$.

IV.4 Pressure and spin density

We now consider the behavior of the system pressure at finite $T$, $P=-F=-V_{\sigma}^{\rm ren}-F_{q}$. In Fig. 7(a), we show $P$ as a function of temperature for various values of $\mu_{\Sigma}$. The pressure is normalized by the Stefan-Boltzmann expression, $P_{\rm SB}=7N_{c}N_{f}\pi^{2}T^{4}/180$, corresponding to the leading-order contribution to the large temperature expansion, shown in Eq. (66). As shown in Fig. 7(a), $P$ becomes independent of $\mu_{\Sigma}$ at large temperatures. At vanishing temperature, $P$ vanishes in the chirally-broken phase, whenever $\mu_{\Sigma}\sqrt{\ell}\lesssim 0.5$ GeV (see discussion in Sec. III.4). Above this threshold value, the zero-point pressure exhibits a quadratic dependence on $\mu_{\Sigma}$, with the thermal contribution becoming dominant when $T^{2}\gtrsim 0.09gf_{\pi}\mu_{\Sigma}\sqrt{\ell}$.

In the case of the spin charge, $Q_{\Sigma}=-\partial F/\partial\mu_{\Sigma}=Q^{\rm ZP}_{\Sigma}+Q^{q}_{\Sigma}$, the zero-point contribution $Q^{\rm ZP}_{\Sigma}$ is shown in Eq. (55), while the thermal quark contribution can be obtained from $F_{q}$, shown in Eq. (60a):

[TABLE]

As discussed in Sec. III.5, $Q_{\Sigma}^{\rm ZP}$ vanishes in the chirally-broken phase, when $\sigma=f_{\pi}$. On the other hand, in the chirally-restored phase,

[TABLE]

becomes negligible, since $\sigma=O(T^{-2})$ is small. Figure 7(b) confirms that after the chiral restoration, $Q_{\Sigma}$ is dominated by the zero-point contribution, $Q_{\Sigma}^{\rm ZP}$, shown in Eq. (55). Before the chiral restoration, $Q_{\Sigma}^{\rm ZP}\simeq 0$ and the buildup before the phase transition is largely due to the quark thermal contribution, $Q_{\Sigma}^{q}$.

V Phase diagram

In this section, we focus on the phase diagram corresponding to the chiral phase transition in the presence of a spin potential, represented in the $(\mu_{\Sigma},T)$ plane.

Panel (a) of Fig. 8 shows the phase diagram corresponding to values of $\ell$ between [math] and $5$. For $\ell=0$, the transition remains of crossover type for any $\mu_{\Sigma}$ and the transition temperature increases with $\mu_{\Sigma}$. For any finite $\ell$, the discussion in Sec. III.4 reveals that, at $T=0$, the system undergoes a first-order phase transition when $\mu_{\Sigma}\sqrt{\ell}\simeq 0.50$ GeV. Since the transition at $\mu_{\Sigma}=0$ is always of crossover type, this implies that for any $\ell>0$, the phase diagram exhibits a critical point (shown by the dark star symbols). Furthermore, the requirement that the system undergoes a phase transition at $T=0$ forces the transition line to bend downwards. Thus, the renormalized meson potential $V_{\sigma}^{\rm ren}$ competes with the effect of the thermal quark contribution, $F_{q}$, as can be seen in the case $\ell=0.5$. When the curvature $\kappa_{\Sigma}^{(2)}$ changes sign, i.e. for $\ell\gtrsim 0.70$, the renormalized meson potential becomes dominant and the transition temperature stops increasing with $\mu_{\Sigma}$, for small values of $\mu_{\Sigma}$. Even when $\kappa_{\Sigma}^{2}$ is positive, our numerical calculations reveal a non-monotonic dependence of the transition temperature on $\mu_{\Sigma}$. The temperature monotonically decreases with $\mu_{\Sigma}$ when $\ell>1.22$, as seen in Fig. 8(a).

For any $0<\ell\lesssim 1.22$, for example $\ell=0.5$, shown in Fig. 8(b), we notice that the competition between $V_{\sigma}^{\rm ren}$ and $F_{q}$ leads to a nonmonotonic behavior of the transition line. On the right side of the plot, a vertical crossing of the transition line at $0.706$ GeV $\lesssim\mu_{\Sigma}\lesssim 0.740$ GeV drives the system from the chirally-restored phase (small $T$) to the chirally-broken phase and then back to the chirally restored phase, via first-order phase transitions. This property is illustrated in panel (c) of Fig. 8, where we can see the behavior of $\sigma/f_{\pi}$ with respect to temperature for the values of $\mu_{\Sigma}$ corresponding to the vertical dotted lines shown in panel (b) of Fig. 8.

Similarly, in the top region, keeping $T$ fixed for the narrow temperature range $T_{c}(0)\simeq 0.15$ GeV $\lesssim T\lesssim 0.17$ GeV and varying $\mu_{\Sigma}$, one drives the system from the chirally-restored phase to the broken phase through a crossover transition, and back to the restored phase through a first-order transition. This behaviour is illustrated in panel (d) of Fig. 8, where $\sigma/f_{\pi}$ is represented with respect to $\mu_{\Sigma}$ for the values of $T$ corresponding to the horizontal dashed lines shown in panel (b) of Fig. 8.

We now discuss the scaling of the transition line as $\ell$ is increased to arbitrarily large values. The transition temperature $T_{c}(0)$ at vanishing spin potential represents a common point on the transition line for all values of $\ell$. Similarly, the transition spin potential at vanishing temperature satisfies the scaling $\mu_{\Sigma}\sqrt{\ell}\simeq 0.5$ GeV. Thus, representing the phase diagram with respect to the $(T,\mu_{\Sigma}\sqrt{\ell})$ axes reveals that the transition line converges as $\ell$ is increased to a unique line, as shown in Fig. 9. When $\ell\gg 1$ but $\ell\mu_{\Sigma}^{2}$ is kept finite, the spin potential becomes vanishingly small and the thermal quark contribution can be evaluated at vanishing $\mu_{\Sigma}$. The only contribution is due to the renormalized meson potential, which contains the quadratic term $\ell\mu_{\Sigma}^{2}$.

Finally, panel (a) of Fig. 10 shows the phase diagram for the value $\ell=2.61$, corresponding to the curvature parameter $\kappa_{\Sigma}^{(2)}\simeq 0.06$ computed in lattice QCD (see discussion in Sec. IV.3). The vertical ($T$) and horizontal ($\mu_{\Sigma}$) axes are scaled by the pseudocritical temperature $T_{c}(0)$ at a vanishing spin potential. At low $\mu_{\Sigma}$, the transition is of crossover type, indicated as a dotted segment. The crossover is separated from the region of first-order phase transition (the solid line) by a critical endpoint (CEP, denoted by the black star) at

[TABLE]

As $\mu_{\Sigma}$ increases further, the line of the first-order phase transition hits the $T=0$ axis at the critical point

[TABLE]

At higher values of the spin potential, $\mu_{\Sigma}>\mu_{\Sigma,{\rm c}}$, the thermodynamic ground state in this model corresponds to a chirally restored phase for all temperatures.

In panel $(b)$ of Fig. 10, we show the behavior of the control parameter $\sigma/f_{\pi}$ with respect to the search parameter function $G=\sqrt{T^{2}+\mu_{\Sigma}^{2}}$, for three values of the angle $\theta=\tan^{-1}(\mu_{\Sigma}/T)$, indicated by the dotted lines in panel (a) of Fig. 10. The angles $\theta\simeq 0.26$ rad and $\theta\simeq 1.05$ rad lie deep in the crossover and first-order regions, respectively, while $\theta\simeq 0.63$ rad is in the first-order region, just after the critical point.

VI Conclusions

In this work, we investigated the impact of a finite spin density on the QCD phase diagram, with spin polarization characterized by the spin potential, $\mu_{\Sigma}$. The spin potential represents the energy difference associated with adding, to a system of relativistic fermions, a particle of one spin orientation relative to the other.222See the discussion on the precise thermodynamical interpretation of the spin potential at the end of Subsec. II.2. It quantifies the preference for a particular spin state and provides a natural thermodynamic variable that incorporates the spin polarization into the thermodynamic description of the system.

To study the effect of the spin polarization on the thermodynamics of strongly interacting matter, we employ the Linear Sigma model coupled to quarks, ${\rm LSM}_{q}$, which describes both quark and mesonic degrees of freedom. We work in the mean-field approximation, in which fluctuations of the bosonic fields around their ground-state values are neglected, and only the contribution arising from the fermion determinant is taken into account. This approach simplifies the treatment of the thermodynamics of the model, while preserving the essential features of quark dynamics.

We demonstrate that the presence of a finite spin potential, $\mu_{\Sigma}\neq 0$, leads to a nontrivial dependence of the energy of a fermionic quasiparticle on its momentum (19). The particularities of the fermionic energy become especially important in the coupling regime where the spin potential dominates over the mass of the fermionic quasiparticle, $\tfrac{1}{2}|\mu_{\Sigma}|>g|\sigma|$. In this domain of parameters — which is especially important in the transition region between the chirally broken and chirally restored phases — the mean-field analysis of the thermodynamic ground state requires a more careful treatment to properly capture the modifications induced by a finite spin density. We discuss these subtle questions in our work.

The nontrivial structure of the energy of the fermionic quasiparticle at a finite spin potential has several important consequences. First of all, one immediately finds that the zero-point (vacuum) free energy generated by a single fermionic loop develops an ultraviolet-divergent contribution that explicitly depends on $\mu_{\Sigma}$. This feature appears to be in sharp contrast with the standard field-theoretical paradigm, in which the ultraviolet divergencies are only associated with vacuum fluctuations of the fields that are independent of the characteristics of matter such as temperature and/or chemical potentials. Therefore, in the standard case, the ultraviolet divergencies — both in the presence or in the absence of matter — can systematically be incorporated through a well-defined loop-by-loop renormalization of the bare couplings of the theory by matter-independent vacuum effects Kapusta and Gale (2006). Therefore, one arrives at an intuitive picture that separates vacuum and matter effects: the thermal and quantum fluctuations of matter generate a finite correction to thermodynamic quantities, while the divergent quantum contributions coming from the vacuum can be systematically isolated and treated within a consistent renormalization scheme.

The explicit dependence of the divergent terms on matter-specific properties challenges the standard renormalization procedure, which is employed to absorb ultraviolet divergences into redefined couplings of the original mode. This observation highlights the importance to use a consistent framework that disentangles medium effects—arising from $\mu_{\Sigma}$-dependent divergent contributions to the total energy—from the genuine divergent terms associated with the zero-point (vacuum) free energy. In this work, we adopt the general strategy of the medium-separation scheme, previously developed for quark matter in the presence of the chiral chemical potential Farias et al. (2016); Azeredo et al. (2024) and the isospin chemical potential Avancini et al. (2019); Lopes et al. (2021); Ayala et al. (2023). In both these cases, the corresponding chemical potentials produce nontrivial modifications of the quasiparticle dispersion relation, analogous to the effects generated by the spin potential in Eq. (19).

The renormalization procedure introduces a free model parameter, $\ell$, which encodes an interplay of vacuum fluctuations with the ones of spin-polarized matter. Although this parameter cannot be determined from field-theoretical considerations, its value can be constrained by matching the curvature of the chiral phase transition predicted within the LSMq to nonperturbative results of lattice QCD simulations Braguta et al. (2025). Due to the well-known sign problem, lattice computations are restricted to the regime of small spin polarizations, where such a comparison is feasible. Once fixed, our approach allows us to make a consistent extrapolation of lattice results beyond this limited domain and, therefore, extend the description of the chiral phase transition to the full $(\mu_{\Sigma},T)$ phase plane.

As a result, we obtain the phase diagram shown in Fig. 10(a). The phase diagram contains a chirally broken phase separated from the chirally restored phase by a transition line. This chiral transition is a crossover in a small domain of values of the spin potential $\mu_{\Sigma}$. For $|\mu_{\Sigma}|\ll T_{c}(0)$, the transition line can be expressed via the series (74), where the curvature of the quadratic term (73) has already been found in the lattice simulations Braguta et al. (2025), while the value of the quartic term, that has not been determined numerically, is predicted by our analytical model in Eq. (76).

The transition becomes of first order at the critical endpoint given in Eq. (79). The transition-line continues at higher spin potentials and intersects the $T=0$ axis at another endpoint (80). At higher spin potentials, the spin-polarized matter resides in the chirally restored phase, regardless of temperature.

Our investigations revealed several non-intuitive properties of the QCD matter at finite spin potential $\mu_{\Sigma}$. First, the spin potential does not introduce a Fermi level while leading to a non-trivial modification to the energy dispersion relation (18). The first-order phase transition encountered at vanishing temperature, $T=0$, is driven by the renormalized mesonic potential. Second, the thermal quark contribution to the thermodynamic potential becomes independent of $\mu_{\Sigma}$ in the high-temperature, chirally-restored phase. Here, the polarization effects on the QCD matter are entirely governed by the renormalized meson potential. Third and finally, the spin polarization saturates with respect to increasing the temperature, being given by a temperature-independent term derived from the renormalized mesonic potential. In the chirally-restored phase, the thermal quark contribution to the spin polarization becomes negligible.

The emergence of the spin degree of freedom as a thermodynamical variable opens a new avenue of research of QCD thermodynamics, including the influence of spin polarization on the phase diagram of QCD, which could potentially include the emergence of new, spin-promoted phases of QCD. With the current setup we have, one possible extension can be to explore the deconfinement phase transition properties, using the Polyakov loop extended LSMq (PLSMq) model.

Note added. While this paper has been in the final stages of preparation, a similar study has emerged in a generalization of a Polyakov–loop enhanced Nambu–Jona-Lasinio model in Ref. Farias and Tavares (2025). Although the predicted phase diagram qualitatively agree with our result in Fig. 10(a), there are also two important differences: contrary to Ref. Farias and Tavares (2025), our approach allows us to match precisely the curvature of the chiral transition to the first-principle lattice data at small polarization which gives us the critical end point (CEP) at substantially lower value of the spin potential (79).

Appendix A Dimensional Regularization

In the context of zero-point renormalization within the standard LSMq model, we introduced the integral $I_{1}$ in Eq. (29). Additionally, while addressing the axial-vector interaction term, we introduced two more integrals, $I_{2}$ and $I_{3}$ in Eqs. (36a),(36b). Below, we present the evaluations of these integrals using dimensional regularization. We also discuss the computation of the finite remainder $\delta I_{\Sigma}$, introduced in Eq. (35).

Under the dimensional regularization prescription, we represent $d$-dimensional momentum integral as

[TABLE]

The angular integral can be computed via

[TABLE]

For future convenience, we note that the integral of even powers of $p^{z}$ can be obtained as

[TABLE]

To cope with divergent terms, the dimension of the spatial integral is replaced by $d\to 3-2\epsilon$, with $\epsilon$ a regularization parameter to be taken to [math] at the end of the calculation. To preserve the dimensionality of the expression, a renormalization energy scale $M$ is employed as a pre-factor $M^{2\epsilon}$. Within the minimal subtraction scheme, one can then express

[TABLE]

The angular integral is evaluated by analytical continuation:

[TABLE]

The next steps involve evaluating the integrals and expanding the results in powers of $\epsilon$. The $1/\epsilon$ contributions represent the divergent terms that must be removed by suitable counterterms.

Applying the dimensional regularization to the integrals $I_{1}$, $I_{2}$ and $I_{3}$, given in Eqs. (29), (36a) and (36b) of the main text, we obtain

[TABLE]

The integrals with respect to $p$ can be performed using the formula:

[TABLE]

leading to

[TABLE]

Taking the limit $\epsilon\to 0$, we employ, for all non-negative integers $n\geqslant 0$, the following relation:

[TABLE]

with $\psi(z)=d\ln\Gamma(z)/dz$ being the polygamma function. Retaining terms up to $O(\epsilon^{0})$, we uncover the results in Eqs. (30) and (37) presented in the main text.

The remainder $\delta I_{\Sigma}$, introduced in Eq. (35), is regular both in the IR and in the UV, and reads:

[TABLE]

Due to the axial symmetry of the integrand, it is convenient to perform the integral with respect to cylindrical coordinates, $d^{3}p\to 2\pi dp^{z}dp_{\perp}p_{\perp}$. We perform the improper integration with respect to $p_{\perp}$ as the limit $\int_{0}^{\infty}dp_{\perp}\,f(p_{\perp})=\lim_{P_{\perp}\to\infty}\int_{0}^{P_{\perp}}dp_{\perp}f(p_{\perp})$. Using the relations

[TABLE]

with $E_{\bf P}^{(s)}=\sqrt{P_{\perp}^{2}+(E_{p_{z}}-s\mu_{\Sigma})^{2}}$ and $E_{P}=\sqrt{P_{\perp}^{2}+E_{p_{z}}^{2}}$, as well as the expansions

[TABLE]

and $E_{P}^{-2n-1}\simeq 0$ for all $n\geq 0$, we arrive at

[TABLE]

When $E_{p_{z}}>\frac{1}{2}|\mu_{\Sigma}|$, the sum over $s=\pm 1/2$ yields $\frac{1}{2}\sum_{s}|E_{p_{z}}-s\mu_{\Sigma}|^{3}=E_{p_{z}}^{3}+\frac{3}{4}E_{p_{z}}\mu_{\Sigma}^{2}$, leading to an exact cancellation of the integrand. Therefore, $\delta I_{\Sigma}$ receives non-vanishing contribution only when $E_{p_{z}}<|\mu_{\Sigma}|$:

[TABLE]

with $p_{\Sigma}=\sqrt{\frac{1}{4}\mu_{\Sigma}^{2}-g^{2}\sigma^{2}}$ and we took into account that $E_{p_{z}}<\frac{1}{2}|\mu_{\Sigma}|$ implies that $\frac{1}{2}|\mu_{\Sigma}|>g|\sigma|$. Performing the above integration leads to the result shown in Eq. (38).

Appendix B Evaluation of the thermal contribution

In this appendix, we discuss the transition from the representation of the thermal part of the thermodynamic potential, $F_{q}$ (57), and its derivative $\partial F_{q}/\partial\sigma$ (59), expressed with respect to the spin-deformed fermion energy $E_{\mathbf{p}}^{(s)}$ defined in Eq. (19), to the more convenient forms in Eqs. (60a) and (60b), respectively. We begin in Sec. B.1 by pointing out that the energy $E_{\mathbf{p}}^{(s)}$ can be degenerate with respect to $\mathbf{p}$, in the case when $s\mu_{\Sigma}>g\sigma$. Sections B.2 and B.3 address separately the cases $g\sigma>s\mu_{\Sigma}$ (case 1) and $g\sigma<s\mu_{\Sigma}$ (case 2).333In the case when $\sigma$ takes a negative value, all instances of $\sigma$ should be replaced by $|\sigma|$.

B.1 Energy branches

Noting the relation (19) between the fermion energy $E_{\mathbf{p}}^{(s)}$ expressed with respect to the spherical coordinates $(p,\theta,\varphi)$, and the momentum magnitude $p$, we derive

[TABLE]

The above relation allows the integration variable to be changed from $p$ to $E_{\mathbf{p}}^{(s)}$ in Eqs. (57) and (59):

[TABLE]

where $\Gamma_{E}$ represents a suitable contour for the energy integral (see below). In the above, we introduced the angular integrals

[TABLE]

where the integration domain $\Gamma_{\theta}$ depends on the value of the energy $E$, as discussed below. Note that now we view $p\to p_{s}(E)$ and $E_{p_{z}}\to E_{p_{z}}^{s}=\sqrt{p_{s}^{2}\cos^{2}\theta+g^{2}\sigma^{2}}$ as functions of the energy $E$ and the spin branch $s=\pm\tfrac{1}{2}$.

Let us now critically analyze the energy structure to get a handle on the energy contour $\Gamma_{E}$. Setting $p=0$ in Eq. (19), the energy evaluates to $E_{0}=|s\mu_{\Sigma}-g\sigma|$. On the other hand, setting $dE_{\mathbf{p}}^{(s)}/dp=0$ at fixed $\theta$ in Eq. (95) implies $E_{p_{z}}=s\mu_{\Sigma}\cos^{2}\theta$. When $s\mu_{\Sigma}>0$, this relation reveals a critical momentum $p_{m}$,

[TABLE]

where the energy attains a minimum value $E_{m}$:

[TABLE]

The particular structures of $E_{0}$ and $E_{m}$ invoke two distinct cases, namely, case 1: $g\sigma>s\mu_{\Sigma}$ and case 2: $g\sigma<s\mu_{\Sigma}$. In case 1, $E_{m}^{2}<0$ and the energy increases monotonically, from $E_{0}=g\sigma-s\mu_{\Sigma}>0$ to infinity. For case 2, $E_{m}^{2}>0$ implies that the energy first decreases from $E_{0}=s\mu_{\Sigma}-g\sigma$ to $E_{m}$, then it increases back to $E_{0}$ and subsequently to infinity. Furthermore, demanding $p_{m}^{2}>0$, the expression under the square root of Eq. (98) must be positive, and this is achieved only when $|\cos\theta|$ exceeds a critical value,

[TABLE]

i.e., when $\cos\theta<-\cos\theta_{c}$ or $\cos\theta>\cos\theta_{c}$.

We now have a complete picture of the behaviour of the phase space in the $(E,\theta)$ coordinates, for case 2: for $|\cos\theta|<\cos\theta_{c}$, the energy is single valued, increasing from $E_{0}=E(p=0)=\tfrac{1}{2}|\mu_{\Sigma}|-g\sigma$ to infinity. When $\cos\theta_{c}<|\cos\theta|\leq 1$, the energy first decreases from $E_{0}$, when $p=0$, to $E_{m}$ when $p=p_{m}$. Note that $E_{m}=E_{0}$ when $|\cos\theta|=\cos\theta_{c}$ ($p_{m}=0$) and $E_{m}=0$ when $\cos\theta=\pm 1$ ($p_{m}=\sqrt{\frac{1}{4}\mu_{\Sigma}^{2}-g^{2}\sigma^{2}}$). From $E_{m}$, both the energy and the momentum magnitude increase monotonically, first to $E_{0}$ and $p_{M}=\sqrt{|\mu_{\Sigma}|(|\mu_{\Sigma}|\cos^{2}\theta-2g|\sigma|)}$, and then up to infinity. This behavior is illustrated in Fig. 11. In the remaining two subsections, we will evaluate both $F_{q}$ and $\partial F_{q}/\partial\sigma$ within the aforementioned two distinct scenarios.

B.2 Case 1: $g\sigma>s\mu_{\Sigma}$

For case 1, the energy contour follows a monotonic increase from $E_{0}^{1}=g\sigma-s\mu_{\Sigma}$ to infinity. Hence, $\Gamma_{E}$ in Eqs. (96) is replaced by the interval $E_{0}^{1}<E<\infty$. Similarly, the contour $\Gamma_{\theta}$ in the angular integrals in Eq. (97) is simply $-1\leq\cos\theta\leq 1$.

We now compute the angular integrals $I^{s}_{\theta,E}$ and $\delta I^{s}_{\theta,E}$. Using the relations $E_{p_{z}}^{s}=\sqrt{p_{s}^{2}\cos^{2}\theta+g^{2}\sigma^{2}}$ and $E^{2}=p_{s}^{2}+g^{2}\sigma^{2}+\frac{1}{4}\mu_{\Sigma}^{2}-2s\mu_{\Sigma}E_{p_{z}}^{s}$, we can derive

[TABLE]

and re-label $p_{s}\to p$ to further simplify $I^{s}_{\theta,E}$ and $\delta I^{s}_{\theta,E}$ as:

[TABLE]

with $E_{p_{z}}^{s}=(p^{2}+g^{2}\sigma^{2}+\frac{1}{4}\mu_{\Sigma}^{2}-E^{2})/s\mu_{\Sigma}$ and

[TABLE]

Since $(p_{1}^{s})^{2}-(p_{0}^{s})^{2}=2s\mu_{\Sigma}(E+s\mu_{\Sigma}-g\sigma)$, $p_{1}^{s}>p_{0}^{s}$ when $s\mu_{\Sigma}>0$ and $p_{1}^{s}<p_{0}^{s}$ if $s\mu_{\Sigma}<0$.

Expressing $\cos\theta$ as

[TABLE]

we arrive at

[TABLE]

Since $p_{s;-}^{2}={\rm sgn}(s\mu_{\Sigma})[p^{2}-(p^{s}_{0})^{2}]$ and $p_{s;+}^{2}=p_{s;-}^{2}+2g\sigma|\mu_{\Sigma}|$, we can switch the integration variable to $y=p_{s;-}^{2}$, leading to

[TABLE]

where $Y_{s}=|\mu_{\Sigma}|(E-g\sigma+s\mu_{\Sigma})$. Using the relations

[TABLE]

with $a_{s}=|\mu_{\Sigma}|g\sigma$, we arrive at

[TABLE]

Replacing $I^{s}_{\theta,E}$ and $\delta I^{s}_{\theta,E}$ in Eqs. (96), we eventually obtain

[TABLE]

Shifting the energy to $E+s\mu_{\Sigma}\to E$ in Eqs. (109) recovers Eqs. (60) in the main text.

B.3 Case 2: $g\sigma<s\mu_{\Sigma}$

As discussed in Sec. B.1, the fermion energy $E_{\mathbf{p}}^{(s)}$ becomes degenerate when $g\sigma<s\mu_{\Sigma}$. Separating $F_{q}=F_{q}^{+}+F_{q}^{-}$, with $F_{q}^{\pm}$ corresponding to $s\mu_{\Sigma}=\pm\tfrac{1}{2}|\mu_{\Sigma}|$, we note that the $F_{q}^{-}$ contribution is compatible with case 1 discussed in Sec. B.2. We therefore focus for the remainder of this section on $F_{q}^{+}$, which we split as follows:

[TABLE]

where the symbols $<$, $>$ and $c$ refer to various integration ranges, as follows:

[TABLE]

where we took into account that the integrands of $I^{s}_{\theta,E}$ and $\delta I^{s}_{\theta,E}$ are even with respect to $\cos\theta$ and therefore the integration range runs only over non-negative values of $\cos\theta$. On the $<$ branch, $\cos\theta<\cos\theta_{c}$, $p_{m}^{2}<0$ and the energy is monotonic and increases from $E_{0}=\tfrac{1}{2}|\mu_{\Sigma}|-g\sigma$ to infinity. The $>$ branch corresponds to the region with $p>p_{M}$ (see end of Sec. B.1), on which the energy again increases monotonically from $E_{0}=\tfrac{1}{2}|\mu_{\Sigma}|-g\sigma$ to $\infty$. The $r=c$ branch spans the critical region, when $\cos\theta>\cos\theta_{c}$ and the energy first decreases from $E_{0}$ to $E_{m}$ and then increases back to $E_{0}$.

To get a handle of what happens within the integration regions shown in Eq. (111), we use Eq. (19) to investigate the relation between the momentum $p_{s}\to\{p_{0},p_{c},p_{1}\}$ and the energies $(E,E_{0})$, corresponding to $\cos\theta\in\{0,\cos\theta_{c},1\}$:

[TABLE]

Let us consider first the $>$ and $<$ branches. Here, $E>E_{0}$ and it can be seen that the above equations allow for two values of $p_{c}$ and $p_{1}$:

[TABLE]

It is easy to check that, when $E>E_{0}$, only the solutions corresponding to the upper (positive) sign are meaningful. In the case of $p_{c}$, we have $p_{c;-}<p_{0}<p_{c;+}$, so that $p_{c;-}$ is not part of the integration domain. In the case of $p_{1}$, the solution originates from $\sqrt{p_{1}^{2}+g^{2}\sigma^{2}}=\frac{1}{2}|\mu_{\Sigma}|\pm E$. Since $E>\frac{1}{2}|\mu_{\Sigma}|-g\sigma$, it is clear that $\frac{1}{2}|\mu_{\Sigma}|-E<g\sigma$, hence $p_{1,-}^{2}<0$ and this is not an acceptable solution. Therefore, the angular integration ranges $\Gamma_{\theta}^{<}=[0,\cos\theta_{c}]$ and $\Gamma_{\theta}^{>}=[\cos\theta_{c},1]$ are converted into the momentum integration ranges $p_{0}\leq p\leq p_{c;+}$ and $p_{c;+}\leq p\leq p_{1;+}$, which can be merged seamlessly. Therefore, the sum $F_{q}^{<}+F_{q}^{>}$ can be computed by merging the angular integration to $0\leq\cos\theta\leq\cos\theta_{c}$, with the momentum $p_{s}(\cos\theta)$ continuously and monotonically covering the range between $p_{0}$ to $p_{1}$, as given in Eqs. (103), with $s\mu_{\Sigma}\to\frac{1}{2}\mu_{\Sigma}$. Since the above discussion applies equally to $\partial F_{q}^{>}/\partial\sigma+\partial F_{q}^{<}/\partial\sigma$, we can directly use the results in Eqs. (109), with the lower end of the energy integration domain changed to $g\sigma-s\mu_{\Sigma}\to\frac{1}{2}|\mu_{\Sigma}|-g\sigma$:

[TABLE]

with $p_{1}^{+}=\sqrt{(E+\frac{1}{2}|\mu_{\Sigma}|)^{2}-g^{2}\sigma^{2}}$.

Moving now to the critical region, we remind that $E$ decreases from $E_{0}$ to $E_{m}$ at $p=p_{m}$, set as a function of $\cos\theta>\cos\theta_{c}$ by Eq. (98), then increases back to $E_{0}$. We now express the integration of the critical region as an integral with respect to $E$ between the minimum value, $E=0$, and the maximum value, $E_{0}=\frac{1}{2}|\mu_{\Sigma}|-g\sigma$, seeking to establish the integration limits for $\cos\theta$, or equivalently, for $p$. It is not difficult to see that the momentum range is set by the two roots $p_{1}^{\pm}$ allowed by Eq. (113):

[TABLE]

which agrees with the functions $p^{s}_{1}$ in Eqs. (103), with $s=\pm 1/2$ and $\mu_{\Sigma}\to|\mu_{\Sigma}|$. We thus have

[TABLE]

with $p^{2}_{\pm}=p^{2}+(\tfrac{1}{2}|\mu_{\Sigma}|\pm g\sigma)^{2}-E^{2}$. Switching now to $y=p^{2}-(p_{1}^{-})^{2}$, we arrive at

[TABLE]

where $y_{\pm}=|\mu_{\Sigma}|(\frac{1}{2}|\mu_{\Sigma}|\pm g\sigma-E)$ and $Y=(p_{1}^{+})^{2}-(p_{1}^{-})^{2}=2|\mu_{\Sigma}|E$. Performing the integration yields

[TABLE]

We now consider the total free energy, $F_{q}=F_{q}^{+}+F_{q}^{-}$, with the notation introduced at the beginning of this subsection. It is convenient to separate this result as $F_{q}=F^{\rm reg}_{q}+F^{c;+}_{q}+F^{c;-}_{q}$, with $F^{\rm reg}_{q}=F_{q}^{<}+F_{q}^{>}+F_{q}^{-}$ having the same form as in Eq. (109a) of Case 1:

[TABLE]

where $p_{1}^{\pm}=\sqrt{(\frac{1}{2}|\mu_{\Sigma}|\pm E)^{2}-g^{2}\sigma^{2}}$.

Splitting similarly $\partial F_{q}/\partial\sigma$, we have

[TABLE]

From Eqs. (119)-(120), we again obtain the Eqs. (60) given in the main text, by shifting the energies as:

[TABLE]

In conclusion, we note that even though the two cases (case 1 in Subsec. B.2 and case 2 in Subsec. B.3) had to be treated differently, suitable manipulations of the corresponding integrals allow the final results to be expressed in a unitary manner, as shown in Eqs. (60) of the main text.

Appendix C Small-mass limit of $\delta F_{q}$

To start this section, we refer back to Eq. (64). Absorbing the prefactors within $\delta F_{q}$ as $\delta F_{q}^{\prime}=-\frac{2\pi^{2}\delta F_{q}}{\mu_{\Sigma}Tg^{2}\sigma^{2}}$ and performing a derivative with respect to $\mu_{\Sigma}$ we can rewrite Eq. (64) as

[TABLE]

where we have used the fact that $s^{2}=1/4$ is independent of $s$, and have therefore factored it out as a common prefactor on the LHS. Using Eq. (62) one can then separate $4T\frac{\partial\delta F_{q}^{\prime}}{\partial\mu_{\Sigma}}$ into two parts, that is:

[TABLE]

with

[TABLE]

In the limit of small $g\sigma$, one can analytically compute these two contributions, which, to leading order in $\mu_{\Sigma}/T$, yield the following:

[TABLE]

Here we have considered the hierarchy $g\sigma<\frac{1}{2}|\mu_{\Sigma}|<T$ and, to evaluate Eq. (125), we made use of the integral Kapusta and Gale (2006):

[TABLE]

Combining the two contributions from Eqs. (125)-(126) using Eq. (122) and integrating with respect to $\mu_{\Sigma}$ one can eventually obtain

[TABLE]

in agreement with Eq. (70).

Acknowledgements.

The authors are grateful to Dr. Nyx Shiva for useful comments on the manuscript. This work was funded by the EU’s NextGenerationEU instrument through the National Recovery and Resilience Plan of Romania - Pillar III-C9-I8, managed by the Ministry of Research, Innovation and Digitization, within the project entitled “Facets of Rotating Quark-Gluon Plasma” (FORQ), contract no. 760079/23.05.2023 code CF 103/15.11.2022.