QCD Equation of State at nonzero baryon density in external magnetic field

This paper is devoted to the study of QCD equation of state in external magnetic field and nonzero baryon density. Our study is carried out by means of lattice simulation with 2+1 dynamical staggered quarks at the physical masses. The simulation is conducted at imaginary baryon chemical potential what allowed us to overcome the sign problem. We expand the pressure in the baryon imaginary chemical potential and study three leading nonzero coefficients in this expansion. These coefficients were calculated for the following values of magnetic field: $eB=0.3$, $0.6$, $1.2$ GeV$^2$ with the lattice sizes $8\times32^3$, $10\times40^3$, $12\times48^3$. Using these data we take continuum limit for the coefficients. Our results indicate considerable enhancement of the expansion coefficients by the magnetic field.


I. INTRODUCTION
Equation of State (EoS) of Quantum Chromodynamics (QCD) plays a fundamental role both from theoretical and practical points of view.From the theoretical perspective EoS contains an important information about thermodynamic properties and QCD phase transitions.On the other hand, from the practical perspective EoS is used for hydrodynamic simulations of heavy-ion collision experiments as well as in different astrophysical applications.
Strong magnetic field might have significant impact on QCD EoS.This can be anticipated if one accounts for dimensional reduction phenomenon realised in magnetic field [34].The idea is that for small magnetic field the quarks live in 3 + 1 dimensions, whereas for strong magnetic field they effectively live in 1 + 1 dimensions.It is reasonable to expect that this phenomenon might lead to strong dependence of EoS on magnetic field.At zero baryon density EoS in external magnetic field was studied on the lattice in [35][36][37].The second-order fluctuations of the baryon number, electric charge and strangeness, which are related to EoS, in external magnetic field were studied in papers [38,39].It is worth to note that perturbatively QCD EoS at nonzero density and magnetic field was studied in Refs.[40,41].
In this paper we are going to study EoS of QCD in external homogeneous magnetic field and nonzero baryon chemical potential.Our study is carried out within lattice simulation with dynamical staggered u-, d-, s-quarks at their physical masses.To avoid the sign problem the simulation is performed at imaginary chemical potentials and at the following setup µ u = µ d = µ, µ s = 0 which approximately corresponds to zero strangeness.We expand the free energy into a series in imaginary chemical potential and restrict our study to three nonzero terms in this expansion.Mainly we focus on the coefficients of this expansion and their dependence on magnetic field.This paper in organized as follows.In the next Section we derive the formulas important for our study.The lattice setup used in the simulation is presented in Section 3. In Section 4 the results of our study are shown.In last Section we discuss the results and draw conclusions.

II. EQUATION OF STATE AT FINITE BARYON DENSITY IN EXTERNAL MAGNETIC FIELD
For a thermodynamic ensemble with a known partition function Z, the pressure p may be found using the general formula where µ u , µ d , µ s are the chemical potentials for u-, dand s-quarks correspondingly.These chemical potentials can be rewritten in terms of the chemical potentials µ B , µ Q , µ S which correspond to the baryon charge B, electric charge Q and strangeness S.
In this paper we restrict our consideration to the following relation between quark chemical potentials: µ u = µ d = µ, µ s = 0, which approximately corresponds to zero strangeness ⟨S⟩ ≃ 0. In this case Eqs.(2) give µ B = 3µ, µ Q = 0, µ S = µ.More accurate tuning of the chemical potential values relevant for modern heavy-ion experiments can be found in Refs.[10,39].For a sufficiently small baryon density the EoS (1) can be expanded in a series in µ B /T : The coefficient c 0 determines QCD EoS at zero baryon density and external magnetic field.This case was studied in papers [35][36][37] and we are not going to consider it.On contrary in this paper we are going to focus on the coefficients c 2 , c 4 , c 6 and their dependence on the magnitude of external magnetic field.The coefficients c 2 , c 4 , c 6 are related to the fluctuations χ BQS ijk of the conserved charges B, Q, S at zero baryon density .
(4) It causes no difficulties to find the formulas which connect the χ BQS ijk to the c 2 , c 4 , c 6 coefficients Lattice simulations at real values of the chemical potential are hampered by the sign problem [42], so we perform the simulations at the imaginary values of the chemical potential, characterized by θ = µ I /T = iµ B /T and expand our results in θ.Using the data we can conduct an analytical continuation to the real values of µ as long as we do not encounter any discontinuity.Since the thermal QCD phase transitions at baryon density and not too large magnetic field are expected to be a crossover, we believe that the analytical continuation procedure is justified in our case.
The partition function itself and the pressure cannot be measured in the lattice simulations.One can only measure pressure derivatives with respect to external parameters.The basic quantity, which we measure on the lattice is the (imaginary) baryon density n I , which is defined as the derivative of the pressure with respect to θ = iµ B /T : By fitting the dependence of the imaginary baryon density n I on the values of the (imaginary) chemical potential θ we can extract the values of the coefficients c 2 , c 4 and c 6 .

III. LATTICE SETUP
We perform simulations in QCD with N f = 2 + 1 staggered fermions at the physical pion mass in the presence of an external magnetic field B and quark chemical potentials µ f .The partition function Z of the system under study has the following form: where det D f [U x,µ , q f B, µ f ] corresponds to the staggered Dirac operator of the quark flavour f .Note that various quark flavours have different quark charges q f and values of the chemical potential µ f , thus we treat all three flavours separately.Each staggered Dirac operator det D f [U x,µ , q f B, µ f ] corresponds to 4 quark tastes and we take the fourth root to have one quark taste.Expression for the D f [U x,µ , q f B, µ f ] has the following form: where am f is the quark mass, η µ = (−1) x1+...+xµ−1 are the standard factors for the staggered fermions, and U x,µ are two-times stout smeared gauge links U x,µ with isotropic smearing parameter ρ = 0.15 [43].The factors u f x,µ ∈ U (1) correspond to the magnetic field B pointing in the z-direction and they are given by the following expressions (N s is the spatial lattice size): On the lattice with the periodic boundary conditions the magnetic flux is quantized.Since different quarks have different quarks charges, we take the minimum absolute value for the quark charge q = |q d | = |q s | = e/3 and the quantization condition for the magnetic field value takes the form It should be noted that the condition n ≪ N 2 s is required in order to avoid lattice artifacts.In our simulations this condition is also fulfilled, moreover the comparison of results obtained at different values of N t suggest rather small lattice artifacts almost for all points used in the simulations.Another restriction comes from the existence of the Roberge-Weiss phase transition [44] in the plane of complex chemical potentials, µ I /T ≤ π, which is also always fulfilled in our simulations.
For the gluon sector we use the tree-level Symanzik improved gauge action: where β = 6/g 2 corresponds to the bare lattice gauge coupling and W n×k x,µν stands for the trace of the rectangular with the size n × k constructed from the gauge links U x,µ , starting from the point x in directions µ and ν.
Bare parameters have been set so as to stay on a line of constant physics [4,5,45], at the isospin symmetric point, m u = m d = m l , a physical strange-to-light mass ratio, m s /m l = 28.15, and a physical pseudo-Goldstone pion mass, m π ≃ 135 MeV.
We perform simulations with three lattice sizes 8×32 3 , 10×40 3 and 12×48 3 , keeping the aspect ratio N s /N t = 4. Previous simulations [46] suggest that such aspect ratio has sufficiently small finite volume effects.By using three different lattice spacings we can study the continuum extrapolation of the studied quantities.Typically, the results for different N t are close to each other and we perform continuum extrapolation using simple quadratic ansatz for the dependence on the lattice spacing a: O(a) = O(0) + Aa 2 .To assess the systematic uncertainty, we do the same procedure keeping only two fine lattice sizes N t = 10 and N t = 12 and take half of the sum as our final result and half of the difference as the systematic uncertainty.

IV. RESULTS OF THE CALCULATION AND DISCUSSION
In Fig. 1 we present the dependence of the density n I on the (imaginary) chemical potential µ I for all three studied values of magnetic field eB = 0.3, 0.6 and 1.2 GeV 2 and for various temperatures in the vicinity of the phase transition.The data shown in Fig. 1 were obtained on the lattice 10 × 40 3 .The density n I for various µ I , eB and T was also calculated on the lattices 8 × 32 3 , 12 × 48 3 .
Fitting our lattice data for the n I by formula (6) we obtained the coefficients c 2 , c 4 , c 6 .In Fig. 2 we plot these coefficients as function of temperature for all values of magnetic field and lattices considered in our study.In addition in Fig. 2 we show continuum limit for the coefficients c 2 , c 4 , c 6 .The numerical values of these coefficients in continuum limit are listed in the Appendix , and the bands in Fig. 2 represent the interpolation of these data by the cubic splines [47].
In Fig. 3 we show the continuum extrapolation for the coefficients c 2 , c 4 , c 6 , but contrary to Fig. 2 we draw each coefficient for all values of magnetic field on one plot.In addition in Fig. 3 we present the results from Ref [38] obtained at the following values of magnetic field eB = 0.314, 0.627, 1.255 GeV 2 and at zero magnetic field, determined from [12].In these works the authors calculated the fluctuations of the conserved charges, from which one can easily reconstruct the coefficients c 2 , c 4 and c 6 using formulae (5).In particular, the results of [38] allow to determine the c 2 only, while using the results of [12] one can find all three coefficients c 2 , c 4 , c 6 considered in this paper.The results of [38] were obtained at slightly larger pion mass m π ≈ 220 MeV and magnetic fields eB = 0.314, 0.627, 1.255 GeV 2 , which are very close to the values used in this work.We see good agreement with our results for the c 2 at high temperatures and some difference between them at T ≤ 140 MeV.We believe, that this discrepancy can be attributed to the higher pion mass in [38].It might also be attributed to the fact that the authors of [38] did not take continuum limit.Now few comments are in order.First let us con-sider continuum extrapolated c 2 coefficient.It seen in Fig. 3 that the c 2 coefficient, i.e. fluctuations, gradually increases with magnetic field.At the largest value of magnetic field eB = 1.2 GeV 2 the c 2 is by an order of magnitude larger than that at zero magnetic field.Moreover magnetic field changes the dependence of the c 2 on temperature.Thus at zero magnetic field the c 2 coefficient is monotonously rising function of temperature, while for large magnetic fields considered in this paper the c 2 clearly develops a peak structure.For larger magnetic fields the peak is more pronounced.In addition the position of the peak shifts to lower temperatures with increasing of the value of eB.We believe that this peak is related to the QCD critical temperature and the shift of the peak position to the left is related to inverse magnetic catalysis phenomenon [20].
Further let us pay attention to the c 4 and c 6 coefficients.They are known to exhibit nontrivial behaviour around the phase transition/crossover.For example, if one looks at the plot for c 4 (middle panel of Fig. 3), then one observes a peak around the transition temperature.At the same time, the dependence of c 6 (lower panel of Fig. 3) is even more complicated, it has positive peak before the transition and a negative dip after the Continuum extrapolated coefficients c2, c4, c6 as functions of temperature.The dashed lines represent the ideal gas approximation for the same coefficients.Various magnetic fields considered in this work are marked by colors as follows eB = 0.3 GeV 2 (green), 0.6 GeV 2 (blue), 1.2 GeV 2 (red).In addition we show the results for the c2 coefficient [38] obtained at the magnetic fields eB = 1.255GeV 2 (red diamonds), eB = 0.627 GeV 2 (blue squares), eB = 0.314 GeV 2 (green circles) and the coefficients c2, c4, c6 at zero magnetic field (brown bands) determined from [12].phase transition.It can clearly seen, that corresponding peaks and dips are shifted to the left with increasing magnetic field, which is in agreement with inverse magnetic catalysis phenomenon.Moreover, their magnitude significantly grows with the magnetic field, while their width decreases.
Using our results for coefficients c 2 , c 4 and c 6 we calculated additional contribution to the pressure coming from nonzero baryon density ∆p = p(µ B ) − p(µ B = 0).In order to perform it we carried out analytical continuation from imaginary to real values of the baryon chemical potential using formula (3).We would like to note here that the problem of analytical continuation is not well defined and becomes unstable at large values of baryon chemical potential.For this reason we consider values of chemical potential µ B /πT ≤ 0.9, within this region the results from analytical continuation are stable against various functional forms for zero magnetic field [10].In Fig. 4 we plot the pressure excess ∆p as a function of temperature for various values of real baryon chemical potential µ B and all magnetic fields considered in this work.To assess the stability of our results, we determined the pressure excess ∆p either using all three coefficients c 2 , c 4 and c 6 or using only two lower coefficients c 2 and c 4 (shown in Fig. 4 by different hatching).It can be clearly seen that the results obtained in two different ways are almost the same, only for the largest studied value of chemical potential µ B /πT = 0.90 there is some difference.It confirms that for the studied values of µ B the analytical continuation is under control.
From the Fig. 4 one can conclude that magnetic field not only enhances the pressure ∆p but also modifies its dependence on temperature and chemical potential.
It is believed that at high temperatures the quarkgluon plasma reveal properties of ideal gas of quarks and gluons.In the ideal gas approximation magnetic field acts only on quarks.The pressure for massless quarks in magnetic field can be written in the following form [24] where we have used the designation μf = µ f /T .Expression (12) allows us to find the coefficients c 2 , c 4 , c 6 for an ideal gas which have the following form It should be noted here that the first term in the coefficient c 2 (13) arises from the lowest Landau level, whereas the second term results from the excited Landau levels.3), as a function of temperature for the values of (real) baryon chemical potential µB/πT = 0.15, 0.30, 0.45, 0.60, 0.75, 0.90 for all magnetic fields used in our study.Different directions of diagonal hatching correspond to the results obtained using all three coefficients c2, c4, c6 or using only two lower coefficients c2, c4 in Eq. ( 3).The dashed lines represent the ideal gas limit.
What concerns the coefficients c 4 , c 6 (14), ( 15) they acquire contribution from the excited Landau levels only.The values of magnetic fields studied in this paper obey the inequality eB > T 2 .Moreover, for the largest magnetic field it takes the form of eB ≫ T 2 .So, in the ideal gas approximation the contribution of the exited Landau levels are suppressed leading to the suppression of the coefficients c 4 and c 6 .However, radiative corrections to these coefficients might change this picture.
In Fig. 2 and Fig. 3 we also show the coefficients c 2 , c 4 , c 6 in the ideal gas approximation.First let us consider the c 2 coefficient.Our continuum limit for the c 2 and ideal gas approximation prediction differ by 20-30% in the region T > 160 MeV for the smallest magnetic field.However, as one increases the value of magnetic field this discrepancy decreases and at the largest magnetic field lattice results are in agreement with the ideal gas approximation.Taking into account the above discussion this result implies that the lowest Landau level plays more and more important role as one increases the magnitude of magnetic field in the temperature range T ≥ 160 MeV.Probably for the largest magnetic field eB = 1.2 GeV 2 the contribution of the lowest Landau level becomes dominant.What concerns the coefficients c 4 , c 6 lattice results are in agreement with the ideal gas approximation in the temperature interval T > 160 MeV.Notice, however, that the uncertainties of our lattice results are quite large in this case.
In Fig. 4 we have shown the excess of pressure ∆p in the ideal gas approximation by the dashed lines.Similarly to Fig. 2 and Fig. 3 lattice results for the ∆p are in agreement with the ideal gas approximation in the temperature interval T > 160 MeV.

V. CONCLUSION
In this paper we conducted the study of QCD equation of state in external magnetic field and nonzero baryon density.Our study is carried out within lattice simulation with 2+1 dynamical staggered quarks at their physical masses.Nonzero baryon density was introduced to the system through baryon chemical potential.To avoid the sign problem the simulation is performed at imaginary chemical potential.We expand the pressure into a series in imaginary chemical potential and focus in our study on the first three nonzero coefficients in this expansion.These coefficients were calculated on the lattices 8 × 32 3 , 10 × 40 3 , 12 × 48 3 for the following values of magnetic field: eB = 0.3, 0.6, 1.2 GeV 2 , and the continuum limit was taken.
Our data indicate that magnetic field modifies the temperature dependence of the coefficients.For instance, large magnetic field gives rise to the peak of the c 2 coefficients which was absent at zero magnetic field.Notice also that for large magnetic fields positions of peaks and dips in the temperature dependence of coefficients are shifted to lower temperatures.It is reasonable to assume that nontrivial behaviour of the coefficients takes place in the crossover region and the shift of this region to lower temperatures is related to the inverse magnetic catalysis phenomenon.
In addition to modifying the temperature dependence, magnetic field enhances the magnitude of the coefficients c 2 , c 4 , c 6 considerably.The coefficients are known to be connected to the fluctuations of the conserved charges, i.e. magnetic field enhances these fluctuations.The origin of this enhancement is not completely clear.It is might be related to the dimensional reduction phenomenon.Indeed, we observe considerable enhancement of the c 2 coefficient in an ideal fermion gas approximation, but this approximation is not sufficient to explain the behaviour of the c 4 and c 6 coefficients, possibly indicating that the enhancement of the fluctuations in the c 4 and c 6 coefficients is of nonperturbative nature.
The authors are grateful to Natalia Kolomoyets for the help at the initial stage of the project.This work 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 Supercomputer "Govorun" of Joint Institute for Nuclear Research.This work was supported by the Russian Science Foundation (project no.23-12-00072).
Appendix: Data tables The continuum limit values of the coefficients c 2 , c 4 , c 6 in the pressure expansion (3) at fixed magnetic fields eB = 0.3, 0.6, 1.2 GeV 2 are listed in Tables I, II, III.Data are presented with full uncertainty, which includes both statistical and systematic errors.

5 − 8 −
FIG. 1.The density nI /T 3 as a function of the imaginary chemical potential µI /T = θ calculated on the lattice 10×40 3 for three values of the magnetic field eB = 0.3, 0.6, 1.2 GeV 2 and several temperatures in the vicinity of the thermal phase transition.

6 FIG. 2 .
FIG.2.The coefficients c2, c4, c6 as functions of temperature for all studied values of magnetic field calculated on the lattices: 8 × 32 3 , 10 × 40 3 , 12 × 48 3 .The orange bands correspond to the continuum extrapolation of the coefficients, their widths reflect full uncertainties, combined from statistical and systematic.The dashed line represents the coefficients in the ideal gas approximation.

TABLE I .
Continuum limit results for the coefficients of the series, Eq. (3), at eB = 0.3 GeV 2 .

TABLE II .
Continuum limit results for the coefficients of the series, Eq. (3), at eB = 0.6 GeV 2 .

TABLE III .
Continuum limit results for the coefficients of the series, Eq. (3), at eB = 1.2 GeV 2 .