INTRODUCTION

Heavy-ion collision experiments make it possible to study the effect of extreme conditions such as high temperatures, strong electromagnetic fields, a high baryon density, and fast rotation on the properties of a strongly interacting quark–gluon plasma. The effect of fast rotation on the quark–gluon plasma is of particular interest. The average vorticity of the formed quark–gluon plasma is estimated from experimental results for the polarization of the \(\Lambda ,\bar {\Lambda }\) hyperons as ω = (9 ± 1) × 1021 s–1 [1]. Such a high angular velocity corresponds to relativistic rotation, which can strongly affect the properties of the quark–gluon plasma.

The effect of relativistic rotation on the properties of QCD is studied in numerous theoretical works. Various models such as the Nambu–Jona-Lasinio model [29] and the hadron resonance gas model [10], as well as the holographic approach to QCD [1114] and other methods [1517], were used. Although interesting results were obtained in these studies, analytical methods of studying QCD have serious disadvantages because QCD is a very complex theory, which cannot be analytically studied without additional assumptions. Additional assumptions lead to uncontrolled systematic errors in the results. Consequently, it is difficult to estimate the reliability of the predictions obtained.

The above problem is absent for the lattice simulation of QCD. This method makes it possible to study the properties of QCD based on the fundamental principles of quantum field theory, controlling the errors of such calculations. We use this method in this work to study gluodynamics. The first study of the properties of rotating QCD was reported in [18]. Later, we performed lattice simulations of the effect of rotation on the thermodynamic properties of gluodynamics [1921] and QCD [22, 23].

The aim of this work is the lattice simulation study of the effect of rotation on the properties of QCD. In particular, we study the effect of rotation on the equation of state of gluodynamics by calculating the first nonzero correction in angular velocity to the free energy of the studied system. It is noteworthy that the equation of state of gluodynamics and QCD, as well as the effect of various extreme conditions on it, was actively studied in lattice simulations (see, e.g., [2428]). The equation of state of rotating gluodynamics was investigated in [29]. In this work, we carry out a similar study by using another computational method.

LATTICE SIMULATION OF ROTATING GLUODYNAMICS

In this section, we briefly describe the lattice simulation study of rotating gluodynamics and derive formulas necessary for calculations. This approach was described in more detail in [19, 20].

We calculate the partition function of the equilibrium system by the Monte Carlo method in the reference frame rotating with the studied “medium” about the z axis. In this reference frame, rotation is manifested in the form of an external gravitational field, specified by the known metric tensor

$${{g}_{{\mu \nu }}} = \left( {\begin{array}{*{20}{c}} {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{array}} \right),$$
(1)

where \(r = \sqrt {{{x}^{2}} + {{y}^{2}}} \) is the distance to the rotation axis.

The partition function of gluodynamics in the external gravitational field can be represented in the continuous space in the form of the integral over the gluon degrees of freedom [1820]:

$$Z = \int DA{\kern 1pt} \exp ( - {{S}_{G}}),$$
(2)

where

$${{S}_{G}} = \frac{1}{{2g_{{YM}}^{2}}}\int {{d}^{4}}x{\kern 1pt} \sqrt {{{g}_{{\text{E}}}}} {\kern 1pt} g_{{\text{E}}}^{{\mu \nu }}g_{{\text{E}}}^{{\alpha \beta }}F_{{\mu \alpha }}^{a}F_{{\nu \beta }}^{a}{\kern 1pt} $$
(3)

is the Euclidian action of the gluon field in the external gravitational field, depending on the Euclidean metric tensor (gE)μν, which can be obtained from Eq. (1) by means of the Wick rotation \(t \to i\tau \). As in the partition function without the gravitational field, the Euclidean time τ varies in the range \(\tau \in (0,\beta )\), and periodic boundary conditions in the Euclidean time \(A_{\mu }^{a}(0,\mathbf{x}) = A_{\mu }^{a}(\beta ,\mathbf{x})\) are imposed on gluon fields. Greek and Latin scripts are used in this work for Lorentz and color indices, respectively.

The substitution of the metric tensor (gE)μν into Eq. (3) gives the action in the form

$${{S}_{G}} = \frac{1}{{2g_{{YM}}^{2}}}\int {{d}^{4}}x{\text{Tr}}{\kern 1pt} [(1 - {{r}^{2}}{{\Omega }^{2}})F_{{xy}}^{a}F_{{xy}}^{a}$$
$$ + \;(1 - {{y}^{2}}{{\Omega }^{2}})F_{{xz}}^{a}F_{{xz}}^{a} + (1 - {{x}^{2}}{{\Omega }^{2}})F_{{yz}}^{a}F_{{yz}}^{a}$$
$$ + \;F_{{x\tau }}^{a}F_{{x\tau }}^{a} + F_{{y\tau }}^{a}F_{{y\tau }}^{a} + F_{{z\tau }}^{a}F_{{z\tau }}^{a}$$
(4)
$$ - \;2iy\Omega (F_{{xy}}^{a}F_{{y\tau }}^{a} + F_{{xz}}^{a}F_{{z\tau }}^{a})$$
$$ + \;2ix\Omega (F_{{yx}}^{a}F_{{x\tau }}^{a} + F_{{yz}}^{a}F_{{z\tau }}^{a}) - 2xy{{\Omega }^{2}}{{F}_{{xz}}}{{F}_{{zy}}}].$$

According to this formula, the action is a complex quantity, which leads to the sign problem. Unfortunately, the direct Monte Carlo simulation of such systems is currently impossible. This problem can be overcome by using a method based on the simulation of the system at an imaginary angular velocity [18, 2029]. This method was used to study the equation of state of rotating gluodynamics in [29]. In this work, we use another approach. It is shown below that the free energy can be expanded in a series in angular velocity. The coefficients of this expansion are operators of gluon fields, which can be calculated by simulating the system without rotation. Thus, the sign problem does not arise in this approach.

The action (4) with the imaginary angular velocity is discretized as in [1820]. We do not present the explicit expression for the lattice action in this work because it is lengthy. As mentioned above, the lattice simulation of the nonrotating system is sufficient to calculate corrections to the free energy in angular velocity. To this end, we used the tree-level improved Symanzik action [30, 31].

Boundary conditions are of particular importance in the simulation of rotating systems. We performed calculations with periodic boundary conditions in the directions τ and z. Periodic, open, or Dirichlet boundary conditions can be imposed in the directions perpendicular to the axis of rotation. Various boundary conditions and their effect on observables were examined in detail in [20], where it was showed that the dependence of various observables on boundary conditions is insignificant because boundary conditions are screened and do not affect the bulk of the studied system. In this work, we use periodic boundary conditions, which are usually used in lattice simulations of gauge field theories (gluodynamics) without r-otation.

The main aim of this work is to study the effect of rotation on the equation of state, which can be determined if the free energy of the studied system F is known. For this reason, we analyze the effect of rotation on the free energy of gluodynamics. We expand the free energy in a series in angular velocity:

$$F = {{F}_{0}} - \frac{{{{\Omega }^{2}}}}{2}{{F}_{2}} + O({{\Omega }^{4}}),$$
(5)

where \({{F}_{0}}\) is the free energy of gluodynamics without rotation and \({{F}_{2}}\) is the next coefficient of the expansion. The free energy of gluodynamics \({{F}_{0}}\) is well studied (see, e.g., [24]), and the coefficient \({{F}_{2}}\) will be calculated below. In this work, we examine only the second coefficient of the expansion because statistical errors of calculations increase rapidly with the size of the system in view of the structure of the studied operators. In addition to \({{F}_{0}}\) and \({{F}_{2}}\), we use the specific quantities \({{f}_{0}} = {{F}_{0}}{\text{/}}V\) and \({{f}_{2}} = {{F}_{2}}{\text{/}}V\), where V is the volume of the studied system.

To determine \({{F}_{2}}\), we represent the action (4) in the form

$${{S}_{G}} = {{S}_{0}} + \Omega {{S}_{1}} + \frac{{{{\Omega }^{2}}}}{2}{{S}_{2}},$$
(6)

where \({{S}_{0}}\) is the action of the gluon fields without rotation and

$$\begin{gathered} {{S}_{1}} = \frac{i}{{g_{{YM}}^{2}}}\int {{d}^{4}}x[xF_{{yx}}^{a}F_{{x\tau }}^{a} + xF_{{yz}}^{a}F_{{z\tau }}^{a} \\ - \;yF_{{xy}}^{a}F_{{y\tau }}^{a} - yF_{{xz}}^{a}F_{{z\tau }}^{a}], \\ \end{gathered} $$
(7)
$$\begin{gathered} {{S}_{2}} = - \frac{1}{{g_{{YM}}^{2}}}\int {{d}^{4}}x[{{r}^{2}}{{(F_{{xy}}^{a})}^{2}} + {{y}^{2}}{{(F_{{xz}}^{a})}^{2}} \\ + \;{{x}^{2}}{{(F_{{yz}}^{a})}^{2}} + 2xyF_{{xz}}^{a}F_{{zy}}^{a}]. \\ \end{gathered} $$
(8)

Using these formulas, the coefficient \({{F}_{2}}\) can be represented in the form

$${{F}_{2}} = T{{\left. {\frac{{{{\partial }^{2}}\log Z}}{{\partial {{\Omega }^{2}}}}} \right|}_{{\Omega = 0}}} = T\left( {{{{\langle \langle {{S}_{2}}\rangle \rangle }}_{T}} - {{{\langle \langle S_{1}^{2}\rangle \rangle }}_{T}}} \right){\kern 1pt} ,$$
(9)

where \({{\langle \langle \mathcal{O}\rangle \rangle }_{T}} = \langle \mathcal{O}\rangle - {{\langle \mathcal{O}\rangle }_{{T = 0}}}\) is the contribution of thermal fluctuations to the mean value of the operator \(\mathcal{O}\), the first term \(T{{\langle \langle {{S}_{2}}\rangle \rangle }_{T}}\) corresponds to the average square of the chromomagnetic field in the studied medium, and the second term \( - T{{\langle \langle S_{1}^{2}\rangle \rangle }_{T}}\) reflects fluctuations of the angular momentum of the gluon field (it is taken into account that \({{\left. {{{{\langle \langle {{S}_{1}}\rangle \rangle }}_{T}}} \right|}_{{\Omega = 0}}} = 0\)).

To calculate the coefficient \({{F}_{2}}\) by Eq. (9), it is sufficient to perform the lattice simulation of gluodynamics without rotation. Consequently, the sign problem is absent in our calculations, as mentioned above. We studied the dependence of the free energy on the angular velocity by another method in [29], where the calculation was carried out by the simulation of gluodynamics rotating at an imaginary angular velocity. The derivative of the partition function with respect to the inverse coupling constant \(1{\text{/}}g_{{YM}}^{2}\) was calculated to extract the free energy. This derivative is related to the mean value of the action of the gluon field \({{\langle \langle {{S}_{G}}\rangle \rangle }_{T}}\) (4). The free energy was calculated by integrating its derivative over the inverse coupling constant.

When studying the equation of state, the corresponding operator at \(T = 0\) is usually subtracted. Similarly, Eq. (9) for the coefficient \({{F}_{2}}\) at a finite temperature is evaluated with subtraction at zero temperature. This procedure physically means that the vacuum of gluodynamics at \(T = 0\) does not rotate.

Formulas (6) and (9) allow us to obtain the expansion of the free energy in angular velocity for continuous theory. Lattice versions of these formulas can be easily derived. To do this, we note that the lattice expression of the action in the rotating reference frame can be represented in the form of expansion (6). The explicit form of the lattice operators \({{S}_{1}}\) and \({{S}_{2}}\) can be found in [20]. We do not present here these lengthy expressions. The lattice formula (9) for the coefficient \({{F}_{2}}\) has the same form, where the lattice expressions are substituted for the operators \({{S}_{1}}\) and \({{S}_{2}}\).

The coefficient \({{F}_{2}}\) is related to the moment of inertia of the system given by the formula

$$I = \frac{M}{\Omega } = - \frac{1}{\Omega }\frac{{\partial F}}{{\partial \Omega }} = {{F}_{2}},$$
(10)

where M is the angular momentum of the system. Further, we use the nonrelativistic expression for the moment of inertia to evaluate the dependence of the coefficient \({{F}_{2}}\) on the volume of the studied system:

$$\begin{gathered} I = \int dV{{r}^{2}}\rho (T,x,\Omega ) = \frac{1}{6}{{\rho }_{0}}(T){{L}_{z}}L_{s}^{4} \\ = \frac{1}{6}{{\rho }_{0}}(T)VL_{s}^{2}, \\ \end{gathered} $$
(11)

where \(\rho (T,x,\Omega )\) and \(V = {{L}_{z}}L_{s}^{2}\) are the mass density and the volume of the studied system, respectively, which is an \({{L}_{s}} \times {{L}_{s}} \times {{L}_{z}}\) parallelepiped. The axis of rotation passes through the center of the parallelepiped and coincides with the z axis. It is taken into account in the second equality in Eq. (11) that the density in the leading approximation is independent of the angular velocity and coordinates; i.e., it is a function \({{\rho }_{0}}(T)\) of only the temperature.

Expression (11) makes it possible to represent expansion (5) in the following form, which is more convenient for the analysis:

$$F = {{F}_{0}}\left( {1 + \frac{1}{2}{{K}_{2}}{{v}^{2}}} \right),$$
(12)

where \({{K}_{2}} = - 4{{F}_{2}}{\text{/}}({{F}_{0}}L_{s}^{2})\) and \(v = \Omega {{L}_{s}}{\text{/}}2\) is the velocity at the center of the side face at a distance of \({{L}_{s}}{\text{/}}2\) from the axis of rotation. The coefficient \({{K}_{2}}\) is used in Eq. (12) instead of the coefficient \({{F}_{2}}\) because it is dimensionless and should be independent of the volume for a sufficiently large system, as follows from Eq. (11). This statement will be verified in the next section. The coefficients \({{F}_{0}}\) and \({{F}_{2}}\) depend on boundary conditions, but our calculations show that this dependence is canceled in the ratio \({{F}_{2}}{\text{/}}{{F}_{0}}\) within the errors.

RESULTS AND DISCUSSION

In this section, we describe the calculation of the coefficient \({{K}_{2}}\) in lattice simulation on \({{N}_{t}} \times {{N}_{z}} \times N_{s}^{2}\) (\({{N}_{s}} = {{N}_{x}} = {{N}_{y}}\)) lattices. For the extrapolation to the c-ontinuous limit \(1{\text{/}}{{N}_{t}} \to 0\) (lattice spacing \(a \to 0\)), we carried out the simulation on \(4 \times 16 \times {{21}^{2}}\), \(5 \times 20 \times {{25}^{2}}\), and \(6 \times 24 \times {{31}^{2}}\) lattices. The lattice sizes were chosen such that the ratios \({{N}_{s}}{\text{/}}{{N}_{t}} = {{L}_{s}}T\) and \({{N}_{z}}{\text{/}}{{N}_{t}} = {{L}_{z}}T\), where \({{L}_{z}} = {{N}_{z}}a\), \({{L}_{s}} = {{N}_{s}}a\), and \(1{\text{/}}T = {{N}_{t}}a\), approximately remain constant with the variation of \({{N}_{t}}\). The size of the lattice in the x and y directions was chosen such that the axis of rotation, when calculating the operators (7) and (8), passes through a site of the lattice and the corresponding lattice operators are symmetric with respect to this axis. As mentioned above, to determine the coefficient \({{F}_{2}}\) at a finite temperature, this coefficient is also calculated at zero temperature. For this subtraction procedure, we used the lattices with the same spatial dimensions but with \({{N}_{t}} = {{N}_{z}}\).

Figure 1 presents the temperature dependence of the ratio \({{f}_{2}}{\text{/}}({{T}^{4}}L_{s}^{2})\) for lattices with different number of discretization steps \({{N}_{t}}\) in the time direction. The temperature in Figs. 1–3 is given in units of the critical temperature Tc of the confinement/deconfinement phase transition in gluodynamics without rotation.

Fig. 1.
figure 1

(Color online) Ratio \({{f}_{2}}{\text{/}}({{T}^{4}}L_{s}^{2})\) for lattices with va-rious time steps \({{N}_{t}}\) versus the temperature in units of the critical temperature Tc of the confinement/deconfinement phase transition in gluodynamics without rotation.

Fig. 2.
figure 2

(Color online) Ratio \({{f}_{2}}{\text{/}}({{T}^{4}}L_{s}^{2})\) for lattices with different volumes in the x (y) and z directions versus the temperature in units of the critical temperature Tc of the confinement/deconfinement phase transition in gluodynamics without rotation.

Fig. 3.
figure 3

(Color online) Coefficient \({{K}_{2}}\) for lattices with various time steps \({{N}_{t}}\) versus the temperature in units of the critical temperature Tc of the confinement/deconfinement phase transition in gluodynamics without rotation. Stars indicate the continuous limit \(1{\text{/}}{{N}_{t}} \to 0\), and the results of the continuous limit are approximated by a spline.

To study the dependence of our results on the size of the system, we calculated the ratio \({{f}_{2}}{\text{/}}({{T}^{4}}L_{s}^{2})\) on \(5 \times 16 \times {{25}^{2}}\), \(5 \times 24 \times {{25}^{2}}\), \(5 \times 20 \times {{21}^{2}}\), and \(5 \times 20 \times \) 312 lattices. Figure 2 shows the temperature dependence of the ratio \({{f}_{2}}{\text{/}}({{T}^{4}}L_{s}^{2})\) for lattices with different spatial volumes in the x (y) and z directions. It is seen that the dependence of the ratio \({{f}_{2}}{\text{/}}({{T}^{4}}L_{s}^{2})\) on the sizes \({{L}_{s}}\) and \({{L}_{z}}\) of the system is weak, which confirms the estimate made in Eq. (11).

To determine the coefficient \({{K}_{2}}\) on the \(4 \times 16 \times {{21}^{2}}\), \(5 \times 20 \times {{25}^{2}}\), and \(6 \times 24 \times {{31}^{2}}\) lattices, we calculated \({{F}_{0}}\) and \({{F}_{2}}\), then, the ratio \({{K}_{2}} = - 4{{F}_{2}}{\text{/}}({{F}_{0}}L_{s}^{2})\) was taken and extrapolated to the continuous limit \(1{\text{/}}{{N}_{t}} \to 0\). The temperature dependence of the coefficient \({{K}_{2}}\) for lattices with different \({{N}_{t}}\) values is presented in Fig. 3, where stars indicate the continuous limit \(1{\text{/}}{{N}_{t}} \to 0\), and the results in the continuous limit are approximated by a spline. Only the temperature range T > Tc is shown in Fig. 3 because \({{F}_{0}}\) at T < Tc is close to zero and the calculation errors increase significantly. For this reason, the plot for the coefficient \({{K}_{2}}\) below the critical temperature becomes noninformative. Figure 3 is the main result of this work.

In this work, the lattice simulation of gluodynamics without rotation was carried out and the first nonzero correction to the free energy in the angular velocity was directly calculated. At the same time, we calculated the coefficient \({{K}_{2}}\) by a different method in [29], where rotating gluodynamics was simulated, the free energy was calculated for different angular velocities, the calculations were completed by the expansion in (imaginary) angular velocity, and the results were analytically continued to real angular velocities \(\Omega \). The results obtained by both methods are in agreement with each other.

As seen in Fig. 3, the coefficient \({{K}_{2}}\) is negative for temperatures T < T * and becomes positive at T > T *, where T * ~ 1.5Tc. According to Eq. (10), this coefficient is related to the moment of inertia of the studied system as \(I = - {{K}_{2}}{{F}_{0}}L_{s}^{2}{\text{/}}4\). Since \({{f}_{0}} = - p < 0\), Fig. 3 indicates that the moment of inertia of gluodynamics is negative up to the temperature T  * ~ 1.5Tc and becomes positive only above this temperature. We suppose that the negative moment of inertia means the thermodynamic instability of the gluon plasma with respect to uniform rotation.Footnote 1 The uniform rotation of the plasma in thermodynamical equilibrium is possible at temperatures above T  *. The thermodynamically stable motion of the gluon plasma at temperatures below T  * is more complex, e.g., rotation at the angular velocity depending on the coordinates. However, additional studies are required to understand the physical mechanisms of this instability.

To reveal the reason for the negative moment of inertia, we consider expression (9) for \({{F}_{2}}\), which is the sum of two operators. The first operator \( \sim {\kern 1pt} {{\langle \langle {{S}_{2}}\rangle \rangle }_{T}}\) is related to the square of the chromomagnetic field in gluodynamics, which appears in an important physical quantity, the gluon condensate \( \sim {\kern 1pt} {{\langle \langle {{(G_{{\mu \nu }}^{a})}^{2}}\rangle \rangle }_{T}} \sim \) \({{\langle \langle {{({\mathbf{E}^{a}})}^{2}}\rangle \rangle }_{T}} + {{\langle \langle {{({\mathbf{H}^{a}})}^{2}}\rangle \rangle }_{T}}\). This operator determines the energy of the vacuum in gluodynamics at \(T = 0\) and scale invariance breaking. With increasing temperature, the magnetic component of the gluon condensate decreases from its value at \(T = 0\), reaches a certain minimum, and begins to increase [24]. The contribution from thermal fluctuations to the mean value of this operator, which appears in Eq. (9), is proportional to \({{\langle \langle {{({\mathbf{H}^{a}})}^{2}}\rangle \rangle }_{T}}\), is negative at 0 < T ≲ 2Tc, and becomes positive only at temperatures above T ~ 2Tc [24]. The second operator \( \sim {\kern 1pt} {{\langle \langle S_{1}^{2}\rangle \rangle }_{T}}\), which is related to fluctuations of the angular momentum, makes a positive contribution to the moment of inertia, thus slightly increasing it. Thus, the negative moment of inertia of gluodynamics at temperatures \(T < T{\text{*}}\) is explained by the contribution from the magnetic component of the scale anomaly.

As follows from Eq. (12) and from the results for \({{K}_{2}}(T)\), rotation reduces the free energy of gluodynamics near the critical temperature of the confinement/deconfinement phase transition. Therefore, it can be expected that this phase transition at \(\Omega \ne 0\) requires an additional “heating” of the system. Furthermore, the higher the angular velocity, the larger the required heating of the system. Consequently, the critical temperature Tc of the confinement/deconfinement phase transition should increase with the angular velocity of rotation. Thus, the result obtained in this work is in qualitative agreement with the results obtained previously in [19, 20].