# On the origin of mixed inhomogeneous phase in vortical gluon plasma

Victor Braguta, Maxim Chernodub, Yana Gershtein, Artem Roenko

► To cite this version:

Victor Braguta, Maxim Chernodub, Yana Gershtein, Artem Roenko. On the origin of mixed inhomogeneous phase in vortical gluon plasma. Journal of High Energy Physics, 2025, 2025 (9), pp.79. 10.1007/JHEP09(2025)079. hal-04819279v2

HAL Id: hal-04819279

https://hal.science/hal-04819279v2

Submitted on 11 Sep 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
# On the origin of mixed inhomogeneous phase in vortical gluon plasma

V. V. Braguta,ª M. N. Chernodub,ᵇ,ᶜ Ya. A. Gershtein,ᵈ A. A. Roenkoª

ªBogoliubov Laboratory of Theoretical Physics, Joint Institute for Nuclear Research, Dubna, 141980 Russia

ᵇUniversité de Tours, Université d'Orléans, CNRS, Institut Denis Poisson, UMR 7013, Tours, 37200, France

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

ᵈMoscow Institute of Physics and Technology, Dolgoprudny, 141700 Russia

E-mail: roenko@theor.jinr.ru

**ABSTRACT:** Recently, lattice simulations of SU(3) Yang-Mills theory revealed that rotating hot gluon matter in thermal equilibrium possesses a novel inhomogeneous phase consisting of the deconfinement phase located in the center region, which is spatially separated from the confinement phase in the periphery. This inhomogeneous two-phase structure is also expected to be produced by vorticity in quark-gluon plasma formed in non-central relativistic heavy-ion collisions. We show that its vortical properties are determined by two types of couplings of the angular velocity to the gluon fields: a linear coupling to the mechanical angular momentum of gluons and a quadratic “magnetovortical” coupling to a chromomagnetic component. We demonstrate numerically that the distinctive inhomogeneous structure of the vortical (quark-)gluon plasma is determined by the latter, while the former plays only a subleading role. We argue that the anisotropy of the gluonic action in the curved co-rotating background can quantitatively explain the remarkable property that the spatial structure of this inhomogeneous phase disobeys the picture based on a straightforward implementation of the Tolman-Ehrenfest law. We also support our findings with Monte Carlo simulations of Yang-Mills plasma at the real-valued angular frequency, which take into account only the magnetic part of the action.

ARXIV ePRINT: 2411.15085
# Contents

<table><tr><td>1 Introduction</td><td>2</td></tr><tr><td>2 Gluons in rotation: Mechanical and magnetic vortical couplings</td><td>3</td></tr><tr><td>2.1 Kinematics of rigid rotation</td><td>3</td></tr><tr><td>2.2 QCD matter in rotating frame</td><td>4</td></tr><tr><td>2.2.1 Fermions</td><td>4</td></tr><tr><td>2.2.2 Gluons</td><td>5</td></tr><tr><td>2.3 Exotic vortical properties of gluonic matter</td><td>7</td></tr><tr><td>3 Lattice simulation of rotating gluodynamics</td><td>9</td></tr><tr><td>3.1 Theoretical background</td><td>9</td></tr><tr><td>3.2 Lattice set up</td><td>10</td></tr><tr><td>3.3 Phase transition and Polyakov loop</td><td>11</td></tr><tr><td>4 Mixed inhomogeneous phase in rotating gluon plasma</td><td>12</td></tr><tr><td>4.1 Observation of mixed inhomogeneous phase</td><td>12</td></tr><tr><td>4.2 Study of local critical temperature</td><td>14</td></tr><tr><td>4.3 The phase diagram for real and imaginary rotation</td><td>19</td></tr><tr><td>5 Decomposition of rotating action</td><td>20</td></tr><tr><td>5.1 Mechanical vs. chromomagnetic vortical effects.</td><td>20</td></tr><tr><td>5.2 Inhomogeneous phases</td><td>22</td></tr><tr><td>6 Local thermalization of rotating gluodynamics</td><td>25</td></tr><tr><td>6.1 Local approximation for inhomogeneous action</td><td>25</td></tr><tr><td>6.2 The results obtained with the local action</td><td>26</td></tr><tr><td>6.3 Inconsistency with the Tolman-Ehrenfest approach</td><td>29</td></tr><tr><td>7 Discussion and conclusion</td><td>31</td></tr><tr><td>A Width of the transition region</td><td>34</td></tr><tr><td>B Calculation of the local critical temperature</td><td>35</td></tr><tr><td>C Analysis of the finite-volume and lattice spacing effects</td><td>36</td></tr><tr><td>D Estimation of systematic uncertainty for the fit coefficients</td><td>37</td></tr></table>
# 1 Introduction

The influence of rapid rotation on the properties of different physical systems proved to be an intriguing and important area of scientific research. There are many examples when rotation plays a significant role in astrophysics [1, 2], condensed matter physics [3, 4], nuclear physics [5], to mention a few. The fastest rotation in Nature is realized in heavy-ion collision experiments that create highly vortical quark-gluon matter [6, 7]. The global polarization of $\Lambda$ and $\bar{\Lambda}$ hyperons obtained by the STAR Collaboration indicate that the vorticity of hot quark-gluon matter in noncentral collisions reaches the enormous value $\sim 10^{22}$ Hz ($\Omega \sim 10$ MeV) [8]. In terms of speed of light $c$, the mechanical vortical motion of plasma acquires significant velocities $v \sim \Omega r \sim (0.1...0.3)c$ at typical distances $r \sim (\text{a few}) \times \text{fm}$ within the expanding fireball.

It is reasonable to assume that the high vorticity of the quark-gluon matter affects its thermodynamic properties and, in particular, modifies the phase diagram of QCD. The last question has been addressed in various papers within different theoretical approaches (see, for instance, [9–27]). Most of the theoretical works predict a decrease in the critical temperatures of deconfinement and chiral transitions in rotating hot QCD matter. On the contrary, some theoretical works point out to quite the opposite behaviour [28–31]. The influence of rotation on the phase transitions in gluodynamics and QCD within lattice simulation has also been studied, following an early work of Ref. [32], in several papers both in pure gluonic plasmas [33–36] as well as in lattice QCD with dynamical quarks [37, 38]. Contrary to the predictions of many theoretical approaches, these first-principles simulations demonstrate that both transition temperatures rise with increasing rotation.

It was also suggested that rotation could lead to a new, inhomogeneous phase in QCD, in which spatially separated hadronic and plasma phases could exist simultaneously [39]. This result, obtained in a low-dimensional confining model, has been subsequently supported by a study in the holography approach [40]. Later, the first-principle lattice simulation conducted in Ref. [41] confirmed that mixed-phase indeed takes place in SU(3) gluodynamics. It was shown that the local critical temperature increases and the mixed state has an inverted phase arrangement compared to Refs. [39, 40]. In particular, the lattice simulations of Ref. [41] have revealed that for the temperatures close to the critical one, there appears a mixed state with deconfinement at the center and confinement close to the boundaries of the studied volume. This observation has also been supported in analytical approaches [42, 43].

In a parallel development, the rotating gluon plasma has numerically been shown to possess a negative moment of inertia in a region of temperatures above the deconfining phase transition $T_c$ up to the “supervortical temperature” $T_s \simeq 1.5T_c$ [44–46]. This phenomenon has been proposed to originate from the negative Barnett effect and the thermal evaporation of the non-perturbative chromomagnetic condensate, which gives, together with the standard angular momentum, a contribution to the mechanical moment of inertia of the gluon plasma [47]. Similar features were also argued to emerge in a relativistic bosonic gas under certain regimes [48, 49].

In this paper, we present a more detailed study of the mixed phase found in Ref. [41]
in an inquiry to shed more light on its origin. In the spirit of Ref. [47], we separate the effects coming from the mechanical moment of thermal gluons from the chromomagnetic contribution and analyze their impact on the formation of the inhomogeneous plasma.

This paper is organized as follows. In Section 2, we discuss the kinematical properties of rigid rotation at the level of the QCD Lagrangian. We stress the non-linear nature of the coupling of the angular velocity to the gluonic fields as opposed to the linear coupling of the rotation to quarks. A special role of the chromomagnetic degrees of freedom is highlighted. In Section 3, we present a lattice framework implemented in the rotating SU(3) gluodynamics and discuss the peculiarities of the rotating system. Section 4 is devoted to the results of the lattice simulations revealing the spatial distribution of the Polyakov loop. This Section also provides details on the calculation of the local critical temperature of the deconfinement transition. In Section 5, we study how the terms in lattice action, which are linear/quadratic in angular velocity, influence the critical temperature and discuss the analytic continuation of our results from the domain of imaginary angular frequencies implemented in the lattice simulations, to the real-valued frequencies, relevant to the realistic plasmas. In Section 6, we show that the metric of the rotating reference frame gives rise to asymmetry in the coupling constants in the gluon action. We demonstrate that it is the peculiar effect of metric that leads to the shift of the critical temperature in the periphery regions compared to the central regions. The last Section 7 is devoted to the conclusions and detailed discussions of the numerical results. Relevant technical details are provided in the Appendix.

# 2 Gluons in rotation: Mechanical and magnetic vortical couplings

## 2.1 Kinematics of rigid rotation

The thermodynamic equilibrium of a rotating fluid is achieved in a state of a rigid rotation when the angular velocity of every fluid element has the same value $\Omega$ about the axis of rotation $\mathbf{n} = \mathbf{\Omega}/\Omega$ with the angular frequency $\Omega = |\mathbf{\Omega}|$. The solid rotation excludes the hydrodynamic friction of the adjacent fluid layers and prevents the generation of heat (entropy) in consistency with the notion of thermodynamic equilibrium [50, 51].

It is convenient to describe the rotating system in a co-rotating reference frame, $x^\mu \equiv (x^0, x^1, x^2, x^3) = (t, x, y, z)$, in which the rigidly rotating system appears static. Due to the symmetry of the problem, we consider a system bounded by a cylinder of the radius $R$ and briefly use the cylindrical coordinates, $x^\mu = (t, r \sin \varphi, r \cos \varphi, z)$, with $\varphi, r$, and $z$ being the azimuthal angle, the radial and longitudinal coordinates, respectively. The system rotates about the $z$ axis, $\mathbf{\Omega} = \Omega \mathbf{e}_z$.

The metric in the co-rotating system can be obtained by a coordinate transformation from the laboratory frame to the co-rotating frame. To this end, we denote by the subscript “lab” the coordinates in the flat Minkowski metric of a laboratory frame, $x_{\text{lab}}^\mu$. The only difference between the laboratory and co-rotating frames appears in the azimuthal coordinate $\varphi = [\varphi_{\text{lab}} - \Omega t]_{2\pi}$, where $[\dots]_{2\pi}$ denote taking a modulo of $2\pi$. The other coordinate components coincide in both frames, $t = t_{\text{lab}}$, $r = r_{\text{lab}}$ and $z = z_{\text{lab}}$. The diffeomorphism
gives us the metric in the co-rotating frame in the Cartesian coordinates:

$$g_{\mu\nu} = \begin{pmatrix} 1 - r^2 \Omega^2 & \Omega y & -\Omega x & 0 \\ \Omega y & -1 & 0 & 0 \\ -\Omega x & 0 & -1 & 0 \\ 0 & 0 & 0 & -1 \end{pmatrix}, \qquad (2.1)$$

where $r = \sqrt{x^2 + y^2}$ is the distance from the rotation axis. The co-rotating coordinates correspond to a non-inertial nature contrary to the inertial laboratory frame. Notice that $g \equiv \det g_{\mu\nu} = -1$.

In a relativistic theory, the solid rotation can only be considered for a physical body lying entirely within a light cylinder. It means that the farthermost point of the body from the rotation axis should possess a linear velocity smaller than the speed of light, $R\Omega < 1$. This requirement prevents the system from breaking the fundamental causality constraint. For the causality-respecting systems, the $g_{00}$ component of the co-rotating metric (2.1) is a non-negative quantity at any point of the rotating body.

## 2.2 QCD matter in rotating frame

### 2.2.1 Fermions

Before going into the discussion of the main subject of this paper, the kinematics of rotating gluons, it is convenient to discuss first how rigid rotation affects the kinematics of fermions. The Lagrangian of a Dirac particle in the co-rotating frame corresponds to the Dirac Lagrangian in the curved spacetime described by the metric (2.1). Without repeating the details of the derivation that can be found elsewhere (for example, in Refs. [9, 10, 52]), we give here the Dirac Lagrangian written the co-rotating coordinates

$$\mathcal{L}_{\psi} = \mathcal{L}_{\psi}^{(0)} + \mathcal{L}_{\psi}^{(1)} \equiv \mathcal{L}_{\psi, \text{lab}} + \mathcal{L}_{\psi, \text{mech}}, \qquad (2.2)$$

which is a linear function of the angular frequency $\Omega$ with $\mathcal{L}_{\psi}^{(n)} \propto \Omega^n$ ($n = 0, 1$). The terms in Lagrangian (2.2) are

$$\mathcal{L}_{\psi, \text{lab}} \equiv \mathcal{L}_{\psi}^{(0)} = \bar{\psi}(i\gamma^{\mu}\partial_{\mu} - m)\psi, \qquad (2.3)$$

$$\mathcal{L}_{\psi, \text{mech}} \equiv \mathcal{L}_{\psi}^{(1)} = \Omega \cdot \mathbf{J}_{\psi}, \qquad (2.4)$$

where $\psi$ is the Dirac 4-spinor and $\bar{\psi} = \psi^{\dagger}\gamma^t$ is its Dirac conjugate. The last term in Eq. (2.4) is linear in the local angular momentum

$$\mathbf{J}_{\psi} \equiv \mathbf{J}_{\psi}(\mathbf{x}) = \psi^{\dagger}(\mathbf{x})\hat{\mathbf{J}}\psi(\mathbf{x}), \qquad (2.5)$$

carried by the fermion field $\psi$. The z-component of the angular momentum operator, which enters the last term in Eq. (2.4), is given by the sum of the angular momentum operator and the spin operator, respectively:

$$\hat{\mathbf{J}}_{\psi,z} = -i(-y\partial_x + x\partial_y) + \frac{1}{2}\Sigma^{xy} \equiv -i\partial_{\varphi} + \frac{1}{2}\begin{pmatrix} \sigma_3 & 0 \\ 0 & \sigma_3 \end{pmatrix}. \qquad (2.6)$$
In the last line, the $\Sigma^{xy}$ component of the spin matrix $\Sigma^{\mu\nu} = (i/2)[\gamma^\mu, \gamma^\nu]$ is written in the Dirac representation with $\sigma_3$ being a Pauli matrix.

The Dirac Lagrangian (2.2) describes a neutral massive spinor field in the co-rotating reference frame equipped with the metric (2.1). The Dirac Lagrangian in the laboratory frame (2.3), $\mathcal{L}_{\text{lab}} \equiv \mathcal{L}(\Omega = 0)$, is given by the very same Eq. (2.2) with angular frequency set to zero, $\Omega = 0$.

It appears that the difference between the Lagrangians in the laboratory and co-rotating frames originates in the single term (2.4) that includes the angular momentum (2.6). The meaning of the term produced by the rotating environment is easy to estimate if one considers the Hamiltonian of the system (2.2):

$$ \hat{H}_{\psi} = \hat{H}_{\psi, \text{lab}} - \Omega \hat{J}_{\psi}. \qquad (2.7) $$

The rotating environment shifts the energy levels of the system via the mechanical coupling between the background angular velocity $\Omega$ and the total angular momentum of the system $\mathbf{J}_{\psi}$. The latter quantity comprises both the orbital and spin contributions. This natural type of mechanical coupling appears in classical mechanics [53] and thermodynamics [54].

For the reasons which will be clear shortly after, we call the last term in the Hamiltonian (2.7) the “mechanical-vortical coupling”. This term is responsible for various effects, including the standard classical mechanical phenomena exhibited by the rotating systems, such as the appearance of the Coriolis force and the centrifugal force [55]. The last statement is in no way trivial because the centrifugal force in mechanical systems is proportional to the square of the angular velocity, $\Omega^2$. In contrast, the coupling of the angular momentum of the system with the background angular velocity in the Hamiltonian (2.7) has a linear nature.

It is evident that the very same considerations also apply to quarks. The Dirac Lagrangian (2.2) will then be modified to include the coupling to gluons via a covariant derivative as well as color and flavor indices of the Dirac spinors. The Hamiltonian for quarks will then be described by the general form given in Eq. (2.7), with the mentioned modifications of its terms.

### 2.2.2 Gluons

Gluon fields in the co-rotating frame are described by SU(3) Yang-Mills theory in the curved spacetime endowed with the metric (2.1). The Lagrangian in the Minkowski spacetime has the following form:

$$ \mathcal{L}_G = -\frac{1}{4g_{YM}^2} g^{\mu\nu} g^{\alpha\beta} F_{\mu\alpha}^a F_{\nu\beta}^a, \qquad (2.8) $$

where $F_{\mu\nu}^a$ is the gluon field stress tensor and $g_{YM}$ is the Yang-Mills coupling constant. Substituting the metric (2.1) into Eq. (2.8) allows us to decompose the Yang-Mills Lagrangian into three terms:

$$ \begin{aligned} \mathcal{L}_G &= \mathcal{L}_G^{(0)} + \mathcal{L}_G^{(1)} + \mathcal{L}_G^{(2)} \\ &\equiv \mathcal{L}_{G, \text{lab}} + \mathcal{L}_{G, \text{mech}} + \mathcal{L}_{G, \text{magn}}, \end{aligned} \qquad (2.9) $$
which, contrary to the fermionic Lagrangian (2.2), is a higher-order, quadratic polynomial of the angular frequency $\Omega$ with $\mathcal{L}_G^{(n)} \propto \Omega^n$ for $n = 0, 1, 2$. Each term in Eq. (2.9) has a certain physical meaning, which we discuss below.

(i) **Laboratory frame term.** The first term in Eq. (2.9) is given by the Yang-Mills Lagrangian in the inertial (non-rotating) laboratory frame:

$$ \mathcal{L}_G^{(0)} = -\frac{1}{4g_{YM}^2} \eta^{\mu\nu} \eta^{\alpha\beta} F_{\mu\alpha}^a F_{\nu\beta}^a, \qquad (2.10) $$

where $\eta^{\mu\nu} = \text{diag}(+1, -1, -1, -1)$ is the flat Minkowski metric. The other two terms in Lagrangian (2.9) encode the effects of the global rotation.

(ii) **Mechanical coupling.** The second contribution to Lagrangian (2.9) is a linear term in the angular frequency $\Omega$. It has precisely the same generic form as the coupling for the fermions (2.4):

$$ \mathcal{L}_{G,\text{mech}} \equiv \mathcal{L}_G^{(1)} = \Omega \cdot \mathbf{J}_G, \qquad (2.11) $$

This term encodes the mechanical coupling between the background angular velocity $\Omega$ and the angular momentum density of gluons, $\mathbf{J}_G$. The gluonic angular momentum density, taken in the non-rotating limit ($\Omega \to 0$), has the following components:

$$ J_{G,i} = \frac{1}{2} \varepsilon_{ijk} J_G^{jk}, \qquad J_G^{ij} = x^i T_G^{jt} - x^j T_G^{it}, \qquad (2.12) $$

where

$$ T_G^{\mu\nu} = -\frac{1}{g_{YM}^2} \left( F^{a,\mu\alpha} F^{a,\nu}_{\alpha} - \frac{1}{4} \eta^{\mu\nu} F^{a,\alpha\beta} F_{\alpha\beta}^a \right), \qquad (2.13) $$

is the gluonic stress-energy tensor in the Minkowski spacetime. The rising and lowering of indices in Eq. (2.13) are performed using the Minkowski metric $\eta^{\mu\nu}$.

The gluonic angular momentum, given in Eqs. (2.12) and (2.13), in laboratory frame has the following familiar form:

$$ \mathbf{J}_G = \frac{1}{g_{YM}^2} \mathbf{r} \times (\mathbf{E}^a \times \mathbf{B}^a), \qquad (2.14) $$

where $\mathbf{S} \equiv \frac{1}{g_{YM}^2}(\mathbf{E}^a \times \mathbf{B}^a)$ is the Poynting vector which represents the directional energy flux for the gluon field. Notice that both chromoelectric, $E_i^a = F_{ti}^a$, and chromomagnetic, $B_i^a = -(1/2)\varepsilon_{ijk}F^{a,jk}$, fields contribute to the gluon angular momentum (2.14).

(iii) **Magnetic coupling.** The gluon Lagrangian (2.9) also has a term quadratic in the angular frequency $\Omega$:

$$ \mathcal{L}_{G,\text{magn}} \equiv \mathcal{L}_G^{(2)} = \frac{1}{2g_{YM}^2} \left[ \Omega^2 (\mathbf{B}^a \cdot \mathbf{r}_{\perp})^2 + r_{\perp}^2 (\mathbf{B}^a \cdot \mathbf{\Omega})^2 \right], \qquad (2.15) $$

where $\mathbf{r}_{\perp} = \mathbf{r} - z\mathbf{e}_z \equiv \mathbf{e}_z \times (\mathbf{r} \times \mathbf{e}_z)$ is a vector which points, along a normal, from the axis of rotation to the point $\mathbf{r} = (x, y, z)$. The quadratic term (2.15) couples the
chromomagnetic field $B^a$ to the background angular velocity $\Omega$. Therefore, we call this term a “magnetovortical coupling”.

One could naively suspect that the magnetovortical term (2.15), being quadratic in $\Omega$, is responsible for the net centrifugal force that occurs in a rotating gluon gas. However, while this term certainly contributes to the mentioned force, it does not lie in its origin. To see this fact, one can note that the centrifugal force does exist for fermions that do not possess the quadratic in $\Omega$ coupling in the fermionic Lagrangian (2.2).

The Lagrangian of gluons in rotating frame is written in terms of rotating gluon fields, and it is a quadratic function in the angular velocity (2.9), highlighting the magnetovortical coupling. However, the Hamiltonian of the rotating gluons has the usual linear form

$$ \hat{H}_G = \hat{H}_{G,\text{lab}} - \Omega \hat{J}_G, \qquad (2.16) $$

as it was explicitly shown in Refs. [38, 56], and the static system in rotating frame corresponds to rotating system in the laboratory frame. The key distinction between gluons and quarks is that the gluons are vector fields, and they transform when we pass from one to another reference frame. As the result, Eqs. (2.9) and (2.16) have different structure, and the magnetic coupling arises explicitly only in the rotating frame. It implies that some degrees of freedom in the rest frame might play a role of the magnetic coupling in rotating frame, but it is a non-trivial task to separate them explicitly in the laboratory frame.

It is also worth noticing that the canonical analyses of the thermodynamics of rotating systems are conventionally formulated in the rotating reference frame, in which the temperature and other thermodynamic parameters (such as the chemical potentials) are specified [54, 57]. A calculation in the laboratory reference frame lies outside of the scope of the present paper since here we concentrate on the conventional analysis in the corotating reference frame. A discussion of the relation between calculations in the two frames in the scope of the lattice gluodynamics can be found in Refs. [38, 56].

## 2.3 Exotic vortical properties of gluonic matter

Our first-principle numerical simulations [44–47] revealed that hot gluonic matter has exotic vortical properties in the phenomenologically important window of temperatures that ranges from slightly lower than the deconfining transition at $T = T_c$ up to supervortical temperature $T_s \simeq 1.5T_c$. The most striking feature of the gluon plasma at these temperatures is that it possesses a negative moment of inertia [47].

In a mechanical sense, the microscopic mechanism of this exotic effect has been proposed to originate from the negative Barnett effect, which implies the negative value of the spin-vortical coupling in rotating hot gluon matter. [47] The negativity of the coupling has an intuitively clear meaning, indicating that in the rotating gluonic fluid, the total spin polarization of vortical gluons and their total orbital momentum should be oriented in opposite directions. Since the orbital momentum points in the direction of the angular velocity, the total angular momentum $J$ stored in the vortical plasma becomes directed opposite to its angular velocity $\Omega$. The latter fact implies the negativity of the moment of inertia, $I < 0$, which is nothing but the coefficient of proportionality between the angular momentum and the angular velocity of the body, $J = I\Omega$.
Initially, the negative moment of inertia has been associated with the instability of the host gluon plasma with respect to the vortical motion [44]. However, it later became clear that the plasma can rotate in the thermodynamic equilibrium while the negative moment of inertia can qualitatively be explained by the negative Barnett effect [47]. The latter property is a particular feature of the strongly interacting vector particles, gluons.

The moment of inertia for a generic statistical system is determined via its free energy F in the co-rotating reference frame [54]:

$$I(T) = \frac{J(T, \Omega)}{\Omega} \Big|_{\Omega \to 0} = -\frac{1}{\Omega} \left( \frac{\partial F(T, \Omega)}{\partial \Omega} \right)_{T \to 0} \quad (2.17)$$

For gluodynamics, the free energy F = F<sub>G</sub> will be defined in the next section.

The moment of inertia of gluons can be represented as a sum of two parts [47]:

$$I_G = \frac{1}{T} \left\langle \left\langle (\mathbf{n} \cdot \mathbf{J}_G^{\text{tot}})^2 \right\rangle \right\rangle_T + \int d^3x \left[ \left\langle \left\langle (\mathbf{B}^a \cdot \mathbf{r}_\perp)^2 \right\rangle \right\rangle_T + \left\langle \left\langle (\mathbf{B}^a \cdot \mathbf{n})^2 \right\rangle \right\rangle_T r_\perp^2 \right], \quad (2.18)$$

where $\mathbf{n} = \Omega/\Omega$ is the axis of rotation, $\langle\langle\mathcal{O}\rangle\rangle_T = \langle\langle\mathcal{O}\rangle\rangle_T - \langle\langle\mathcal{O}\rangle\rangle_{T=0}$ is the thermal part of the expectation value of an operator $\mathcal{O}$, and $\mathbf{J}_G^{\text{tot}}$ is the total gluonic angular momentum (2.14) integrated over the volume of the system (a rigorous definition is given in [47]).

The first term in Eq. (2.18) appears due to the linear coupling (2.11) of the angular gluonic degrees of freedom to the angular velocity. This term has a form of a susceptibility which represents a conventional mechanical contribution to the moment of inertia. As noted in Ref. [47], the mechanical term gives always a non-negative contribution to the moment of inertia because $\langle\langle (\mathbf{n} \cdot \mathbf{J}_G^{\text{tot}})^2 \rangle\rangle_T \ge 0$.

The presence of the last term in Eq. (2.18) is the consequence of the quadratic coupling (2.15) of gluons to the angular velocity, which appears in the gluonic Lagrangian in the rotating reference frame (2.9). While the linear coupling is commonly encountered across a wide range of statistical systems, the quadratic coupling has an exotic nature, which can be traced back to Ref. [32].

At the level of quantum chromodynamics, the negativity of the moment of inertia can be associated with the thermal melting of the chromomagnetic condensate [44] which provides a negative contribution to the second term in Eq. (2.18). The gluon condensate is a nonperturbative quantity proportional to the expectation value of the gluon field strength squared, $\langle\langle \text{Tr } F_{\mu\nu}F^{\mu\nu} \rangle\rangle > 0$. As we have already mentioned, the magnetovortical term (2.15) couples the bilinear of the chromomagnetic field to (a square of) the angular velocity and, thus, provides a contribution of the chromomagnetic condensate to the moment of inertia. The first-principle numerical calculations in Ref. [47] demonstrated that this contribution takes a negative value in a relevant temperature window. Thus, the $I < 0$ effect has been attributed to the thermal melting of the gluon condensate.<fnref n="1" /> The present paper aims to

<footnote n="1">¹The gluon condensate is a gauge-invariant quantity. Notice that at zero temperature, it also possesses the Lorentz invariance, implying that its division into chromomagnetic and chromoelectric components has no Lorentz-invariant sense. Indeed, the condensate is described by the Lorentz invariant quantity Tr $F_{\mu\nu}F^{\mu\nu} \equiv (\mathbf{B}^a)^2 - (\mathbf{E}^a)^2$ and not by the electric field $\mathbf{E}^a$ or the magnetic field $\mathbf{B}^a$ separately. Still, we use the generally accepted term “chromomagnetic condensate” because, in the nonperturbative ground state, one has $\langle\langle \text{Tr } F_{\mu\nu}F^{\mu\nu} \rangle\rangle > 0$, which is interpreted as the prevalence of the chromomagnetic condensate over the chromoelectric one [58].</footnote>
elucidate further the expected crucial role of the magnetovortical term (2.15) in the exotic
vortical properties of hot gluon matter.

3 Lattice simulation of rotating gluodynamics

**3.1 Theoretical background**

We study the rotating matter of hot gluons using first-principle numerical simulations of the lattice Yang-Mills theory.<fnref n="2" /> The matter resides in thermodynamic equilibrium in the non-inertial co-rotating reference frame, which is described by metric (2.1). The partition function in the Euclidean spacetime, suitable for Monte Carlo simulations can be written as a path integral over the gluon fields

$$
Z \equiv e^{-F_G/T} = \int D A e^{-S_G[A]}. \qquad (3.1)
$$

where $F_G$ is the free energy of gluons. The $S_G$ in Eq. (3.1) is the Euclidean action of
rotating Yang-Mills theory:

$$
S_G = \frac{1}{4g_{YM}^2} \int d^4 x \sqrt{g_E} g_E^{\mu\nu} g_E^{\alpha\beta} F_{\mu\alpha}^a F_{\nu\beta}^a, \quad (3.2)
$$

where $F_{\mu\nu}^a$ is the field-strength tensor of the gluon field. The Euclidean metric $g_{\mu\nu}^E$ is
related to the Minkowski metric (2.1) via the Wick rotation $t \rightarrow -i\tau$. The imaginary-
time direction $\tau$ is compactified to a circle of the length of the inverse temperature $1/T_0$
with periodic boundary conditions for the gluon fields. The explicit form of the Euclidean
action (3.2) is as follows:

$$
S_G = \frac{1}{2g_{YM}^2} \int d^4 x \left[
\begin{aligned}
& F_{x\tau}^a F_{x\tau}^a + F_{y\tau}^a F_{y\tau}^a + F_{z\tau}^a F_{z\tau}^a + \\
& \qquad + (1-x^2\Omega^2) F_{yz}^a F_{yz}^a + (1-y^2\Omega^2) F_{xz}^a F_{xz}^a + (1-r^2\Omega^2) F_{xy}^a F_{xy}^a - \\
& \qquad - 2yi\Omega (F_{xy}^a F_{y\tau}^a + F_{xz}^a F_{z\tau}^a) + 2xi\Omega (F_{yx}^a F_{x\tau}^a + F_{yz}^a F_{z\tau}^a) - 2xy\Omega^2 F_{xz}^a F_{zy}^a
\end{aligned}
\right]. \quad (3.3)
$$

For later use, we rewrite this action in cylindrical coordinates $dx^\mu = (d\tau, dr, r d\varphi, dz)$ where
it has a quite simple form:

$$
S = \frac{1}{2g_{YM}^2} \int d^4 x \left[ \begin{aligned}
& F_{r\tau}^a F_{r\tau}^a + F_{\varphi\tau}^a F_{\varphi\tau}^a + F_{z\tau}^a F_{z\tau}^a + F_{rz}^a F_{rz}^a + \\
& \qquad + (1-r^2\Omega^2) (F_{\varphi z}^a F_{\varphi z}^a + F_{r\varphi}^a F_{r\varphi}^a) + 2ri\Omega (F_{\varphi r}^a F_{r\tau}^a + F_{\varphi z}^a F_{z\tau}^a)
\end{aligned} \right], \quad (3.4)
$$

where $\hat{\varphi}$ denotes the polar direction with $d\hat{\varphi} = rd\varphi$.

Notice that the thermodynamic equilibrium of a system of particles subjected to an ex-
ternal static gravitational field has, in general, an inhomogeneous temperature that encodes
the redshift of the thermal wavelength by a gravitational force. The local temperature is

<footnote n="2"><sup>2</sup>We follow the approach of Refs. [34, 37] to which we refer an interested reader for technical details.</footnote>
given by the Ehrenfest-Tolman (TE) law: $T(r)\sqrt{g_{00}} = \text{const}$ [59, 60]. For a rigidly rotating system with the metric (2.1), the TE law gives $T(r)\sqrt{1 - r^2\Omega^2} = T_0$ [34], where $T_0 = T(r = 0)$ corresponds to the temperature at the axis of rotation, which is determined by the extension $1/T_0$ of the lattice in the temporal direction. For notational simplicity, we denote hereafter by $T$ the value of temperature precisely at the rotation axis, where the metric-induced effects vanish unless stated otherwise.

## 3.2 Lattice set up

In this work, we discretize rotating terms in the action (3.3) following Refs. [32, 34] and use the tree-level improved Symanzik gauge action for the terms without rotation [61, 62]. The explicit lattice expression for action (3.3) can be found in Ref. [44].

We perform the simulations on the lattices with the spacetime geometry $N_t \times N_z \times N_s^2$ ($N_x = N_y = N_s$). The rotational axis is located at the center of the $xy$ plane. In our simulations, the transverse spatial extension of the lattice $N_s$ is an odd number, so the rotation axis passes through the sites of our lattice. It is worth mentioning that the temporal extension $N_t$ is related to the temperature $T$ as $T = T_0 = 1/(aN_t)$, where $a$ is the lattice spacing. The size of the lattice along the axis of rotation $N_z$ is arbitrary, and one can take an infinite volume limit in $z$-direction. On the contrary, one has to restrict the size of the lattice in the transverse plane so that the velocity at any point of the considered volume is smaller than the speed of light. Taking into account the rectangular geometry of the lattices, the velocity $v$ at the middle point of the boundary plane is limited as $v = \Omega R < 1/\sqrt{2}$, where $R = a(N_s - 1)/2$.

Boundary conditions are particularly important in the simulations. In our study, we impose periodic boundary conditions in $\tau$- and $z$-directions. What concerns $x$- and $y$-directions, one can impose open, periodic, and Dirichlet boundary conditions. These boundary conditions were studied in paper [34]. It was found that though there is some dependence on boundary conditions, this dependence doesn't play a significant role for bulk observables because of screening. Nevertheless, in order to reduce the influence of boundary on our results, we perform simulations on the lattices with large $N_s$: $N_s/N_t \gtrsim 24$. In this paper mostly, we are going to use open boundary conditions (OBC). In addition to studying the systematic uncertainty of our results, in some cases, periodic boundary conditions (PBC) will be employed.

Action (3.3), as well as its lattice version, are complex-valued functions. As a result, Monte Carlo methods cannot be applied for the direct calculation of the partition function of this system because of the sign problem. In order to overcome this difficulty, we carry out a lattice study of rotating gluodynamics at imaginary values of angular frequency $\Omega_I$,

$$ \Omega_I = -i\Omega. \qquad (3.5) $$

We also introduce the imaginary velocity at the boundary $v_I$:

$$ v_I = \Omega_I R = \Omega_I a \frac{N_s - 1}{2}. \qquad (3.6) $$

The results obtained in this way can be analytically continued to real values of angular frequency. Note that this continuation works only for the bounded systems subjected to
the causality condition [63], which can be formulated for the imaginary angular frequencies on the square-shaped cylinder as $v_I < 1/\sqrt{2}$. The causality condition is satisfied in our lattice simulations.

### 3.3 Phase transition and Polyakov loop

In the absence of global vorticity, Yang-Mills theory has two homogeneous phases: the low-temperature confinement phase and the high-temperature deconfinement phase. In order to distinguish these phases, the expectation value of the Polyakov loop is commonly used

$$L(\mathbf{r}) = \mathrm{Tr} \mathcal{P} \exp \left( \oint_0^{1/T} d\tau A_4(\tau, \mathbf{r}) \right), \qquad (3.7)$$

where $\mathcal{P}$ is the path-ordering operator. On the lattice, it has the following definition:

$$L(\mathbf{r}) = \mathrm{Tr} \left[ \prod_{\tau=0}^{N_t-1} U_4(\tau, \mathbf{r}) \right], \qquad (3.8)$$

where $U_4(\tau, \mathbf{r})$ is the link variable in the temporal direction.

The Polyakov loop can be parameterized by the free energy of static quark $F_Q$ as $\langle L \rangle = e^{-F_Q/T}$. In the confinement phase, infinite energy is needed to create a free quark, whereas, in the deconfinement, the free energy of a static quark is finite. As a result, for static hot gluon matter, the Polyakov loop takes a zero value in the confinement phase, and this quantity is globally nonzero in the deconfinement phase.

The analytical form of the action (3.3) implies that rotating gluodynamics is homogeneous in z-direction and inhomogeneous in the transverse directions. Therefore, the homogeneity of the gluon configurations is lost in the vortical plasmas. However, we believe that the local expectation value of the Polyakov loop remains an order parameter also for rotating gluodynamics, and this local order parameter can be used to distinguish spatially separated confinement and deconfinement phases.

To study the spatial distribution of the Polyakov loop, we introduce the local Polyakov loop in x, y-plane as

$$L(x, y) = \frac{1}{N_z} \sum_{z=0}^{N_z-1} L(x, y, z), \qquad (3.9)$$

where $L(x, y, z) = L(\mathbf{r})$ is defined by Eq. (3.7) and calculate its ensemble average. Using the local Polyakov loop (3.9), one can introduce local Polyakov loop susceptibility,

$$\chi_L(r) = \langle |L(r)|^2 \rangle - \langle |L(r)| \rangle^2, \qquad (3.10)$$

which can also be employed to detect the spatial confinement/deconfinement transition.

Lattice calculations are performed at finite $N_z$, which leads to high uncertainty in the determination of the position of the confinement/deconfinement transition. To reduce the associated statistical error, we calculated the mentioned quantities within a thin cylinder ($r - \delta r/2, r + \delta r/2$). We justified our approach by demonstrating numerically that the finiteness of $\delta r$ brings only a minor systematic error to the estimation of the critical temperature, as we demonstrate below.
<figure><img src="image_2.png" /><figcaption>**Figure 1.** (top) The distribution of the local Polyakov loop in x, y-plane for lattice of size 5 × 30 × 181² with open boundary conditions at the fixed imaginary velocity at the boundary $v_I^2 = (\Omega_I R)^2 = 0.16$ and different on-axis temperatures $T/T_{c0} = 0.80, 0.90, 0.95, 1.05$. (bottom) The Polyakov loop at the x axis ($y=0$). The vertical lines mark the phase boundaries with shaded uncertainties. The violet (blue) data points correspond to periodic (open) boundary conditions. Movies showing the evolution of inhomogeneity as temperature rises are available as ancillary files (see Supplementary material).</figcaption></figure>

In this work, we focus on the confinement–deconfinement phase transition, characterized via the Polyakov loop and its susceptibility. The deconfining transition point is conventionally identified either by the inflection point of the Polyakov loop —marking its change from near-zero to finite values—or by the peak of the susceptibility of the Polyakov loop. These points coincide in the thermodynamic limit for a genuine phase transition. In pure SU(3) Yang–Mills theory, where the transition is weakly first-order, the location of the transition point can be determined directly from the bare (unrenormalized) Polyakov loop and its susceptibility.

# 4 Mixed inhomogeneous phase in rotating gluon plasma

## 4.1 Observation of mixed inhomogeneous phase

Taking into account the explicit coordinate dependence of the action (3.3), one can expect that various local observables are inhomogeneous quantities in the transverse x, y-plane, which is normal to the rotation axis. In this subsection, we study the spatial distribution of the Polyakov loop in the transverse plane following the first study of the inhomogeneous phases in vortical gluodynamics [41]. In this section, we present an improved and detailed analysis of the results obtained earlier in Ref. [41].

We fix the imaginary boundary velocity at the value $v_I^2 = (\Omega_I R)^2 = 0.4$ and vary the on-axis temperature $T$ of our system. The distributions of the local Polyakov loop (3.9) are shown in Fig. 1, where the temperature is given in terms of the critical temperature $T_{c0}$ of a non-rotating gluon matter. The results presented in this subsection were obtained on the lattice $5 \times 30 \times 181^2$. Notice that in this subsection, we only average the Polyakov loop over the $z$-direction.
<figure><img src="image_3.png" /><figcaption>**Figure 2.** The same as in Fig. 1 but for the temperature $T/T_{c0} = 0.95$ and different imaginary velocities at the boundary, $v_I^2 = (\Omega_I R)^2 = 0.04, 0.12, 0.24, 0.48$. Movies on the phase evolution with increasing angular velocity are available as ancillary files (see Supplementary material).</figcaption></figure>

At low temperatures (see the left panel in Fig. 1 with $T/T_{c0} = 0.80$), the system is in the confinement phase, and the local Polyakov loop fluctuates close to zero in the whole space. If we heat the system, it remains in this state up to a certain temperature when an inhomogeneous two-phase structure appears. Then, imaginary rotation results in the deconfinement phase at the periphery of the system, while the region near the rotation axis stays in the confinement phase. As the (on-axis) temperature increases, the radius of the inner region shrinks (see two central panels in Fig. 1 with $T/T_{c0} = 0.90, 0.95$). Finally, at the temperature above the confinement/deconfinement transition of non-rotating gluodynamics, the whole volume is in the deconfinement phase (see the right panel in Fig. 1 with $T/T_{c0} = 1.05$).

In our second numerical calculation, we fix the on-axis temperature at the value $T/T_{c0} = 0.95$, which is slightly below the critical temperature of the non-rotating system, and vary the imaginary boundary velocity $v_I$. The results are presented in Fig. 2. From these figures, one sees that for a slow rotation, $v_I^2 = 0.04$, the whole system is in the confinement phase. For larger velocities $v_I^2 = 0.12, 0.24$, there appears to be a mixed inhomogeneous phase with deconfinement at the periphery and confinement in the center. This picture also holds for the largest imaginary velocity $v_I^2 = 0.48$, which is quite close to the causality limit $v_I^2 < 0.5$. Interestingly, the system in the center stays in the confinement phase even for such a fast rotation.

To summarize, our results indicate that the imaginary rotation of gluodynamics gives rise to the appearance of a mixed-phase state. In this state, the confinement phase is located in the center region, whereas the deconfinement phase is located in the periphery of the studied volume. For imaginary rotation, this mixed-phase state takes place for on-axis temperatures smaller than the critical temperature of the non-rotating system $T_{c0}$. For larger temperatures, $T > T_{c0}$, the mixed-phase state does not emerge, and the whole imaginary-rotating system resides in the deconfinement phase with a minor dependence of the local Polyakov loop on the radius (see the right panel in Fig. 1).
The lower rows of Figs. 1 and 2 show the dependence of our results on boundary conditions. These plots show that boundary conditions affect only small regions adjacent to boundaries. The results for PBC and OBC agree in bulk. The movies showing the evolution of the Polyakov loop distribution due to the variation of simulation parameters for both PBC and OBC are attached as ancillary files to this article (see Supplementary material).

We mention that our results do not depend on the size of the lattice in the $z$ direction. Notice also that as one can see from Fig. 1, the regions with different phases are well spatially separated, and the transition region between these phases has a regular circle form, despite the square geometry of the lattice. Thus, cylindrical symmetry is restored in the bulk of our system for our lattice parameters. The restoration of the rotation symmetry is a natural consequence of the property that we work close to the continuum limit of the lattice Yang-Mills theory.

It is important to discuss the width of the transition region between phases. In our study, we have found that this width can be estimated as ~3–5 fm, and it weakly depends on the parameters of the simulations (see Appendix A). So, one can conclude that in order to observe the mixed inhomogeneous phase in rotating gluon plasma, one needs lattices of rather large sizes in the transverse $xy$ directions.

## 4.2 Study of local critical temperature

The mixed-phase picture assumes that the phase transition at different points of the rotating system occurs at different values of the (on-axis) temperature. To characterize this phenomenon quantitatively, we introduce the local (pseudo)critical temperature $T_c(r)$. This quantity is defined as the temperature on the rotation axis, for which the system, with angular velocity $\Omega_I$, undergoes the confinement/deconfinement phase transition at a distance $r$. This function determines the position of the transition region. Indeed, for a given temperature $T = T_c(r)$, one has the confinement phase at distances smaller than $r$ and the deconfinement phase at larger distances.

In order to find the local critical temperature, we calculate the expectation value of the Polyakov loop $L(r)$ and its susceptibility $\chi_L(r)$. The critical temperature $T_c(r)$ can be found from the peak position in the susceptibility $\chi_L(r)$. Notice, however, that to reach acceptable accuracy, one needs either lattices with large $N_z$ extent or quite large statistics. As was mentioned above, to overcome this problem, we average the Polyakov loop within a thin cylinder of the width $\delta r$. In addition, to reduce the influence of boundary conditions, we discard $\delta b$ layers adjacent to the boundary and do not account for them in the calculation. Technical details of the fitting procedure can be found in Appendix B.

In order to study how our results depend on the simulation parameters, we calculated the critical temperature $T_c(r)$ for the fixed value of the imaginary velocity $v_I^2 = (\Omega R)^2 = 0.16$. The simulation was carried out on the lattice $5 \times 30 \times 181^2$ with periodic/open boundary conditions, the width of ignored boundary layer $\delta b \cdot T = 5$, and the following widths of the circular layer: $\delta r \cdot T = 1, 3, 5$. The results are presented in Fig. 3.

One can see that the results for different $\delta r$ are in good agreement with each type of boundary condition. The smaller width of the averaging layer $\delta r$ allows us to compute the
<figure><img src="image_4.png" /><figcaption>**Figure 3.** The local critical temperature (shown in units of the critical temperature $T_{c0}$ of nonrotating Yang-Mills theory) as a function of radial distance for the periodic (PBC) and open (OBC) boundary conditions, several values of the averaging width $\delta r \cdot T = 1, 3, 5$ and the width of the ignored boundary layer $\delta b \cdot T = 5$ at the fixed velocity of rotation $v_I^2 = (\Omega_I R)^2 = 0.16$. The dotted lines show the fit of the data points by function (4.1).</figcaption></figure>

local critical temperature closer to the rotation axis. However, as can be anticipated, the uncertainty of the calculation with smaller $\delta r$ is larger.

Although the results shown in Fig. 3 were obtained with relatively large width of the ignored boundary layer $\delta b \cdot T = 5$, one can still distinguish some difference between the results for open and periodic boundary conditions at $r/R \sim 1$. However, this difference is quite small, i.e., we can state that boundary effects are under control in our simulations. We analyzed the results for several thicknesses of the boundary layer $\delta b$ and found that for $\delta b \cdot T \gtrsim 3$, boundary effects do not play a significant role.

To fit the data for the dependence of the critical temperature on the distance $r$ from the axis or rotation, we used the following quartic formula

$$ \frac{T_c(r)}{T_{c0}} = C_0 - C_2 \left(\frac{r}{R}\right)^2 + C_4 \left(\frac{r}{R}\right)^4. \qquad (4.1) $$

We found that the quality of the fit normalized by the number of the degrees of freedom (n.d.f.) used in the fit are $\chi^2$/n.d.f. $\sim 1$ for $\delta r \cdot T \ge 3$ and $\chi^2$/n.d.f. $\sim 2-3$ for $\delta r \cdot T = 1$, indicating that the fit (4.1) describes our data reasonably well. The results for the best fit are shown in Fig. 3 by dotted lines.

In addition, we tried a quadratic fit, i.e., by setting the quartic constant to zero, $C_4 = 0$, to describe our data. This approximation appears to be inconsistent with our results if the data is taken in the whole range in $r$. However, the parabolic fit can still be used if one restricts the region of the fit by $r/R \lesssim 0.5$, as it was done in Ref. [41]. In this case, the fit quality is good $\chi^2$/n.d.f. $\sim 1-3$, and the results of the $C_0, C_2$ are in agreement with that obtained through the fitting by the formula (4.1) in the whole available region.
<figure><img src="image_5.png" /><figcaption>**Figure 4.** The local critical temperature as a function of radial distance *r*, given in units of the transverse size of the system, for various imaginary velocities $v_I^2$. The results were obtained on the lattice $5 \times 30 \times 181^2$ with OBC. Dotted lines represent the fit of the data by the equation (4.1).</figcaption></figure>

Further, we consider the dependence of the coefficients $C_i$, $i = 0, 2, 4$ on the lattice spacing $a$ and transverse size of the lattice used in this simulation. We found that the $C_0$ and $C_2$ coefficients depend weakly on these parameters for the OBC and PBC, while the $C_4$ shows certain dependence on particularities of the lattice setup. In Appendix C, we present the relevant details and the results of this analysis. Since our results weakly depend on the boundary conditions, hereafter, we are going to use open boundary conditions, which are more physically motivated than the periodic ones.

Let us now study how the coefficients in the formula (4.1) depend on the velocity of rotation. In Fig. 4 we show the local temperature $T_c(r)$ as a function of radius for various imaginary velocities $v_I^2 = 0.04, \dots, 0.48$. These results were obtained for the averaging width $\delta r \cdot T = 3$ and the ignored boundary layer $\delta b \cdot T = 5$ on the lattice $5 \times 30 \times 181^2$.

The results shown in Fig. 4 are qualitatively similar to those in Fig. 3. However, from Fig. 4, one sees that the behavior of the $T_c(r)$ strongly depends on the rotation velocity $v_I^2$. In particular, the larger the imaginary velocity $v_I^2$, the stronger the dependence on the distance of the local critical temperature $T_c(r)$.

In order to determine the coefficients $C_0, C_2, C_4$, the data in Fig. 4 were fitted by formula (4.1). The quality of these fits is $\chi^2$/n.d.f. $\lesssim 2$. We consider the layer width $\delta r$ as a source of systematic uncertainties and average the values of best-fit coefficients for different $\delta r$ to improve the quality of the data points (see the details in Appendix D). The resulting coefficients $C_0, C_2, C_4$ are shown in Fig. 5.

Our data for the $C_0, C_2, C_4$ can be well described by the linear functions of $v_I^2$:

$$C_0(v_I^2) = a_0 - \kappa_0 v_I^2, \qquad (4.2a)$$

$$C_2(v_I^2) = a_2 + \kappa_2 v_I^2, \qquad (4.2b)$$

$$C_4(v_I^2) = a_4 + \kappa_4 v_I^2. \qquad (4.2c)$$
<figure><img src="image_6.png" /><figcaption>**Figure 5.** The coefficients $C_0, C_2, C_4$ that encode spatial inhomogeneities of the local critical temperature (4.1) as functions of $v_I^2$. Results are obtained for open boundary conditions on the lattice $5 \times 30 \times 181^2$. Dotted lines show linear fits (4.2) of the data. The results of the quadratic fit from Ref. [41] are also presented. The uncertainties include both statistical and systematic contributions.</figcaption></figure>

from which one can determine the parameters $a_0, a_2, a_4$ and $\kappa_0, \kappa_2, \kappa_4$.

It is clear from formulas (4.1) and (4.2) that the coefficients $a_0, a_2, a_4$, which do not depend on the imaginary angular frequency $\Omega_I$, describe the critical temperature of gluodynamics in the absence of rotation:

$$ \frac{T_c(r)}{T_{c0}} = a_0 - a_2 \left(\frac{r}{R}\right)^2 + a_4 \left(\frac{r}{R}\right)^4, \quad [\Omega = 0]. \qquad (4.3) $$

In non-rotating gluon plasma, the inhomogeneities in the critical temperature may only appear due to the presence of the boundaries. Therefore, for a sufficiently large volume, we expect that the boundary effects in bulk are small, indicating that in the thermodynamic limit, $a_0 \to 1$, $a_2 \to 0$, and $a_4 \to 0$. The deviations of the parameters $a_0, a_2, a_4$ from these values would signal the presence of finite volume effects. In our calculation, we found that the parameters $a_2$ and $a_4$ are indeed zero within a small numerical uncertainty while the deviation of the leading coefficient $a_0$ from unity is smaller than 1% (cf. Fig. 5). These properties undoubtedly demonstrate that finite volume effects are small and under control in our study.

Given a moderate dependence of the critical temperature on the boundary conditions, it is important to determine the properties of the vortical hot gluonic matter at lattices that approach the continuum limit. To find the continuum limit for the coefficients $\kappa_0, \kappa_2, \kappa_4$ in addition to the lattice $5 \times 30 \times 181^2$ we repeated our calculations on the lattices
<figure><img src="image_7.png" /><figcaption>**Figure 6.** The continuum limit extrapolation for the fit coefficients $\kappa_i$ in Eq. (4.2). In addition, we show the results from our first study [41] (the blue markers).</figcaption></figure>

4 × 24 × 145², 6 × 36 × 217² that have different size in lattice units as well as different lattice spacings. The results found in this way are presented in Fig. 6.

The continuum limit for the $\kappa_0$ appears to be zero within the small numerical uncertainties. Taking into account the formula (4.2a), this property implies that the coefficient $C_0$ is unity up to the finite volume effects. Noticing that the coefficient $C_0$ determines the ratio $T_c(r = 0)/T_{c0}$ for rotating gluodynamics we conclude that the critical deconfining temperature on the rotating axis $T_c(r = 0)$ for rotating gluon matter coincides with the critical temperature of the same system without rotation. This result should be confronted with perturbative estimations of Ref. [15] where such dependence has been found.

Now let us consider the vortical coefficients $\kappa_2$ and $\kappa_4$ that appear in Eqs. (4.2b) and (4.2c), respectively. The continuum limit for these coefficients, obtained through the quartic fit of the data, gives us the following values:

$$ \kappa_2 = 1.051(29), \quad \kappa_4 = 0.300(34). \qquad (4.4) $$

It is worth comparing this result to the value

$$ \kappa_2 = 0.902(33), \qquad (4.5) $$

which has been obtained in Ref. [41] through the quadratic ($C_4 = 0$) fit of the data in the bulk with $a_2 = 0$. We also found that the rotation-dependent correction $\kappa_0$ to the on-axis critical temperature $C_0$, given by Eq. (4.2a), has a vanishingly small value, $\kappa_0 = 0.008(6)$.
To summarize the results of this section, the critical temperature for rotating gluodynamics can be well described by the following formula<fnref n="3" />:

$$ \frac{T_c(r\Omega_I)}{T_{c0}} = 1 - \kappa_2(\Omega_I r)^2 + \kappa_4(\Omega_I r)^2 \left(\frac{r}{R}\right)^2. \qquad (4.6) $$

Our data suggests that the coefficient $\kappa_2$ is universal in bulk, i.e., it does not depend on the lattice parameters such as the transverse lattice size and the type of boundary conditions. On the contrary, the quartic coefficient $\kappa_4$ depends on these parameters of the system.

Equation (4.6) can be rewritten as follows:

$$ \frac{T_c(r, \Omega_I)}{T_{c0}} = 1 - (\Omega_I r)^2 \left( \kappa_2 - \kappa_4 \left( \frac{r}{R} \right)^2 \right). \qquad (4.7) $$

In this form, the quartic $\kappa_4$-term resembles a finite volume correction to the $\kappa_2$ coefficient. In addition to the observed $\kappa_4$-term, the $r^4$-dependence of the critical temperature can be parameterized by the term $\sim (\Omega_I r)^4$. In our simulations, the dominant contribution results from the $\kappa_4$ term, and, with our numerical accuracy, we are not able to distinguish it from the suspected $\sim (\Omega_I r)^4$ term. However, one can expect that there still exists a contribution of the form $\sim (\Omega_I r)^4$ to the critical temperature, which will be important for sufficiently large angular velocity. We will discuss this $(\Omega_I r)^4$ contribution in Section 6 within the approximation of local thermalization.

## 4.3 The phase diagram for real and imaginary rotation

We have already discussed that for imaginary rotation, the mixed inhomogeneous phase can be observed for the on-axis temperatures below the critical temperature of non-rotating gluodynamics: $T = T_{c0} - \Delta T$ with $\Delta T > 0$. In this case, the confinement/deconfinement phases are localized to the central/periphery regions of the studied volume. The critical distance $r$, where the transition from the confinement to the deconfinement takes place, can be found from the equation: $T_{c0} - \Delta T = T_c(r, \Omega_I)$, where the $T_c(r, \Omega_I)$ is given by (4.6) and $\Omega_I$ and $\Delta T$ serve as the parameters of this equation.

To carry out the analytical continuation of the phase diagram from imaginary to real rotation, we make the following substitution:

$$ \Omega_I^2 \rightarrow -\Omega^2. \qquad (4.8) $$

Following the analytical continuation (4.8) applied to Eq. (4.6), one finds that for the rotation with a real angular velocity, the deconfinement phase is localized in the central region, while the confinement phase is in the periphery of the system. As we mentioned in Ref. [41] and also will discuss below, this is an exotic phenomenon that is totally unexpected from generic considerations based on the Tolman-Ehrenfest law. The critical distance from the rotation axis, where the spatial deconfinement/confinement transition appears, is the solution of the following equation: $T_{c0} + \Delta T = T_c(r, \Omega)$. This procedure gives us the phase

<footnote n="3">³Previous lattice studies [33–35, 37] have found a dependence of the critical temperature $T_c$ on $\Omega_I$ without accounting for a dependence of the phase structure on the distance from the center of rotation. These findings should be interpreted as bulk-averaged critical temperatures.</footnote>
<figure><img src="image_8.png" /></figure>

**Figure 7.** The critical radius $r$ of the spatial boundary between the confining and deconfining phases for, simultaneously, the imaginary frequencies $\Omega_I$ (on-axis temperature $T = T_{c0} - \Delta T$) and the real frequencies $\Omega$ (on-axis temperature $T = T_{c0} + \Delta T$) for various temperature offsets $\Delta T > 0$. The shading corresponds to a degree of uncertainty around the central curve. The size of the system is $R = 13.5$ fm, which is typical for our lattice simulations. The regions with the horizontal hatching represent the diagram from Ref. [41], obtained with the help of a quadratic fit [i.e., with $\kappa_4 \equiv 0$ in Eq. (4.6)].

diagram, which is shown in Fig. 7.<fnref n="4" /> Note that Eq. (4.6) is a quadratic function in $\Omega_I$. Therefore, for a given temperature offset $\Delta T$, the diagram for the real angular velocity at the temperature $T_{c0} + \Delta T$ has the same shape as for the imaginary rotation at the symmetric values of temperature $T_{c0} - \Delta T$.

In Fig. 7, we also plot the diagram using the results of quadratic fit (i.e. formula (4.6) without $\kappa_4$-term) from Ref. [41]. The minor difference between the quadratic and quartic results appears for large velocities of rotation, mostly near the boundaries $r/R \gtrsim 0.8$. Nevertheless, these results are in agreement with each other within the statistical uncertainties.

# 5 Decomposition of rotating action

## 5.1 Mechanical vs. chromomagnetic vortical effects.

It is worth asking what is the origin of the peculiar behaviour of gluon plasma subjected to rotation? More concretely, which terms in the action are responsible for the statistical properties of gluonic gas? To answer this question below, we decompose the action into mechanical and magnetic terms, similar to Eq. (2.9), and analyze their effects associated with rotation separately.

<footnote n="4">⁴In constructing the diagram in Fig. 7, we include a minor deviation of $a_0$ from unity due to finite volume effects to uncertainty.</footnote>
Similarly to Lagrangian (2.9), the action of rotating gluodynamics (3.4) is a quadratic function in the imaginary angular frequency $\Omega_I$:

$$S(\Omega_I) = S_0 + S_1\Omega_I + S_2\Omega_I^2, \qquad (5.1)$$

where $S_0$ is a well-known action in gluodynamics without rotation. The contributions $S_1$ and $S_2$ in the cylindrical coordinates can be written as

$$S_1 \equiv S_{\text{mech}} = -\frac{1}{g_{YM}^2} \int d^4 x \, r \left[ F_{\phi r}^a F_{r \tau}^a + F_{\phi z}^a F_{z \tau}^a \right], \qquad (5.2)$$

$$S_2 \equiv S_{\text{magn}} = \frac{1}{2g_{YM}^2} \int d^4 x \, r^2 \left[ (F_{\phi z}^a)^2 + (F_{r \phi}^a)^2 \right]. \qquad (5.3)$$

The terms $S_1$ and $S_2$ have a clear physical meaning [45]. The first term, $S_1$, is related to the total angular momentum of gluons and, according to the discussion in Section 2, represents a linear “mechanical” coupling with the angular velocity. The second term, $S_2$, is given by the squares of two components of chromomagnetic field $(B_r^a)^2$ and $(B_z^a)^2$. This term represents a quadratic “chromomagnetic” coupling with the angular velocity, which stresses the highlighted role of the chromomagnetic condensate in the rotational motion of the gluon plasma.

In papers [44–47], the moment of inertia of gluon plasma was studied. It was found that the mechanical coupling $S_1$ and the chromomagnetic coupling $S_2$ have opposite influences on the moment of inertia. In particular, the term that couples to the angular velocity linearly, $S_1$, gives a positive contribution to the moment of inertia. On the contrary, due to the specific behavior of magnetic gluon condensate, the contribution of the term which couples to the angular frequency quadratically, $S_2$, is negative up to the temperatures $\sim 2T_c$. It was also suggested that these operators pull the critical temperature in opposite directions: for imaginary rotation, the linear mechanical coupling to rotation $S_1$ leads to an increase of the critical temperature in agreement with our intuition, while the quadratic chromomagnetic coupling $S_2$ tends to decrease the critical temperature.

In this section, we are going to study the roles of the $S_1$ and $S_2$ operators in the mixed inhomogeneous phase of rotating gluodynamics. To do this we introduce the factors $\lambda_1$ and $\lambda_2$ to the action (5.1):

$$S_G = S_0 + \lambda_1 S_1 \Omega_I + \lambda_2 S_2 \Omega_I^2. \qquad (5.4)$$

These factors, which take the values $\lambda_1 = 0, 1$ and $\lambda_2 = 0, \pm 1$, allow us to switch on/off the $S_1$ and $S_2$ operators in the simulation. In particular, the following regimes are possible:

* Im1) $\lambda_1 = 1, \lambda_2 = 0$ (the linear coupling only): The imaginary rotation acts only through the linear term $S_1$;

* Im2) $\lambda_1 = 0, \lambda_2 = 1$ (the quadratic coupling only): The imaginary rotation acts only through the quadratic term $S_2$;

* Im12) $\lambda_1 = 1, \lambda_2 = 1$ (the full action): This case, which corresponds to the imaginary rotation of the whole action, does not require additional study since it has been considered in Section 4.
<figure><img src="image_9.png" /></figure>

**Figure 8.** The distribution of the local Polyakov loop in x, y-plane for the lattice size 5 × 30 × 181² with OBC at the fixed absolute value of rotational velocity |v²| = |ΩR|² = 0.16 and various regimes corresponding, from the left plot to the right plot to the linear imaginary coupling (Im1), the quadratic imaginary coupling (Im2), the full action in rotation (Im12), and the real angular velocity acting only the quadratic-coupling term (Re2), respectively. On each panel, the temperatures are chosen in such a way that the mixed phase occurs. Movies on the phase evolution with increasing temperature for different regimes are available as ancillary files (see Supplementary material).

Note that in the case of Im2, there is, actually, no sign problem. Therefore, for this choice of the factors $\lambda_i$, we can conduct lattice simulation with real angular velocity. In order to study this case, we introduce one more regime

Re2) $\lambda_1 = 0$, $\lambda_2 = -1$ (real rotation for the quadratic coupling): The rotation acts only through the quadratic term $S_2$, but due to the negative value of the coupling, the angular velocity takes a real value: $\Omega^2 = -\lambda_2\Omega_I^2 > 0$.

We stress that the separation of the action into three different parts (5.4) does not violate the gauge invariance of the action since all three contributions, $S_0$, $S_1$ and $S_2$ are invariant under SU(3) gauge transformations.

## 5.2 Inhomogeneous phases

The distributions of the local Polyakov loop for different regimes Im1, Im2, Im12, and Re2 are shown in Fig. 8. The simulations were carried out at the fixed value of the velocity $v_I^2 = 0.16$ for the Im1, Im2, Im12 and $v_I^2 = -0.16$ for the Re2 on the lattice 5 × 30 × 181² with OBC. We have adjusted the temperatures on the central axis (indicated in the plots) to observe the mixed-phase structures. The movies with the phase evolution in different regimes for both OBC and PBC are attached as ancillary files (see Supplementary material).

**Linear coupling at imaginary rotation.** In the first regime, Im1, the mixed-phase structure can be observed for temperatures $T \gtrsim T_{c0}$, which are slightly above the critical temperature of non-rotating gluodynamics. This mixed phase arises only in a narrow temperature range because the effects of linear coupling are relatively weak (in the left panel of Fig. 8, we show the local Polyakov loop distribution for the temperature $T/T_{c0} = 1.014$, for which the mixed phase is clearly visible). However, in this case, the deconfinement phase is localized in the center region, whereas the confinement phase arises at the periphery of the system. As one increases the temperature at the rotation axis, the deconfinement phase
<figure><img src="image_10.png" /><figcaption>**Figure 9.** The effect of the quadratic-only coupling $S_2$ of the vorticity to the gluon action: The distribution of the local Polyakov loop as a function of $x$ coordinate for lattice size $5 \times 30 \times 181^2$, open (OBC) and periodic (PBC) boundary conditions for the Im2 and Re2 regimes. Data are shown, respectively, for $|v^2| = 0.16$ at $T/T_{c0} = 0.95$ with imaginary angular velocity (Im2) and $T/T_{c0} = 1.05$ with real angular velocity (Re2).</figcaption></figure>

captures larger and larger volumes. Finally, at some temperature, the whole system transits to the deconfinement phase; for the temperatures, $T \lesssim T_{c0}$, the mixed-phase structure does not appear. Thus, the mixed-phase structure of the Im1-regime has properties opposite to that observed in Section 4, where the imaginary rotation with the full action (the case Im12) has been studied.

**Quadratic coupling at imaginary rotation.** A different picture is observed for the Im2-regime, where the rotation couples only to the chromomagnetic component of the gluon field. In this case, the mixed-phase structure has properties qualitatively similar to that of rotating gluodynamics (Im12-regime studied in Section 4): the confinement in the central region and the deconfinement at the periphery. However, quantitatively, the Im2-regime differs a bit from the regime Im12 (see two central panels in Fig. 8, where the phase structures are shown for these regimes at the same value of the (on-axis) temperature $T/T_{c0} = 0.95$).

**Quadratic coupling at real rotation.** The Re2-regime is designed to realize real rotation for the $S_2$ term. For a fixed $|v_I|$, the gluon matter has properties opposite to those of the Im2-regime: the deconfinement arises in the central region, and the confinement is located at the periphery. In addition, in the Re2-regime, the mixed phase appears for the temperatures $T > T_{c0}$, which differ from the mixed-phase domain in the Im2-regime, $T < T_{c0}$. In the right panel of Fig. 8 the distribution of the local Polyakov loop is shown for $T/T_{c0} = 1.05$.

To demonstrate the opposite behavior of the regimes Im2 and Re2, in Fig. 9, we show the local Polyakov loop as a function of $x$ coordinate for both regimes and OBC/PBC.
<figure><img src="image_11.png" /><figcaption>**Figure 10.** The local critical temperature as a function of radial distance for different regimes of rotation. The results from the lattice 5 × 30 × 181² with open (OBC) and periodic (PBC) boundary conditions at fixed rotation velocity $|v^2| = |\Omega R|^2 = 0.16$. Dotted lines represent the fit of the data by the function (4.1).</figcaption></figure>

The results were obtained for the positive (negative) shift of the on-axis temperature $\Delta T = 0.05 T_{c0}$ and $|v_I^2| = 0.16$, i.e. $T/T_{c0} = 0.95$, $v_I^2 = 0.16$ and $T/T_{c0} = 1.05$, $v_I^2 = -0.16$. One sees that the confinement phase translates into the deconfinement phase after analytic continuation and vice versa. The boundary between the phases remains at the same position. Notice also that the spatial boundary conditions affect the results only in the regions close to the boundary within a finite-temperature correlation length.

**Phase boundaries.** Using the methods developed in the previous section, we compute the local critical temperatures for all regimes considered in this section. The results of this calculation are shown in Fig. 10. It is seen from this plot that the rotation in regime Im1 increases the critical temperature of gluodynamics, whereas in regime Im2, it decreases due to the rotation. This result confirms our suggestion proposed in paper [44] that the term with total gluon angular momentum $S_1$ increases, while the chromomagnetic field contribution $S_2$ decreases the critical temperature. It is also seen from Fig. 10 that the radial dependence of the $T_c(r)$ for the Im1-regime is very weak. For this reason, the local critical temperatures of the regimes Im2 and Im12, which differ by the accounting of $S_1$ term, are close to each other. The next points to be discussed are the Im2- and Re2-local critical temperatures. Our results for the critical temperature indicate that within the uncertainty, the results for regimes Im2 and Re2 are connected through the analytical continuation procedure $\Omega^2 \leftrightarrow -\Omega_I^2$ in Eq. (4.6).<fnref n="5" />

Briefly summarizing this section, we conclude that the standard linear coupling of the

<footnote n="5">⁵This agreement is excellent up to the large radial distances $r/R \sim 0.8$; the minor discrepancy appears only in the vicinity of the boundary. It may indicate that the next-order terms may be needed to reconstruct the boundary effects fully after analytic continuation.</footnote>
vorticity to the angular momentum of the gluons, represented by the action $S_1$, Eq. (5.2), leads to an increase of the critical temperature of the deconfinement phase transition as the imaginary angular frequency increases. This effect is anticipated from earlier analytical studies [39] and, at the same time, it is opposite to the first-principle results of lattice simulations of the full gluonic action [33, 34]. The reason why the imaginary rotation leads to the decrease of the critical temperature in fully-rotated gluodynamics appears to be related to the dominance of the quadratic coupling of the vorticity to the chromomagnetic component of the gluon fields, which is characterized by the action $S_2$, Eq. (5.3). This result also agrees with our earlier finding that the mechanical properties of the rotating gluon plasma are predominantly influenced by the chromomagnetic condensate [44, 47].

# 6 Local thermalization of rotating gluodynamics

From a simple kinematic point of view, the influence of the rotation on the gluonic fields has two major effects: Yang-Mills action (3.4) becomes both inhomogeneous and anisotropic. In details,

(i) Inhomogeneity: the gluonic action develops an explicit dependence of the coefficients in front of various gluonic operators on the distance from the rotation axis $r$;

(ii) Anisotropy: The rotation affects chromoelectric and chromomagnetic components of gluonic fields differently.

In the previous sections, we studied both effects simultaneously. In this section, we will concentrate mostly on the impact of anisotropy on the gluonic action in the rotating reference frame.

## 6.1 Local approximation for inhomogeneous action

The results presented in the previous sections were obtained on the lattices with quite large extensions in the transverse directions, which complicates the simulation considerably. The study of rotating gluodynamics on large lattices is necessary since we put the axis of rotation to the center of the lattice and thus effectively half the lattice sizes in the $x$- and $y$-directions. Furthermore, the phenomena studied in this paper cannot be observed on small lattices because of the rather wide confinement/deconfinement transition region and boundary effects. Another property which complicates the study of rotating gluodynamics is an explicit dependence of its action (3.3), (3.4) on the coordinate $r$, which complicates theoretical study of this system.

In order to simplify the situation from both numerical and theoretical perspectives, we introduce the following approximation. Let us consider the system at distances $r \sim r_0 \gg \zeta$ from the rotation action, where $\zeta$ is the correlation length of gluodynamics. It is clear that in this approximation, the coefficients of the gluon fields operators in action (3.4), which depend on the product $\Omega r$, remain approximately the same if one varies the distance by few $\zeta$ in the vicinity of $r_0$. Therefore, it is reasonable to assume that in this region, thermodynamic parameters of the system correspond to the action (3.4) with the
coefficients fixed to their values at a distance $r = r_0$. One can expect that this idea works for small distances $r_0 \lesssim \zeta$ as well. This property can be seen as follows: if one considers slow rotation $\Omega_I \zeta \ll 1$, which is always valid in our simulations, the coefficients preceding various terms in action (3.4) do not change considerably. One could expect a prevailing contribution coming from the linear $\Omega_I r$ term over the terms that are quadratic in $\Omega_I r$. However, according to the results of the previous section, this linear term does not play an important role for the spatial layers located at small distances from the rotation axis.

In this section, we are going to study rotating gluodynamics through the simulation of the action (3.4) with the radius $r$ fixed at some value $r_0$. In this case, the action does not depend on the coordinate explicitly, and the system becomes homogeneous. Below, we will call this approximation as *local thermalization*. Actually, one can fix the coefficients at any point on the circle $r = r_0$. In our simulation, we choose the point $x = r_0$, $y = 0$. The local velocity at this point is $u_I = \Omega_I r_0$, and the homogeneous local action to be simulated can be written as

$$S_G = \frac{1}{2g_{YM}^2} \int d^4 x \left[ F_{x\tau}^a F_{x\tau}^a + F_{y\tau}^a F_{y\tau}^a + F_{z\tau}^a F_{z\tau}^a + F_{xz}^a F_{xz}^a + \right. \\ \left. + (1+u_I^2) F_{yz}^a F_{yz}^a + (1+u_I^2) F_{xy}^a F_{xy}^a - 2u_I (F_{yx}^a F_{x\tau}^a + F_{yz}^a F_{z\tau}^a) \right]. \quad (6.1)$$

Notice that in this approximation, the observables depend on the $u_I$, and they are equal to the ones in non-rotating gluodynamics at $u_I = 0$.

Lattice simulation with the action (6.1) can be carried out using standard lattice methods with periodic boundary conditions in all directions. We discretize the local action (6.1) in the same way as the action of full rotating system (3.3) (see Section 3.2). We determine the critical lattice coupling from the peak of the Polyakov loop susceptibility. Then, using the $\beta$ function of non-rotating gluodynamics, we calculate the critical temperature $T_c$, which can be interpreted as a critical temperature on the axis of rotation.

## 6.2 The results obtained with the local action

First, we carry out the study of the local thermalization approximation on the lattice $5 \times 20^3$. Note that the action (6.1) supports the decomposition using Eq. (5.4) and the Im1-, Im2-, Im12-, Re2- regimes proposed in the previous section can be realized. In Fig. 11, we show the results for the (on-axis) critical temperature as a function of the $u^2 = -u_I^2$ calculated for different regimes. For comparison, the results from Fig. 10, which were computed on the lattice $5 \times 30 \times 181^2$ with the action (3.3), are also shown. It is seen from Fig. 11 that the results for the actions (5.4) and (6.1) are in good agreement with each other in the shown region $|u^2| \lesssim 0.2$. This result confirms that for sufficiently small velocity $u^2$, the local thermalization approximation works within the uncertainty of the calculation for all the regimes.

In the studied region of rotational velocity $|u^2| \lesssim 0.5$ the data are well described by the polynomial:

$$\frac{T_c(u)}{T_{c0}} = 1 + k_2 u^2 + k_4 u^4. \qquad (6.2)$$
<figure><img src="image_12.png" /><figcaption>**Figure 11.** The critical (on-axis) temperature as the function of the local velocity $u^2 = -u_I^2$ for various regimes: Im1, Im2, Im12, Re2. The data for open (OBC) and periodic (PBC) boundary conditions are calculated in Section 5 on the lattice $5 \times 30 \times 181^2$. The data obtained within the local thermalization approximation are calculated on the lattice $5 \times 20^3$. The dashed curves represent the fitting of the data by the formula (6.2).</figcaption></figure>

In Eq. (6.2), we denote the vortical curvature $\kappa$ by $k$ to distinguish the methods. Here, a few comments are in order. First, at $u_I^2 = 0$, the system under study coincides with non-rotating gluodynamics; therefore, $T_c(u=0)/T_{c0} = 1$. Next, although in the equations (4.6) and (6.2) the $\kappa_2$ and $k_2$ are similar, the terms $\kappa_4$ and $k_4$ have different physical meanings. The former results from boundary contribution to the $\kappa_2$ coefficient, while the latter describes the influence of rotation on the critical temperature. We would like to reiterate that the data from Section 4 does not allow us to find the term analogous to the $k_4$. However, in the local thermalization approximation, which is free from boundary effects, the $k_4$ term can be studied.

The simulation of rotating gluodynamics in the Re2-regime allows one to carry out a lattice study with real angular velocity. Notice that this is not the rotation of the full action (6.1), since in this regime, we disregard the linear in $\Omega$ term. Therefore, it is rather difficult to transfer the conclusions from regime Im2 to physical regime Im12, especially in the region of large velocities $u^2 \sim 1$. However, the Re2-regime is straightforwardly related to the Im2-regime through the procedure of analytical continuation. Thus, one can inspect this procedure for different analytic functions.

In Fig. 12, we show the critical temperatures as a function of $u^2$ for these regimes. The region $u^2 < 0$ corresponds to the rotation with imaginary angular velocity realized in the Im2-regime, while the region $u^2 > 0$ represents a real rotation in the Re2-regime.

To fit the data shown in Fig. 12 in the region $u^2 < 0$, we tried to implement the formula (6.2). Then, we continue the result, Eq. (6.2), to the domain of real velocity $u^2 > 0$ (the green hatched region in Fig. 12) and compare it with lattice data for regime Re2. It turns out that the analytic continuation can describe the data only up to $u^2 \lesssim 0.2$.
<figure><img src="image_13.png" /><figcaption>**Figure 12.** The critical temperature for the Im2- and Re2-regimes as a function of $u^2$ (the red open circles). The region $u^2 < 0$ corresponds to the rotation with imaginary angular velocity realized in the Im2-regime. The region $u^2 > 0$ represents the real rotation in the Re2-regime. The calculation was conducted on the lattice $5 \times 20^3$. The hatched regions represent an analytic continuation of the results from region $u^2 < 0$ to region $u^2 > 0$, based on the functions (6.2) and (6.3). The red dashed line shows the best fit of all the data by the rational function (6.3).</figcaption></figure>

The addition of the next order $k_6 u^6$ term does not improve the situation considerably.

In our analysis, we found that the rational function with a pole structure,

$$ \frac{T_c(u)}{T_{c0}} = \frac{1 + c_2 u^2}{1 - b_2 u^2}, \qquad (6.3) $$

can predict the data for $u^2 > 0$ from the results for $u^2 < 0$ reasonably well (see blue hatched region in Fig. 12). We also note that this rational function describes all data from regimes Im2 and Re2 simultaneously (see red dashed line in Fig. 12), whereas the polynomial function (6.2) is applicable only for each subset of data $u^2 \ge 0$ separately. Interestingly, the fit of all data gives $b_2 \approx 0.9$, i.e., the function has a pole at $u^2 \ge 1$.

In order to study how the local thermalization approximation works in the continuum limit, we calculated the critical temperature for the local action (6.1) on the lattices $4 \times 16^3$, $5 \times 20^3$, $6 \times 24^3$ and $8 \times 32^3$. The results are presented in Fig. 13. In addition, we plot the critical temperatures calculated with the full action (3.3) on the lattices $4 \times 24 \times 145^2$, $5 \times 30 \times 181^2$ and $6 \times 36 \times 217^2$ in Section 4. For consistency, we consider the region $u_I^2 < 0.5$ for local action. It is seen from Fig. 13 that the results are in agreement with each other. Further, we fit the data in Fig. 13 by the function (6.2) and (6.3) for each $N_t$ and extrapolate the coefficients to the continuum limit. Our results for the coefficients $k_2, k_4$ are

$$ k_2 = 0.869(31), \quad k_4 = 0.388(53), \qquad (6.4) $$

and the coefficients $c_2, b_2$ take the values

$$ c_2 = 0.206(66), \quad b_2 = 0.694(101). \qquad (6.5) $$
<figure><img src="image_14.png" /><figcaption>Figure 13. The (on-axis) critical temperature, calculated for the system with local action (6.1) (filled points), as a function of the $u_I$. The dotted (dashed) lines represent the best fit of the data by formula (6.2) (by function (6.3)). The data are obtained on the lattices $4 \times 16^3$, $5 \times 20^3$, $6 \times 24^3$, $8 \times 32^3$. In addition, we show the critical temperatures calculated with the full action (3.3) (empty points) on the lattices $4 \times 24 \times 145^2$, $5 \times 30 \times 181^2$, $6 \times 36 \times 217^2$ in Section 4.</figcaption></figure>

Although the function with singularity better describes the data in regimes Im2 and Re2, we do not have any strong theoretical reasons to prefer one of these options for full action. Hence, the difference between the functions (6.2) and (6.3) for $u^2 > 0$ may be considered as a systematic uncertainty of the analytic continuation procedure.

## 6.3 Inconsistency with the Tolman-Ehrenfest approach

According to the Tolman-Ehrenfest (TE) law [59, 60] in thermodynamic equilibrium in the rotating reference frame, the local temperature is given by the formula

$$T(r) = \frac{T_0}{\sqrt{1 - (\Omega r)^2}} = \frac{T_0}{\sqrt{1 + (\Omega_I r)^2}}. \qquad (6.6)$$

The first equality corresponds to the real rotation, while the second one corresponds to the imaginary rotation.

The TE law (6.6) suggests that real rotation effectively heats the system outside of the rotation axis, which naturally leads to the inhomogeneous phase in the rotating plasma in thermal equilibrium. Indeed, the phase transition at some distance $r$ from the rotation axis is expected to occur when the effective local temperature $T(r)$ is equal to the critical temperature for the non-rotating system $T_{c0}$. This assumption leads to the following prediction for local critical temperature:

$$\frac{T_c(u)}{T_{c0}} = \sqrt{1 - u^2}. \qquad (6.7)$$

As a result, one can expect that for real rotation, deconfinement appears at the periphery and confinement at the center. This mixed-phase structure, described by Eq. (6.7), was
suggested in Ref. [39] on the basis of a calculation in a low-dimensional confining model. Nevertheless, this theoretically motivated phase structure is in contradiction with the numerical simulations of this paper. Thus, a naive implementation of the TE law, which leads to formula (6.7), does not describe correctly the properties of the rotating gluodynamics found in our first-principle numerical simulations. We conclude that the thermodynamics of rotating plasma is governed not only by the TE law but also by other factors that we discuss later.

In addition to the observed inconsistency with the TE prediction, there are many theoretical papers where it was shown that real rotation decreases the bulk-averaged critical temperature (see, for instance, [9, 13, 25]). In view of these theoretical works, it is interesting why, within lattice simulation, the average over the volume critical temperature is increased by real rotation.

To answer this question, take a look at the action (6.1). We have already found that the linear in the imaginary frequency $\Omega_I$ term does not affect significantly the critical temperature, and, therefore, it can be disregarded. The action (6.1) without the linear term can be written in the following form

$$S_G = \int d^4 x \left[ \beta \left( (F_{x\tau}^a)^2 + (F_{y\tau}^a)^2 + (F_{z\tau}^a)^2 + (F_{xz}^a)^2 \right) + \tilde{\beta} \left( (F_{yz}^a)^2 + (F_{xy}^a)^2 \right) \right], \quad (6.8)$$

where the effect of rotation is encoded in the effective couplings

$$\tilde{\beta} = (1 - \Omega^2 r_0^2)\beta \equiv (1 + \Omega_I^2 r_0^2)\beta, \quad \beta = \frac{1}{2g_{YM}^2}. \qquad (6.9)$$

The coupling $\tilde{\beta} = \tilde{\beta}(r_0)$ depends on the distance $r_0$ from the rotation axis, thus encoding the effect of the curved metric in the co-rotating coordinates.

Thus, the external gravitational field (2.1), associated with the rotational motion, generates asymmetry in the coupling constants $\tilde{\beta}$ and $\beta$ of different components of the field-strength tensor squared $(F_{\mu\nu})^2$ that enter the action (6.8). The thermodynamic effect of the asymmetry in couplings (6.9) can be assessed in the lattice gauge theory, which reveals that the anisotropy generates a rather peculiar effect: it leads to a change in the lattice spacing in the spatial and temporal directions [64], and consequently, can affect the critical temperature of the deconfinement transition. While the discussed gluon action (6.8) differs from the one considered in Ref. [64], the remarkable effect of the action anisotropy can indeed be seen from our numerical results on the critical temperature shown in Fig. 12. In an analytical form, this effect is excellently reproduced by the rational function (6.3) with the fitting parameters (6.5).

Figure 12 indicates that the larger the asymmetry $\tilde{\beta}/\beta$, the smaller the critical temperature. Taking into account that, for imaginary rotation, $\Omega_I$, the asymmetry of the couplings (6.9) increases with the rise in the distance $r_0$ from the rotation action, $\tilde{\beta}/\beta = 1 + \Omega_I^2 r_0^2$, we arrive to the immediate conclusion that the periphery of the system transits to the deconfinement phase earlier than the central regions, and the average critical temperature decreases with imaginary rotation. These arguments remain true for real rotation, but in this case, they work in the opposite direction: for real rotation $\tilde{\beta}/\beta < 1$,
the confinement/deconfinement phases take place in the periphery/center regions, and the
averaged critical temperature is increased.

It should be noted here that the asymmetry generated by the external gravitational
field is related to the magnetovortical term (5.3) in the gluon action. In the absence of this
term, the mixed phase structure would be in a qualitative agreement with the TE prediction
(see the plot for Im1-regime in Fig. 8). Notice also that the same asymmetry terms, which
are associated with the chromomagnetic fields, give rise to the negative moment of inertia
measured in papers [44–46].

Thus, we conclude that the Tolman-Ehrenfest law is not the only factor that deter-
mines the thermodynamics of the rotating gluon plasma. The additional factor is the
external gravity that generates the asymmetry in the gluon action, which, in turn, in-
fluences the dynamics of gluon matter. This mechanism cannot be accounted for by the
simple Tolman-Ehrenfest formula (6.6). Therefore, the Tolman-Ehrenfest effect, if applied
straightforwardly, would lead to a qualitatively wrong prediction for the behavior of the
inhomogeneous phase of gluon plasma. In this sense, the Tolman-Ehrenfest picture is
“violated” in the rotating gluon plasma.

7 Discussion and conclusion

Our work focuses on the lattice investigation of the confinement/deconfinement phase tran-
sition in rigidly rotating gluon plasma. To incorporate the rotation in the system, we em-
ployed a reference frame that is co-rotating with the gluon plasma, effectively describing
the rotational effects by an external gravitational field with a suitable metric. Since direct
lattice simulations of the full gluonic action with nonzero angular velocity are hindered
by the sign problem, we performed numerical simulations using an imaginary angular ve-
locity. The results at real angular velocities are subsequently recovered via an analytic
continuation.

We have found that there appears to be a mixed-phase state for certain simulation
parameters, thus confirming and refining our earlier studies [41]. For *imaginary* rotation in
this state, the confinement phase is located in the center region, whereas the deconfinement
phase is localized in the periphery of the studied volume. The mixed-phase state can be
realized for on-axis temperatures smaller than the critical temperature of the non-rotating
gluodynamics $T_{c0}$. For the higher temperatures $T > T_{c0}$, the mixed inhomogeneous phase
does not appear.

To describe the inhomogeneous phase, we introduced the local (pseudo)critical tem-
perature $T_c(r)$. This is the temperature on the rotation axis, for which the system, with
imaginary angular velocity $\Omega_I$, undergoes the confinement/deconfinement phase transition
at the distance $r$. From the construction, it is clear that for a given temperature $T = T_c(r)$,
one has the confinement phase at distances smaller than $r$ and the deconfinement phase at
larger distances.

We calculated the local critical temperature for various parameters of the simulations
from the peak of the local Polyakov loop susceptibility. Disregarding the finite volume
effects, which do not play an important role in the bulk of the studied volume, our results
can be well described by the formula

$$ \frac{T_c(r)}{T_{c0}} = 1 - \kappa_2 (\Omega_I r)^2, \quad [\text{imaginary rotation}], \qquad (7.1) $$

where $r$ is the distance from the boundary to the rotation axis. The vortical curvature $\kappa_2$ weakly depends on the parameters of the simulations, in particular, on boundary conditions, and seems to be a universal quantity. The results imply that the local critical temperature on the axis of rotation $T_c(r = 0)$ within the uncertainty coincides with the critical temperature of non-rotating gluodynamics $T_{c0}$. For the system with open boundary conditions, we also measured the sub-leading term, $\sim \kappa_4 r^4$, (see Eq. (4.6) with the coefficients (4.4)), but we expect that the value of coefficient $\kappa_4$ strongly depends on the boundary conditions and it gives a minor contribution in bulk.

Formula (7.1) allowed us to carry out analytical continuation from imaginary to real rotation through the substitution $\Omega_I^2 \to -\Omega^2$:

$$ \frac{T_c(r)}{T_{c0}} = 1 + \kappa_2 (\Omega r)^2, \quad [\text{real rotation}], \qquad (7.2) $$

which explicitly indicates that the local critical temperature growths increase both in the distance from the rotation axis $r$ and the angular velocity $\Omega$. The coefficient $\kappa_2$ in Eqs. (7.1) and (7.2) is given in (4.5). Using Eq. (7.2), we concluded that the mixed-phase state can be realized for real rotation. However, in this case, the inhomogeneous phase appears at temperatures above the $T_{c0}$, and the deconfinement phase is localized in the central region while the confinement phase appears to be in the periphery of the system. Confirming the earlier arguments put forward in Ref. [44], the radial dependence of the local critical temperature (7.2) does not demonstrate the behavior suggested by the Tolman-Ehrenfest law in the rotating gluon matter.

Although in our paper, we studied gluodynamics, one can expect that our qualitative results will not change if the dynamical fermions are accounted for in the simulation. This conclusion can be drawn from the results of paper [37] where it was shown that in rotating QCD, the gluon sector of the theory plays the major role, while the contribution of dynamical fermions is not significant.

The action of rotating gluodynamics is a polynomial of the second power of angular velocity. In view of this observation, the gluon operators that contribute to the linear and quadratic in the $\Omega_I$ terms can be studied separately, i.e., one can switch on/off different operators and study the action in this regime. These operators have clear physical meaning. The linear operator is related to the total angular momentum of gluons. The quadratic one is related to the squares of the chromomagnetic fields. It is also interesting to note that the quadratic term is responsible for the negative values of the moment of inertia of gluon plasma [44–47] below the supervortical temperature $T_s \simeq 1.5T_{c0}$.

Taking into account these interesting properties of the linear and quadratic operators, we have studied their influence on the mixed-phase state and local critical temperature. We found that the linear term generates the inhomogeneous phase transition in rotating gluodynamics, but its properties are opposite to that of the full theory. In particular, for
imaginary rotation, the linear operator gives rise to the deconfinement in the center and the confinement in the periphery region. The local critical temperature increases with $\Omega_I r$. The properties of the quadratic operator are similar to that of the full theory. We also found that the linear term plays a minor role in the dynamics of the full system, and the properties of the mixed-phase state are determined by the quadratic magnetovortical term.

Lattice simulation of rotating gluodynamics without linear term is not spoiled by the sign problem. Therefore, one can conduct the lattice simulations for both real and imaginary angular velocity and study how analytic continuation works. Our results for the critical temperature and the distribution of the Polyakov loop suggest that within a small numerical uncertainty, the simulations with imaginary and real angular velocity are connected through the analytical continuation procedure $\Omega^2 \leftrightarrow -\Omega_I^2$.

In the last section of our work, we proposed and studied the local thermalization hypothesis. The idea is that for sufficiently slow rotation, the coefficients of the operators in the gluon action, which depend on the $\Omega_I r$, remain approximately the same if one varies the distance by a few correlation lengths of gluodynamics. Therefore, instead of the direct simulation of the whole inhomogeneous system, one can study homogeneous gluon action with the anisotropic couplings fixed at some distance. This approximation has valuable advantages over the direct simulation of rotating gluodynamics. In particular, there is no need to perform lattice simulations at huge spatial volumes. In addition, in the homogeneous case, there is no systematic uncertainty caused by the boundary effects.

We have conducted a lattice study of the local thermalization hypothesis and found that the results are in good agreement with that of the direct simulation of rotating gluodynamics. Furthermore, in the local thermalization approximation, it is possible to observe the $\sim (\Omega_I r)^4$ contribution to the local critical temperature, which is not visible in the direct simulations because of the boundary effects. The results for local critical temperature in this approximation are well described by Eq. (6.2) with the continuum limit coefficients (6.4).

Finally, the local thermalization approximation allows us to show that the external gravitational field of the rotating reference frame gives rise to asymmetry in the coupling constants of the gluon action. In turn, this leads to the shift of the critical temperature in the periphery regions as compared to the central regions of the rotating system and, as a result, to the appearance of the mixed inhomogeneous confinement-deconfinement phase.

**Summary of our results.** For the convenience of the reader, we briefly summarize the results of this paper:

1. The critical temperature of the deconfining transition on the rotating axis $T_c(r=0)$ for rotating gluonic matter coincides with the critical temperature of non-rotating gluonic matter $T_{c0}$. In other words, rotation produces no effect on the critical deconfining temperature on the rotating axis, thus affirming our earlier finding of Ref. [41].

2. The local critical temperature $T_c(r)$ in the rotating system is mainly determined by the local velocity $\Omega_I r$ (see Section 4.2). This allowed us to introduce the approximation of local thermalization, which gave consistent results for $T_c(r)$ (see Section 6.2).
3. In agreement with Ref. [41], we confirm the existence of a novel inhomogeneous phase in the rotating Yang-Mills theory. Gluon matter rotating with a real-valued angular velocity resides in the deconfinement phase close to the axis of rotation surrounded by the confinement phase at larger radii. At a given temperature $T$ at the rotation axis, the deconfinement phase extends up to the radius $r$ given by the solution of the equation $T = T_c(r)$.

4. The emergence of this exotic inhomogeneous phase appears to be a result of the quadratic magnetovortical coupling of the angular velocity to the magnetic compo- nent of the gluonic field. The standard linear coupling to the mechanical angular momentum plays a subleading role [cf. Eqs. (2.9) and (5.1)].

5. The gluon matter with the magnetovortical coupling is not affected by the sign prob- lem, so it can be simulated with real-valued angular velocity. We confirm the crucial role of the magnetovortical coupling in the formation of the new inhomogeneous phase. We also demonstrate the validity of the analytical continuation of our results from imaginary to real-valued angular frequencies.

6. We show that the source of the disagreement of our first-principle result with a naive prediction suggested by the Tolman-Ehrenfest law is the result of the magnetovortical coupling of rotation to the magnetic gluonic condensate. This term is produced by the anisotropy of the gluonic action in curved spacetime in a co-rotating non-inertial reference frame of rotating plasma.

## Acknowledgments

The work of VVB, YaAG and AAR has been carried out using computing resources of the Federal collective usage center Complex for Simulation and Data Processing for Mega- science Facilities at NRC "Kurchatov Institute", http://ckp.nrcki.ru/ and the Supercom- puter "Govorun" of Joint Institute for Nuclear Research. The work of VVB and AAR, which consisted of the lattice calculation of the observables used in the paper and in- terpretation of the data, was supported by the Russian Science Foundation (project no. 23-12-00072). The work of MNC was funded by the EU's NextGenerationEU instrument through the National Recovery and Resilience Plan of Romania - Pillar III-C9-I8, man- aged 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. The authors are grateful to Victor Ambruş, Matteo Buzzegoli, Andrey Kotov and Dmitrii Sychev for useful discussions.

## A Width of the transition region

To estimate the spatial width of the transition region, we consider the radial dependence of the local Polyakov loop $L(r)$, defined by Eq. (3.9). In the vicinity of the transition, we
<figure><img src="image_15.png" /><figcaption>Figure 14. The spatial width $2\Delta r$ of the transition region as a function of the velocity of rotation for a fixed temperature $T/T_{c0} = 0.95$ (left), and as a function of temperature for a fixed velocity of rotation $v_I^2 = 0.24$ (right). The results are obtained on the lattices $4 \times 24 \times 145^2$, $5 \times 30 \times 181^2$ and $6 \times 36 \times 217^2$.</figcaption></figure>

fit the data by the following function:

$$L(r) = A_0 + A_1 \tanh \left( \frac{r - r_0}{\Delta r} \right), \qquad (A.1)$$

where $\Delta r$ is a half-width of the spatial transition. Note that here, we do not introduce any additional averaging over cylindrical layers and we do use values of the local Polyakov loop calculated in all spatial points in x, y-plane. The transition width $2\Delta r$ is approximately $\sim 3-5$ fm with a mild dependence on the rotation velocity and the temperature (see Fig. 14).

## B Calculation of the local critical temperature

For a given circular layer of the width $\delta r$ with an average radius $r$, we calculate a local Polyakov loop susceptibility (3.10) for different values of the (on-axis) temperature. The resulting data are fitted by the Gaussian function of a general form

$$\chi_L(T) = c_0 + c_1 \exp \left( -\frac{(T - T_c)^2}{2 c_2^2} \right) + c_3 \tanh \left( \frac{T - T_c}{c_4} \right), \qquad (B.1)$$

what gives us an estimation of the corresponding local critical temperature $T_c$. The role of the last term in Eq. (B.1) is to make the fitting function compatible with different asymptotic values below and above the critical temperature $T_c$.

To estimate the systematic uncertainty, we fit the data using different subsets of data points (at least 12 data points are always used for a 6-parameter fit). As the final result for the local critical temperature, we use the fit coefficient $T_c$ from the subset, which provides the best-fit quality ($\chi^2$/n.d.f. is usually about $\sim 1-3$). An example of this procedure is illustrated in Fig. 15.
<figure><img src="image_16.png" /><figcaption>**Figure 15.** The temperature dependence of the susceptibility of the Polyakov loop for a given circular layer $44 \le r/a < 56$ of width $\delta r \cdot T = 3$ with boundary offset $\delta b \cdot T = 4$. The results of Gaussian fit (B.1) for different subsets of data points are shown by red and orange lines. The best fit with $\chi^2$/n.d.f. = 1.422 is shown by the red line. Orange lines represent the fit results, which are used to estimate systematic uncertainty. The vertical shaded region denotes the uncertainty of the local critical temperature for this circular layer, which includes both statistical and systematic contributions. The data are obtained on the lattice $4 \times 24 \times 145^2$ with OBC.</figcaption></figure>

## C Analysis of the finite-volume and lattice spacing effects

The results for local critical temperature $T_c(r)$, which were presented in Section 4.2, are well described in the whole range of available data by the quartic function

$$ \frac{T_c(r)}{T_{c0}} = C_0 - C_2 \left(\frac{r}{R}\right)^2 + C_4 \left(\frac{r}{R}\right)^4. \qquad (C.1) $$

Note that the first two terms in Eq. (C.1) are not sufficient to describe all data points. However, the quadratic fitting function (i.e., Eq. (C.1) with $C_4 = 0$) can be used if we cut the fit range at $(r/R)^{(\max)} \sim 0.5$. In this case, fit quality becomes appropriate $\chi^2$/n.d.f. $\sim 1-3$. The quadratic fit was also used in Ref. [41].

To study the effects of finite lattice spacing, we check the results for the lattices with $N_t = 4, 5, 6$ and the same aspect ratios for fixed boundary velocity $v_I^2 = 0.16$. The comparison of the fit coefficients for the lattices $4 \times 24 \times 145^2, 5 \times 30 \times 181^2, 6 \times 36 \times 217^2$ with open and periodic boundary conditions are presented in Fig. 16. The uncertainties include both statistical and systematic contributions, as it is described in Appendix D.

From Fig. 16, one can conclude that the results for two higher coefficients $C_0, C_2$ are almost independent of boundary conditions, while the coefficient $C_4$ are slightly differ for PBC and OBC. Note that the results for the coefficients $C_0$ and $C_2$ from the quartic and quadratic fit are in reasonable agreement with each other. The results are mildly dependent on the lattice spacing.

The results for the local critical temperature, presented above, were calculated on the lattices with the transversal aspect ratio $N_s/N_t = 36$. To study the dependence of the
<figure><img src="image_17.png" /><figcaption>**Figure 16.** The coefficients $C_0, C_2, C_4$ as a function of the $N_t$. The results are obtained on the lattices $4 \times 24 \times 145^2$, $5 \times 30 \times 181^2$ and $6 \times 36 \times 217^2$ at the fixed rotating velocity $v_I^2 = 0.16$ for open and periodic boundary conditions.</figcaption></figure>

result on the spatial size of the system, we calculate the local critical temperatures using the lattices with $N_s/N_t = 24, 48$ and the same value of boundary velocity $v_I^2 = 0.16$. A comparison of the fit coefficients for the lattices $4 \times 24 \times 97^2, 4 \times 24 \times 145^2$, and $4 \times 24 \times 193^2$ are presented in Fig. 17.

From Fig. 17, it follows that the coefficient $C_0$, which determines the local critical temperature at the axis of rotation, tends to unity with increasing of the system size $R \propto N_s$. It means that a small deviation of $T_c(0)$ from $T_{c0}$ may be associated with the effects of finite system radius. The second coefficient, $C_2$, is almost independent of $N_s$ and boundary conditions within the uncertainty, while the next coefficient, $C_4$, is different for open and periodic boundary conditions. This coefficient, $C_4$, is larger in the case of open boundary conditions, but it decreases with the growth of the lattice size. As the size of the lattice increases, the results for OBC and PBC approach each other. This result means that the effects of the boundary conditions also become smaller as the system radius increases.

## D Estimation of systematic uncertainty for the fit coefficients

We check how the fit coefficients $C_0, C_2$, and $C_4$ depend on the averaging width, and no significant dependence has been found. A minor dependence for the coefficient $C_0$ is observed (this dependence is more or less profound depending on the lattice size and the rotation velocity). This coefficient becomes closer to unity with the decrease of the layer
<figure><img src="image_18.png" /><figcaption>**Figure 17.** The coefficients $C_0, C_2, C_4$ as a function of transversal lattice size $N_s/N_t$. The results are obtained on the lattices $4 \times 24 \times 97^2$, $4 \times 24 \times 145^2$, $4 \times 24 \times 193^2$ with open and periodic boundary conditions at the fixed rotating velocity $v_I^2 = 0.16$.</figcaption></figure>

width $\delta r$. At the same time, other coefficients are not sensitive to the averaging width. So, we do not take $\delta r \to 0$ limit of the coefficients and consider the choice of the circular layer width $\delta r$ as a source of additional systematic error. We average the fit coefficients $C_i$ overall considered values $\delta r \cdot T = 1, 2, 3, 4, 5$ and use this value as the result for a given lattice and rotational velocity $v_I$.

Note that the neighboring circular layers may partially intersect with each other. To check the stability of the results with respect to the distance between neighboring layers, we fit the data from layers shifted by different steps of 2 and $N_t$ lattice sites along the radius and add these results to the averaging.

Another source of the statistical ambiguity is associated with the value of the boundary skip $\delta b$ (for the quartic fit) or fit range $x = r/R < x_{max}$ (for the quadratic fit). We check how the fit coefficients vary with minor changes in this parameter ($x_{max} = 0.4, 0.5$ or $\delta b = 4, 5$) and also include this variation into the systematic uncertainty.

So, the final results for the fit coefficient include all systematic uncertainties discussed above, as well as a statistical one. An example of this calculation is shown in Fig. 18 for the lattice $4 \times 24 \times 145^2$ with OBC at $v_I^2 = 0.16$.

## References

[1] A.L. Watts et al., *Colloquium: Measuring the neutron star equation of state using x-ray timing*, *Rev. Mod. Phys.* **88** (2016) 021001 [1602.01081].
<figure><img src="image_19.png" /><figcaption>**Figure 18.** The coefficients $C_0, C_2, C_4$ of the best fit of the local critical temperature by the function (4.1), and the fit quality $\chi^2$/n.d.f., obtained with different widths of circular layer $\delta r \cdot T$ (w), boundary skips $\delta b \cdot T$ (s), and distance between the inner radius of used neighboring layers (d). The horizontal lines with hatched uncertainty show the final results. The data are obtained on the lattice $4 \times 24 \times 145^2$ with OBC at $v_I^2 = 0.16$.</figcaption></figure>

[2] I.A. Grenier and A.K. Harding, *Gamma-ray pulsars: a gold mine*, *Comptes Rendus Physique* **16** (2015) 641 [1509.08823].

[3] G. Basar, D.E. Kharzeev and H.-U. Yee, *Triangle anomaly in Weyl semimetals*, *Phys. Rev. B* **89** (2014) 035142 [1305.6338].

[4] K. Landsteiner, *Anomalous transport of Weyl fermions in Weyl semimetals*, *Phys. Rev. B* **89** (2014) 075124 [1306.4932].

[5] M. Arabgol and T. Sleator, *Observation of the nuclear barnett effect*, *Phys. Rev. Lett.* **122** (2019) 177202.

[6] Y. Jiang, Z.-W. Lin and J. Liao, *Rotating quark-gluon plasma in relativistic heavy ion collisions*, *Phys. Rev. C* **94** (2016) 044910 [1602.06580].

[7] M. Baznat, K. Gudima, A. Sorin and O. Teryaev, *Helicity separation in Heavy-Ion Collisions*, *Phys. Rev. C* **88** (2013) 061901 [1301.7003].

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

[9] M.N. Chernodub and S. Gongyo, *Interacting fermions in rotation: chiral symmetry restoration, moment of inertia and thermodynamics*, *JHEP* **01** (2017) 136 [1611.02598].

[10] Y. Jiang and J. Liao, *Pairing Phase Transitions of Matter under Rotation*, *Phys. Rev. Lett.* **117** (2016) 192302 [1606.03808].
[11] 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].

[12] X. Wang, M. Wei, Z. Li and M. Huang, Quark matter under rotation in the NJL model with vector interaction, Phys. Rev. D 99 (2019) 016018 [1808.01931].

[13] X. Chen, L. Zhang, D. Li, D. Hou and M. Huang, Gluodynamics and deconfinement phase transition under rotation from holography, JHEP 07 (2021) 132 [2020.14478].

[14] A.A. Golubtsova and N.S. Tsegelnik, Probing the holographic model of $N=4$ SYM rotating quark-gluon plasma, Phys. Rev. D 107 (2023) 106017 [2211.11722].

[15] S. Chen, K. Fukushima and Y. Shimada, Perturbative Confinement in Thermal Yang-Mills Theories Induced by Imaginary Angular Velocity, Phys. Rev. Lett. 129 (2022) 242002 [2207.12665].

[16] P. Singha, V.E. Ambrus and M.N. Chernodub, Inhibition of the splitting of the chiral and deconfinement transition due to rotation in QCD: The phase diagram of the linear sigma model coupled to Polyakov loops, Phys. Rev. D 110 (2024) 094053 [2407.07828].

[17] Y. Jiang, Rotating SU(2) gluon matter and deconfinement at finite temperature, Phys. Lett. B 853 (2024) 138655 [2312.06166].

[18] 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].

[19] H.-L. Chen, Z.-B. Zhu and X.-G. Huang, Quark-meson model under rotation: A functional renormalization group study, Phys. Rev. D 108 (2023) 054006 [2306.08362].

[20] Y.-Q. Zhao, S. He, D. Hou, L. Li and Z. Li, Phase diagram of holographic thermal dense QCD matter with rotation, JHEP 04 (2023) 115 [2212.14662].

[21] G. Yadav, Deconfinement temperature of rotating QGP at intermediate coupling from M-theory, Phys. Lett. B 841 (2023) 137925 [2203.11959].

[22] N.R.F. Braga, L.F. Faulhaber and O.C. Junqueira, Confinement-deconfinement temperature for a rotating quark-gluon plasma, Phys. Rev. D 105 (2022) 106003 [2201.05581].

[23] S.M.A.T. Mehr and F. Taghinavaz, Chiral phase transition of a dense, magnetized and rotating quark matter, Annals Phys. 454 (2023) 169357 [2201.05398].

[24] N. Sadooghi, S.M.A. Tabatabaee Mehr and F. Taghinavaz, Inverse magnetorotational catalysis and the phase diagram of a rotating hot and magnetized quark matter, Phys. Rev. D 104 (2021) 116022 [2108.12760].

[25] 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].

[26] 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].

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

[28] Y. Jiang, Chiral vortical catalysis, Eur. Phys. J. C 82 (2022) 949 [2108.09622].

[29] K. Mameda and K. Takizawa, Deconfinement transition in the revolving bag model, Phys. Lett. B 847 (2023) 138317 [2308.07310].
[30] Y. Chen, X. Chen, D. Li and M. Huang, *Deconfinement and chiral restoration phase transition under rotation from holography in an anisotropic gravitational background*, Phys. Rev. D **111** (2025) 046006 [2405.06386].

[31] 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 under rotation*, Phys. Rev. D **109** (2024) 116017 [2402.16595].

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

[33] 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*, JETP Lett. **112** (2020) 6.

[34] 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].

[35] V. Braguta, A.Y. Kotov, D. Kuznedelev and A. Roenko, *Lattice study of the confinement/deconfinement transition in rotating gluodynamics*, PoS LATTICE2021 (2022) 125 [2110.12302].

[36] 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 fast imaginary rotation*, Phys. Rev. D **107** (2023) 114502 [2209.15534].

[37] 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].

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

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

[40] N.R.F. Braga and O.C. Junqueira, *Inhomogeneity of a rotating quark-gluon plasma from holography*, Phys. Lett. B **848** (2024) 138330 [2306.08653].

[41] 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].

[42] S. Chen, K. Fukushima and Y. Shimada, *Inhomogeneous confinement and chiral symmetry breaking induced by imaginary angular velocity*, Phys. Lett. B **859** (2024) 139107 [2404.00965].

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

[44] 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].

[45] 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.

[46] V.V. Braguta, M.N. Chernodub, I.E. Kudrov, A.A. Roenko and D.A. Sychev, *Moment of*
inertia and supervortical temperature of gluon plasma, PoS LATTICE2023 (2024) 181 [2311.03947].

[47] 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 chromomagnetic condensate*, Phys. Rev. D 110 (2024) 014511 [2310.16036].

[48] E. Siri and N. Sadooghi, *Thermodynamic properties of a relativistic Bose gas under rigid rotation*, Phys. Rev. D 110 (2024) 036016 [2405.09481].

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

[50] W. Israel, *Nonstationary irreversible thermodynamics: A Causal relativistic theory*, Annals Phys. 100 (1976) 310.

[51] C. Cercignani and G. Kremer, *The relativistic Boltzmann equation: theory and applications*, Birkhäuser Verlag, Basel, Switzerland (2002).

[52] V.E. Ambruş and E. Winstanley, *Rotating fermions inside a cylindrical boundary*, Phys. Rev. D 93 (2016) 104014 [1512.05239].

[53] L.D. Landau and E.M. Lifshitz, *Mechanics*, Butterworth-Heinemann, Oxford, England, 3 ed. (1982).

[54] L.D. Landau and E.M. Lifshitz, *Statistical Physics*, Butterworth-Heinemann, Oxford, England, 3 ed. (1996).

[55] M. Matsuo, J. Ieda and S. Maekawa, *Mechanical generation of spin current*, Frontiers in Physics 3 (2015) 54.

[56] S. Wang, J.-X. Chen, D. Hou and H.-C. Ren, *Strong Coupling Expansion of Gluodynamics on a Lattice under Rotation*, 2505.15487.

[57] L.D. Landau and E.M. Lifshitz, *Statistical Physics, Part 2*, Butterworth-Heinemann, Oxford, England, 3 ed. (1996).

[58] G.K. Savvidy, *Infrared Instability of the Vacuum State of Gauge Theories and Asymptotic Freedom*, Phys. Lett. B 71 (1977) 133.

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

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

[61] G. Curci, P. Menotti and G. Paffuti, *Symanzik's Improved Lagrangian for Lattice Gauge Theory*, Phys. Lett. B 130 (1983) 205.

[62] M. Luscher and P. Weisz, *Computation of the Action for On-Shell Improved Lattice Gauge Theories at Weak Coupling*, Phys. Lett. B 158 (1985) 250.

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

[64] F. Karsch, *SU(N) Gauge Theory Couplings on Asymmetric Lattices*, Nucl. Phys. B 205 (1982) 285.