# Linear sigma model with quarks and Polyakov loop in rotation: phase diagrams, Tolman-Ehrenfest law and mechanical properties

Pracheta Singha, Sergiu Busuioc, Victor Ambruş, Maxim Chernodub

► To cite this version:

Pracheta Singha, Sergiu Busuioc, Victor Ambruş, Maxim Chernodub. Linear sigma model with quarks and Polyakov loop in rotation: phase diagrams, Tolman-Ehrenfest law and mechanical properties. Phys.Rev.D, 2025, 112 (9), pp.094031. 10.1103/knn8-sv3k. hal-05019110

HAL Id: hal-05019110

https://hal.science/hal-05019110v1

Submitted on 20 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Linear sigma model with quarks and Polyakov loop
in rotation: phase diagrams, Tolman-Ehrenfest law
and mechanical properties

Pracheta Singha,ª Sergiu Busuioc,ª Victor E. Ambruș<fnref n="1" />ª and Maxim N. Chernodub²b,ª

ªDepartment of Physics, West University of Timișoara,
Bd. Vasile Pârvan 4, Timișoara 300223, Romania

<sup>b</sup>Institut Denis Poisson, CNRS - UMR 7013, Université de Tours, Tours 37200, France

E-mail: pracheta.singha@e-uvt.ro, sergiu.busuioc@e-uvt.ro, victor.ambrus@e-uvt.ro, maxim.chernodub@univ-tours.fr

ABSTRACT: We study the effect of rotation on the confining and chiral properties of QCD using the Polyakov-enhanced linear sigma model coupled to quarks. Working in the homogeneous approximation, we obtain the phase diagram at finite temperature, baryon density and angular frequency, taking into account the causality constraint enforced by the spectral boundary conditions at a cylindrical surface. We explicitly address various limits with respect to system size, angular frequency and chemical potential. We demonstrate that, in this model, the critical temperatures of both transitions diminish in response to the increasing rotation, being in contradiction with the first-principle lattice results. In the limit of large volume, the thermodynamics of the model is consistent with the Tolman-Ehrenfest law. We also compute the mechanical characteristics of rotating plasma such as the moment of inertia and the $K_4$ shape coefficient.

<footnote n="1">¹Corresponding author.</footnote>

<footnote n="2">²Corresponding author.</footnote>
# Contents

1 Introduction 2

2 Polyakov-loop extended linear sigma model 5

2.1 Meson Lagrangian 6

2.2 Polyakov loop Lagrangian 7

2.3 Quark Lagrangian 9

2.4 Grand canonical potential 9

2.5 Saddle point approximation 12

2.6 Thermodynamics of the model 14

2.7 Vanishing temperature limit 18

3 Rotating PLSM<sub>q</sub> model 19

3.1 Corotating frame 19

3.2 Grand canonical potential 20

3.3 Saddle point approximation 22

3.4 Rotating LSM<sub>q</sub> model with boundary 23

3.5 Vanishing temperature limit 24

4 Phase Diagram 25

4.1 Static system: Effect of boundary 27

4.2 System under rotation 29

5 Comparison to Tolman-Ehrenfest law 33

5.1 Kinetic theory formulation of the rotating system 34

5.2 Vanishing chemical potential: LSM<sub>q</sub> model 35

5.3 Vanishing chemical potential: PLSM<sub>q</sub> model 36

5.4 Finite-density matter at vanishing temperature 39

5.5 Tolman-Ehrenfest predictions vs. numerical results 40

6 Mechanical properties of rotating plasma 42

6.1 Moment of inertia 42

6.2 Shape coefficients 45

7 Conclusions 50

A The Polyakov potential 52

A.1 Polynomial Polyakov potential 52

A.2 Choice of $T_0$ 54

A.2.1 Constant $T_0$ 55

A.2.2 Chemical-potential-dependent $T_0$ 55
# 1 Introduction

The experimental observation of highly vortical quark-gluon plasmas by the RHIC collaboration has ignited a strong interest in the theoretical community to explore and understand the properties of this most vortical fluid ever created [1, 2]. The vortical motion of plasma has been addressed from the point of view of hydrodynamic and transport-based models, which were reviewed, along with the relevant experimental results, in Refs. [3, 4]. On general grounds, one can expect that the presence of vorticity in the system can also substantially affect the thermodynamic properties and the phase structure of the quark-gluon plasma, which will be the main topic of our paper.

Due to the strong non-perturbative nature of quark-gluon plasma, its dynamics cannot be reliably explored theoretically using conventional, perturbative methods. This fact leaves us with two options: we can obtain information from the first-principle approaches based on numerical Monte Carlo techniques (i.e., lattice QCD) [5–13], or we can use effective theoretical models that describe the infrared properties of QCD. Unfortunately, numerical methods often cannot provide us with the relevant information about the mechanisms underlying the non-perturbative effects. However, numerical data can serve as a valuable source of knowledge that can be used to constrain effective theoretical models, consequently allowing us to acquire information about the effects that lie behind the physical phenomena in question. So far, most of these theoretical models fail to describe the lattice data under rotation.

The rotation of the quark-gluon plasma (QGP) implies several puzzles. The earlier numerical simulations that focused on the global behavior of rotating gluon matter have revealed that the critical temperature of the deconfining phase transition is a rising function of the angular velocity (vorticity) [6, 7]. This result challenges the naive prediction, based on the Tolman-Ehrenfest (TE) law [14, 15], that rotation makes the outer layers of the rotating system hotter, and thus rotation should decrease the critical temperature of the phase transition. Analytical predictions coming from the standard effective field theories, such as the Nambu-Jona-Lasinio (NJL) model [16, 17], the linear sigma model coupled to quarks ($LSM_q$) [18] and other effective models [19, 20] appear to confirm the aforementioned TE effect [21, 22], being, thus, in contradiction with the lattice data [6, 7].

One possible explanation for the “inverse Tolman-Ehrenfest effect” revealed by lattice data is related to the melting of the gluonic condensate that alters dramatically the moment of inertia of the gluon plasma [10]. Such an effect cannot be explained by standard effective models of QCD, however qualitative agreement with the lattice data can still be captured in certain analytical approaches in which the model couplings are specially-selected functions of the rotational velocity to match the lattice data [21, 23–25] (see also Ref. [26] for evidence from perturbative QCD).

In this paper, we aim to study the phase diagram of strongly-interacting matter under rotation. Without attempting to resolve the discrepancy between model calculations and lattice QCD results, we perform this study by considering a rigidly rotating plasma in global thermal equilibrium, while rigorously enforcing the causality constraint by imposing boundary conditions, such that the local speed of corotating observers does not exceed the
speed of light. While heavy-ion collisions produce plasma with a space- and time-dependent vorticity [3], rigid rotation corresponds to a state of global thermodynamic equilibrium, thus serving as a useful probe to study the response of the quark-gluon plasma to vorticity, also in numerical (lattice) approaches [5–12].

We employ the Polyakov-enhanced [27] linear sigma model [18] coupled to quarks ($PLSM_q$) as an effective model that captures both the deconfinement and the chiral phase transitions (see Ref. [28] for a review). In this model, the potential for the mesonic self-interactions includes a small chiral symmetry breaking term that leads to a finite constitutive quark mass in the vacuum state. The Polyakov sector, modelled as a temperature-dependent effective potential, is related to a non-vanishing vacuum expectation value of the time component of the gluon four-potential. One drawback of this model is that neither the meson, nor the Polyakov potentials depend directly either on the kinematic state of the system, or on the presence of the boundary. Instead, the rotation and the boundary have an indirect effect on the meson and Polyakov loop expectation values through their coupling to fermionic loops. Although the inclusion of rotational effects either through the meson self-interactions [29] or at the level of the gluonic sector [26] may improve agreement with lattice data, we do not pursue this avenue in this work.

While it is known that rotation produces radial inhomogeneities in the quark-gluon plasma [12, 13, 30], we will use a simplified approach, treating both the chiral condensate $\sigma$ and the expectation values $L$ and $L^*$ of the Polyakov loop and its conjugate as homogeneous quantities. This approach has been adopted in various investigations of the rotating plasma (see, however, the study that incorporates inhomogeneities of rotating fields [30, 31]; as well as a detailed comparison of homogeneous and inhomogeneous approaches in a complementary setup in Ref. [22]). On the other hand, lattice data reported in Ref. [13] suggests that rotation has a sub-percent effect on the QCD plasma on the rotation axis, while drastically affecting the system phase away from the rotation axis. Despite the fact that the effective model approach does not agree with the latest lattice results, we believe that it is essential to understand down to the fine details how these models respond to both boundary and rotational effects.

The particular flavor of $PLSM_q$ model that we consider involves a logarithmic potential for the Polyakov loop, as introduced in Ref. [32] (see also Ref. [28] for a review). With respect to the polynomial ansatz [33], the logarithmic potential has the advantage of automatically limiting the values of the Polyakov loop and its conjugate between the physically-accessible values $0 \le L, L^* < 1$. Moreover, we take, for simplicity, the deconfinement temperature $T_0$ at the value 0.27 GeV, corresponding to the phase transition temperature in the $SU(3)$ pure Yang-Mills theory. In principle, this temperature should be adjusted based on the number $N_f$ of fermionic flavours and in the presence of finite quark chemical potential $\mu$ (see Ref. [27] for details). Our numerical experiments indicate that, when $T_0 \ne 0.27$ GeV, deconfinement occurs before chiral symmetry restoration, which contradicts lattice data. Moreover, the polynomial potential [33] leads to an unbounded growth of the expectation values of the Polyakov loop $L$ and its conjugate, $L^*$, which is exacerbated in the presence of rotation. Therefore, we discuss these extensions only briefly, in the appendix to this work, and instead focus predominantly on the logarithmic model
with fixed deconfinement temperature $T_0 = 0.27$ GeV.

As already mentioned, we investigate the effects of rotation on the phase diagram by considering a state of rigid rotation with constant angular velocity $\Omega$. To enforce the causality constraint of the resulting state, we enclose the system in a cylindrical boundary of radius $R$, such that $\Omega R \le 1$. On the boundary of the system, we follow Refs. [34, 35] and employ the spectral boundary conditions [36] (see Ref. [37] for the implementation within a cylindrical surface) to constrain the allowed fermionic states. This leads to a quantization of the transverse momentum by means of the roots of the Bessel functions, thus requiring numerical methods for the computation of all fermionic expectation values. The spectral boundary conditions have the advantage of preserving chiral symmetry [35], unlike the MIT bag model [38–41], which explicitly breaks this symmetry.

We performed a thorough investigation of the phase structure of strongly-interacting matter extensively over the parameter space of the model. While allowing the angular velocity $\Omega$ to vary between 0 (static system) and $\Omega R = 1$ (causal limit), we considered system sizes from small (where the fermionic contributions are suppressed) to very large (thermodynamic limit) and studied the impact of both rotation and boundary effects on the phase diagram. In the limit of small volume, where fermionic contributions are suppressed, the Polyakov sector is dominant and the deconfinement transition proceeds as for the unbounded, pure-glue system. For large system sizes, it is expected that the thermodynamics of the system should be compatible with the classical Tolman-Ehrenfest (TE) law.

We formulate a quantitative implementation of the famous TE law by considering the thermodynamic limit of large volume (large cylinder radius $R$). To retain the causality constraint, we keep $|\Omega R| \le 1$ fixed, thus leading to decreasing $\Omega$ as $R$ is increased. Therefore, while the boundary effects become negligible, so do quantum effects, which are known to be proportional to $\Omega^2$ [42, 43]. Using the local thermal equilibrium approximation for the density operator, introduced in Refs. [42, 44, 45], we consider the angular velocity dependence of the pseudicritical chiral restoration temperature $T_{pc}$ at vanishing chemical potential, as well as the critical chemical potential $\mu_c$ at vanishing temperature.

In the unbounded, non-rotating system, the deconfinement and chiral restoration temperatures are approximately the same both in QCD, and in the PLSM$_q$ model that we consider. Recently, Ref. [46] reported that the chiral and deconfining transition split under the influence of rotation. On the other hand, the existing numerical calculations imply that the window for this effect in the phase diagram is rather narrow [11]. Our recent effective model calculation found that the splitting is promoted by the finite-volume effects only and is absent for large enough volumes [47]. We now revisit this intriguing phenomenon, by considering both the limit of small volume (discussed above), as well as that of fast rotation. In the latter case, the TE analogy provides unexpected evidence for a rotation-induced split occurring in the thermodynamic limit $R \to \infty$, as the causal limit is approached ($\Omega R \to 1$). Notice that, in the thermodynamic limit, the fastest possible rotation is necessarily slow in terms of its physical amplitude, as $\Omega \to 0$ when $R \to \infty$ in order to keep the product $\Omega R$ finite.

Another intriguing feature of the rotating quark-gluon plasma is related to its mechanical properties, namely to the moment of inertia $I$. Lattice calculations [10, 11] indicate
that $I$ becomes negative in a so-called “supervortical” range of temperatures. We therefore compute this quantity within the PLSM$_q$ model and compare with the predictions of the TE law. Furthermore, we consider the shape coefficients $K_{2n}$, in particular, $K_2$ – related to the moment of inertia at vanishing angular velocity – and $K_4$ – corresponding to shape deformations induced by rotation (see Refs. [10, 48]). In the limit of both $T \to \infty$ (when the system becomes asymptotically free) and $R \to \infty$ (when boundary effects become negligible), these coefficients approach the values of 2! and 4!, respectively, as expected [48, 49] for a free (non-interacting) conformal gas. At intermediate temperatures, we study the behavior of these coefficients in the context of the TE law.

The structure of the paper is as follows. In Section 2, we review the PLSM$_q$ in the absence of rotation, in the thermodynamic limit, thereby introducing our notation. We describe the meson, quark and Polyakov loop sectors of the model, the construction of the grand canonical potential, the mean-field approximation and the system thermodynamics. We also discuss the analytical and numerical strategy for determining the first-order, second-order or cross-over transitions corresponding to deconfinement or chiral restoration. The vanishing-temperature limit is particularly discussed, as in this limit, PLSM$_q$ reduces to the LSM$_q$ without the Polyakov loop variable.

In Section 3, we construct the bounded PLSM$_q$ system under rotation, satisfying spectral boundary conditions on the cylindrical boundary.

Section 4 starts with the phase diagram of the bounded, non-rotating PLSM$_q$ model in Subsec. 4.1, where we treat the limits of small, intermediate and large system sizes. Subsection 4.2 discusses the effect of rotation on the phase diagram of the system.

Section 5 presents the treatment of the rotating system based on the Tolman-Ehrenfest law. We consider the angular velocity dependence of the pseudocritical temperature at vanishing chemical potential in both the LSM$_q$ (Subsec. 5.2) and the PLSM$_q$ (Subsec. 5.3) models, as well as of the critical chemical potential at vanishing temperature (Subsec. 5.4).

The mechanical properties (moment of inertia, shape coefficients) of the rotating plasma are studied in Section 6, where we also discuss the comparison with the TE prediction.

Our conclusions are presented in Sec. 7. The paper also includes Appendix A, which discusses subtleties of the choice of the potential for the Polyakov loop, including the shape of the potential (logarithmic vs. polynomial), as well as the dependence of the deconfinement temperature on both the number of quark flavors and on the quark chemical potential.

## 2 Polyakov-loop extended linear sigma model

In this section, we review the basic features of the Polyakov-loop-extended linear sigma model coupled to quarks in a uniform, non-rotating background. The Lagrangian of the model reads as follows:

$$ \mathcal{L} = \mathcal{L}_M + \mathcal{L}_L + \mathcal{L}_q. \qquad (2.1) $$

We address the meson Lagrangian, $\mathcal{L}_M$, in Subsec. 2.1. The Polyakov Lagrangian $\mathcal{L}_L(L)$ is described in Subsec. 2.2. Finally, the quark Lagrangian $\mathcal{L}_q$, including the quark-meson and quark-gluon interaction terms, is discussed in Subsec. 2.3. Subsection 2.4 presents the
computation of the free energy. The saddle point equations giving the mean-field values of the meson fields and of the Polyakov loops $L$ and $L^*$ are discussed in Subsec. 2.5. Subsection 2.6 discusses the thermodynamics of the model, while its vanishing temperature limit is presented in Subsec. 2.7.

## 2.1 Meson Lagrangian

The mesonic Lagrangian $\mathcal{L}_M = T_M - V_M$ is composed of kinetic and potential parts,

$$ T_M = \frac{1}{2} (\partial_\mu \sigma \partial^\mu \sigma + \partial_\mu \vec{\pi} \cdot \partial^\mu \vec{\pi}), \quad V_M = \frac{\lambda}{4} (\sigma^2 + \vec{\pi}^2 - v^2)^2 - h\sigma, \qquad (2.2) $$

where $\sigma$ and $\vec{\pi} = (\pi^1, \pi^2, \pi^3)$ are the scalar singlet and pseudoscalar triplet meson fields, respectively.

In this paper, we work in the mean-field approximation, by which the meson fields are replaced by their (classical) expectation values, thus ignoring their fluctuating (quantum) contributions. The mean-field values $\sigma$ and $\vec{\pi}$ minimize the grand potential, thereby ensuring the thermodynamic stability of the solution. In this paper, we also ignore space-time gradients of the meson fields, such that the kinetic part $T_M$ becomes thermodynamically irrelevant. The effects of the spatial variations of the sigma field, albeit in a simpler model without the Polyakov loop, have been explored in the recent study of Ref. [22]. However, for simplicity, we restrict ourselves to the homogeneous approximation.

Denoting by $\sigma_0$ and $\vec{\pi}_0$ the mean-field vacuum expectation values of the mesons, the minimization of $V_M$ gives the following equations:

$$ \lambda \vec{\pi}_0 (\sigma_0^2 + \vec{\pi}_0^2 - v^2) = 0, \quad \lambda \sigma_0 (\sigma_0^2 + \vec{\pi}_0^2 - v^2) - h = 0. \qquad (2.3) $$

The symmetry-breaking term $h$ in the second equation leads to a vanishing pion condensate, $\vec{\pi}_0 = 0$, while the scalar condensate $\sigma_0$ is obtained as a solution of the following cubic equation:

$$ \lambda \sigma_0 (\sigma_0^2 - v^2) - h = 0. \qquad (2.4) $$

In principle, Eqs. (2.3) receive corrections from fermion loops, which renormalize the coupling constants. Such corrections can be treated self-consistently beyond the mean-field approximation [50], leading to the same vacuum structure. For brevity, we do not discuss these details in the present paper.

<table><caption>**Table 1.** Parameters of the meson Lagrangian $\mathcal{L}_M$ in Eq. (2.2).</caption><thead><tr><th colspan="4">Input parameters</th></tr></thead><tbody><tr><td>σ<sub>0</sub> = f<sub>π</sub></td><td>m<sub>π</sub> = √λ(σ<sub>0</sub><sup>2</sup> - v<sup>2</sup>)</td><td>m<sub>σ</sub> = √λ(3σ<sub>0</sub><sup>2</sup> - v<sup>2</sup>)</td><td>m<sub>q</sub> = gσ<sub>0</sub></td></tr><tr><td>0.093 GeV</td><td>0.138 GeV</td><td>0.600 GeV</td><td>0.307 GeV</td></tr></tbody><thead><tr><th colspan="4">Model parameters</th></tr></thead><tbody><tr><td>h = f<sub>π</sub>m<sub>π</sub><sup>2</sup></td><td>λ = (m<sub>σ</sub><sup>2</sup> - m<sub>π</sub><sup>2</sup>)/(2f<sub>π</sub><sup>2</sup>)</td><td>v<sup>2</sup> = f<sub>π</sub><sup>2</sup> - m<sub>π</sub><sup>2</sup>/λ</td><td>g = m<sub>q</sub>/σ<sub>0</sub></td></tr><tr><td>0.00177 (GeV)<sup>3</sup></td><td>19.7</td><td>0.00768 (GeV)<sup>2</sup></td><td>3.3</td></tr></tbody></table>
<table><caption>**Table 2.** Parameters for the Polyakov Lagrangian given in Eq. (2.7).</caption><thead><tr><th colspan="4">Polyakov potential parameters</th></tr><tr><th>a<sub>0</sub></th><th>a<sub>1</sub></th><th>a<sub>2</sub></th><th>b<sub>3</sub></th></tr></thead><tbody><tr><td>3.51</td><td>-2.47</td><td>15.2</td><td>-1.75</td></tr></tbody></table>

The model parameters $\lambda$, $v$ and $h$ are fixed by the vacuum properties of the model, namely the expectation value $\sigma_0$ of the sigma meson and the masses of the $\sigma$ and $\pi$ mesons:

$$ \sigma_0 = f_{\pi}, \quad m_{\sigma}^2 = \frac{\partial^2 V_M}{\partial \sigma^2} = \lambda(3\sigma_0^2 - v^2), \quad m_{\pi}^2 = \frac{\partial^2 V_M}{\partial \pi^2} = \lambda(\sigma_0^2 - v^2), \quad (2.5) $$

where $f_{\pi} = 0.093$ GeV is the pion decay constant, while $m_{\sigma} = 0.600$ GeV is the mass of the sigma meson [51] and $m_{\pi} = 0.138$ GeV is the pion mass [52]. The latter two relations in Eq. (2.5) fix the parameters $\lambda$ and $v$, while the explicit symmetry breaking term $h$ is determined by Eq. (2.3), as follows:

$$ \lambda = \frac{m_{\sigma}^2 - m_{\pi}^2}{2f_{\pi}^2}, \quad v^2 = f_{\pi}^2 - \frac{m_{\pi}^2}{\lambda}, \quad h = \lambda\sigma_0(\sigma_0^2 + \pi_0^2 - v^2) = f_{\pi}m_{\pi}^2. \quad (2.6) $$

The values of the model parameters are summarized in Table 1.

## 2.2 Polyakov loop Lagrangian

The Polyakov Lagrangian $\mathcal{L}_L$ models the (non-linear) thermal properties of the gluon sector of QCD and contains only the potential term,

$$ \frac{\mathcal{L}_L}{T^4} = -\frac{V_L}{T^4} = \frac{1}{2}a(T)L^*L - b(T)\ln\mathcal{H}, \quad \mathcal{H} = 1 - 6L^*L + 4(L^{*3} + L^3) - 3(L^*L)^2, \quad (2.7) $$

where $\mathcal{H}$ is the SU(3) Haar measure associated with the gluonic path integral, while $a(T)$ and $b(T)$ are temperature-dependent functions, parametrized as

$$ a(T) = a_0 + a_1 \left(\frac{T_0}{T}\right) + a_2 \left(\frac{T_0}{T}\right)^2, \quad b(T) = b_3 \left(\frac{T_0}{T}\right)^3. \quad (2.8) $$

The values of the constant parameters $a_0$, $a_1$, $a_2$ and $b_3$ are summarized in Table 2. For simplicity, we set $T_0 = 0.27$ GeV, corresponding to the transition temperature of the pure-glue SU(3) theory. In Eq. (2.7), we considered the logarithmic model for the Polyakov potential, which will be used in the main text of this work. We dedicate Appendix A to a discussion of other models for the Polyakov potential, as well as of different values and models for $T_0$.

The degrees of freedom appearing in $\mathcal{L}_L$ are the color-traced Polyakov loops, $L$ and $L^*$, defined as:

$$ L = \frac{1}{3} \mathrm{Tr}_c \Phi(x), \quad \Phi(x) = \mathcal{P} \exp \left[ -i \int_0^\beta d\tau \mathcal{A}_4(x, \tau) \right], \quad (2.9) $$

where $\mathcal{A}_4 \equiv \mathcal{A}_{\tau} = i\mathcal{A}_0$ is the temporal component of the gluon field, expressed with respect to the imaginary time. In the mean-field approximation, $\mathcal{A}_4$ is a constant, classical field.
In the Polyakov gauge [53], $\mathcal{A}_4 = \frac{1}{2}\mathcal{A}_4^{(3)}\lambda_3 + \frac{1}{2}\mathcal{A}_4^{(8)}\lambda_8$ is diagonal, with $\lambda_3 = \text{diag}(1, -1, 0)$ and $\lambda_8 = \frac{1}{\sqrt{3}}\text{diag}(1, 1, -2)$ being the diagonal Gell-Mann matrices representing the Abelian (Cartan) subalgebra of the SU(3) Lie algebra. The diagonal entries of $\mathcal{A}_4$, corresponding to its eigenvalues, can be taken as

$$ \mathcal{A}_4 = \text{diag}(\phi, \phi', -\phi - \phi'), \qquad (2.10) $$

where the scalar parameters $\phi$ and $\phi'$ are related to the vector components of the gluonic fields $\mathcal{A}_4^{(3)}$ and $\mathcal{A}_4^{(8)}$ via

$$ \phi = \frac{1}{2}\mathcal{A}_4^{(3)} + \frac{1}{2\sqrt{3}}\mathcal{A}_4^{(8)}, \quad \phi' = -\frac{1}{2}\mathcal{A}_4^{(3)} + \frac{1}{2\sqrt{3}}\mathcal{A}_4^{(8)}, \quad -\phi - \phi' = -\frac{1}{\sqrt{3}}\mathcal{A}_4^{(8)}. \qquad (2.11) $$

With the above notations, the Polyakov loops $L$ and $L^*$ are

$$ L = \frac{1}{3}(e^{-i\beta\phi} + e^{-i\beta\phi'} + e^{i\beta(\phi+\phi')}), \quad L^* = \frac{1}{3}(e^{i\beta\phi} + e^{i\beta\phi'} + e^{-i\beta(\phi+\phi')}), \qquad (2.12) $$

In the effective model description, $L$ and $L^*$ are treated as two independent, real fields, chosen by imposing the minimization of the grand potential. In the absence of the quark degrees of freedom, $L$ and $L^*$ satisfy the following equations:

$$ \frac{1}{T^4} \frac{\partial V_L}{\partial L} = -\frac{1}{2} a L^* - \frac{6b}{\mathcal{H}} (L^* - 2L^2 + LL^{*2}) = 0, \qquad (2.13a) $$

$$ \frac{1}{T^4} \frac{\partial V_L}{\partial L^*} = -\frac{1}{2} a L - \frac{6b}{\mathcal{H}} (L - 2L^{*2} + L^* L^2) = 0. \qquad (2.13b) $$

Imposing $L = L^*$, the above equations allow for three solutions: $L = L^* = 0$ and

$$ L^{\pm} = \frac{1}{3} \pm \frac{2}{3} \sqrt{1 + \frac{9b(T)}{a(T)}}. \qquad (2.14) $$

In general, $L = 0$ and $L = L^+$ correspond to the minima of the effective potential, while $L = L^-$ gives a local maximum. At small $T$, the global minimum is attained for $L = L^* = 0$, corresponding to the confined phase, while $L^\pm$ are complex-valued quantities, corresponding to the deconfinement phase. When $T \to \infty$, the system becomes deconfined and $L, L^*$ asymptotically approach the value 1, as follows:

$$ L, L^* \to 1 + \frac{3b_3}{a_0} \left(\frac{T_0}{T}\right)^3 + O(T^{-4}) \quad [\text{as } T \to \infty]. \qquad (2.15) $$

Imposing the Stefan-Boltzmann limit for the gluon pressure, $P_g = -V_L = \mathcal{L}_L$, fixes the $a_0$ parameter:

$$ \lim_{T \to \infty} P_g = \frac{1}{2} a_0 T^4 = \nu_g \frac{\pi^2 T^4}{90}, \qquad (2.16) $$

where $\nu_g = 2(N_c^2 - 1) = 16$ is the gluon degeneracy factor.

The model also consistently imposes a first-order phase transition when $T = T_0$ that takes place in a pure glue SU(3) gauge theory. This property is achieved by demanding that $V_L = 0$ when $T = T_0$ and $L = L^+$, such that

$$ \ln[(1 - L^+)^3(1 + 3L^+)] = \frac{a_c L^{+2}}{2b_3}, \quad a_c = a(T_0) = a_0 + a_1 + a_2. \qquad (2.17) $$
Using the relation $L^+ = \frac{1}{3} + \frac{2}{3}\sqrt{1 + 9b_3/a_c}$, we get $b_3 = -0.107817(a_0 + a_1 + a_2)$. This fixes $b_3$ once $a_1$ and $a_2$ are known. These latter two parameters are found by fitting $P_g$ to the lattice data. At the first-order phase transition, the Polyakov loop jumps to the value $L^{+,c} \simeq 0.448$.

## 2.3 Quark Lagrangian

The quark Lagrangian $\mathcal{L}_q$ reads

$$ \mathcal{L}_q = \bar{\psi} [iD\!\!\!/ + \mu\gamma^0 - g(\sigma + i\gamma_5\vec{\tau} \cdot \vec{\pi})] \psi, \qquad (2.18) $$

where the Feynman slash denotes contraction with the Dirac matrices, $A = \gamma^\mu A_\mu$, $D_\mu = \partial_\mu - iA_\mu$ is the covariant derivative, $A_\mu = \frac{1}{2}A_\mu^a t_a$ is the gluon potential and $\mu$ is the quark chemical potential. The strength $g$ of the coupling to the meson fields is determined by fixing the effective mass of the fermions in the vacuum to $M_q = 0.307$ GeV:

$$ g = \frac{M_q}{\sigma_0} = \frac{M_q}{f_\pi} \simeq 3.3, \qquad (2.19) $$

as shown in Table 1.

In this work, we consider the $N_f = 2$ flavor model, such that $\psi = (u,d)^T$. Each of the components of the field $\psi$ is a spinor that also carries an $SU(3)$ color index $a = 1,2,3$ (not explicitly shown). The gamma matrices $\gamma^\mu = (\gamma^0, \gamma)$ are taken in the Dirac-Pauli representation,

$$ \gamma^0 = \begin{pmatrix} I_2 & 0 \\ 0 & -I_2 \end{pmatrix}, \qquad \gamma = \begin{pmatrix} 0 & \boldsymbol{\sigma} \\ -\boldsymbol{\sigma} & 0 \end{pmatrix}, \qquad (2.20) $$

where $\boldsymbol{\sigma}$ are the Pauli matrices,

$$ \sigma^1 = \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix}, \quad \sigma^2 = \begin{pmatrix} 0 & -i \\ i & 0 \end{pmatrix}, \quad \sigma^3 = \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix}. \qquad (2.21) $$

## 2.4 Grand canonical potential

The partition function of PLSM$_q$ follows from its Lagrangian (2.1)

$$ \mathcal{Z} = e^{-\beta V_{3d}(V_M + V_L)} \mathcal{Z}_q, \quad \mathcal{Z}_q = \int [i d\psi^\dagger] [d\psi] e^{-S_E}, \quad S_E = \int_0^\beta d\tau \int_V d^3 x \mathcal{L}_E, \\ \mathcal{L}_E = -\bar{\psi} [\gamma^0 (-\partial_\tau - iA_4 + \mu) + i\gamma \cdot \nabla - g\sigma] \psi, \qquad (2.22) $$

where we omitted an irrelevant overall normalization factor. The kinetic parts of the meson and gluon sectors were neglected, while the potentials $V_M$ and $V_L$ are given in Eqs. (2.2) and (2.7), respectively. The volume-averaged grand canonical potential of the model, $F = -(T \ln \mathcal{Z})/V$, is given by

$$ F = V_M + V_L + F_q, \quad F_q = -\frac{T}{V} \ln \mathcal{Z}_q. \qquad (2.23) $$
Note that $F$ has units of [energy]/[length]$^3$.

The fermionic part $\mathcal{Z}_q$ of the partition function can be evaluated by considering a decomposition $\psi(X) = \sum_j \hat{b}_j \psi_j(X)$ of the Dirac spinor with respect to the plane-wave basis $\psi_j(X)$, satisfying the following eigenvalue equations:

$$(-i\partial_{\tau}, \sigma_f^3, A_4, h, \mathbf{P})\psi_j = (w_j, f_j, A_{4:j}, \lambda_j, \mathbf{p}_j)\psi_j, \quad (2.24)$$

with $w_j = 2\pi T(n_j + \frac{1}{2})$ the Matsubara frequencies; the flavor $f_j = \pm 1$ for the u/d quarks; $A_{4:j} \equiv A_{4:j_c}$ the eigenvalues of the gluon field shown in Eq. (2.10), $(A_{4:1}, A_{4:2}, A_{4:3}) = (\phi, \phi', -\phi - \phi')$; $\lambda_j = \pm 1/2$ the eigenvalues of the helicity operator $h = \mathbf{S} \cdot \mathbf{P}/|\mathbf{P}|$; and $\mathbf{P} = -i\nabla$ the momentum operator. In addition, we require $\psi_j$ to be eigenmodes of the Hamiltonian operator, $H = iA_4 - \mu + 2\gamma^5 \mathbf{S} \cdot \mathbf{P} + g\sigma\gamma^0$:

$$H \psi_j = \zeta_j \mathcal{E}_j \psi_j, \quad \mathcal{E}_j = E_j - \zeta_j(\mu - iA_{4:j}), \quad E_j = \sqrt{\mathbf{p}_j^2 + g^2\sigma^2}, \quad (2.25)$$

where $\zeta_j = \pm 1$ distinguishes between the particle and anti-particle branches. In the above discussions, we used $j \equiv (\zeta_j, w_j, f_j, c_j, \lambda_j, \mathbf{p}_j)$ to collectively denote the eigenvalues characterizing the mode $\psi_j$.

The modes described above have the following form:

$$\psi_j(X) = \frac{e^{iw_j\tau + i\mathbf{p}_j \cdot \mathbf{x}}}{(2\pi)^3} \begin{pmatrix} \frac{1}{2}(1+f_j) \\ \frac{1}{2}(1-f_j) \end{pmatrix} \otimes \begin{pmatrix} \delta_{c_j,1} \\ \delta_{c_j,2} \\ \delta_{c_j,3} \end{pmatrix} \otimes \frac{1}{\sqrt{2}} \begin{pmatrix} \mathfrak{E}_j \\ 2\zeta_j \lambda_j \mathfrak{E}_j^\zeta \end{pmatrix} \otimes \xi_j, \quad (2.26)$$

where $\mathfrak{E}_j = \sqrt{1 + \zeta_j \frac{g\sigma}{E_j}}$ and $\xi_j \equiv \xi_{\lambda_j}(\mathbf{p}_j)$ satisfying $\boldsymbol{\sigma} \cdot \mathbf{p} \xi_j(\mathbf{p}) = 2\lambda_j \xi_j(\mathbf{p})$ given by

$$\xi_{1/2}(\mathbf{p}) = \begin{pmatrix} e^{-\frac{i}{2}\varphi_p \cos \frac{\theta_p}{2}} \\ e^{\frac{i}{2}\varphi_p \sin \frac{\theta_p}{2}} \end{pmatrix}, \quad \xi_{-1/2}(\mathbf{p}) = \begin{pmatrix} -e^{-\frac{i}{2}\varphi_p \sin \frac{\theta_p}{2}} \\ e^{\frac{i}{2}\varphi_p \cos \frac{\theta_p}{2}} \end{pmatrix} \quad (2.27)$$

It can be easily checked that the modes are orthonormal:

$$\int_0^\beta d\tau \int d^3x \psi_j^\dagger \psi_{j'} = \beta\delta(j, j'), \quad \delta(j, j') = \delta_{\zeta_j, \zeta_{j'}} \delta_{n_j, n_{j'}} \delta_{f_j, f_{j'}} \delta_{c_j, c_{j'}} \delta_{\lambda_j, \lambda_{j'}} \delta^3(\mathbf{p}_j - \mathbf{p}_{j'}). \quad (2.28)$$

Considering now that the spatial volume $V$ is a fictitious rectangular box and the momentum modes are discretized via $p_i = 2\pi n_i/L_i$, with $n_i$ being integer numbers and $L_i$ being the length of side $i$, the fermionic partition function can be written as

$$\mathcal{Z}_q = \int \prod_j (i d\hat{b}_j^\dagger d\hat{b}_j) e^{-S}, \quad S = -\sum_{j,j'} i \hat{b}_j^\dagger D(j, j') \hat{b}_{j'}, \quad (2.29)$$

with $D(j, j')$ being a diagonal matrix:

$$D(j, j') = \beta(-w_j + i\zeta_j \mathcal{E}_j)\delta(j, j'). \quad (2.30)$$

Using the integration formula $\int \prod_j [i d\hat{b}_j^\dagger d\hat{b}_j] e^{i\hat{b}^\dagger \cdot D \cdot \hat{b}} = \det D$ over the Grassmann variables $i\hat{b}_j^\dagger$ and $\hat{b}_j$, the partition function becomes simply

$$\mathcal{Z}_q = \prod_j \left[ \int i d\hat{b}_j^\dagger d\hat{b}_j e^{i \hat{b}_j^\dagger [\beta(i \zeta_j \mathcal{E}_j - w_j)] \hat{b}_j} \right] = \prod_j [\beta(i \zeta_j \mathcal{E}_j - w_j)]. \quad (2.31)$$
The fermion contribution to the volume-averaged grand potential, $F_q = -\frac{T}{V} \ln \mathcal{Z}_q$, becomes

$$F_q = -2N_f T \sum_{n=-\infty}^{\infty} \sum_{c=1}^{3} \sum_{\zeta=\pm 1} \int \frac{d^3 p}{(2\pi)^3} \ln[\beta(i\zeta\mathcal{E}_\zeta^c - w_n)], \quad (2.32)$$

where now we introduced the notation

$$\mathcal{E}_\zeta^c = \mathcal{E}_\zeta + i\zeta A_{4;c}, \quad \mathcal{E}_\zeta = E - \zeta\mu, \quad E = \sqrt{\mathbf{p}^2 + g^2\sigma^2}. \quad (2.33)$$

Averaging over positive and negative $n$ leads to

$$F_q = -2N_f T \sum_{n=-\infty}^{\infty} \sum_{c=1}^{3} \sum_{\zeta=\pm 1} \int \frac{d^3 p}{(2\pi)^3} \frac{1}{2} \ln[\beta^2 w_n^2 + (\beta \mathcal{E}_\zeta^c)^2]. \quad (2.34)$$

The summation over $n$ can be performed using the method described in Ref. [54], leading to:

$$\frac{1}{2} \sum_{n=-\infty}^{\infty} \ln[(2n+1)^2\pi^2 + (\beta \mathcal{E}_{\zeta}^c)^2] = \frac{1}{2} \beta \mathcal{E}_{\zeta}^c + \ln(1 + e^{-\beta \mathcal{E}_{\zeta}^c}) + \mathcal{N}, \quad (2.35)$$

where $\mathcal{N} = \frac{1}{2} - \ln(1+e) + \frac{1}{2} \sum_{n=-\infty}^{\infty} 1$ is an (infinite) constant which is independent of the thermodynamic parameters of the state and can therefore be discarded. The two terms above give rise to the zero-point (vacuum) energy and the thermal contribution, such that $F_q = F_q^{\text{vac}} + F_q^{\beta}$, with [55, 56]:

$$F_q^{\text{vac}} = -2N_c N_f \int \frac{d^3 p}{(2\pi)^3} E, \quad F_q^{\beta} = -2N_f T \sum_{\zeta=\pm 1} \int \frac{d^3 p}{(2\pi)^3} \sum_c \ln \left(1 + e^{-\beta \mathcal{E}_\zeta^c}\right). \quad (2.36)$$

In the vacuum term, we took into account that $\sum_{\zeta,c} \mathcal{E}_\zeta^c = 2N_c E$, since $\sum_c A_{4;c} = \phi + \phi' - (\phi + \phi') = 0$. In the linear sigma model, the (regularized) vacuum contribution can be employed to renormalize the parameters of the meson Lagrangian and can therefore be safely neglected. In this work, we follow Ref. [55] and ignore such vacuum contributions, therefore tacitly assuming that $F_q^{\text{vac}} = 0$. Since now $F_q = F_q^{\beta}$, we will omit the $\beta$ superscript and instead identify the fermionic free energy with its thermal contribution.

The sum over the colour index can be performed as follows:

$$F_{\zeta} \equiv \frac{1}{3} \sum_{c=1}^{3} \ln(1 + e^{-\beta \mathcal{E}_{\zeta}^{c}}) = \frac{1}{3} \ln(1 + 3L_{\zeta} e^{-\beta \mathcal{E}_{\zeta}} + 3L_{-\zeta} e^{-2\beta \mathcal{E}_{\zeta}} + e^{-3\beta \mathcal{E}_{\zeta}}), \quad (2.37)$$

where we introduced the notation $L_{\pm} = \frac{1}{3} \sum_c e^{\mp i\beta A_{4;c}}$, with $L_{+} \equiv L$ and $L_{-} \equiv L^*$. Finally, the fermion free energy becomes

$$F_q = -2N_f N_c T \sum_{\zeta=\pm 1} \int \frac{d^3 p}{(2\pi)^3} F_{\zeta}. \quad (2.38)$$
## 2.5 Saddle point approximation

The thermodynamic and phase space structure of the model are guided by the mean-field values: $\Phi = \Phi_{\text{mf}}$, $\sigma = \sigma_{\text{mf}}$ and $\pi = \pi_{\text{mf}} = 0$. We find these mean-field values by imposing the minimization of the total free energy at thermodynamic equilibrium by solving the following equations:

$$ \frac{\partial F}{\partial \sigma} = \frac{\partial F}{\partial L} = \frac{\partial F}{\partial L^*} = 0, \qquad (2.39) $$

where we remember that within the scope of the Polyakov loop-enhanced models, the two conjugate fields $L$ and $L^*$ are treated as real but independent numbers.

The derivatives appearing in Eq. (2.39) can be expressed with the help of several distribution functions, which we define below:

$$ f_{\pm} \equiv -\frac{\partial F_{\pm}}{\partial(\beta\mathcal{E}_{\pm})} = \frac{L_{\pm}e^{-\beta\mathcal{E}_{\pm}} + 2L_{\mp}e^{-2\beta\mathcal{E}_{\pm}} + e^{-3\beta\mathcal{E}_{\pm}}}{1 + 3L_{\pm}e^{-\beta\mathcal{E}_{\pm}} + 3L_{\mp}e^{-2\beta\mathcal{E}_{\pm}} + e^{-3\beta\mathcal{E}_{\pm}}}, \qquad (2.40a) $$

$$ f_{\pm, \pm} \equiv \frac{\partial F_{\pm}}{\partial L_{\pm}} = \frac{e^{-\beta\mathcal{E}_{\pm}}}{1 + 3L_{\pm}e^{-\beta\mathcal{E}_{\pm}} + 3L_{\mp}e^{-2\beta\mathcal{E}_{\pm}} + e^{-3\beta\mathcal{E}_{\pm}}}, \qquad (2.40b) $$

$$ f_{\pm, \mp} \equiv \frac{\partial F_{\pm}}{\partial L_{\mp}} = \frac{e^{-2\beta\mathcal{E}_{\pm}}}{1 + 3L_{\pm}e^{-\beta\mathcal{E}_{\pm}} + 3L_{\mp}e^{-2\beta\mathcal{E}_{\pm}} + e^{-3\beta\mathcal{E}_{\pm}}}, \qquad (2.40c) $$

where the upper (lower) sign corresponds to the (anti)particle sector. With the above definitions in hand, we obtain

$$ \frac{\partial F}{\partial \sigma} = \frac{\partial V_M}{\partial \sigma} + g\langle\bar{\psi}\psi\rangle, \qquad \frac{\partial F}{\partial L_{\pm}} = \frac{\partial V_L}{\partial L_{\pm}} - 2N_f N_c T \sum_{\zeta=\pm 1} \int \frac{d^3 p}{(2\pi)^3} f_{\zeta, \pm}, \qquad (2.41) $$

where we identified the fermion condensate in the first equation,

$$ \langle \bar{\psi} \psi \rangle = 2N_f N_c g \sigma \sum_{\zeta=\pm 1} \int \frac{d^3 p}{(2\pi)^3 E} f_{\zeta}, \qquad (2.42) $$

while the derivative of the mesonic and of the Polyakov loop potentials are given by

$$ \frac{\partial V_M}{\partial \sigma} = \lambda(\sigma^2 - v^2)\sigma - h, \qquad (2.43a) $$

$$ \frac{\partial V_L}{\partial L_{\pm}} = -\frac{6T^4 b(T)}{\mathcal{H}}(L_{\mp} - 2L_{\pm}^2 + L_{\mp}^2 L_{\pm}) - \frac{T^4}{2}a(T)L_{\mp}. \qquad (2.43b) $$

where $\mathcal{H} = 1 - 6L^*L + 4(L^{*3} + L^3) - 3(L^*L)^2$ was introduced in Eq. (2.7).

As the terms coming from Eqs. (2.43) involve cubic powers of the order parameters, the saddle-point equations (2.39) will typically give multiple solutions. As a rule of thumb, the physically acceptable solutions must satisfy the constraints:

$$ 0 < \sigma \le f_{\pi}, \qquad 0 \le L, L^* < 1. \qquad (2.44) $$

The first constraint arises on physical grounds as we expect that both thermal fluctuations and finite-density environment restore the chiral symmetry ($\sigma \to 0$) instead of enhancing its
breaking ($\sigma \to f_\pi$). Similarly, in the vacuum, the Polyakov loop expectation value vanishes. As the medium becomes more energetic (denser, at high $\mu$; or hotter, at large $T$), the $\sigma$ condensate melts and the system becomes deconfined, i.e. $L, L^* \to 1$. Note that $L$ and $L^*$ cannot exceed 1 as a result of the mathematical definition of the Polyakov loop (2.12). Values outside the intervals indicated in Eq. (2.44) cannot represent physically meaningful systems.

When Eqs. (2.39) allow for multiple solutions satisfying Eqs. (2.44), one must distinguish between unstable and (meta)stable solutions. This selection can be done at the level of the volume-averaged grand potential, $F = V_M + V_L + F_q$, where the fermionic contribution can be written in terms of the pressure $P_q = -\frac{i}{6}\bar{\psi}\gamma\cdot\nabla\psi$,

$$F_q = -P_q, \quad P_q = 2N_f N_c \sum_{\zeta=\pm 1} \int \frac{d^3 p}{(2\pi)^3} \frac{p^2}{3E} f_\zeta. \qquad (2.45)$$

The unstable solutions correspond to points of local maxima or saddle points. Points of local (but not lowest) minima correspond to metastable states. The global (lowest-free-energy) minimum represents the true stable ground state of the system. Under quantum and thermal fluctuations, it is reasonable to expect that the quantum system will always find its way from any of the metastable states to the true ground state. Typical examples of the dependence of the total grand potential, $F = V_M + V_L + F_q$, on $\sigma$ can be found in Fig. 1 of Ref. [47].

In order to find the solutions of the saddle point equations, we write Eqs. (2.39) in a compact form as

$$\frac{\partial F}{\partial p_i} = 0, \quad p_i \in \{\sigma, L, L^*\}. \qquad (2.46)$$

Equation (2.46) is satisfied when the control parameters take their mean-field values, $p_i \to p_i^{\text{mf}}$. Considering now some values $\tilde{p}_i$ of the control parameters, which are in a small vicinity of their mean-field values, we can set up the Newton-Raphson algorithm by expanding $(\partial F / \partial p_i)$ in a Taylor series around the mean-field set of values $p_i^{\text{mf}}$:

$$\left. \frac{\partial F}{\partial p_i} \right|_{p=p_{\text{mf}}} = \left. \frac{\partial F}{\partial p_i} \right|_{p=\tilde{p}} + \left. \frac{\partial^2 F}{\partial p_i \partial p_j} \right|_{p=\tilde{p}} (p_j^{\text{mf}} - \tilde{p}_j) + \dots \qquad (2.47)$$

We can find $p_i^{\text{mf}}$ by solving

$$M_{ij}(\tilde{p})(p_j^{\text{mf}} - \tilde{p}_j) = -\left.\frac{\partial F}{\partial \tilde{p}_i}\right|_{\tilde{p}_i}, \quad M_{ij}(\tilde{p}) = \left.\frac{\partial^2 F}{\partial p_i \partial p_j}\right|_{\tilde{p}_i}, \qquad (2.48)$$

where $M_{ij}$ is the 3 × 3 Jacobian matrix, which we split into the potential and fermionic parts as follows:

$$M_{ij} = M_{ij}^{V} + M_{ij}^{q}, \quad M_{ij}^{V} = \frac{\partial^2 (V_M + V_L)}{\partial p_i \partial p_j}, \quad M_{ij}^{q} = \frac{\partial^2 F_q}{\partial p_i \partial p_j}. \qquad (2.49)$$

The contributions $M_{\sigma\sigma}^V$ and $M_{\sigma\pm}^V$ are simply

$$M_{\sigma\sigma}^{V} = \lambda(3\sigma^2 - v^2), \quad M_{\sigma\pm}^{V} = M_{\pm\sigma}^{V} = 0. \qquad (2.50)$$
The double derivatives with respect to L and L* are given by:

$$ 
\begin{aligned}
\mathcal{M}_{\pm\pm}^{V} &= \frac{6T^4 b(T)}{\mathcal{H}} (4L_{\pm} - L_{\mp}^2) - \frac{1}{T^4 b(T)} \left( \frac{\partial V_L}{\partial L_{\pm}} + \frac{T^4}{2} a(T) L_{\mp} \right)^2, \\
\mathcal{M}_{\pm\mp}^{V} &= \frac{T^4}{2} a(T) + \frac{6T^4 b(T)}{\mathcal{H}^2} [1 + 2L_+L_- + 21L_+^2L_-^2 - 4(2 + L_+L_-)(L_+^3 + L_-^3)].
\end{aligned} 
\quad (2.51) $$

The elements of the fermionic contribution $\mathcal{M}_{ij}^q$ read

$$ 
\begin{aligned}
\mathcal{M}_{\sigma\sigma}^q &= 2N_f N_c g^2 \sum_{\zeta=\pm 1} \int \frac{d^3 p}{(2\pi)^3 E} \left( \frac{\mathbf{p}^2}{E^2} f_{\zeta} + \frac{g^2 \sigma^2}{E T} f'_{\zeta} \right), \\
\mathcal{M}_{\sigma\pm}^q &= 2N_f N_c g^2 \sigma \sum_{\zeta=\pm 1} \int \frac{d^3 p}{(2\pi)^3 E} f_{\zeta,\pm} (1 + \delta_{\zeta,\mp} - 3f_{\zeta}), \\
\mathcal{M}_{\pm\pm}^q &= 2N_f N_c^2 T \sum_{\zeta=\pm 1} \int \frac{d^3 p}{(2\pi)^3} f_{\zeta,\pm}^2, \quad \mathcal{M}_{\pm\mp}^q = 2N_f N_c^2 T \sum_{\zeta=\pm 1} \int \frac{d^3 p}{(2\pi)^3} f_{\zeta,\pm} f_{\zeta,-},
\end{aligned} 
\quad (2.52) $$

where $\delta_{\zeta,\mp}$ is the Kronecker symbol and

$$ f'_{\zeta} \equiv \frac{\partial f_{\zeta}}{\partial(\beta\mathcal{E}_{\zeta})} = -3f_{\zeta}(1-f_{\zeta}) + 2L_{+}f_{\zeta,+} + 2L_{-}f_{\zeta,-}. \quad (2.53) $$

While the solution of Eq. (2.48) can be obtained by a straightforward inversion of the 3 × 3 matrix $\mathcal{M}_{ij}$, using the standard Newton iteration $p_i^{(n+1)} = p_i^{(n)} - [\mathcal{M}^{(n)}]_{ij}^{-1} \frac{\partial F}{\partial p_i^{(n)}},$ the stability of this method is guaranteed only when the initial guess is (sufficiently) close to the convergence, mean-field values $p_i^{\text{mf}}$. In our numerical code, we instead employ the Broyden method [57, 58], which adjusts the numerical step implied by Eq. (2.48) to increase the overall stability. Another advantage of this routine is that it provides a numerical approximation to the Jacobian matrix $\mathcal{M}_{ij}$ for the cases where obtaining the Jacobian analytically becomes cumbersome.

In our numerical analysis, we employed the Jacobian in a two-fold way. First, we have used the analytically calculated Jacobian to obtain the control parameters as a function of temperature and chemical potential (and also as a function of radius and angular velocity for the bounded case to be considered later). In addition, we used the numerical Jacobian to obtain the phase diagrams.

## 2.6 Thermodynamics of the model

In the absence of a thermodynamic discontinuity, the system goes continuously from one phase to the other, undergoing a so-called crossover phase transition. Pinpointing an exact position of the pseudo-critical parameters of the crossover transition has a certain degree of arbitrariness, which depends on the choice of the observables used to identify the crossover. In this paper, we employ the phase-space gradients of the control parameters, defined as

$$ dp_{\text{mf}} = \left(\frac{\partial p}{\partial T}\right)_{\text{mf}} dT + \left(\frac{\partial p}{\partial \mu}\right)_{\text{mf}} d\mu, \quad (2.54) $$
where $p \in \{\sigma, L, L^*\}$ is one of the control parameters and the “mf” subscript reminds that the variations are taken along the mean-field, thermodynamic trajectories of the control parameters imposed by the saddle-point equations (2.39). We compute these derivatives following Ref. [59] by noting that due to the constraint in Eq. (2.46), the free energy $F$ becomes an implicit function only of the thermodynamic parameters $T$ and $\mu$ since the control parameters $p_i \to p_i^{\text{mf}} \equiv p_i^{\text{mf}}(T, \mu)$ are replaced by their mean-field values. Considering now that $F(T, \mu) \to F(T, \mu; p_i)$ is an explicit function of the thermodynamic parameters $(T, \mu)$ and of the control parameters $p_i$, we can take the derivative of Eq. (2.46) with respect to a thermodynamic parameter $G_a$ (say, $G_1 = T$ and $G_2 = \mu$), leading to

$$ M_{ij} \frac{\partial p_j}{\partial G_a} + \frac{\partial^2 F}{\partial G_a \partial p_i} = 0, \qquad (2.55) $$

where the Jacobian matrix $M_{ij}$ was introduced in Eq. (2.48). Therefore, the derivatives of the control parameters with respect to the thermodynamic parameters can be obtained by inverting the $M_{ij}$ matrix,

$$ \frac{\partial p_i}{\partial G_a} = -M_{ij}^{-1} \frac{\partial^2 F}{\partial G_a \partial p_j}. \qquad (2.56) $$

The terms appearing on the right-hand side of the above equation can be computed explicitly starting from Eq. (2.41):

$$ \begin{aligned} \frac{\partial^2 F}{\partial T \partial \sigma} &= -\frac{2N_f N_c g^2 \sigma}{T^2} \sum_{\zeta=\pm 1} \int \frac{d^3 p}{(2\pi)^3 E} \mathcal{E}_\zeta f'_\zeta, \\ \frac{\partial^2 F}{\partial \mu \partial \sigma} &= -\frac{2N_f N_c g^2 \sigma}{T} \sum_{\zeta=\pm 1} \zeta \int \frac{d^3 p}{(2\pi)^3 E} f_\zeta, \\ \frac{\partial^2 F}{\partial T \partial L_\pm} &= \frac{\partial^2 V_L}{\partial T \partial L_\pm} - 2N_f N_c \sum_{\zeta=\pm 1} \int \frac{d^3 p}{(2\pi)^3} f_{\zeta, \pm} [1 + \beta \mathcal{E}_\zeta (1 + \delta_{\zeta, \mp} - 3f_\zeta)], \\ \frac{\partial^2 F}{\partial \mu \partial L_\pm} &= -2N_f N_c \sum_{\zeta=\pm 1} \zeta \int \frac{d^3 p}{(2\pi)^3} f_{\zeta, \pm} (1 + \delta_{\zeta, \mp} - 3f_\zeta), \end{aligned} \qquad (2.57) $$

where the second derivatives of the Polyakov potential read:

$$ \begin{aligned} \frac{\partial^2 V_L}{\partial T \partial L_\pm} &= -\frac{6}{\mathcal{H}}(L_\mp - 2L_\pm^2 + L_\mp^2 L_\pm) \frac{\partial(T^4 b(T))}{\partial T} - \frac{L_\mp}{2} \frac{\partial(T^4 a(T))}{\partial T}, \\ \frac{\partial(T^4 b(T))}{\partial T} &= b_3 T_0^3, \qquad \frac{\partial(T^4 a(T))}{\partial T} = 4T^3 a_0 + 3T_0 T^2 a_1 + 2T_0^2 T a_2. \end{aligned} \qquad (2.58) $$

We now study the phase diagram of the (P)LSM models discussed so far. To help identify the transition line, we run a scan over the $(T, \mu)$ plane along diagonal lines passing through the origin, $T = \mu = 0$, which makes an angle $\theta$ to the vertical axis. Introducing $G = \sqrt{T^2 + \mu^2}$, the temperature and chemical potential along this diagonal line search have the expressions $(T, \mu) = G(\cos \theta, \sin \theta)$, while the slope is computed as

$$ \frac{\partial p_i}{\partial G} = \cos \theta \frac{\partial p_i}{\partial T} + \sin \theta \frac{\partial p_i}{\partial \mu}, \qquad (2.59) $$
<figure><img src="image_2.png" /><figcaption>**Figure 1.** (a) Chiral and deconfinement phase diagram for the unbounded, nonrotating $LSM_q$ and $PLSM_q$ models. The critical points are denoted by the star "*". (b) The splitting between chiral and deconfinement phase diagram in $T-\mu$ plane. (c) A comparison between the numerically-calculated gap (shown with dashed blue lines and open squares) and the analytically-calculated slope (shown with colored dotted lines and filled circles) for the chiral phase diagram, as a function of $G = (T^2 + \mu^2)^{1/2}$, with increasing values of $\theta = \tan^{-1}(\mu/T)$ denoted by the color gradient from yellow to dark-green colors. (d) The dependence of the Polyakov loop and its slope on the temperature $T = G \cos \theta$ for two different values $\theta_1 = 0.58$ rad and $\theta_2 = 1.29$ rad, marked in the phase diagram of (b). The slopes are scaled by the constant factors $N_1 = 0.2$ and $N_2 = 0.005$, respectively, for presentation purposes. The corresponding chiral transition temperatures are marked by the black dotted lines. See text for more details.</figcaption></figure>

with $p_i \in \{\sigma, L, L^*\}$.

In Fig. 1(a), we present the phase diagram obtained by such a diagonal search for our 3 order parameters: $\sigma/f_\pi$, $L$ and $L^*$. For all three cases, the open symbols signify the pseudocritical point of a crossover transition, while the solid symbols mark the 1st order phase transition. The chiral phase diagram for the $LSM_q$ model is also presented. Comparing the $PLSM_q$ and $LSM_q$ phase diagrams shows that, due to the coupling with the background Polyakov loop, the chiral critical point moves towards higher temperature and lower chemical potential (denoted by the green star "*" in Fig. 1(a)).

Let us discuss the criteria we used to mark a first-order phase transition in our numerical analysis. We have already mentioned that, to pinpoint the crossover temperature for different order parameters, one can calculate the slope along diagonal lines of search, as given in Eq. (2.59), using Eq.(2.56). Along with this, one can define a numerical derivative,
using the second-ordered centered scheme, for the control parameters ($\sigma, L, L^*$) with respect to $T$ and $\mu$ (and $\Omega$, for the rotating case). We refer to this numerical derivative as “gap” to distinguish it from the “slope” calculated Eq. (2.59). In the crossover region, the two evaluations agree with each other up to $O(\delta G^2)$ and the ratio between “gap” and “slope” will approach 1 as the step $\delta G$ is decreased. However, in the first-order region, the slope evaluated on either side of the discontinuity is finite, but the “gap” goes to infinity as we refine around the region of discontinuity. Thus, the “gap” to “slope” ratio tends to $\infty$ when the transition is of first order. Around the critical point, both the slope and the gap tend to infinity. In our analysis, we have used the heuristic criterion that when “gap” / “slope” > 1.3, the transition is marked as first order. We marked the critical point on the curve where the ratio between the gap and the slope first exceeds the value 1.3.

Our calculation results for the slope (the dashed lines with the filled circles and the color gradient) and gap (the solid blue lines and the empty squares) for the order parameter $\sigma/f_\pi$, computed in the PLSM$_q$ model, are shown in Fig. 1(c). Here, we considered 6 values of $\theta(=\tan^{-1}(\mu/T))$ and we represented our results with respect to the parameter $G = \sqrt{T^2 + \mu^2}$ of our diagonal line search. The color gradient in the slope from yellow ($\theta=0$) to dark-green ($\theta=0.722 \text{ rad}$) corresponds to increasing values of $\theta$. For each value of $\theta$, the comparison between the gap and the slope shows that they both have a distinct peak structure at the same value of $G$ and the magnitudes of the two peaks agree for smaller values of $\theta$. We mark these points as points of crossover transition in the $T-\mu$ chiral phase diagram shown in Figs. 1(a) and 1(b). As $\theta$ increases, it can be seen that, although the peak positions of the gap and slope still agree, their magnitudes differ. While the peak of the slope has a finite value, the value of the gap tends to infinity. Such a feature corresponds to a first-order phase transition.

Note that for the chiral phase transition, in Fig. 1(c) the peak position for the gap and slope always agree. This situation changes as we focus on the deconfinement sector. In Fig. 1(d), we plotted the expectation value of the traced Polyakov loop ($L$) and the scaled slope, $N\partial L/\partial G$, where $N$ is a constant scaling factor, as a function of temperature for two particular values of $\theta$: a smaller one, $\theta_1 = 0.58 \text{ rad}$, where the chiral and deconfinement phase transitions are simultaneous, and a larger $\theta_2 = 1.29 \text{ rad}$, where we see a split between these transitions. These particular values of $\theta$ are highlighted on the phase diagram in Fig. 1(b). The slopes are appropriately scaled, $N_1 = 0.2$ and $N_2 = 0.005$, for presentation purposes.

In Fig. 1(d), we also showed the chiral first-order phase transition temperatures for these two values of $\theta$. Since the fields are coupled to each other, the gap in all control parameters ($\sigma, L$ and $L^*$) will simultaneously exhibit a peak at the point of discontinuity in the system, i.e. when $T = T_c^{(\sigma)}$. For $\theta = 0.58 \text{ rad}$, the slope of $L$ also shows a peak there and the value of $L$ changes significantly over the chiral phase transition. In this instance, we see that the chiral and deconfinement phase transitions are simultaneous. However, for $\theta = 1.29 \text{ rad}$, one can see that, although there is a small jump in $L$ around $T = T_c^{(\sigma)}$, the value of $L$ after the jump is significantly smaller, $L \approx 0.03$. Moreover, the slope only shows a small ridge at $T = T_c^{(\sigma)}$, having a much more prominent peak at a larger temperature, $T = T_c^{(L)} \approx 0.15 \text{ GeV}$. After the small jump at $T_c^{(\sigma)}$, we can observe a smooth increment
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Double-peak structure of the slope $\partial L/\partial G$ of the Polyakov loop with respect to $G = \sqrt{T^2 + \mu^2}$, computed on diagonal trajectories under various angles $\theta$. The deconfinement transition point is taken where the slope has the highest peak: the left peak for $\theta = 0.88$ and 0.92 rad; the right peak for $\theta = 0.96$, 0.99 and 1.03 rad.</figcaption></figure>

in $L$, and around the peak of the slope (at $T = T_c^{(L)}$), $L \approx 0.3$. We argue that the value of $L$ as small as 0.03 cannot be associated with the deconfinement phase transition. Rather, a confinement-deconfinement crossover can be marked at the position of the slope peak, shown as the intersection point between the green dashed line and blue open circles in Fig. 1(b), where $L \approx 0.3$.

There is a splitting between the chiral and deconfinement phase diagram for higher values of the chemical potential and, while the chiral transition shows a first-order behaviour, the deconfinement transition becomes of crossover type. This phenomenon, resembling a first-order phase transition, can be understood by looking at the slope $\partial L/\partial G$, introduced in Eq. (2.59). The pseudocritical parameters for $L$ are defined as those parameters ($T_{pc}$, $\mu_{pc}$) for which $\partial L/\partial G$, computed along the line with $(T, \mu) = G(\cos \theta, \sin \theta)$ at fixed $\theta$, attains a maximum value. As shown in Fig. 2, $\partial L/\partial G$ has a double peak structure around the point where the deconfinement and chiral restoration transition temperatures start to disagree, $T_{pc}(L) \neq T_{pc}(\sigma)$. For small values of the chemical potential $\mu$, the first peak, developed at a lower temperature at the point of the chiral restoration, dominates over the second peak. As $\mu$ increases, the magnitude of the second peak also increases, eventually becoming dominant, leading to a split between the deconfinement and chiral restoration transitions.

## 2.7 Vanishing temperature limit

At vanishing temperature, the distribution functions become step functions:

$$f_{\pm}\bigg|_{T \to 0} = \theta(\pm\mu - E), \quad f_{\pm,+}\bigg|_{T \to 0} = f_{\pm,-}\bigg|_{T \to 0} = 0. \qquad (2.60)$$

Since both functions $f_{\zeta,L}$ and $f_{\zeta,L*}$ vanish, the equations $\partial F/\partial L = 0$ and $\partial F/\partial L* = 0$ are automatically satisfied for any value of $L$ and $L*$. We therefore focus on the $\sigma$ sector, where we expect a first-order phase transition. As a consequence, the Polyakov loop variable does
not play any role at vanishing temperature and the model reduces to the linear sigma model coupled to quarks only (LSM_q).

We must therefore evaluate the fermion condensate [cf. Eq. (2.42)] to enforce $\partial F/\partial \sigma = 0$, and calculate the pressure [cf. Eq. (2.45)], to ensure the minimization of the free energy. Both of these quantities can be assessed analytically:

$$ \begin{aligned} \langle \bar{\psi} \psi \rangle &= \frac{N_f N_c g \sigma}{\pi^2} \int_0^{p_f} \frac{dp}{E} \frac{dp^2}{E} = \frac{N_f N_c g \sigma}{2\pi^2} \left[ p_f |\mu| - (g\sigma)^2 \operatorname{artanh} \frac{p_f}{|\mu|} \right] \theta(|\mu| - g\sigma), \\ P_q &= \frac{N_f N_c}{3\pi^2} \int_0^{p_f} \frac{dp}{E} \frac{dp^4}{E} = \frac{N_f N_c}{24\pi^2} \left[ |\mu| p_f (2\mu^2 - 5g^2\sigma^2) + 3(g\sigma)^4 \operatorname{artanh} \frac{p_f}{|\mu|} \right] \theta(|\mu| - g\sigma), \end{aligned} \quad (2.61) $$

where $p_f = \sqrt{\mu^2 - g^2\sigma^2}$ is the Fermi momentum, and we employed the integration of the phase space $d^3p = p^2dpd\Omega$ with the $d\Omega$ integral trivially giving $4\pi$. The above results are valid only when the Fermi energy exceeds the effective mass, i.e., $|\mu| > g\sigma$, as implied by the step function appearing on the right-hand sides of Eqs. (2.61).

## 3 Rotating PLSM_q model

To facilitate the analysis of the rotating PLSM_q model, we will employ corotating coordinates, defined by the angular velocity $\Omega$ entering the grand canonical ensemble. To avoid superluminal corotating velocities, the system must be enclosed in a cylindrical box, where suitable boundary conditions must be enforced. In Subsec. 3.1, we briefly present the model Lagrangian in corotating coordinates. In Subsec. 3.2, we derive the grand potential for the bounded, rotating system, obeying spectral boundary conditions. In Subsec. 3.3, we discuss the saddle point equations giving the mean-field values for the mesons and Polyakov loops. Finally, Subsections 3.4 and 3.5 discuss the LSM_q limit and the vanishing temperature limit, respectively.

### 3.1 Corotating frame

We now switch to corotating coordinates, defined with respect to the static cylindrical coordinates (denoted with the subscript M) by

$$ t = t_M, \quad \varphi = \varphi_M - \Omega t_M, \quad (3.1) $$

where $\Omega$ represents the angular velocity of the corotating observer. The line element $ds^2 = g_{\mu\nu}dx^\mu dx^\nu$ expressed with respect to the corotating coordinates reads

$$ ds^2 = (1 - \rho^2 \Omega^2) dt^2 - 2 \rho^2 \Omega dt d\varphi - d\rho^2 - \rho^2 d\varphi^2 - dz^2. \quad (3.2) $$

The metric can be brought to Minkowski form by introducing the non-holonomic, orthonormal tetrad (vielbein) vector field, $e_\alpha = e_\alpha^\mu \partial_\mu$, and the associated one-forms $\omega^\alpha = \omega_\mu^\alpha dx^\mu$, defined by

$$ \begin{aligned} e_{\hat{t}} &= \partial_t + \Omega(y\partial_x - x\partial_y), & e_{\hat{x}} &= \partial_x, & e_{\hat{y}} &= \partial_y, & e_{\hat{z}} &= \partial_z, \\ \omega^{\hat{t}} &= dt, & \omega^{\hat{x}} &= dx - \Omega y dt, & \omega^{\hat{y}} &= dy + \Omega x dt, & \omega^{\hat{z}} &= dz, \end{aligned} \quad (3.3) $$
such that $g_{\mu\nu} = \eta_{\hat{\alpha}\hat{\beta}}\omega_{\mu}^{\hat{\alpha}}\omega_{\nu}^{\hat{\beta}}$ and $\omega_{\mu}^{\hat{\alpha}}e_{\hat{\beta}}^{\mu} = \delta_{\hat{\beta}}^{\hat{\alpha}}$.

The Cartan coefficients associated with the tetrad are defined by the commutation relation $[e_{\hat{\alpha}}, e_{\hat{\beta}}] = c_{\hat{\alpha}\hat{\beta}}e_{\hat{\sigma}}$. Taking into account the commutation relations $[e_{\hat{t}}, e_{\hat{x}}] = \Omega e_{\hat{y}}$ and $[e_{\hat{t}}, e_{\hat{y}}] = -\Omega e_{\hat{x}}$, the non-vanishing Cartan coefficients are $c_{\hat{t}\hat{x}}^{\hat{y}} = c_{\hat{y}\hat{t}}^{\hat{x}} = \Omega$. The connection coefficients $\Gamma_{\hat{\rho}\hat{\sigma}\hat{\alpha}} = \frac{1}{2}(c_{\hat{\rho}\hat{\sigma}\hat{\alpha}} + c_{\hat{\rho}\hat{\alpha}\hat{\sigma}} - c_{\hat{\sigma}\hat{\alpha}\hat{\rho}})$ vanish identically, except for

$$ \Gamma_{\hat{x}\hat{y}\hat{t}} = \Omega. \qquad (3.4) $$

The quark Lagrangian (2.18) can be written in curvilinear Minkowski coordinates as

$$ \mathcal{L}_q = \bar{\psi} (i \not{D} + i \not{\slashed}{\mathcal{D}} + \mu \gamma^0 - g\sigma) \psi, \qquad (3.5) $$

where $\Gamma_{\mu} = -\frac{i}{2}\omega_{\mu}^{\hat{\alpha}}\Gamma_{\hat{\beta}\hat{\gamma}\hat{\alpha}}S^{\hat{\beta}\hat{\gamma}}$ is the spin connection coefficient and $S^{\hat{\beta}\hat{\gamma}} = \frac{1}{4}[\gamma^{\hat{\beta}}, \gamma^{\hat{\gamma}}]$ represent the spin part of the generators of the Lorentz transformations. We take the gamma matrices $\gamma^i$ and $\gamma^i$ with respect to the tetrad vector fields in the Dirac representation, shown in Eq. (2.20). With the above convention, $\Gamma_{\hat{t}} = 0$ and $\Gamma_{\hat{t}} = -i\Omega S^z$ and Eq. (3.5) becomes

$$ \mathcal{L}_q = \bar{\psi}[\gamma^0(i\partial_t + A_0 + \mu + \Omega J^z) + i\gamma \cdot \nabla - g\sigma]\psi. \qquad (3.6) $$

Switching now to imaginary time, $\tau = it$, the quark contribution to the partition function reads

$$ \mathcal{Z}_q = \int [i d\psi^\dagger] [d\psi] \exp(-S_E), \qquad S_E = \int_0^\beta d\tau \int_V d^3x \mathcal{L}_E, \qquad (3.7) $$

where $V$ represents the system volume (a cylinder of radius $R$ and infinite vertical extent), while the Euclidean Lagrangian reads

$$ \mathcal{L}_E = -\bar{\psi} [\gamma^0 (-\partial_\tau + \Omega J^z - iA_4 + \mu) + i\gamma \cdot \nabla - g\sigma] \psi. \qquad (3.8) $$

## 3.2 Grand canonical potential

Let us now compute the partition function $\mathcal{Z}_q$ and the associated thermodynamic potential, $F_q = -\frac{T}{V} \ln \mathcal{Z}_q$. As outlined in Sec. 2.4, we seek a mode decomposition of the spinor $\psi$ that diagonalizes the Euclidean Lagrangian. Since the momentum operator $\mathbf{P}$ does not commute with the angular momentum operator $J^z$ appearing in Eq. (3.8), we must employ cylindrical modes that are instead eigenvectors of $J^z$ and $P^z = -i\partial_z$. Thus, Eq. (2.24) is replaced by

$$ (H, -i\partial_{\tau}, \sigma_f^3, A_4, h, J^z, P^z)\psi_j = (\zeta_j \tilde{\mathcal{E}}_j, w_j, f_j, A_{4:j}, \lambda_j, m_j, k_j)\psi_j, \qquad (3.9) $$

where we included the Hamiltonian $H = iA_4 - \mu - \Omega J^z + 2\gamma^5\mathbf{S} \cdot \mathbf{P} + g\sigma\gamma^0$, whose eigenvalues satisfy

$$ \tilde{\mathcal{E}}_j = \mathcal{E}_j - \Omega m_j, \qquad \mathcal{E}_j = E_j - \zeta_j(\mu - iA_{4:j}), \qquad E_j = \sqrt{\mathbf{P}_j^2 + g^2\sigma^2}. \qquad (3.10) $$

The modes described above have the following form:

$$ \psi_j(X) = C_j \frac{e^{i w_j \tau + i k_j z}}{2\pi} \left(\frac{1}{2}(1+f_j)\right) \otimes \left(\frac{\delta_{c_j,1}}{\delta_{c_j,2}}\right) \otimes \left(\frac{1}{\sqrt{2}} \left(2\zeta_j \lambda_j \mathfrak{E}_j^\zeta\right)\right) \otimes \phi_j(\rho, \varphi), \qquad (3.11) $$
where $C_j$ is a normalization constant, $\mathfrak{E}_j = \sqrt{1 + \zeta_j \frac{g\sigma}{E_j}}$ was introduced in Eq. (2.26) and $\phi_j$ satisfying $(-i\partial_\varphi + S^z)\phi_j = m_j\phi_j$ and $-i\sigma \cdot \nabla\phi_j = 2\lambda_j p_j\phi_j$ is given by

$$ \phi_j(\rho, \varphi) = \frac{1}{\sqrt{2}} \left( \begin{array}{c} \mathbf{p}_j^+ e^{i(m_j - \frac{1}{2})\varphi} J_{m_j - \frac{1}{2}}(q_j \rho) \\ 2i\lambda_j \mathbf{p}_j^- e^{i(m_j + \frac{1}{2})\varphi} J_{m_j + \frac{1}{2}}(q_j \rho) \end{array} \right), \qquad \mathbf{p}_j^{\pm} = \sqrt{1 \pm \frac{2\lambda_j k_j}{p_j}}. \qquad (3.12) $$

We now wish to impose orthogonality of the modes $\psi_j$ with respect to the integration over a cylinder of radius $R$. As discussed in Ref. [37], this can be ensured in a number of ways. In the spectral boundary conditions model, the transverse momentum $q_j$ is quantized according to

$$ q_j R = \begin{cases} \xi_{m_j - \frac{1}{2}, l_j}, & m_j > 0, \\ \xi_{-m_j - \frac{1}{2}, l_j}, & m_j < 0, \end{cases} \qquad (3.13) $$

with $\xi_{nl}$ being the l'th nonzero root of the Bessel function $J_n$, i.e. $J_n(\xi_{nl}) = 0$ ($n$ is a non-negative integer). We thus impose the normalization $\int_0^\beta d\tau \int_V d^3x \psi_j^\dagger \psi_{j'} = \beta\delta(j, j')$, following Eq. (2.28), which gives the normalization constant $C_j$ as

$$ C_j = \frac{\sqrt{2}}{R|J_{|m_j|+\frac{1}{2}}(\xi_{|m_j|-\frac{1}{2},l_j})|} \qquad (3.14) $$

The delta function reads

$$ \delta(j, j') = \delta_{\zeta_j \zeta_{j'}} \delta_{n_j n_{j'}} \delta_{f_j f_{j'}} \delta_{c_j c_{j'}} \delta_{\lambda_j \lambda_{j'}} \delta_{m_j m_{j'}} \delta_{l_j l_{j'}} \delta(k_j - k_{j'}). \qquad (3.15) $$

Taking again the decomposition $\psi(X) = \sum_j \hat{b}_j \psi_j(X)$ using the cylindrical modes introduced in Eq. (3.11), the partition function looks formally identical to the one in the static case, with $\mathcal{E}_j$ replaced by $\tilde{\mathcal{E}}_j$ in Eqs. (2.29)–(2.31). Repeating for definiteness Eq. (2.31), we have $\tilde{Z}_q = \prod_j [\beta(i\zeta_j \tilde{\mathcal{E}}_j - w_j)]$. Upon taking the logarithm to express the grand potential, $\tilde{F}_q = -\frac{T}{V} \ln \tilde{Z}_q$, the continuous momentum integration in Eq. (2.32) is replaced using the following rule:

$$ \int \frac{d^3 p}{(2\pi)^3} \rightarrow \frac{1}{\pi R^2} \sum_b \int \frac{dk}{2\pi}, \qquad \sum_b = \sum_{m=-\infty}^{\infty} \sum_{l=1}^{\infty}, \qquad (3.16) $$

keeping in mind that $m = \pm \frac{1}{2}, \pm \frac{3}{2}, \dots$ runs over all odd half-integer values. This leads to the following modification of Eq. (2.36):

$$ \tilde{F}_q = -\frac{2N_f N_c T}{\pi R^2} \sum_{\zeta=\pm 1} \sum_b \int_{-\infty}^{\infty} \frac{dk}{2\pi} \tilde{F}_{\zeta}, \qquad (3.17) $$

where $\tilde{F}_\zeta$ is $F_\zeta$ from Eq. (2.37), with $\mathcal{E}_\zeta$ replaced by $\tilde{\mathcal{E}}_\zeta$:

$$ \tilde{F}_{\zeta} \equiv \frac{1}{3} \ln(1 + 3L_{\zeta} e^{-\beta \tilde{\mathcal{E}}_{\zeta}} + 3L_{-\zeta} e^{-2\beta \tilde{\mathcal{E}}_{\zeta}} + e^{-3\beta \tilde{\mathcal{E}}_{\zeta}}). \qquad (3.18) $$

Note that for the numerical analysis, one can further simplify the expression given in Eq. (3.16) by replacing the sum over the odd half-integers $m$ with a sum over non-negative
integers, $\sum_{m=-\infty}^{\infty} f_m = \sum_{n=0}^{\infty}[f_{n+\frac{1}{2}} + f_{-n-\frac{1}{2}}]$, with the convention $n = |m| - \frac{1}{2}$. The corotating effective energies are modified to

$$ \tilde{\mathcal{E}}_{m>0,l} \rightarrow \tilde{E}_{n+\frac{1}{2},l} - \zeta\mu, \qquad \tilde{\mathcal{E}}_{m<0,l} \rightarrow \bar{E}_{n+\frac{1}{2},l} - \zeta\mu, \qquad (3.19) $$

where we introduced

$$ \tilde{E}_{n+\frac{1}{2},l} = E_{n+\frac{1}{2},l} - \Omega \left(n + \frac{1}{2}\right), \qquad \bar{E}_{n+\frac{1}{2},l} = E_{n+\frac{1}{2},l} + \Omega \left(n + \frac{1}{2}\right), \qquad (3.20) $$

while $E_{n+\frac{1}{2},l} = \sqrt{k^2 + R^{-2}\xi_{nl}^2 + g^2\sigma^2}$. In arriving at Eq. (3.19), we took into account the fact that $\tilde{E}_{-ml} = E_{ml}$, due to the property (3.13).

## 3.3 Saddle point approximation

As in the unbounded case discussed in Sec. 2.5, the control parameters of the model are $p \in \{\sigma, L, L^*\}$. Their expectation values at finite $T$, $\mu$ and $\Omega$ can be obtained by imposing the minimization of the free energy via Eqs. (2.39), with $F$ replaced by $\tilde{F} = V_M + V_L + \tilde{F}_q$, where the meson ($V_M$) and Polyakov ($V_L$) contributions remain unchanged. The derivatives $\partial\tilde{F}/\partial p_i$ remain as given in Eqs. (2.41)-(2.43b), with the replacement of the integration measure as shown in Eq. (3.16), as well as the replacement of the distributions $f_\zeta$ and $f_{\zeta;\pm}$ by their rotating counterparts:

$$ \tilde{f}_{\pm} \equiv -\frac{\partial \tilde{F}_{\pm}}{\partial(\beta \tilde{\mathcal{E}}_{\pm})} = \frac{L_{\pm}e^{-\beta \tilde{\mathcal{E}}_{\pm}} + 2L_{\mp}e^{-2\beta \tilde{\mathcal{E}}_{\pm}} + e^{-3\beta \tilde{\mathcal{E}}_{\pm}}}{1 + 3L_{\pm}e^{-\beta \tilde{\mathcal{E}}_{\pm}} + 3L_{\mp}e^{-2\beta \tilde{\mathcal{E}}_{\pm}} + e^{-3\beta \tilde{\mathcal{E}}_{\pm}}}, \qquad (3.21a) $$

$$ \tilde{f}_{\pm, \pm} \equiv \frac{\partial \tilde{F}_{\pm}}{\partial L_{\pm}} = \frac{e^{-\beta \tilde{\mathcal{E}}_{\pm}}}{1 + 3L_{\pm}e^{-\beta \tilde{\mathcal{E}}_{\pm}} + 3L_{\mp}e^{-2\beta \tilde{\mathcal{E}}_{\pm}} + e^{-3\beta \tilde{\mathcal{E}}_{\pm}}}, \qquad (3.21b) $$

$$ \tilde{f}_{\pm, \mp} \equiv \frac{\partial \tilde{F}_{\pm}}{\partial L_{\mp}} = \frac{e^{-2\beta \tilde{\mathcal{E}}_{\pm}}}{1 + 3L_{\pm}e^{-\beta \tilde{\mathcal{E}}_{\pm}} + 3L_{\mp}e^{-2\beta \tilde{\mathcal{E}}_{\pm}} + e^{-3\beta \tilde{\mathcal{E}}_{\pm}}}, \qquad (3.21c) $$

$$ \tilde{f}'_{\zeta} \equiv \frac{\partial \tilde{f}_{\zeta}}{\partial(\beta \tilde{\mathcal{E}}_{\zeta})} = -3\tilde{f}_{\zeta}(1 - \tilde{f}_{\zeta}) + 2L_{+}\tilde{f}_{\zeta,+} + 2L_{-}\tilde{f}_{\zeta,-}. \qquad (3.21d) $$

As in the unbounded case, we solve the saddle point equations using Eq. (2.48), with the matrix $M_{ij}$ replaced by $\tilde{M}_{ij}$, constructed as discussed above: the meson and Polyakov parts remain the same, as given in Eqs. (2.50) and (2.51), while the fermion contributions in Eq. (2.52) are modified by replacing the momentum integration via Eq. (3.16) and by replacing $(f_\zeta, f_{\zeta,\pm}, f'_\zeta)$ by $(\tilde{f}_\zeta, \tilde{f}_{\zeta,\pm}, \tilde{f}'_\zeta)$.

When discussing the pseudocritical parameters corresponding to a crossover transition, we define the gradient of the control parameters $p \in \{\sigma, L, L^*\}$ by extending Eq. (2.54) to include the effect of a finite angular velocity $\Omega$:

$$ dp_{\text{mf}} = \left(\frac{\partial p}{\partial T}\right)_{\text{mf}} dT + \left(\frac{\partial p}{\partial \mu}\right)_{\text{mf}} d\mu + \left(\frac{\partial p}{\partial \Omega}\right)_{\text{mf}} d\Omega. \qquad (3.22) $$

The evaluation of the derivatives $\partial p_i/\partial G_a$ of the control parameters with respect to the thermodynamic parameters $A_a \in \{T, \mu, \Omega\}$ can be performed as shown in Eq. (2.56), where
the derivatives of the free energy appearing on the right-hand side of this equation with respect to T and µ are identical to those given in Eqs. (2.57), with ($F, f_ζ, f_{ζ,±}, \tilde{\mathcal{E}}_ζ$) replaced by ($\tilde{F}, \tilde{f}_ζ, \tilde{f}_{ζ,±}, \tilde{\mathcal{E}}_ζ$) and the integration measure modified according to Eq. (3.16). In addition, we compute the derivatives involving the rotation angular velocity:

$$ 
\begin{aligned}
\frac{\partial^2 \tilde{F}}{\partial \Omega \partial \sigma} &= -\frac{N_f N_c g^2 \sigma}{\pi^2 R^2 T} \sum_{b, \zeta = \pm 1} m \int_{-\infty}^{\infty} \frac{dk}{E} f'_{\zeta}, \\
\frac{\partial^2 \tilde{F}}{\partial \Omega \partial L_{\pm}} &= -\frac{N_f N_c}{\pi^2 R^2} \sum_{b, \zeta = \pm 1} m \int_{-\infty}^{\infty} dk \tilde{f}_{\zeta, \pm} (1 + \delta_{\zeta, \mp} - 3\tilde{f}_{\zeta}),
\end{aligned}
\quad (3.23) $$

where we remind the reader that $m = \pm \frac{1}{2}, \pm \frac{3}{2}, \dots$ takes odd half-integer values, while $\delta_{\zeta, \zeta'}$ is the Kronecker symbol.

### 3.4 Rotating LSM_q model with boundary

We now discuss the limit $L, L^* \to 1$, when the system is deconfined and the PLSM_q model reduces to the LSM_q model. Here, the only control parameter is $\sigma$, and the corresponding volume-averaged grand potential for the quark sector can be obtained by putting $L = L^* = 1$ in Eq. (3.17):

$$ \tilde{F}_q = -\frac{N_f N_c T}{\pi^2 R^2} \sum_{b, \zeta = \pm 1} \int_{-\infty}^{\infty} dk \, \tilde{F}_{\zeta}, \quad \tilde{F}_{\zeta} = \ln(1 + e^{-\beta \tilde{\mathcal{E}}_{\zeta}}). \quad (3.24) $$

The expectation value of the $\sigma$ meson, $\sigma_{mf}$, can be found by imposing the minimization of the grand potential,

$$ \frac{\partial \tilde{F}}{\partial \sigma} = \frac{\partial \tilde{V}_{\mathcal{M}}}{\partial \sigma} + \frac{\partial \tilde{F}_q}{\partial \sigma} = 0, \qquad \frac{\partial \tilde{F}_q}{\partial \sigma} = \frac{N_f N_c g^2 \sigma}{\pi^2 R^2} \sum_{b, \zeta = \pm 1} \int_{-\infty}^{\infty} \frac{dk}{E} \tilde{f}_{\zeta}, \quad (3.25) $$

where $\partial V_{\mathcal{M}}/\partial \sigma$ is given in Eq. (2.43a) and

$$ \tilde{f}_{\zeta} = -\frac{\partial \tilde{F}_{\zeta}}{\partial (\beta \tilde{\mathcal{E}}_{\zeta})} = \frac{1}{e^{\beta \tilde{\mathcal{E}}_{\zeta}} + 1}, \qquad \tilde{f}'_{\zeta} = -\tilde{f}_{\zeta}(1 - \tilde{f}_{\zeta}). \quad (3.26) $$

Eq. (3.25) is solved numerically using Eq. (2.48) and the following contribution from the fermionic part,

$$ \tilde{\mathcal{M}}_{\sigma\sigma}^q = \frac{\partial^2 \tilde{F}_q}{\partial \sigma^2} = \frac{N_f N_c g^2}{\pi^2 R^2} \sum_{b, \zeta = \pm 1} \int_{-\infty}^{\infty} \frac{dk}{E} \left\{ \frac{\mathbf{p}^2}{E^2} \tilde{f}_{\zeta} + \frac{g^2 \sigma^2}{E T} \tilde{f}'_{\zeta} \right\}, \quad (3.27) $$

with $\mathbf{p}^2 = \mathbf{p}_{kml}^2 = k^2 + R^{-2}\xi^2_{|m|-\frac{1}{2},l}$. The mesonic contribution $\tilde{\mathcal{M}}_{\sigma\sigma}^V = \mathcal{M}_{\sigma\sigma}^V$ is identical to that of the unbounded case, given in Eq. (2.51). When searching for the pseudocritical parameters of a cross-over phase transition, the slopes $\partial\sigma/\partial G_a$, with $G_a \in \{T, \mu, \Omega\}$, can be evaluated using:

$$ \left( \frac{\partial^2 \tilde{F}}{\partial T \partial \sigma}, \frac{\partial^2 \tilde{F}}{\partial \mu \partial \sigma}, \frac{\partial^2 \tilde{F}}{\partial \Omega \partial \sigma} \right) = - \frac{N_f N_c g^2 \sigma}{\pi^2 R^2 T} \sum_{b, \zeta = \pm 1} \int_{-\infty}^{\infty} \frac{dk}{E} (\beta \tilde{\mathcal{E}}_{\zeta}, \zeta, m) \tilde{f}'_{\zeta}. \quad (3.28) $$
<figure><img src="image_4.png" /><figcaption>Figure 3. The order parameters (a) $\sigma/f_\pi$ and (b) $L$, and their associated analytically-calculated slope (shown with colored dashed lines) as functions of $G = (T^2 + \mu^2)^{1/2}$, with $\tan\theta \simeq 54.5$ and $R = 1$ fm, in the absence of rotation. The slopes are scaled by the constant factors (a) $N_\sigma = 0.1$ and (b) $N_L = 10$, respectively.</figcaption></figure>

### 3.5 Vanishing temperature limit

In the $T \to 0$ limit, the distributions $\tilde{f}_\zeta$ and $\tilde{f}_{\zeta, \pm}$ reduce to

$$ \tilde{f}_{\zeta} \to \theta(\zeta\mu - \tilde{\mathcal{E}}_{\zeta}), \quad \tilde{f}_{\zeta,+}, \tilde{f}_{\zeta,-} \to 0. \qquad (3.29) $$

The Polyakov sector completely decouples from the quark sector and the resulting model becomes equivalent to the LSM$_q$ model without the Polyakov extension. The thermodynamic potential $\tilde{F}_{T=0} = V_M + \tilde{F}_q^{T=0}$ receives no contribution from the Polyakov sector, while the fermionic contribution reads

$$ \begin{aligned} \tilde{F}_q^{T=0} &= -\frac{2N_f N_c}{\pi^2 R^2} \sum_{b, \zeta = \pm 1} \int_0^\infty \frac{dk \, k^2}{E_{ml}(k)} \theta[\zeta\mu + \Omega m - E_{ml}(k)] && (3.30) \\ &= -\frac{2N_f N_c}{\pi^2 R^2} \sum_{b, \zeta = \pm 1} \left( \frac{1}{2} E_{ml}^f k_{ml}^f - M_{ml}^2 \operatorname{artanh} \frac{k_{ml}^f}{E_{ml}^f + M_{ml}} \right) \theta(\zeta\mu + \Omega m - M_{ml}), \end{aligned} $$

where $E_{ml}(k) = \sqrt{k^2 + M_{ml}^2}$ is the Minkowski energy and $M_{ml} = \sqrt{q_{ml}^2 + g^2\sigma^2}$ is the transverse mass. For notational brevity, we introduced the Fermi energy $E_{ml}^f = |\zeta\mu + \Omega m|$ and the Fermi momentum, $k_{ml}^f = \sqrt{(E_{ml}^f)^2 - M_{ml}^2}$. The saddle point equation becomes

$$ \frac{\partial \tilde{F}_q^{T=0}}{\partial \sigma} = \frac{4g^2 \sigma N_f N_c}{\pi^2 R^2} \sum_{b, \zeta = \pm 1} \operatorname{artanh} \frac{k_{ml}^f}{E_{ml}^f + M_{ml}} \theta(\zeta\mu + \Omega m - M_{ml}). \qquad (3.31) $$
For the evaluation of the slope, $\partial\sigma/\partial\mu$, we require:

$$ \frac{\partial^2 \tilde{F}_q^{T=0}}{\partial \sigma^2} = -\frac{2g^2 N_f N_c}{\pi^2 R^2} \sum_{b, \zeta = \pm 1} \left( \frac{g^2 \sigma^2 E_{ml}^f}{M_{ml}^2 k_{ml}^f} - 2 \operatorname{artanh} \frac{k_{ml}^f}{E_{ml}^f + M_{ml}} \right) \theta(\zeta\mu + \Omega m - M_{ml}), \\ \\ \frac{\partial^2 \tilde{F}_q^{T=0}}{\partial \mu \partial \sigma} = \frac{2g^2 \sigma N_f N_c}{\pi^2 R^2} \sum_{b, \zeta = \pm 1} \frac{1}{k_{ml}^f} \theta(\zeta\mu + \Omega m - M_{ml}). \qquad (3.32) $$

The presence of $k_{ml}^f$ in the denominators of some of the terms appearing in Eq. (3.32) indicates that, at small temperatures, the slope of $\sigma$ will diverge whenever modes for which $k_{ml}^f$ vanishes are supported by the system. In particular, the step function $\theta(\zeta\mu+\Omega m-M_{ml})$ “fires” precisely when $k_{ml}^f = 0$, such that when $\mu$ is continuously increased, leading to the inclusion of a new mode in the b summation, the contribution of this mode will be arbitrarily large. In practice, this is unlikely to happen, as the mean-field value of $\sigma$ is constantly adjusted, eventually leading to finite $k_{ml}^f$ at the point where the grand potential takes its minimum. Signatures of this effect are still visible, as the slope can still develop peaks whenever a new root is allowed within the Fermi level. These oscillations are the analogues of the Shubnikov-de Haas effect arising due to the periodic modulation of the density of states at the Fermi level in a conducting material subjected to a background, external magnetic field [60].

In panel (a) of Fig. 3, we demonstrate this effect at small, but non-vanishing temperature, for a small system of $R = 1$ fm, by considering the variation of $\sigma$ and $\partial\sigma/\partial G$ along a diagonal phase-space line with $(T, \mu) = G(\cos\theta, \sin\theta)$, at a fixed value of $\theta$ (we took $\tan\theta = 54.5$). While our analysis in this subsection is focused mainly on the behavior of $\sigma$ and its slope, panel (b) of Fig. 3 indicates that the peak structure is also inherited in the slope $\partial L/\partial G$ of the Polyakov loop. Moreover, our results show that following the development of one peak, the slope can become negative, leading to a non-monotonic dependence of $L$ on the control parameter, $G$.

The above discussion, together with Fig. 3, illustrates the conceptual difficulties in pinpointing the parameters of the deconfinement phase transition. At low temperatures, the first-order chiral transition unambiguously occurs at some value of $G$. At this point, the Polyakov loop exhibits only a modest jump, such that the system remains in a confined state. As $G$ increases, the Polyakov loop exhibits a generally increasing trend, but its nonmonotonicity makes it difficult to identify the point of highest peak in its slope. In practice, our algorithm picks the parameters of the deconfinement transition by looking for the closest peak in the vicinity of $L_{th} = 0.3$.

# 4 Phase Diagram

We now discuss the phase transition properties of strongly-interacting matter enclosed in a cylindrical boundary, with or without rotation. For brevity, we will only consider the logarithmic potential (2.7) for the Polyakov loop in what follows.

We first remark that the effects of the boundary and/or of the rotation manifest themselves only through the fermionic contributions to the saddle-point equations. On one
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Variation of $\sigma/f_\pi$ (curves starting at 1) and $L, L^*$ ($L^*$ is shown with gray dotted lines for $\mu = 0.25$ GeV) for a static ($\Omega R = 0$) cylindrical system of various radii for (a) $\mu = 0$ GeV and (b) $\mu = 0.25$ GeV.</figcaption></figure>

hand, we anticipate that a reduced volume (low radius $R$) will lead to a suppression of the fermionic contributions, due to the additional mass gap introduced by the transverse momentum quantization. Thus, the presence of the boundary produces an effect of delaying chiral restoration. Deconfinement also takes place at higher temperatures compared to the case without boundaries, up to the point where the temperature approaches that of the pure glue system, $T_0$. Since the Polyakov potential $V_L$ is not affected by either boundary or rotation, the deconfinement transition, in the absence of fermionic contributions, occurs as in the pure glue system described in Sec. 2.2, namely as a first-order transition at $T = T_0$. When the radius increases, the gap due to the boundary becomes smaller and we can expect to recover the thermodynamics of the unbounded system.

On the other hand, rotation has the effect of amplifying the thermal fermionic contributions, as suggested by the Tolman-Ehrenfest (TE) law [14, 15]. The strength of the effect is characterized by the velocity of the system at the boundary, $v_R = \Omega R$, or, more precisely, by its Lorentz factor $\Gamma_R = (1 - \Omega^2 R^2)^{-1/2}$. Compared to the static system, the chiral phase transition in the rotating medium should occur at lower temperatures, as indicated in Eq. (5.33), discussed in the following section of this paper. While Eq. (5.33) predicts that the chiral restoration temperature $T_\sigma^{\Omega R}$ vanishes when $\Omega R \to 1$, this is not the case at any finite radius $R < \infty$. The finite mass gap induced by the boundary is sufficient to render the system causal and all expectation values finite, even in the causality limit $\Omega R \to 1$. Thus, at finite $R$, quantum effects have a sizable influence on the system thermodynamics, especially at low $R$. When the radius of the system is pushed to infinity, keeping $\Omega R$ fixed, we can expect these quantum effects to become subleading and to recover the Tolman-Ehrenfest predictions.

We address the nonrotating case in Subsec. 4.1, with an emphasis on the suppression effects of the boundary and the restoration of the unbounded system thermodynamics as $R \to \infty$. The rotating system is discussed in Subsec. 4.2, where we explore the influence of the rigid rotation on the chiral and deconfinement phase transitions in cylindrical systems with radii of a few fermi that match sizes of typical realistic systems.
## 4.1 Static system: Effect of boundary

We begin the analysis of the effect of the boundary with a set of plots, shown in Fig. 5, presenting the temperature dependence of $\sigma$ (panels a,b), $L$ and $L^*$ (panels c,d) for different values of the radius $R$ and two different values of the chemical potential: $\mu = 0$ (a,c) and $\mu = 0.25$ GeV (b,d). For vanishing chemical potential, all control parameters $\sigma$, $L$ and $L^*$ exhibit a smooth variation with respect to $T$. For $\mu = 0.25$ GeV, the first-order phase transition expected from the unbounded case (shown with red lines, marked as “Unb.” and discussed in Sec. 2) gradually softens into a crossover transition as the system size $R$ is decreased. The softening effect becomes observable for $R \lesssim 3$ fm. The transition point, be it of first-order or crossover type, moves to higher temperatures as $R$ is decreased. This behavior agrees with the previously reported results with the bag boundary conditions in Ref. [39].

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Variation of $\sigma/f_\pi$ (a,b) and $L, L^*$ (c,d) for fixed values of the chemical potential, $\mu = 0$ (a,c) and $\mu = 0.25$ GeV (b,d) for a static ($\Omega R = 0$) cylindrical systems of various radii $R$. In panels (c) and (d), $L$ is shown using filled symbols and $L^*$ is shown with solid lines.</figcaption></figure>

We now investigate the restoration of the first-order deconfinement phase transition in the limit where the fermionic contributions are suppressed due to the small system size. For simplicity, we focus on the case of vanishing chemical potential. Figure 6 confirms that as $R < 1$ fm, both the deconfinement and the chiral restoration transitions are inhibited. As $R \to 0$, the Polyakov sector decouples from the fermionic sector and the deconfinement
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Temperature dependence of $\sigma/f_\pi$ (the dotted lines with filled symbols), $L$ (the empty symbols) and $L^*$ (the solid lines) for a static system at vanishing chemical potential, for small system sizes. The pure glue line corresponds to the minimum of the Polyakov-loop potential determined by Eq. (2.13), in the absence of fermions.</figcaption></figure>

transition takes place as predicted by the pure glue theory discussed in Sec. 2.2, namely $L = L^* = 0$ for $T < T_0$, while for $T > T_0$, $L = L^* = L_+$ introduced in Eq. (2.14). This solution is shown with a red dashed line in Fig. 6 and it can be seen that it approximates very well the results obtained at $R = 0.2$ fm. As $R$ is increased, the transition temperature decreases, but the deconfinement transition remains of first order, until a critical radius is reached ($R_c \simeq 0.45$ fm for $T_0 = 0.27$ GeV). When the radius is further increased, the deconfinement transition becomes of crossover type. Remarkably, the chiral restoration transition takes place at significantly larger temperatures. The pseudocritical temperatures $T_\sigma$ for the chiral restoration transition, obtained for $R = 0.2, 0.3, 0.45$ and $0.7$ fm, are 0.476, 0.377, 0.310 and 0.262 GeV, respectively.

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Phase diagram in the $T-\mu$ plane for (a) the chiral restoration and (b) deconfinement at zero angular velocity for different radii $R$. The solid and dashed line segments mark the positions of, respectively, the first-order phase transitions and crossovers separated by the endpoints (the black-filled circles). The black lines in panel (b) indicate the chiral transition lines.</figcaption></figure>

We now investigate the effect of the boundary on the phase diagram of the system.
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Phase diagram in the T-μ plane with the radius R = 0.4 fm for (a) the chiral (purple) and deconfinement (green and blue for L and L*) transitions in the PLSM_q model and the chiral transition (orange) in the LSM_q model at zero angular velocity and (b) the deconfinement transition for various angular velocities in the PLSM_q model. The dashed lines denote a crossover transition and the solid lines indicate a first-order phase transition, while the critical points separating these two regimes are indicated by the black-filled circles. Note that, in panel (a), the deconfinement transition lines for L and L* coincide with each other and they both have two endpoints. Panel (b) shows a smaller range of chemical potential. The low-μ critical point touches the vertical axis when ΩR ≈ 0.5345 and disappears at higher angular velocities.</figcaption></figure>

Figure 7 shows the phase diagram in the T-μ plane indicating separately (a) the chiral symmetry restoration and (b) the deconfinement transition lines. Starting from R = 1 fm, the transition lines for the bounded system approach gradually the transition of the unbounded system (indicated as R = ∞) as the radius R increases. Generally, the critical point separating the crossover and first-order phase transition regions travels towards higher temperatures and smaller chemical potentials as the system size is reduced.

As can be noticed from Figs. 5 and 6, as the radius of the system R decreases, the chiral pseudocritical temperature $T_σ$ increases more rapidly than the deconfinement temperature $T_L$, leading to a visible split between the two transitions. In Ref. [47], we have shown that this splitting between the transitions is an effect of the boundaries that become significant only in small transverse volumes. The splitting becomes insignificant when R is increased.

## 4.2 System under rotation

We now turn on rotation. Figures 9 and 10 show the effect of rotation on the temperature dependence of the σ expectation value for various system sizes (R = 1, 2, 3 and 5 fm), at μ = 0 and μ = 0.25 GeV. It can be seen that for all these parameters, rotation promotes the chiral symmetry restoration, leading to a decrease of σ as ΩR is increased at fixed R, μ and T. This result was reported in the recent literature [61], albeit for a system without boundaries. Furthermore, panels (b)–(d) of Fig. 10 indicate that increasing ΩR pushes the system towards a first-order phase transition, similar to the effect of increasing the radius R. In the case of a system of radius R = 5 fm at μ = 0.25 GeV, chiral symmetry is restored at any temperature when ΩR ≳ 0.9. Similarly, we can expect that increasing ΩR leads to
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Normalized scalar condensate $\sigma/f_\pi$, at $\mu = 0$, as a function of temperature $T$, for different angular velocities and various system sizes: (a) $R = 1$, (b) 2, (c) 3 and (d) 5 fm.</figcaption></figure>

an increase of the expectations of the traced Polyakov loop $L$ and its conjugate, $L^*$, leading to a decrease of the deconfinement transition temperature.

We now consider the effect of rotation on the phase diagram of the model. Considering the system size $R$ and the rotation parameter $\Omega$ fixed, we constructed the phase diagram by considering the variations of the control parameters $\sigma$, $L$ and $L^*$ with respect to $T$ and $\mu$. Our results are summarized in Fig. 11 for systems of sizes $R = 1$ fm (a), 2 fm (b) and 5 fm (c) fm. As expected, as $\Omega R$ is increased, the transition line moves inwards towards lower temperatures and chemical potentials, in a monotonous way, i.e. without crossings. It can be seen that the system with $\Omega R$ behaves regularly, having a perfectly well-defined phase diagram, even in the causal limit $\Omega R = 1$.

Rotation maintains the phase structure, the system exhibiting a crossover transition at small $\mu$ and a first-order phase transition at large $\mu$, separated by a critical point. The trajectory of the critical point in the $(T, \mu)$ plane as $\Omega R$ is varied at fixed $R$ is shown in Fig. 11(d), for various values of $R$. At small $R$, the critical point travels almost on rectangular segments: first towards decreasing temperatures at fixed chemical potential; then towards lower chemical potentials at fixed $T$. At intermediate values of $R$, this trajectory becomes more diagonal, until eventually, at $R = 10$ fm, it becomes an almost vertical line: the chemical potential at the critical point seems to become independent of rotation. At the same time, as $R$ is increased, the critical point for the non-rotating case ($\Omega R = 0$)
<figure><img src="image_11.png" /><figcaption>Figure 10. Same as Fig. 9 for a non-vanishing chemical potential $\mu = 0.25$ GeV.</figcaption></figure>

approaches that of the static, unbounded system, shown as a red dot in panel (d) of Fig. 11.

The phase diagram for the deconfinement transition is shown in Fig. 12. As already mentioned in Figs. 7(b) and 1(a,b), the deconfinement transition is delayed with respect to the chiral one when the chemical potential exceeds some threshold value. This is confirmed also in the case of rotation: at low values of $\mu$, the deconfinement transition is either simultaneous with the chiral one, or it takes place earlier. At large $\mu$, the deconfinement transition takes place after the chiral one, being a crossover, even when the chiral transition is of first-order type.

Coming back to panels (b)–(d) of Fig. 10, we have seen that increasing $\Omega R$ can turn a crossover transition into a first-order one. To systematically study this effect, we constructed in Fig. 13 the phase diagram in the $T-\Omega R$ plane, computed by considering variations of the control parameters $\sigma$, $L$ and $L^*$ with respect to $T$ and $\Omega R$, at various values of $\mu$, for (a) $R=2$ (b) 3 and (c) 5 fm. As the parameter range of $\Omega R$ is restricted between 0 and 1, the transition line has a finite extent also with respect to temperature, extending down to the transition temperature when $\Omega R = 1$. When the system undergoes a crossover transition, we compute the transition point by maximizing the slope along a line in the $T-\Omega R$ plane, namely

$$ \frac{dp}{dG}(\theta) = \cos\theta \frac{\partial p}{\partial T} + \alpha \sin\theta \frac{\partial p}{\partial \Omega R}, \qquad (4.1) $$

where $G = \sqrt{T^2 + (\Omega R/\alpha)^2}$, with the constant $\alpha$ required for dimensional reasons. For
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Phase diagram for the chiral transition in the T-μ plane, at fixed angular velocities and various sizes: (a) R = 1; (b) R = 2; and (c) R = 5 fm. (d) The trajectory of the critical point in the T-μ plane as a function of ΩR, for various radii R, expressed in fm. The critical point corresponding to the static, unbounded system is indicated by a red dot at the top left of panel (d).</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 12.** Phase diagram for the deconfinement transition, shown with colored lines, computed in the T-μ plane at selected angular velocities, for various system sizes: (a) R = 2 fm and (b) R = 5 fm. The black lines represent the chiral transition lines shown also in Figs. 11(b) and (c), respectively.</figcaption></figure>

definiteness, we take it as $\alpha = 1 \text{ GeV}^{-1}$.
<figure><img src="image_14.png" /><figcaption>**Figure 13.** Phase diagram in the $T-\Omega R$ plane for the chiral transition, computed at selected values of the chemical potential $\mu$, for systems of various sizes: (a) $R=1$; (b) $R=3$; (c) $R=5$ and (d) $R=10$ fm.</figcaption></figure>

We now consider the trajectory of the critical point in the $T-\Omega R$ plane. It can be seen in Fig. 13 that the critical point moves towards the left (towards lower $\Omega R$) as $\mu$ increases. Fig. 13(c) shows that, when $R=5$ fm, the critical point in the $T-\Omega R$ phase diagram exists for a very narrow range of chemical potentials, when $0.16 \text{ GeV} \lesssim \mu \lesssim 0.2 \text{ GeV}$. For values of $\mu$ smaller than this range, the transition is of crossover type for all values of $\Omega R$; while for larger $\mu$, the phase transition becomes of first order for all values of $\Omega R$. This property can be understood by going back to Fig. 11(d), where we have shown the trajectory of the critical point as a function of $\Omega R$ for different values of $R$ in the $T-\mu$ plane. We observe that, as $R$ increases, the trajectories become almost parallel to the temperature axis, thus the range of chemical potential over which the critical point moves as a function of $\Omega R$, at fixed radius $R$, decreases as $R$ is increased.

# 5 Comparison to Tolman-Ehrenfest law

The present section aims to complete our study by clarifying to what extent the Tolman-Ehrenfest picture agrees with the $LSM_q$ and $PLSM_q$ results. Namely, we will discuss the Tolman-Ehrenfest prediction for a rigidly rotating system and present a comparative study
to investigate the range of applicability of such classical ideas to study the phase structure of the (P)LSMq models under rotation.

To this end, we will take a step back and calculate the chiral transition temperature under rotation within our effective model framework, under the assumption of validity of the TE law for the local thermodynamic state for a system under rotation. We begin in Subsec. 5.1 with the formulation of the rotation (P)LSMq model within kinetic theory, under the assumption of the validity of the TE law. We then apply our discussions to the LSMq and PLSMq models for vanishing chemical potential in Subsections 5.2 and 5.3. We discuss the vanishing temperature case in Subsec. 5.4. The numerical validation of the TE predictions in Subsec. 5.5.

## 5.1 Kinetic theory formulation of the rotating system

As is well known, a rigidly-rotating state is described by the following temperature four-vector [62]:

$$ \beta^{\mu} \partial_{\mu} = \frac{1}{T} (\partial_{t} + \Omega \partial_{\varphi}). \qquad (5.1) $$

Writing $\beta^\mu = \beta(x)u^\mu$, we can extract the local temperature $T_\rho = 1/\sqrt{\beta^\mu \beta_\mu}$ and the four-velocity,

$$ T_{\rho} = \Gamma_{\rho} T, \qquad u^{\mu} \partial_{\mu} = \Gamma_{\rho} (\partial_{t} + \Omega \partial_{\varphi}), \qquad \Gamma_{\rho} = (1 - \rho^{2} \Omega^{2})^{-1/2}, \qquad (5.2) $$

where $\rho = \sqrt{x^2 + y^2}$ is the distance to the rotation axis, taken as the z axis for convenience.

At the level of statistical mechanics, rigid rotation can be modeled by replacing the density operator $\hat{\rho} = e^{-\beta(\hat{H}-\mu\hat{Q})}$ by $\hat{\rho} = e^{-\beta(\hat{H}-\Omega\hat{J}^z-\mu\hat{Q})}$, where $\hat{H}$ is the Hamiltonian, $\hat{Q}$ is the charge and $\hat{J}^z$ is the total angular momentum along the rotation axis. For a classical particle with momentum $p^\mu$, the orbital angular momentum along the z axis is $L^z = -yp^x + xp^y$, such that in kinetic theory, implementing rigid rotation amounts to changing the Boltzmann factor as follows:

$$ e^{-\beta \mathcal{E}_{\zeta}} \rightarrow e^{-\beta(p^0 + \Omega y p^x - \Omega x p^y - \zeta \mu)} = e^{-\beta_{\rho}(u_{\mu} p^{\mu} - \zeta \mu_{\rho})}, \qquad (5.3) $$

where $\beta_\rho = 1/T_\rho$ is the local temperature, $\mu_\rho = \mu\Gamma_\rho$ is the local chemical potential (such that $\beta_\rho\mu_\rho = \beta\mu = \text{const}$ [62]) and the four-vector $u^\mu$ is the velocity corresponding to rigid rotation. The above procedure is equivalent to the local thermal equilibrium approximation for the density operator $\hat{\rho}$, discussed in Refs. [42, 44, 45]. Applying this prescription at the level of the free energy, we have

$$ F_q = -\frac{2N_f N_c T}{V} \int_V d^3x \sum_{\zeta=\pm 1} \int \frac{d^3p}{(2\pi)^3} F_\zeta|_{E \rightarrow u_\mu p^\mu / \Gamma_\rho}, \qquad (5.4) $$

where V represents a cylindrical volume of radius R, while $F_\zeta$ is given in Eq. (2.37), with the replacement $\beta\mathcal{E}_\zeta \rightarrow \beta_\rho(u_\mu p^\mu - \zeta\mu_\rho)$.

We now perform a Lorentz transformation to the local rest frame, such that $u^\mu \rightarrow \delta_0^\mu$. The integration measure $d^3p$ is not Lorentz covariant, however, $d^3p/E$ is. Writing $u_\mu p^\mu =$
$g_{\mu\nu}\Lambda^\mu_\alpha\Lambda^\nu_\rho u^\alpha p^\rho$, where $\Lambda^\mu_\alpha u^\alpha = \delta^\mu_0$ is the Lorentz transformation to the rest frame, we have

$$F_q = -\frac{2N_f N_c T}{V} \sum_{\zeta=\pm 1} \int d^3 x \int \frac{d^3 p}{(2\pi)^3 E} p^\mu \Lambda_\mu^0 F_\zeta^\rho, \quad (5.5)$$

where now $F_\zeta^\rho$ is given in Eq. (2.37), with $\beta\mathcal{E}_\zeta \to \beta_\rho(E - \zeta\mu_\rho)$. Substituting $p^\mu\Lambda_\mu^0 = \Gamma_\rho(E + \mathbf{v} \cdot \mathbf{p})$, the part involving the local three-velocity $\mathbf{v} = \mathbf{u}/\Gamma_\rho$ can be seen to make a vanishing contribution, being odd with respect to $\mathbf{p} \to -\mathbf{p}$, such that

$$F_q = -\frac{4N_f N_c}{R^2} \int_0^R d\rho \, \rho T_\rho \sum_{\zeta=\pm 1} \int \frac{d^3 p}{(2\pi)^3} F_\zeta^\rho, \quad (5.6)$$

with $T_\rho = \Gamma_\rho T$. Introducing the integration variable $x = p/T_\rho$, we have

$$F_q = -\frac{4N_f N_c}{R^2} \int_0^R d\rho \, \rho T_\rho^4 \sum_{\zeta=\pm 1} \int_0^\infty \frac{dx \, x^2}{2\pi^2} F_\zeta^\rho, \quad (5.7)$$

where $F_\zeta^\rho = \frac{1}{3}\ln(1 + 3L_\zeta e^{-y_\zeta^\rho} + 3L_{-\zeta} e^{-2y_\zeta^\rho} + e^{-3y_\zeta^\rho})$ and

$$\mathcal{Y}_\zeta^\rho = y_\rho - \zeta\beta\mu, \quad y_\rho = \frac{E}{T_\rho} = \sqrt{x^2 + z_\rho^2}, \quad z_\rho = \frac{g\sigma}{T_\rho}. \quad (5.8)$$

We now apply the same prescription to the fermion condensate and to $\partial F_q/\partial L_\pm$, keeping in mind that both the condensate $\sigma$ and the Polyakov-loop values $L_\pm$ are coordinate-independent quantities (note that $L_+ \equiv L$ and $L_- \equiv L^*$). The saddle point Eqs. (2.41) become

$$\frac{\partial V_M}{\partial \sigma} + \frac{2N_f N_c g^2 \sigma}{\pi^2 R^2} \int_0^R d\rho \, \rho T_\rho^2 \sum_{\zeta=\pm 1} \int_0^\infty \frac{dx \, x^2}{y_\rho} f_\zeta(\mathcal{Y}_\rho) = 0, \quad (5.9a)$$

$$\frac{\partial V_L}{\partial L_\pm} - \frac{2N_f N_c}{\pi^2 R^2} \int_0^R d\rho \, \rho T_\rho^4 \sum_{\zeta=\pm 1} \int_0^\infty dx \, x^2 f_{\zeta, \pm}(\mathcal{Y}_\rho) = 0. \quad (5.9b)$$

The purpose of this section is to study the effect of rotation, as captured by the TE law, on the phase transition of the (P)LSM$_q$ models. To simplify the analysis, we will consider only the limiting cases of vanishing quark chemical potential ($\mu = 0$), discussed in Subsections 5.2 and 5.3 for the LSM$_q$ and PLSM$_q$ models, respectively, and vanishing temperature, considered in Subsec. 5.4.

## 5.2 Vanishing chemical potential: LSM$_q$ model

We now consider the LSM$_q$ model at a vanishing chemical potential. The omission of the Polyakov-loop degrees of freedom is equivalent to setting $L = L^* = 1$ and taking only the saddle-point equation for the $\sigma$ meson. Substituting the meson potential (2.2) into Eq. (5.9a) gives us:

$$\lambda(\sigma^2 - v^2)\sigma - h + g\langle\bar{\psi}\psi\rangle_R = 0, \quad \langle\bar{\psi}\psi\rangle_R = \frac{4N_f N_c g\sigma}{\pi^2 R^2} \int_0^R d\rho \, \rho T_\rho^2 \int \frac{dx \, x^2}{y_\rho(e^{y_\rho} + 1)}, \quad (5.10)$$
where we took into account that $f_ζ(\mathcal{Y}_ρ) = (e^{y_ρ} + 1)^{-1}$ is the distribution for a neutral Fermi gas. The above form of the saddle point equation serves to indicate that the difference between the rotating and non-rotating systems comes from the evaluation of the fermion condensate. At large temperatures, or equivalently, small $z_ρ = gσ/T_ρ$, we have $y_ρ \simeq x$ and the fermion condensate can be approximated by

$$ \langle \bar{\psi} \psi \rangle_R \approx \frac{N_f N_c g \sigma T^2}{\pi^2} \langle \Gamma_\rho^2 \rangle_R I_\sigma, \quad I_\sigma = 2 \int_0^\infty \frac{dx \, x}{e^x + 1} = \frac{\pi^2}{6}, \qquad (5.11) $$

where we neglected terms of order $O[(gσ/T)^3]$. In this approximation, the average over the cylindrical volume affects only the Lorentz factor and gives

$$ \langle \Gamma_{\rho}^{2} \rangle_{R} = \frac{2}{R^{2}\Omega^{2}} \ln \Gamma_{R}, \quad \Gamma_{R} = (1 - R^{2}\Omega^{2})^{-1/2}. \qquad (5.12) $$

Coming back to the saddle point equation, we replace $σ \to σ_c$ and the volume average in Eq. (5.10), thus obtaining

$$ \lambda(\sigma_c^2 - v^2)\sigma_c - h + g^2\sigma(T_\sigma^{\Omega R})^2\langle\Gamma_\rho^2\rangle_R = 0, \qquad (5.13) $$

where we used $N_f = 2$ and $N_c = 3$ to simplify the expression. From the form of the above equation, it is reasonable to expect that the conditions for the crossover transition are met at a value $σ_c$ similar to that encountered in the static case (see Subsec. 2.5), when the effective quark mass is about half of its constituent mass:

$$ m_{q;c} = g\sigma_c \approx 0.45gf_\pi \approx 0.138 \text{ GeV}. \qquad (5.14) $$

Applying the same approximations to the static case, $σ_c$ satisfies the equation

$$ \lambda(\sigma_c^2 - v^2)\sigma_c - h + g^2\sigma_c(T_\sigma^{\text{nr}})^2 = 0, \qquad (5.15) $$

with $T_\sigma^{\text{nr}} = 0.147$ GeV being the pseudocritical temperature for the non-rotating (“nr”) system, as shown in the previous sections. Even though the ratio $m_{q;c}/T_\sigma^{\text{nr}}$ is not negligible, we will see that our approximations are still viable. Comparing Eqs. (5.13) and (5.15), we conclude that in the rotating case, the pseudocritical temperature $T_\sigma^{\Omega R}$ in a rotating system is linked to the pseudocritical temperature $T_\sigma^{\text{nr}}$ in the static case via

$$ T_{\sigma}^{\Omega R} = \frac{T_{\sigma}^{\text{nr}}}{\sqrt{\langle \Gamma_{\rho}^{2} \rangle_{R}}} = \frac{R\Omega T_{\sigma}^{\text{nr}}}{\sqrt{2 \ln \Gamma_{R}}}. \qquad (5.16) $$

The above function is shown with the dotted black line in Fig. 14.

## 5.3 Vanishing chemical potential: PLSM$_q$ model

Proceeding now to the PLSM$_q$ model, we consider the case of vanishing chemical potential, when $L = L^*$ and the distributions $f_ζ$ and $f_{ζ;\pm}$, introduced in Eqs. (2.40), become

$$ f_{\zeta}(y) = \frac{Le^{-y} + 2Le^{-2y} + e^{-3y}}{(e^{-y} + 1)(e^{-y} + e^{i\vartheta})(e^{-y} + e^{-i\vartheta})}, $$

$$ \sum_{\zeta=\pm 1} f_{\zeta;\pm}(y) = \frac{e^{-y}}{(e^{-y} + e^{i\vartheta})(e^{-y} + e^{-i\vartheta})}, \qquad (5.17) $$
<figure><img src="image_15.png" /><figcaption>**Figure 14.** Tolman-Ehrenfest prediction in the (P)LSM<sub>q</sub> models. Panel (a) shows, at vanishing chemical potential, $\mu = 0$, the pseudocritical temperatures $T_\sigma$ in the LSM<sub>q</sub> (dashed line) and PLSM<sub>q</sub> (dot-dashed line) models. The Polyakov loop, calculated in the PLSM<sub>q</sub> model at the transition point, $L_c \equiv L(T = T_c)$, is shown by the solid line. Panel (b) shows the critical chemical potential $\mu_c$ and the sigma meson $\sigma_c/f_\pi$ measured after chiral symmetry restoration, for a system rotating with angular velocity $\Omega$ at vanishing temperature, $T = 0$, in both models. All quantities in both panels are shown for systems confined in a fictitious cylinder of radius $R$, as a function of the Lorentz factor $\Gamma_R = 1/(1 - R^2\Omega^2)^{-1/2}$ of a co-rotating observer sitting on the system edge.</figcaption></figure>

where $y \equiv y_\rho$ is introduced in Eq. (5.8) and the phase $\vartheta$ satisfies

$$ (e^{-y} + e^{i\vartheta})(e^{-y} + e^{-i\vartheta}) = e^{-2y} + (3L-1)e^{-y} + 1, \qquad (5.18) $$

such that

$$ \sin \vartheta = \frac{\sqrt{3}}{2} \sqrt{(1-L)(1+3L)}, \qquad \cos \vartheta = \frac{3L-1}{2}. \qquad (5.19) $$

The angle $\vartheta$ decreases from $2\pi/3$ when $L=0$ (confined state) towards 0 as $L \to 1$ (deconfined state). We can thus split the fraction appearing in Eq. (5.17) as

$$ f_{\varsigma}(y) = \frac{1}{3} \left( \frac{1}{e^y + 1} + \frac{1}{e^{y+i\vartheta} + 1} + \frac{1}{e^{y-i\vartheta} + 1} \right), $$

$$ \sum_{\varsigma=\pm 1} f_{\varsigma;\pm}(y) = \frac{i}{2 \sin \vartheta} \left( \frac{1}{e^{y+i\vartheta} + 1} - \frac{1}{e^{y-i\vartheta} + 1} \right). \qquad (5.20) $$

Focusing now on the saddle point equation for the $\sigma$ meson and repeating the steps discussed in Sec. 5.2, we obtain:

$$ \langle \bar{\psi} \psi \rangle_R = \frac{N_f N_c g \sigma T^2}{\pi^2} \langle \Gamma_{\rho}^2 \rangle_R I_{\sigma}, \qquad I_{\sigma} = \sum_{\varsigma=\pm 1} \int_0^{\infty} dx x \frac{dx^2}{y_{\rho}} f_{\varsigma}(y_{\rho}). \qquad (5.21) $$

The term $I_\sigma$ can be evaluated in the limit of negligible mass as

$$ I_{\sigma} \approx \frac{2}{3} \int_{0}^{\infty} dx \, x \left( \frac{1}{e^x + 1} + \frac{1}{e^{x+i\vartheta} + 1} + \frac{1}{e^{x-i\vartheta} + 1} \right) = \frac{\pi^2}{18} - \frac{2}{3} \left[ \mathrm{Li}_2(-e^{i\vartheta}) + \mathrm{Li}_2(-e^{-i\vartheta}) \right], \qquad (5.22) $$
where $Li_n(z) = \sum_{k=1}^{\infty} z^k/k^n$ is the polylogarithm function, satisfying

$$ \int_0^\infty \frac{dx \, x^n}{e^{x+a} + 1} = -n! \mathrm{Li}_{n+1}(-e^{-a}), \qquad (5.23) $$

Using the property $Li_2(-e^a) + Li_2(-e^{-a}) = -\frac{\pi^2}{6} - \frac{a^2}{2}$, with $a = i\vartheta$, we find

$$ I_{\sigma} = \frac{\pi^2}{6} \left( 1 - \frac{2\vartheta^2}{\pi^2} \right). \qquad (5.24) $$

Compared to the LSM model, the saddle point equation for $\sigma$ (5.13) gets modified to

$$ \lambda(\sigma_c^2 - v^2)\sigma_c - h + g^2\sigma_c(T_{\sigma}^{\Omega R})^2\langle\Gamma_{\rho}^2\rangle_R \left(1 - \frac{2\vartheta_{\Omega R}^2}{\pi^2}\right) = 0, \qquad (5.25) $$

where the angle $\vartheta_{\Omega R}$ sets the value of the Polyakov loop (5.19) in the rotating system. The above equation establishes a connection between the pseudocritical temperature and Polyakov loop in the PLSM model, compared to the pseudocritical temperature of the LSM model:

$$ (T_{\sigma}^{\mathrm{nr}})_{\mathrm{PLSM}} = \frac{(T_{\sigma}^{\mathrm{nr}})_{\mathrm{LSM}}}{\sqrt{1 - 2\vartheta_{\mathrm{nr}}^2/\pi^2}}. \qquad (5.26) $$

Knowing that $T_{\sigma}^{\mathrm{nr}} = 0.147$ GeV in the LSM model and $T_{\sigma}^{\mathrm{nr}} = 0.212$ GeV in the Polyakov-loop enhanced model, we can estimate

$$ \vartheta_{\mathrm{nr}} = \frac{\pi}{\sqrt{2}} \sqrt{1 - \left(\frac{(T_{\sigma}^{\mathrm{nr}})_{\mathrm{LSM}}}{(T_{\sigma}^{\mathrm{nr}})_{\mathrm{PLSM}}}\right)^2} \approx 1.6 \, \mathrm{rad} \Rightarrow L_{\mathrm{nr}} \approx 0.313 \qquad (5.27) $$

The value estimated above differs from the actual values $L_{\mathrm{nr}} \approx 0.289$ by less than 10%. This discrepancy can be expected to arise due to the major difference between the ratio $m_q/T_{\mathrm{pc}}$ of the effective quark mass $m_q = g\sigma f_{\pi}$ and the pseudocritical temperature $T_{\mathrm{pc}}$, as $\sigma_c$ has roughly the same value in both the LSM and PLSM models, while $T_{\mathrm{pc}}$ exhibits a significant variation. In any case, it becomes clear that the estimation of the pseudocritical temperature $T_{\sigma}$ for the chiral phase transition depends crucially on estimating the Polyakov loop $L$, for which we must employ Eq. (5.9b).

For the Polyakov loop part [Eq. (5.9b)], we can compute the first mass correction by writing $f_{\zeta}(y_{\rho}) = f_{\zeta}(x) + (df_{\zeta}/dy_{\rho})_{y_{\rho}=x} \times (y_{\rho} - x) + \dots$, with $y_{\rho} - x \approx z_{\rho}^2/2x$, i.e.:

$$ \int dx \, x^2 f_{\zeta}(y_{\rho}) = \int dx \left[ x^2 f_{\zeta}(x) - \frac{z_{\rho}^2}{2} f_{\zeta}(x) + O(z_{\rho}^4) \right]. \qquad (5.28) $$

Thus, Eq. (5.9b) becomes

$$ \frac{\partial V_L}{\partial L_{\pm}} - \frac{N_f N_c T^4}{\pi^2} \left( \langle \Gamma_{\rho}^4 \rangle_R I_{L;2} - \frac{g^2 \sigma^2}{2T^2} \langle \Gamma_{\rho}^2 \rangle_R I_{L;0} \right) = 0, \qquad (5.29) $$

where

$$ \begin{aligned} I_{L;2} &= \sum_{\zeta=\pm 1} \int_0^\infty dx \, x^2 f_{\zeta;\pm} = -\frac{i}{\sin \vartheta} \left[ \mathrm{Li}_3(-e^{-i\vartheta}) - \mathrm{Li}_3(-e^{i\vartheta}) \right], \\ I_{L;0} &= \sum_{\zeta=\pm 1} \int_0^\infty dx \, f_{\zeta;\pm} = -\frac{i}{2 \sin \vartheta} \left[ \mathrm{Li}_1(-e^{-i\vartheta}) - \mathrm{Li}_1(-e^{i\vartheta}) \right]. \end{aligned} \qquad (5.30) $$
Using the relations

$$ \mathrm{Li}_3(-e^{-a}) - \mathrm{Li}_3(-e^a) = \frac{\pi^2 a}{6} + \frac{a^3}{6}, \qquad \mathrm{Li}_1(-e^{-a}) - \mathrm{Li}_3(-e^a) = a, \qquad (5.31) $$

we find

$$ I_{L;2} = \frac{\pi^2 \vartheta}{6 \sin \vartheta} \left( 1 - \frac{\vartheta^2}{\pi^2} \right), \qquad I_{L;0} = \frac{\vartheta}{2 \sin \vartheta}. \qquad (5.32) $$

In the PLSM<sub>q</sub> model, both the transition temperature T<sub>σ</sub><sup>OR</sup> and the Polyakov loop L must be determined simultaneously and self-consistently. This amounts to solving the coupled system of equations

$$ T_{\sigma}^{\Omega R} = \frac{T_{\sigma}^{\mathrm{nr}} \sqrt{1 - 2\vartheta_{\mathrm{nr}}^2 / \pi^2}}{\sqrt{\langle \Gamma_{\rho}^2 \rangle_R \sqrt{1 - 2\vartheta_{\Omega R}^2 / \pi^2}}}, $$

$$ \left. \frac{\partial V_L}{\partial L_{\pm}} \right|_{\substack{L=L^*=L_{\Omega R} \\ T=T_{\sigma}^{\Omega R}}} = \frac{N_f \vartheta_{\Omega R} (T_{\sigma}^{\Omega R})^4}{2 \sin \vartheta_{\Omega R}} \left[ \left(1 - \frac{\vartheta_{\Omega R}^2}{\pi^2}\right) \langle \Gamma_{\rho}^4 \rangle_R - \frac{3g^2 \sigma_c^2}{4\pi^2 (T_{\sigma}^{\Omega R})^2} \langle \Gamma_{\rho}^2 \rangle_R \right], \qquad (5.33) $$

where ⟨Γ<sub>ρ</sub><sup>2</sup>⟩<sub>R</sub> = 2 ln Γ<sub>R</sub>/(R<sup>2</sup>Ω<sup>2</sup>) and ⟨Γ<sub>ρ</sub><sup>4</sup>⟩<sub>R</sub> = Γ<sub>R</sub><sup>2</sup>. We employ as input parameters the transition temperature T<sub>σ</sub><sup>nr</sup> = 0.212 GeV in the non-rotating case and the sigma meson σ<sub>c</sub> = 0.45f<sub>π</sub> at the crossover transition. We determine the non-rotating Polyakov loop parameter self-consistently, giving L<sub>nr</sub> = 0.283, in good agreement with the numerical results obtained in Sec. 2. The derivative ∂V<sub>L</sub>/∂L<sub>±</sub> can be found in Eq. (2.43b).

Figure 14(a) shows, in addition to the pseudocritical temperature T<sub>σ</sub> in the LSM<sub>q</sub>, the dependence of the pseudocritical chiral restoration temperature T<sub>σ</sub> and of the Polyakov loop L computed at T = T<sub>σ</sub> in the PLSM<sub>q</sub>, thus allowing us to make a comparison between the predictions of these two models below. The quantities are shown as functions of the Lorentz factor Γ<sub>R</sub> = (1 − R<sup>2</sup>Ω<sup>2</sup>)<sup>−1/2</sup> according to this Tolman-Ehrenfest prediction. When Γ<sub>R</sub> ≲ 10, the decrease of T<sub>σ</sub> with Γ<sub>R</sub> is visibly slower in the PLSM<sub>q</sub> model than in the LSM<sub>q</sub> model. This slowing down occurs because L increases from its value L<sub>nr</sub> without rotation towards 1, while the angle ϑ<sub>ΩR</sub> decreases towards 0. It is interesting to remark that, when Γ<sub>R</sub> increases above ≃ 5, the Polyakov loop L<sub>c</sub> at the chiral transition becomes increasingly large, signaling that the system is already deconfined when chiral restoration takes place. The decrease of the deconfinement temperature below the chiral restoration one for fast rotation is consistent with the results reported in Fig. 6 of Ref. [47] in the case when ΩR = 1.

## 5.4 Finite-density matter at vanishing temperature

As already discussed in Sections 2.7 and 3.5, the Polyakov sector plays no role when T=0. Therefore, both PLSM<sub>q</sub> and LSM<sub>q</sub> models give the same predictions at vanishing temperature. In this case, the local fermion condensate can be evaluated as it was done in Eq. (2.61). Hence, the averaged fermion condensate is

$$ \langle \bar{\psi} \psi \rangle_R = \frac{N_f N_c g \sigma}{\pi^2 R^2} \int_0^R d\rho \, \rho \left[ p_f^\rho |\mu_\rho| - (g\sigma)^2 \mathrm{artanh} \frac{p_f^\rho}{|\mu_\rho|} \right] \theta(|\mu_\rho| - g\sigma), \qquad (5.34) $$
where $\mu_\rho = \Gamma_\rho\mu$ and $p_f^\rho = \sqrt{\mu_\rho^2 - g^2\sigma^2}$ is the local Fermi momentum. Switching the integration variable to $\zeta = (g\sigma/\mu_\rho)^2$ and using $d\zeta = -2(g\sigma\Omega/\mu)^2 d\rho d\rho$ leads to

$$ \langle \bar{\psi}\psi \rangle_R = \frac{N_f N_c \mu^2 g \sigma}{2\pi^2 \Omega^2 R^2} \theta(1 - \zeta_R) \int_{\zeta_R}^{\zeta_{\max}} d\zeta \left( \frac{1}{\zeta} \sqrt{1 - \zeta} - \operatorname{artanh}\sqrt{1 - \zeta} \right), \quad (5.35) $$

with

$$ \zeta_{\max} = \begin{cases} g^2 \sigma^2 / \mu^2, & \mu > g\sigma, \\ 1, & \text{otherwise.} \end{cases} \quad (5.36) $$

The above integral can be performed analytically, leading to the following saddle point equation:

$$ \lambda(\sigma^2 - v^2)\sigma - h + \theta(1 - \zeta_R) \frac{N_f N_c \mu^2 g^2 \sigma}{2\pi^2 \Omega^2 R^2} \left( 3\sqrt{1 - \zeta} - (2 + \zeta)\operatorname{artanh}\sqrt{1 - \zeta} \right)_{\zeta_R}^{\zeta_{\max}} = 0. \quad (5.37) $$

As the system undergoes a first-order phase transition, the above equation typically allows for three simultaneous roots. The physical solution must minimize the total free energy $F = \frac{\lambda}{4}(\sigma^2 - v^2) - h\sigma + F_q$, where the contribution from the quark sector $F_q$ can be computed by inserting the pressure $P_q = -F_q(\rho)$ given in Eq. (2.61) into Eq. (5.6):

$$ \begin{aligned} F_q &= -\frac{N_f N_c}{12\pi^2 R^2} \int_0^R d\rho \, \rho \left[ |\mu_\rho| p_f^\rho (2\mu_\rho^2 - 5g^2\sigma^2) + 3(g\sigma)^4 \operatorname{artanh} \frac{p_f^\rho}{|\mu_\rho|} \right] \theta(|\mu_\rho| - g\sigma) \\ &= -\theta(1 - \zeta_R) \frac{N_f N_c \mu^2 g^2 \sigma^2}{24\pi^2 \Omega^2 R^2} \left[ 3(4 + \zeta)\operatorname{artanh}\sqrt{1 - \zeta} - \frac{\sqrt{1 - \zeta}}{\zeta}(2 + 13\zeta) \right]_{\zeta_R}^{\zeta_{\max}}. \end{aligned} \quad (5.38) $$

Figure 14(b) shows the dependence of the critical chemical potential $\mu_c$ at which the chiral symmetry gets restored along with the mean-field value of the scalar condensate $\sigma_c$, computed above the phase transition, with respect to the Lorentz factor $\Gamma_R = (1 - R^2\Omega^2)^{-1/2}$. It can be seen that the scalar condensate $\sigma_c$ increases as $\Omega R$ becomes larger according to this Tolman-Ehrenfest prediction.

## 5.5 Tolman-Ehrenfest predictions vs. numerical results

We are now in a position to address the validity of the Tolman-Ehrenfest (TE) law. As highlighted in Sec. 4.1, boundary effects can have a significant influence when the dimensionless product $RT_c$ is of order unity. Therefore, we can expect the TE law to become applicable in the limit where boundary effects become negligible, i.e., in the limit of large system size, $R \to \infty$.

When the system size increases, $R \to \infty$, a meaningful comparison with the equivalent system introduced in Sec. 5.1 can be drawn by enforcing that the product $\Omega R$ is kept constant. This implies that the angular velocity decreases, $\Omega \to 0$. Consequently, the quantum corrections to the thermal expectation values in the rigidly-rotating system, which are of order $O(\Omega^2)$ (see, e.g., Refs. [42, 63]), become negligible.

Figure 15 shows —with the horizontal dashed lines— the TE prediction for the $LSM_q$ (panel a) and the PLSM_q (panel b) models, obtained by solving Eqs. (5.13) and (5.33),
<figure><img src="image_16.png" /><figcaption>**Figure 15.** The pseudocritical temperature $T_σ$ of the chiral crossover transition for (a) LSM$_q$ and (b) PLSM$_q$, as a function of the system radius $R$. The dotted black lines represent the Tolman-Ehrenfest prediction (5.33).</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 16.** The critical chemical potential $\mu_c$ for the chiral transition, shown as a function of system size $R$, at vanishing temperature $T = 0$, for various values of $\Omega R$. The dashed lines correspond to the TE prediction obtained in Subsec. 5.4. The results are valid for both LSM$_q$ and PLSM$_q$.</figcaption></figure>

respectively. It can be seen that the numerical results for the pseudocritical temperature $T_σ$ of the chiral symmetry restoration transition at $\mu = 0$ and fixed $\Omega R$, indeed converge, in both models, to the TE predictions as the system size $R$ is increased.

Finally, Fig. 16 shows the critical chemical potential $\mu_c$, computed at vanishing temperature and for various values of $\Omega R$, as a function of the system size $R$. In this limit, both models give the same results. It can be seen that $\mu_c$ approaches the TE prediction (shown with horizontal dashed lines) as $R$ increases.
<figure><img src="image_18.png" /><figcaption>Figure 17. Scaled moment of inertia, $\tilde{IR}^2$, as a function of $(1 - \Omega R)^{-1}$, for systems of various radii, (a) $R = 1$ fm; (b) $R = 3$ fm; (c) $R = 5$ fm; and (d) $R = 10$ fm at vanishing chemical potential $\mu = 0$. The horizontal arrows indicate the $\Omega R \to 1$ limits, while the dashed lines indicate the small-$\Omega R$ limit in Eq. (6.8). In panel (d), the TE prediction for $T = 0.16$ GeV is shown with black dashed lines for comparison, see text for details.</figcaption></figure>

# 6 Mechanical properties of rotating plasma

## 6.1 Moment of inertia

The mechanical response of the system at finite angular velocity $\Omega$ can be analyzed in terms of the moment of inertial associated with the volume-averaged angular momentum $\tilde{L}$, defined as

$$ \tilde{L} = -\frac{\partial \tilde{F}}{\partial \Omega} = \tilde{I}(\Omega)\Omega, \qquad (6.1) $$

where $\tilde{F}$ is the volume-averaged grand potential of the rotating system. In our notations, we use tildas “~” to denote the quantities computed on the basis of the eigenmode series in the rotating cylinder. These results will later be compared with the prediction of the Tolman-Ehrenfest law, which are used without tildas.

Notice that, in Eq. (6.1), the quantity $\tilde{I}$ corresponds to the volume-averaged moment of inertia, that is the moment of inertia divided by the volume of the cylinder. Conveniently,
<figure><img src="image_19.png" /><figcaption>Figure 18. Same as Fig. 17 for a system with finite quark chemical potential $\mu = 0.25$ GeV.</figcaption></figure>

$\tilde{I}$ does not contain the factor associated with the length of the cylinder in the longitudinal direction.

Since both the meson ($V_M$) and the Polyakov ($V_L$) potentials do not explicitly depend on the angular frequency $\Omega$, the angular momentum of the system receives contributions only from the quark part $\tilde{F}_q$.

Therefore, the definition of the angular momentum (6.1) in particular and the mechanical response of the system to the state of rotation in general involves only the density of the quark grand potential, $\tilde{F}_q$, where the dependence on the angular frequency $\Omega$ enters explicitly via Eqs. (3.17)–(3.19).

In the last section, we observed the rotation-induced effects on the phase transition properties of the system. One can then naturally expect the signature of the phase transition to be reflected in the behaviour of the volume-averaged moment of inertia, $\tilde{I}$. We can evaluate $\tilde{I}$ as a function of $T$, $\mu$ and $\Omega$ from its definition:

$$ \tilde{I}(\Omega) = -\frac{1}{\Omega} \frac{\partial \tilde{F}_q}{\partial \Omega} = \frac{2N_f N_c}{\pi^2 R^2 \Omega} \sum_{b, \zeta = \pm 1} m \int_0^\infty dk \, \tilde{f}_\zeta, \quad (6.2) $$

where $\tilde{f}_\zeta$ denotes the Fermi-Dirac distribution for the rotating state (3.26) and the sum over the collective quantum number “b” is defined in Eq. (3.16).
Figures 17 ($\mu$ = 0) and 18 ($\mu$ = 0.25 GeV) show the dependence of the dimensionless moment of inertia $\tilde{IR}^2$ with respect to $(1 - \Omega R)^{-1}$, for various system sizes [$R$ = 1, 3, 5 and 10 fm, shown in the panels (a)-(d)], and for a selection of temperatures (shown as different lines with points), chosen such that the system undergoes the chiral symmetry restoration as $\Omega R$ is increased (see Figs. 11 and 13). There is a monotonic increase of $\tilde{IR}^2$ with both $T$ and $\Omega R$ and $\tilde{IR}^2$ varies over several orders of magnitude as the causal limit is approached (shown with horizontal arrows). In the case of $\mu = 0$ (Fig. 17), the phase transition appears as an inflection point, beyond which the increase of $\tilde{IR}^2$ with respect to $(1 - \Omega R)^{-1}$ becomes milder. At $\mu = 0.25$ GeV (Fig. 18), $\tilde{IR}^2$ exhibits a jump to a larger value when the system undergoes a first-order phase transition, while in panel (a), corresponding to a small system of size $R = 1$ fm, the transition is of crossover type and the situation is similar to that seen in Fig. 17(a). These two figures also show the moment of inertia in the causal limit, $\Omega R = 1$, as horizontal arrows. In the case when $R$ is small, the difference between the causal limit and the value corresponding to the last visible point, where $(1 - \Omega R)^{-1} = 100$, corresponds to a factor of 2–10. At $R = 10$ fm, this factor already increases to about $10^2$.

To gain more insight into the large-$\Omega R$ behavior of the volume-averaged moment of inertia $\tilde{I}$, let us compute this quantity from the TE approximation for the quark potential, given in Eq. (5.7). We use the notation $I$ to denote the volume-averaged moment of inertia in the TE approach.

Taking into account the relation,

$$ \frac{\partial F_{\zeta}^{\rho}}{\partial \Omega} = \frac{\partial(\beta_{\rho} E)}{\partial \Omega} \times \frac{\partial F_{\zeta}^{\rho}}{\partial(\beta_{\rho} E)} = \rho^2 \Omega \Gamma_{\rho}^2 \beta_{\rho} E f_{\zeta}^{\rho}, \qquad (6.3) $$

we arrive at

$$ IR^2 = \frac{4N_f N_c T^4}{2\pi^2} \int_0^R d\rho \, \rho^3 \Gamma_{\rho}^6 \int_0^{\infty} \frac{dx \, x^2}{y_{\rho}} \left( \frac{4x^2}{3} + z_{\rho}^2 \right) f_{\zeta}^{\rho}. \qquad (6.4) $$

At large $\Gamma_R = (1 - R^2\Omega^2)^{-1/2}$, the magnitude of the quantity $z_\rho$ becomes negligible and we can estimate the above quantity in the limit of a chirally-restored and deconfined ($L = 1$) system. In this case, the coordinate dependence is contained in the $\rho^3\Gamma_\rho^6$ factor, which can be integrated straightforwardly:

$$ IR^2 \approx \frac{4N_f N_c T^4 R^4}{3\pi^2} \Gamma_R^4 \int_0^\infty \frac{dx \, x^3}{e^x + 1} = \frac{7\pi^2 N_f N_c}{90} (R\Gamma_R T)^4. \qquad (6.5) $$

If the system is not deconfined, then

$$ IR^2 = \frac{7\pi^2 N_f N_c}{90} (R\Gamma_R T)^4 \left( 1 - \frac{20\vartheta^2}{7\pi^2} + \frac{10\vartheta^4}{7\pi^4} \right), \qquad (6.6) $$

where the angle $\vartheta$ is related to the mean-field value of the Polyakov loop through Eq. (5.19).
## 6.2 Shape coefficients

For slow rotation, one can expand the volume-averaged quark grand potential in a power series with respect to the velocity at the boundary $v_R = \Omega R$ as follows [48]:

$$ \tilde{F}_q = \tilde{F}_q^{(0)} \sum_{n=0}^{\infty} \frac{v_R^{2n}}{(2n)!} \tilde{K}_{2n}, \quad \tilde{K}_{2n} = \frac{1}{\tilde{F}_q^{(0)}} \left. \frac{\partial^{2n} \tilde{F}_q}{\partial v_R^{2n}} \right|_{\Omega=0}, \qquad (6.7) $$

with $\tilde{F}_q^{(0)}$ being the quark grand potential in the absence of rotation. In the thermodynamic ground state, this system respects the invariance under the time-reversal transformation, $t \to -t$ which also flips the direction of rotation, $\Omega \to -\Omega$. This invariance, together with analyticity, implies that the quark potential (6.7) is an even function of the angular frequency $\Omega$.

At small $\Omega R$, the volume-averaged moment of inertia can be expressed in terms of the shape coefficients $K_{2n}$ via

$$ \tilde{I} = -R^2 \tilde{F}_q^{(0)} \left[ \tilde{K}_2 + \frac{1}{6} \Omega^2 R^2 \tilde{K}_4 + O(\Omega^4 R^4) \right], \qquad (6.8) $$

where we used Eq. (6.7) together with the definition of $\tilde{I}$ in Eq. (6.2). We remind that $\tilde{F}_q^{(0)} < 0$ in Eq. (6.8).

The terms in the series (6.8) can be intuitively understood. The coefficient $K_2$ can be associated with a normalized moment of inertia in the limit of vanishing rotation, $\tilde{I}_{\Omega=0} = -R^2 \tilde{F}_q^{(0)} K_2$. The coefficient $K_4$ determines the first nontrivial rotation-dependent contribution to the moment of inertia that arises due to rotation. One can understand the appearance of a nonzero $K_4$ coefficient as a result of the inertial (centrifugal) forces that arise in the rotating system: the rigid rotation tends to pull the rotating matter (for example, an element of a fluid) outwards, thus increasing the density of matter at the periphery of rotation and making the moment of inertia larger. The leading contribution of this effect to the moment of inertia is represented by a non-vanishing $K_4$ coefficient in Eq. (6.8). In a free rotating system, the shape coefficients are $K_{2n} = (2n)!$ [48]. The limit of slow rotation (6.8) is shown with the dashed gray lines in Figs. 17 and 18, being in agreement with the results obtained through the direct evaluation of Eq. (6.2).

In the (P)LSM$_q$ model, the coefficients $K_{2n}$ can be obtained by differentiating the quark potential, given in Eq. (3.17), with respect to $v_R = \Omega R$. This derivative is understood to be taken with respect to the angular velocity, as we are looking for the system's response to the increase of angular velocity, $\Omega$, while keeping the system size $R$ constant. Since $\tilde{F}_q$ depends on $\Omega$ only through the argument $\tilde{\mathcal{E}}_\zeta = \mathcal{E}_\zeta - \Omega m$, these derivatives are simply

$$ \tilde{F}_q^{(0)} \tilde{K}_{2n} = \frac{2N_f N_c}{\pi^2 R^{2n+2} T^{2n-1}} \sum_{b, \zeta = \pm 1} m^{2n} \int_0^\infty dk \, f_\zeta^{(2n-1)}, \qquad (6.9) $$

where $f_\zeta^{(n)} \equiv \partial^n f_\zeta / \partial(\beta \mathcal{E}_\zeta)^n$ is a $n$-th derivative, with $n$ a non-negative integer. Note that the above expression is evaluated at vanishing angular velocity. Using Eq. (2.53) to compute $f'_\zeta$, as well as the following relations,

$$ f'_{\zeta\zeta'} = \frac{\partial f_{\zeta\zeta'}}{\partial(\beta \mathcal{E}_\zeta)} = -f_{\zeta\zeta'}(2 - \delta_{\zeta\zeta'} - 3f_\zeta), \qquad (6.10) $$
we obtain

$$ \tilde{F}_q^{(0)} \tilde{K}_2 = \frac{2N_f N_c}{\pi^2 R^4 T} \sum_{b, \zeta} m^2 \int_0^\infty dk f'_\zeta, \qquad (6.11) $$

$$ \tilde{F}_q^{(0)} \tilde{K}_4 = \frac{2N_f N_c}{\pi^2 R^6 T^3} \sum_{b, \zeta = \pm 1} m^4 \int_0^\infty dk \left[ 9f'_\zeta(1 - 3f_\zeta + 3f_\zeta^2 + f'_\zeta) + 2 \sum_{\zeta' = \pm} L_{\zeta'\zeta} f_{\zeta\zeta'}(2 - \delta_{\zeta\zeta'} - 3f_\zeta)(5 - \delta_{\zeta\zeta'} - 9f_\zeta) \right]. \qquad (6.12) $$

The coefficients $K_2$ and $K_4$, determined by Eqs. (6.11) and (6.12), are presented in Figs. 19, 20, 21 and 22 below.

Figure 19 shows that $K_2$ and $K_4$ approach the free, non-interacting limits $K_2 = 2!$ and $K_2 = 4!$ [48] when the temperature increases and the radius becomes large, $R \to \infty$, so that the system reaches the thermodynamic limit. As the temperature is decreased, the coefficients $K_2$ and $K_4$ decrease until the system reaches the chiral phase transition, after which they start increasing. For $\mu = 0$ (upper panels) and $\mu = 0.25$ GeV (lower panels), both $K_2$ and $K_4$ diverge as $T \to 0$ and the system goes to the chirally broken, confined phase.

To gain further insight into the coefficients $K_2$ and $K_4$, we compute them, as well as the quark potential, using the TE approximation, by differentiating $F_q$ in Eq. (5.7) with respect to $\Omega R$:

$$ \begin{aligned} F_q^{(0)} &= -\frac{N_f N_c T^4}{3\pi^2} \sum_{\zeta=\pm 1} \int_0^\infty \frac{dx}{y} \frac{x^4}{y} f_\zeta, \\ F_q^{(0)} K_2 &= -\frac{2N_f N_c T^4}{3\pi^2} \sum_{\zeta=\pm 1} \int_0^\infty \frac{dx}{y} \left(x^2 + \frac{3z^2}{4}\right) f_\zeta, \\ F_q^{(0)} K_4 &= -\frac{N_f N_c T^4}{\pi^2} \sum_{\zeta=\pm 1} \int_0^\infty \frac{dx}{y} (8x^2y^2 + z^4) f_\zeta, \end{aligned} \qquad (6.13) $$

where $y = \sqrt{x^2 + z^2}$ and $z = g\sigma/T$ [see Eq. (5.8) regarding the notation].

In the high-temperature limit, where $z \to 0$ and $L, L^* \to 1$, it is not difficult to see that $K_2 \to 2!$ and $K_4 \to 4!$, as expected [48]. This property is clearly seen in Fig. 20, where these free-field expressions are gradually reached in the large volume limit. The volume-averaged quark potential in the static limit becomes:

$$ F_q^{(0)} = -\frac{7\pi^2 N_f N_c T^4}{180} \left(1 - \frac{20\vartheta^2}{7\pi^2} + \frac{10\vartheta^4}{7\pi^4}\right). \qquad (6.14) $$

Taking into account the first mass correction in the limit of finite but small mass, $z \ll 1$, and writing $F_q^{(0)} = F_{q;0}^{(0)} + z^2 F_{q;2}^{(0)} + O(z^4)$, we have $K_2 F_q^{(0)} = 2! F_{q;0}^{(0)} + z^2 F_{q;2}^{(0)} + O(z^4)$ and $K_4 F_q^{(0)} = 4! F_{q;0}^{(0)} + 8z^2 F_{q;2}^{(0)} + O(z^4)$. Explicitly, one gets:

$$ F_{q;2}^{(0)} = \frac{N_f N_c T^4}{12} \left(1 - \frac{2\vartheta^2}{\pi^2}\right). \qquad (6.15) $$
The first two terms in the first two $K_{2n}$ coefficients are:

$$K_2 = 2! + \frac{15}{7\pi^2} \left(\frac{g\sigma}{T}\right)^2 \Phi(\vartheta), \quad K_4 = 4! + \frac{240}{7\pi^2} \left(\frac{g\sigma}{T}\right)^2 \Phi(\vartheta), \qquad (6.16)$$

where the higher-order corrections are not shown. The function

$$\Phi(\vartheta) = \frac{1 - \frac{2\vartheta^2}{\pi^2}}{1 - \frac{20\vartheta^2}{7\pi^2} + \frac{10\vartheta^4}{7\pi^4}}, \qquad (6.17)$$

is expressed via the variable $\vartheta$ that encodes the value of the expectation value of the Polyakov loop (5.19). Thus, the chiral symmetry breaking, encoded in the condensate $\sigma$, and the effect of confinement of color affect both shape coefficients (6.16).

The angle $\vartheta$ vanishes, $\vartheta = 0$, as the system reaches the asymptotic deconfined state at $L \to 1$ and $\vartheta = 2\pi/3$ in the confined state ($L = 0$). The function (6.17) is a rising, positive-valued function of $\vartheta$ with $1 \le \Phi \le 9$ in this segment of values. Therefore, the chiral symmetry breaking (characterized by a nonzero condensate $\sigma$) and confinement of color (implying nonzero $\vartheta$) imply that both $K_2$ and $K_4$, presented in Eq. (6.16), are larger than their free field values. Both effects are readily understood at the intuitive level, since they imply the mass gap generation in the medium, which enhances its inertial mechanical response.

In the opposite limit of vanishing temperature, $z \gg 1$, the shape coefficients $\tilde{K}_2$ and $\tilde{K}_4$ can be derived from the grand potential of a rotating system at vanishing temperature, given in Eq. (3.30):

$$\tilde{F}_q^{T=0}(0) \tilde{K}_2^{T=0} = -\frac{2N_f N_c}{\pi^2 R^4} \sum_b m^2 \frac{E_{ml}^f}{k_{ml}^f} \theta(|\mu| - M_{ml}), \qquad (6.18a)$$

$$\tilde{F}_q^{T=0}(0) \tilde{K}_4^{T=0} = -\frac{6N_f N_c}{\pi^2 R^6} \sum_b m^4 \frac{M_{ml}^2 E_{ml}^f}{(k_{ml}^f)^5} \theta(|\mu| - M_{ml}), \qquad (6.18b)$$

where $E_{ml}^f = \mu$, $M_{ml} = \sqrt{q_{ml}^2 + g^2\sigma^2}$ and $k_{ml}^f = \sqrt{\mu^2 - M_{ml}^2}$ represent the $\Omega = 0$ limit of the same quantities introduced in Sec. 3.5. From the TE expressions in Eqs. (6.13), we have

$$F_q^{T=0}(0) K_2 = -\frac{N_f N_c}{6\pi^2} p_f^3 E_f \theta(\mu - g\sigma), \quad F_q^{T=0}(0) K_4 = -\frac{N_f N_c}{\pi^2} p_f E_f (p_f^2 + E_f^2) \theta(\mu - g\sigma), \qquad (6.19)$$

where $E_f = |\mu|$ and $p_f = \sqrt{E_f^2 - g^2\sigma^2}$ are the Fermi energy and momentum, respectively.

Figure 21 shows the behaviour of the shape coefficients as functions of temperature at finite chemical potential and finite radius of the system, $R = 5$ fm. At low temperatures, both coefficients demonstrate substantial oscillations in their values as the chemical potential varies in a rather small interval. This seemingly irregular behaviour disappears quickly as temperature increases.

The zero-temperature oscillations can clearly be seen in Fig. 22, where both coefficients, computed from Eqs. (6.18), are shown as functions of the chemical potential $\mu$ for various
<figure><img src="image_20.png" /><figcaption>**Figure 19.** Scaled shape coefficients $2!/\tilde{K}_2$ (left) and $4!/\tilde{K}_4$ (right), shown as functions of temperature for $\mu = 0$ (top) and $\mu = 0.25$ GeV (bottom), for various system sizes $R = 5, 10$ and $40$ fm, along with the TE prediction (6.13), shown using black lines.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Figure 20.** Shape coefficients $\tilde{K}_2$ and $\tilde{K}_4$, as functions of system radius $R$, for high temperature $T = 1$ GeV and vanishing chemical potential $\mu = 0.0$. The dotted lines show the high-temperature thermodynamic limit for these coefficients where the system becomes free, $K_{2n}^{\text{free}} = (2n)!$.</figcaption></figure>

radii $R$ of the system. They irregularly oscillate around the corresponding TE values (shown by the horizontal dashed lines). The reason of the oscillations is rooted in an analogue of the Shubnikov-de Haas effect, that arises due to the modulation of the density of states
<figure><img src="image_22.png" /><figcaption>**Figure 21.** Scaled shape coefficients, $2!/\tilde{K}_2$ (a) and $4!/\tilde{K}_4$ (b), as functions of temperature $T$, for fixed system size ($R = 5$ fm), computed at various chemical potentials $\mu \in \{0.332, 0.334, 0.336, 0.338\}$.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Figure 22.** Scaled shape coefficients, $2!/\tilde{K}_2$ and $4!/\tilde{K}_4$, as functions of the chemical potential $\mu$ at $T = 0$ for various system sizes $R = 1, 3, 5$ and $10$ fm. The horizontal lines show the TE expressions in the thermodynamic limit where the coefficients correspond to the free fields, $K_2 = 2!$ and $K_4 = 4!$. The vertical lines mark the $R$-dependent critical chemical potential above which the system resides in the chirally-restored phase.</figcaption></figure>

at the Fermi level [60]. In our system, the irregular oscillations arise when the chemical potential crosses the next transverse level determined by the transverse mass eigenvalues $M_{ml}$ and associated with the radial wavefunctions in the cylinder. Mathematically, the irregular behaviour appears due to the presence of the Heaviside functions $\theta(|\mu| - M_{ml})$, that emerge as a low-temperature limit of the Fermi-Dirac distribution in Eq. (6.18), as well as due to the negative powers of the longitudinal momentum $k_{ml}$. As the temperature increases, the sharp Fermi surfaces with $|\mu| = M_{ml}$ erode, and the oscillations disappear bringing the system back to the smooth behaviour of Fig. 21.
# 7 Conclusions

In this paper, we studied the effects of the uniform steady rotation on the phase structure of the quark-gluon plasma. We work in the effective model approach represented by the linear sigma model coupled to quarks and to the Polyakov loop ($PLSM_q$). The attractive features of this model are its renormalizability and the presence of the coupling between chiral and confining degrees of freedom. We fixed the coupling constants of the model to their vacuum values, without allowing for heuristic modification of the model parameters that would fine-tune effects of rotation or temperature. We put the system in a cylinder of a finite radius $R$, imposing the spectral boundary conditions on the quark degrees of freedom, in order to protect the system from breaking causality. We adopted the uniform approximation to the ground state, described by the scalar condensate that leads to the breaking of the chiral symmetry, and the Polyakov loop expectation value that encodes the confinement of color.

Our main focus was to study both the boundary and the rotation effects on the phase diagram of the $PLSM_q$ model. We explored a wide parameter range, covering regimes corresponding to different temperatures, various radii of the rotating plasma and different values of the quark chemical potential. We found that, in the studied range of parameters, the model predicts that the rotation decreases bothe (pseudo)critical temperatures of the transitions corresponding to both color deconfinement and chiral symmetry breaking.

For small system sizes, we showed that the suppression of the fermionic sector allows the deconfinement transition to proceed as in the pure-glue case, at the temperature $T = T_0 = 0.27$ GeV. This is a consequence of the fact that the Polyakov sector is not influenced directly by the presence of the boundary. Since the Polyakov potential is insensitive also to the chemical potential $\mu$, as we increased $\mu$ and the fermionic sector became more important, we showed that the first-order pure-glue transition turned into a crossover transition. Thus, we found a small-$\mu$ critical point occuring for small enough systems, separating the first-order and crossover regimes for the deconfinement phase transition. On the other hand, we showed that the chiral transition takes place at a larger temperature, eventually merging into the deconfinement transition when $\mu$ is sufficiently increased. We also found that rotation inhibits the small-$\mu$ endpoint appearing in small systems.

At finite volume, we found that rotation promotes both chiral restoration and deconfinement, leading to a melting of the $\sigma$ meson condensate and an enhancement of the Polyakov loop expectation values. We performed an extensive investigation of the effects of both boundary and rotation on the phase diagram, from small to large systems, while varying the angular velocity between 0 and the causal limit $\Omega R = 1$. In all cases, we found that increasing $\Omega R$ promotes both the deconfinement and the chiral restoration transitions, leading to smaller (pseudo)critical values of the temperature and chemical potential at the transition point. We considered the trajectory of the critical point for the chiral restoration transition as a function of $\Omega R$, for various system sizes. As the system size is increased and $R \to \infty$, the critical point for the static (non-rotating) case smoothly tends to that for the unbounded system. Surprisingly, when $R$ becomes large, the chemical potential corresponding to the critical point becomes independent of $\Omega$ and, as $\Omega R$ is increased, the critical
point travels in the $T-\mu$ plane on an almost vertical line, towards lower temperatures.

In the large system size, we considered a comparison of the numerical QFT results with the Tolman-Ehrenfest prediction of the uniformly rotating systems, by which the local temperature increases as the distance from the axis of rotation increases. At vanishing chemical potential, we were able to find an analytical relation between the pseudocritical chiral restoration temperatures for the rotating and non-rotating systems, validated by comparison to the numerical results in the large volume limit. Within this setup, we were able to reveal that the expectation value of the Polyakov loop $L$ at the pseudocritical temperature of the chiral transition increases to 1 towards the causal limit, as $\Omega R \to 1$. Thus, in the causal limit, the system is already deconfined when the chiral restoration takes place. We also considered within our TE model the dependence of the critical chemical potential in the limit of vanishing temperature. All of the aforementioned results were supported by the numerical data obtained for large system sizes, where boundary effects become subleading.

While the TE picture is consistent with the PLSM$_q$ results for all studied parameter ranges, it is conflicting with lattice results for the purely gluonic plasma [6, 7, 64], as well as for the theory with dynamical fermions [8, 12]. We stress that at high temperatures, where the QGP becomes weakly coupled, the TE behavior should be restored, with only a minor non-perturbative pressure contribution from magnetic gluons [65]. Thus, some of our findings may hold far above the deconfinement transition.

In addition, we analyzed the mechanical properties of rotating quark-gluon matter within PLSM$_q$. We found that the dimensionless moment of inertia experiences a substantial increase at the phase transition line and a strong enhancement deeper in the deconfined phase. These properties can be attributed to the liberation of the degrees of freedom above the deconfinement temperature and a rapid increase of the energy density as temperature increases in the deconfinement phase. The Tolman-Ehrenfest law predicts a rapid rise of the form $I \sim \Gamma_R^4$, with $\Gamma_R = (1 - \Omega^2 R^2)^{-1/2}$ being the Lorentz factor of the boundary, which is confirmed by the numerical results as the velocity of the boundary $\Omega R$ approaches the causal limit.

We also studied the shape coefficients $K_2$ and $K_4$ of the rotating medium. The first quantity, $K_2$, corresponds to the ratio of the moment of inertia to the free energy of the system at vanishing rotation. The second quantity, $K_4$ effectively encodes the change of the shape of the rotating plasma with increasing rotation, due to the mass re-distribution within the system. Due to this particular normalization, these dimensionless quantities characterize the relative significance of the mechanical moment of inertia and the shape of the rotating plasma, factoring out the number of degrees of freedom. In the large-temperature limit, both dimensionless coefficients approach their free-field values, $K_2 = 2!$ and $K_4 = 4!$, in the large-volume limit $R \to \infty$. While, in the deconfinement phase, these quantities are slowly-varying functions of temperature, they experience a knee-like structure at the transition point and rise significantly as the temperature decreases. The latter fact implies that the number of the thermal degrees of freedom (represented by the magnitude of the free energy) decreases faster than these mechanical characteristics of the rotating plasma. Expectedly, as the radial size of the plasma increases, both these quantities
approach the TE limit of the PLSM<sub>q</sub> model in both phases. The latter limit includes also the effects of the fermionic mass and a finite mean-field value of the Polyakov loop that are nonvanishing at both sides of the (pseudo-)critical transition.

Interestingly, we showed analytically that the TE law implies that both the chiral symmetry breaking and the quark confinement lead to the increase of the K<sub>2</sub> and K<sub>4</sub> coefficients in the zero-density plasma (μ = 0). On the contrary, at vanishing temperature, the varying quark chemical potential leads to quasi-periodic oscillations of the mechanical coefficients K<sub>2</sub> and K<sub>4</sub>, that appear as a result of the modulation of the density of states at decrete thresholds of the nested energy levels that resemble the Shubnikov-de Haas oscillations [60].

**Acknowledgments.** We thank P. Aasha for fruitful discussions. 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.

## A The Polyakov potential

In this appendix, we briefly review the polynomial model for the effective Polyakov potential V<sub>L</sub>, proposed in Ref. [33]. Compared to the logarithmic potential [32] considered in the main text, the qualitative outcomes from the polynomial potential are mostly similar, since both models are constructed based on the Z<sub>3</sub> symmetry properties of SU(3) pure gauge theory, their parameters being fine-tuned to reproduce lattice pure gauge thermodynamics. We also discuss the extension allowing for the transition temperature T<sub>0</sub> appearing in the Polyakov potential to become a function of the number of fermion flavors and of the fermion chemical potential, as described in Ref. [27].

In Subsec. A.1, we describe the structure and features of the polynomial potential in the pure-glue theory. In Subsec. A.2, we discuss the models for the transition temperature T<sub>0</sub>.

### A.1 Polynomial Polyakov potential

In the polynomial parametrization, the potential can be written as the following,

$$ \frac{V_L^{\text{P.}}}{T^4} = -\frac{b_2(T)}{2}L^*L - \frac{b_3}{6}[(L^*)^3 + L^3] + \frac{b_4}{4}(L^*L)^2, \qquad (A.1) $$

where b<sub>3</sub> and b<sub>4</sub> are constants, while b<sub>2</sub>(T) is parametrized

$$ b_2(T) = a_0 + a_1 \left(\frac{T_0}{T}\right) + a_2 \left(\frac{T_0}{T}\right)^2 + a_3 \left(\frac{T_0}{T}\right)^3. \qquad (A.2) $$

Imposing the Stefan-Boltzmann limit (assuming L, L* → 1 when T → ∞) implies

$$ \frac{a_0}{2} + \frac{b_3}{3} - \frac{b_4}{4} = \frac{8\pi^2}{45}. \qquad (A.3) $$
Demanding that $V_L$ is minimized at a fixed temperature with respect to both $L$ and $L^*$ gives $L = L^*$ and

$$-\frac{L}{2}[b_2(T) + b_3L - b_4L^2] = 0. \qquad (A.4)$$

Imposing that $L = 1 - \delta L$, with $\delta L \to 0$ when $T \to \infty$, we find another constraint involving $a_0, b_3$ and $b_4$, as well as the large-temperature behaviour of $L$:

$$a_0 + b_3 - b_4 = 0, \qquad L \approx 1 - \frac{a_1}{b_3 - 2b_4} \frac{T_0}{T}. \qquad (A.5)$$

As opposed to the case of the logarithmic potential, $L$ approaches 1 at a slower rate. The approach is made from below if $a_0/(b_3 - 2b_4) > 0$, which is true for the parametrization that we will consider here (see below).

We can see that Eq. (A.4) gives $L = 0$ as a root, meaning that the confined state always corresponds to a local minimum of the potential. When $T > T_-$, a second minimum develops. This can be found by setting to 0 the term between the square brackets:

$$L^{\pm} = \frac{b_3}{2b_4} \left( 1 \pm \sqrt{1 + \frac{4b_2(T)b_4}{b_3^2}} \right). \qquad (A.6)$$

The temperature $T_-$ corresponds to the solution

$$b_2(T_-) = -\frac{b_3^2}{4b_4}, \qquad (A.7)$$

in which case $L^+ = L^- = b_3/2b_4$. As a side remark, since $0 < L^- \le L^+ < 1$, we must have $b_3 < 2b_4$ and then $a_1 < 0$ such that $L$ approaches 1 from below. Substituting now $T = T_-$ and $L = L^* = L^\pm$ into Eq. (A.1) gives

$$V_L^{\text{P}}(T_-, L^{\pm}) = \frac{b_3^4}{192b_4^3} > 0, \qquad (A.8)$$

such that at $T = T_-$, the solution $L = 0$ still corresponds to the global minimum of the Polyakov potential. Demanding now that the first-order phase transition occurs as $T = T_0$ and denoting $b_c \equiv b(T_0) = a_0 + a_1 + a_2 + a_3$, we impose $V_L^{\text{P}}(T_0, L^+) = 0$, such that

$$-\frac{L^{+2}}{12}(3b_c + b_3L^+) = 0. \qquad (A.9)$$

Substituting the solution for $L^+$, we arrive at

$$b_c = a_0 + a_1 + a_2 + a_3 = -\frac{2b_3^2}{9b_4}, \qquad (A.10)$$

which gives the value $L^{+,c} = 2b_3/3b_4$ immediately after the phase transition. The remaining coefficients, namely $a_1, a_2$ and $a_3$, are found by fitting the lattice data. Following Refs. [27, 28, 33],

$$a_1 = -1.95, \quad a_2 = 2.625, \quad a_3 = -7.44. \qquad (A.11)$$
This leads to the following system of equations:

$$ \frac{a_0}{2} + \frac{b_3}{3} - \frac{b_4}{4} = \frac{8\pi^2}{45}, \quad a_0 + b_3 - b_4 = 0, \quad a_0 + \frac{2b_3^2}{9b_4} = -a_1 - a_2 - a_3 \approx 6.765, \quad (A.12) $$

which is solved by $a_0 \approx 6.75$ and $b_3 \approx 0.82$ and $b_4 \approx 7.56$. This also implies a mere $L^{+,c} \approx 0.0722$ after the phase transition. For $T_0$, we take the pure transition temperature, $T_0 = 0.27$ GeV, which is the same as in the case of the logarithmic potential.

One of the major drawbacks of the polynomial potential is that, under the coupling to fermions, the expectation value of Polyakov loop $L$ goes beyond 1 as $T \to \infty$ [27], as can be seen from Fig.23(a), which will effectively imply an unphysical negative free energy for the heavy quark. On the other hand, the logarithmic term in the logarithmic potential for the Polyakov loop ensures that $L$ never exceeds 1, even in the presence of fermions, approaching it asymptotically as the temperature increases [66].

The situation worsens as we turn on the rotation. As can be seen in Fig.23(b), the expectation value of the Polyakov loop ($L_c$) at the chiral symmetry restoration temperature ($T_σ$), estimated according to the Tolmann-Ehrenfest prediction in Eq. (5.33), increases sharply with increasing $Γ_R$ for the polynomial potential, while for the logarithmic potential, the Polyakov loop $L$ asymptotes to 1 as $Γ_R \to \infty$. We thus conclude that the polynomial potential for the Polyakov loop is not suitable for the study of systems undergoing fast rotation.

<figure><img src="image_24.png" /><figcaption>**Figure 23.** (a) Expectation value of the Polyakov loop as a function of temperature at zero chemical potential for the cases of the logarithmic and polynomial Polyakov potentials. (b) The expectation value of the Polyakov loop at the chiral transition temperature, as a function of $Γ_R$, according to the Tolmann-Ehrenfest prediction.</figcaption></figure>

## A.2 Choice of $T_0$

So far, for both the logarithmic and the polynomial potentials, we have considered $T_0 = 0.27$ GeV, corresponding to the confinement-deconfinement transition, measured on the lattice for the SU(3) pure gauge theory. However, in presence of dynamical quarks, one can expect this transition temperature to be modified. In this section, we will discuss the
model proposed in Ref. [27], starting with an $N_f$-dependent, but $\mu$-independent transition temperature in Subsec. A.2.1, and incorporating the $\mu$-dependence in Subsec. A.2.2.

### A.2.1 Constant $T_0$

First, we will consider a constant $T_0 = 0.208$ GeV, as suggested in Ref. [27] for the case of $N_f = 2$ quark flavors. In Fig. 24, we can see that, with $T_0 = 0.208$ GeV, the temperatures for the chiral and confinement transitions decrease for both the logarithmic and the polynomial potentials, compared to the $T_0 = 0.270$ GeV case. However, there is an unphysical splitting between the chiral and confinement transitions in the crossover region for the logarithmic potential at $T_0 = 0.208$ GeV. This behaviour, with the chiral temperature being higher than the deconfinement temperature, is not expected for the non-rotating system. On the other hand, the phase diagram at $T = 0.270$ GeV agrees qualitatively with the phenomenological expectations and, therefore, it is used in the main text of this paper.

For the polynomial potential, no such unphysical splitting can be seen, for either $T_0 = 0.270$ GeV or $0.208$ GeV. However, we still do not consider the polynomial potential in our numerical analysis, due to the limitation that $L$ and $L^*$ are free to increase without bounds.

<figure><img src="image_25.png" /><figcaption>Figure 24. The chiral and confinement phase diagram for two different values of $T_0 = 0.270$ GeV and $T_0 = 0.208$ GeV for logarithmic Polyakov potential (a) and polynomial Polyakov potential (b). The first order transition is denoted by a solid line in both figures.</figcaption></figure>

### A.2.2 Chemical-potential-dependent $T_0$

Lastly, we will consider the chemical-potential-dependent $T_0$, following Ref. [27]. For brevity, we follow Ref. [27] and consider only the polynomial potential. In the presence of quark degrees of freedom, $T_0$ becomes a function of the quark chemical potential $\mu$ and quark number of flavors $N_f$. This modification can be obtained by computing the transition temperature via [27]

$$T_0 = T_{\tau} e^{-1/(\alpha_0 b)}, \qquad (A.13)$$

where $\alpha_0 = 0.304$ is the QCD coupling at the $\tau$ scale given by $T_{\tau} = 1.770$ GeV. The dependence on $\mu$ can be incorporated through the coefficient $b$ [27]. At one-loop level, one
can write [50, 67]

$$b(\mu) = b_0 - b_{\mu} \frac{\mu^2}{(\hat{\gamma}T_{\tau})^2}, \quad b_0 = \frac{1}{6\pi}(11N_c - 2N_f), \quad b_{\mu} \approx \frac{16}{\pi}N_f, \qquad (A.14)$$

For the $N_f = 2$ model considered in this paper, $T_0$ evaluates to 0.208 GeV at $\mu = 0$. The parameter $\hat{\gamma}$ is in principle free, however for definiteness we take it as $\hat{\gamma} = 1$.

The model in Eqs. (A.13)–(A.14) has a discontinuity in the critical temperature $T_0$ when $b(\mu) = 0$. As $\mu \to \mu_\infty = T_\tau \sqrt{b_0/b_\mu}$ from below, $T_0 \to 0$. When $\mu$ approaches $\mu_\infty$ from above, $T_0 \to \infty$. This behavior is clearly unphysical and it can be traced to the validity of Eq. (A.14) only at small $\mu$. However, inducing the phase transition in a bounded system requires a larger $\mu$ and/or $T$ than in the unbounded system, as the boundary generically suppresses fermionic expectation values. Therefore, the study of strongly-interacting matter in the presence of a boundary requires a model that gives a regular behavior of $T_0$ for any value of $\mu$. One may extend the model presented in Eq. (A.14) by replacing $b(\mu)$ by its first-order Padé approximant,

$$b(\mu) = b_0 \left[ 1 + \frac{b_{\mu}}{b_0} \left( \frac{\mu}{\hat{\gamma}T_{\tau}} \right)^2 \right]^{-1}, \qquad (A.15)$$

which has the advantage of giving a finite critical temperature $T_0$ for any value of $\mu$, which tends asymptotically to 0 as $\mu \to \infty$.

In the case when $T_0$ depends on the chemical potential, the last relation in Eq. (2.57) is modified to:

$$\frac{\partial^2 F}{\partial \mu \partial L_{\pm}} = \frac{\partial^2 V_L}{\partial \mu \partial L_{\pm}} - 6N_f \sum_{\zeta=\pm 1} \zeta \int \frac{d^3 p}{(2\pi)^3} f_{\zeta, \pm} (1 + \delta_{\zeta, \mp} - 3f_{\zeta}). \qquad (A.16)$$

The second derivative $\partial^2 V_L / \partial \mu \partial L_\pm$ of the Polyakov potential in the logarithmic model (2.7) reads:

$$\frac{\partial^2 V_L}{\partial \mu \partial L_{\pm}} = -T^4 \frac{\partial T_0}{\partial \mu} \left[ \frac{6}{\mathcal{H}} (L_{\mp} - 2L_{\pm}^2 + L_{\mp}^2 L_{\pm}) \frac{\partial b}{\partial T_0} + \frac{L_{\mp}}{2} \frac{\partial a}{\partial T_0} \right], \\ \frac{\partial b}{\partial T_0} = 3b_3 \frac{T_0^2}{T^3}, \quad \frac{\partial a}{\partial T_0} = \frac{a_1}{T} + \frac{2a_2 T_0}{T^2}. \qquad (A.17)$$

In the case of the polynomial potential, we have

$$\frac{\partial^2 V_L}{\partial \mu \partial L_{\pm}} = -T^4 \frac{\partial T_0}{\partial \mu} \frac{L_{\mp}}{2} \frac{\partial b_2}{\partial T_0}, \quad \frac{\partial b_2}{\partial T_0} = \frac{a_1}{T} + \frac{2a_2 T_0}{T^2} + \frac{3a_3 T_0^2}{T^3}. \qquad (A.18)$$

The derivative of $T_0$ with respect to $\mu$ can be evaluated starting from Eq. (A.13),

$$\frac{\partial T_0}{\partial \mu} = \frac{T_0}{\alpha_0 b^2(\mu)} \frac{\partial b(\mu)}{\partial \mu}. \qquad (A.19)$$

In the case of the potential proposed in Refs. [27, 50, 67], the derivative of $b(\mu)$ with respect to $\mu$ evaluates to

$$\frac{\partial b(\mu)}{\partial \mu} = - \frac{2\mu b_{\mu}}{(\hat{\gamma}T_{\tau})^2}, \qquad (A.20)$$
while when $b(\mu)$ is represented using the Padé approximant in Eq. (A.15), we have

$$ \frac{\partial b(\mu)}{\partial \mu} = - \frac{2\mu b_{\mu}}{b_0^2 T_{\tau}^2} b^2(\mu). \qquad (A.21) $$

We implemented the above model in our numerical code and used the results reported in Ref. [27] to validate our implementation. Figure 25 confirms that our algorithm is able to correctly capture the temperature dependence of both the scaled scalar condensate, $\sigma/f_{\pi}$, and of the Polyakov loop variables, $L$ and $L^*$, for various values of the chemical potential $\mu$.

In Sec. 2.6, we reported a split between the deconfinement and chiral phase diagrams, for the usual logarithmic Polyakov potential with constant $T_0 = 0.27$ GeV at higher values of the chemical potential. Figure 26 shows that, with these modifications, the coincidence of the two transitions persists for a longer range in the chemical potential.

We would like to point out that the split between the chiral and deconfinement phase transition temperature at high chemical potential can also be observed with the polynomial Polyakov potential used in Ref. [27]. In Fig. 25 we present a comparison between the results we obtained and the results reported in Ref. [27] for $\sigma/f_{\pi}$ and $L, L^*$ with the potential and parameters taken from Ref. [27]. We can see an excellent agreement between the two works. In Fig. 25(b) $L$ and $L^*$ are scaled by, $N = 1.1$, as also discussed in Ref. [27]. However, we observe a split in the phase diagram in Fig.26(a), between the two transitions, unlike what is reported in Ref. [27]. Here we have used the algorithm discussed in Sec.2.6 to obtain the phase diagram.

<figure><img src="image_26.png" /><figcaption>Figure 25. Temperature dependence of (a) scaled $\sigma$ and (b) scaled $L$ and $L^*$ for different values of chemical potential with polynomial Polyakov potential as discussed in Ref. [27]</figcaption></figure>

We thus conclude that, in order to incorporate the effects of dynamical quarks on the effective Polyakov potential, a detailed investigation is necessary.

References

[1] STAR collaboration, Global polarization measurement in Au+Au collisions, Phys. Rev. C 76 (2007) 024915 [0705.1691].
<figure><img src="image_27.png" /><figcaption>**Figure 26.** Chiral and deconfinement phase diagram with chemical potential dependent $T_0$ in the polynomial Polyakov potential with $b(\mu)$ (a) as in Eq. (A.14) and (b) as in Eq. (A.15). (c) The temperature dependence of the Polyakov loop and its slope for two different values of $\theta$, marked in the phase diagram of (a). (d) Double-peak structure of the slope $\partial L/\partial G$ of the Polyakov loop with respect to $G = \sqrt{T^2 + \mu^2}$, computed on diagonal trajectories under various angles $\theta$.</figcaption></figure>

[2] STAR collaboration,
Global $\Lambda$ hyperon polarization in nuclear collisions: evidence for the most vortical fluid,
*Nature* **548** (2017) 62 [1701.06657].

[3] X.-G. Huang, J. Liao, Q. Wang and X.-L. Xia,
Vorticity and Spin Polarization in Heavy Ion Collisions: Transport Models, *Lect. Notes Phys.* **987** (2021) 281 [2010.08937].

[4] F. Becattini, J. Liao and M. Lisa,
Strongly Interacting Matter Under Rotation: An Introduction, *Lect. Notes Phys.* **987** (2021) 1 [2102.00933].

[5] A. Yamamoto and Y. Hirono, *Lattice QCD in rotating frames*, *Phys. Rev. Lett.* **111** (2013) 081601 [1303.6292].

[6] V. V. Braguta, A. Y. Kotov, D. D. Kuznedelev and A. A. Roenko,
*Study of the Confinement/Deconfinement Phase Transition in Rotating Lattice SU(3) Gluodynamics*,
*Pisma Zh. Eksp. Teor. Fiz.* **112** (2020) 9.

[7] V. V. Braguta, A. Y. Kotov, D. D. Kuznedelev and A. A. Roenko,
*Influence of relativistic rotation on the confinement-deconfinement transition in gluodynamics*,
Phys. Rev. D 103 (2021) 094515 [2102.05084].

[8] V. V. Braguta, A. Kotov, A. Roenko and D. Sychev, Thermal phase transitions in rotating QCD with dynamical quarks, PoS LATTICE2022 (2023) 190 [2212.03224].

[9] V. V. Braguta, I. E. Kudrov, A. A. Roenko, D. A. Sychev and M. N. Chernodub, Lattice Study of the Equation of State of a Rotating Gluon Plasma, JETP Lett. 117 (2023) 639.

[10] V. V. Braguta, M. N. Chernodub, A. A. Roenko and D. A. Sychev, Negative moment of inertia and rotational instability of gluon plasma, Phys. Lett. B 852 (2024) 138604 [2303.03147].

[11] J.-C. Yang and X.-G. Huang, QCD on Rotating Lattice with Staggered Fermions, 2307.05755.

[12] V. V. Braguta, M. N. Chernodub and A. A. Roenko, New mixed inhomogeneous phase in vortical gluon plasma: First-principle results from rotating SU(3) lattice gauge theory, Phys. Lett. B 855 (2024) 138783 [2312.13994].

[13] V. V. Braguta, M. N. Chernodub, Y. A. Gershtein and A. A. Roenko, On the origin of mixed inhomogeneous phase in vortical gluon plasma, 2411.15085.

[14] R. Tolman and P. Ehrenfest, Temperature Equilibrium in a Static Gravitational Field, Phys. Rev. 36 (1930) 1791.

[15] R. C. Tolman, On the Weight of Heat and Thermal Equilibrium in General Relativity, Phys. Rev. 35 (1930) 904.

[16] Y. Nambu and G. Jona-Lasinio, Dynamical model of elementary particles based on an analogy with superconductivity. II., Phys. Rev. 124 (1961) 246.

[17] Y. Nambu and G. Jona-Lasinio, Dynamical Model of Elementary Particles Based on an Analogy with Superconductivity. 1., Phys. Rev. 122 (1961) 345.

[18] M. Gell-Mann and M. Lévy, The axial vector current in beta decay, Il Nuovo Cimento 16 (1960) 705-726.

[19] Y. Fujimoto, K. Fukushima and Y. Hidaka, Deconfining Phase Boundary of Rapidly Rotating Hot and Dense Matter and Analysis of Moment of Inertia, Phys. Lett. B 816 (2021) 136184 [2101.09173].

[20] Y. Chen, D. Li and M. Huang, Inhomogeneous chiral condensation under rotation in the holographic QCD, Phys. Rev. D 106 (2022) 106002 [2208.05668].

[21] Y. Jiang, Inhomogeneous SU(2) gluon matter under rotation, Phys. Rev. D 110 (2024) 054047 [2406.03311].

[22] S. Morales-Tejera, V. E. Ambrus and M. N. Chernodub, Firewall boundaries and mixed phases of rotating quark matter in linear sigma model, 2502.19087.

[23] Y. Jiang, Chiral vortical catalysis, Eur. Phys. J. C 82 (2022) 949 [2108.09622].
[24] F. Sun, J. Shao, R. Wen, K. Xu and M. Huang, Chiral phase transition and spin alignment of vector mesons in the polarized-Polyakov-loop Nambu-Jona-Lasinio model, Phys. Rev. D 109 (2024) 116017 [2402.16595].

[25] R. M. Nunes, R. L. S. Farias, W. R. Tavares and V. S. Timóteo, Chiral vortical catalysis constrained by LQCD simulations, 2412.14541.

[26] S. Chen, K. Fukushima and Y. Shimada, Inhomogeneous confinement and chiral symmetry breaking induced by imaginary angular velocity, 2404.00965.

[27] B.-J. Schaefer, J. M. Pawlowski and J. Wambach, The Phase Structure of the Polyakov-Quark-Meson Model, Phys. Rev. D 76 (2007) 074023 [0704.3234].

[28] K. Fukushima and V. Skokov, Polyakov loop modeling for hot QCD, Prog. Part. Nucl. Phys. 96 (2017) 154 [1705.00718].

[29] E. Siri and N. Sadooghi, Bose-Einstein condensation in a rigidly rotating relativistic boson gas, Phys. Rev. D 111 (2025) 036011 [2411.12581].

[30] M. N. Chernodub, Inhomogeneous confining-deconfining phases in rotating plasmas, Phys. Rev. D 103 (2021) 054027 [2012.04924].

[31] M. N. Chernodub, V. A. Goy and A. V. Molochkov, Inhomogeneity of a rotating gluon plasma and the Tolman-Ehrenfest law in imaginary time: Lattice results for the Polyakov loop, Phys. Rev. D 107 (2023) 114502 [2209.15534].

[32] S. Roessner, C. Ratti and W. Weise, Polyakov loop, diquarks and the two-flavour phase diagram, Phys. Rev. D 75 (2007) 034007 [hep-ph/0609281].

[33] C. Ratti, M. A. Thaler and W. Weise, Phases of QCD: Lattice thermodynamics and a field theoretical model, Phys. Rev. D 73 (2006) 014019 [hep-ph/0506234].

[34] S. Ebihara, K. Fukushima and K. Mameda, Boundary effects and gapped dispersion in rotating fermionic matter, Phys. Lett. B 764 (2017) 94 [1608.00336].

[35] Z. Zhang, C. Shi, X.-T. He, X. Luo and H.-S. Zong, Chiral phase transition inside a rotating cylinder within the Nambu-Jona-Lasinio model, Phys. Rev. D 102 (2020) 114023 [2012.01017].

[36] M. Hortacsu, K. D. Rothe and B. Schroer, Zero Energy Eigenstates for the Dirac Boundary Problem, Nucl. Phys. B 171 (1980) 530.

[37] V. E. Ambrus and E. Winstanley, Rotating fermions inside a cylindrical boundary, Phys. Rev. D 93 (2016) 104014 [1512.05239].

[38] A. Chodos, R. L. Jaffe, K. Johnson, C. B. Thorn and V. F. Weisskopf, A New Extended Model of Hadrons, Phys. Rev. D 9 (1974) 3471.

[39] M. N. Chernodub and S. Gongyo, Interacting fermions in rotation: chiral symmetry restoration, moment of inertia and thermodynamics, JHEP 01 (2017) 136 [1611.02598].
[40] M. N. Chernodub and S. Gongyo,
Effects of rotation and boundaries on chiral symmetry breaking of relativistic fermions,
*Phys. Rev. D* **95** (2017) 096006 [1702.08266].

[41] M. N. Chernodub and S. Gongyo,
Edge states and thermodynamics of rotating relativistic fermions under magnetic field, *Phys. Rev. D* **96** (2017) 096014 [1706.08448].

[42] F. Becattini and E. Grossi,
Quantum corrections to the stress-energy tensor in thermodynamic equilibrium with acceleration,
*Phys. Rev. D* **92** (2015) 045037 [1505.07760].

[43] V. E. Ambrus, Quantum non-equilibrium effects in rigidly-rotating thermal states, *Phys. Lett. B* **771** (2017) 151 [1704.02933].

[44] F. Becattini, L. Bucciantini, E. Grossi and L. Tinti,
Local thermodynamical equilibrium and the beta frame for a quantum relativistic fluid, *Eur. Phys. J. C* **75** (2015) 191 [1403.6265].

[45] M. Buzzegoli and F. Becattini,
General thermodynamic equilibrium with axial chemical potential for the free Dirac field,
*JHEP* **12** (2018) 002 [1807.02071].

[46] F. Sun, K. Xu and M. Huang,
Splitting of chiral and deconfinement phase transitions induced by rotation, *Phys. Rev. D* **108** (2023) 096007 [2307.14402].

[47] P. Singha, V. E. Ambrus and M. N. Chernodub,
Inhibition of splitting of the chiral and deconfinement transition due to rotation in QCD: the phase diagram of
*JHEP* **12** (2018) 07828.

[48] V. E. Ambrus and M. N. Chernodub,
Rigidly rotating scalar fields: Between real divergence and imaginary fractalization, *Phys. Rev. D* **108** (2023) 085016 [2304.05998].

[49] T. Pătuleanu, A. D. Fodor, V. E. Ambrus and C. Crucean,
Dirac fermions under imaginary rotation, *Phys. Rev. D* **75** (2007) 014014.

[50] T. K. Herbst, J. M. Pawlowski and B.-J. Schaefer,
The phase structure of the Polyakov–quark–meson model beyond mean field, *Phys. Lett. B* **696** (2011) 58 [1008.0081].

[51] PARTICLE DATA GROUP collaboration, *Review of particle physics*, *Phys. Rev. D* **110** (2024) 030001.

[52] M. Daum, R. Frosch and P. R. Kettle, *The charged and neutral pion masses revisited*, Phys. Lett. B **796** (2019) 11.

[53] G. 't Hooft,
*Topology of the Gauge Condition and New Confinement Phases in Nonabelian Gauge Theories*,
*Nucl. Phys. B* **190** (1981) 455.

[54] J. I. Kapusta and C. Gale, *Finite-Temperature Field Theory: Principles and Applications*,
*Cambridge Monographs on Mathematical Physics*. Cambridge University Press, 2 ed., 2006,
10.1017/CBO9780511535130.

[55] O. Scavenius, A. Mocsy, I. N. Mishustin and D. H. Rischke,
Chiral phase transition within effective models with constituent quarks, *Phys. Rev. C* **64** (2001) 045202 [nucl-th/0007030].

[56] C. Ratti, S. Roessner and W. Weise, Quark number susceptibilities: Lattice QCD versus PNJL model, *Phys. Lett. B* **649** (2007) 57 [hep-ph/0701091].

[57] C. G. Broyden, A class of methods for solving nonlinear simultaneous equations, *Math. Comput.* **19** (1965) 577.

[58] W. H. Press, S. A. Teukolsky, W. T. Vetterling and B. P. Flannery, *Numerical Recipes in C*. Cambridge University Press, Cambridge, USA, second ed., 1992.

[59] M. N. Chernodub and V. E. Ambrus, *Phase diagram of helically imbalanced QCD matter*, *Phys. Rev. D* **103** (2021) 094015 [2005.03575].

[60] L. Schubnikow and W. J. de Haas, *Neue erscheinungen bei der widerstandsänderung von wismuthkristallen im magnetfeld bei der temperatur von flüssigem wasserstoff (i)*, *Proceedings of the Royal Netherlands Academy of Arts and Sciences* **33** (1930) 130.

[61] F. Sun, K. Xu and M. Huang, *Quarkyonic phase induced by Rotation*, **2307.14402**.

[62] C. Cercignani and G. M. Kremer, *The Relativistic Boltzmann Equation: Theory and Applications*, Springer, 2002.

[63] V. E. Ambrus and E. Winstanley, *Rotating quantum states*, *Phys. Lett. B* **734** (2014) 296 [1401.6388].

[64] V. V. Braguta, M. N. Chernodub, I. E. Kudrov, A. A. Roenko and D. A. Sychev, *Negative Barnett effect, negative moment of inertia of the gluon plasma, and thermal evaporation of the chromo*- Phys. Rev. D **110** (2024) 014511 [2310.16036].

[65] K. Kajantie, M. Laine, K. Rummukainen and Y. Schroder, *The Pressure of hot QCD up to g$^6 \ln(1/g)$*, *Phys. Rev. D* **67** (2003) 105008 [hep-ph/0211321].

[66] S. K. Ghosh, T. K. Mukherjee, M. G. Mustafa and R. Ray, PNJL model with a Van der Monde term, *Phys. Rev. D* **77** (2008) 094024 [0710.2790].

[67] T. K. Herbst, J. M. Pawlowski and B.-J. Schaefer, *Phase structure and thermodynamics of QCD*, *Phys. Rev. D* **88** (2013) 014007 [1302.1426].