Dense and magnetized QCD from imaginary chemical potential

It is widely believed that strongly interacting matter existed as a Quark-Gluon Plasma (QGP) during the first microseconds of the Universe, where quarks and gluons were subjected to high temperatures and densities. During that epoch, the so-called electroweak phase transition may have generated magnetic fields comparable with the typical energy scale of the strong interactions [1]. In the present Universe, this strongly interacting phase of matter is believed to exist inside the core of neutron stars, some of which also show strong magnetic fields. The study of such extreme environments in the lab has been pioneered by experiments of Heavy-Ion Collisions (HICs). These experiments also account for the strongest magnetic fields ever produced on Earth. In peripheral collisions, for instance, the motion of charged particles generates field lines that reinforce each other at the center of the collision. The sizeable impact of such fields on experimentally measurable observables has been recently seen by the STAR Collaboration [2]. In phenomenological studies, theoretical estimates indicate that magnetic fields reach $\sqrt{eB}\sim 0.1\text{ GeV}$ – for typical Au+Au collisions – at RHIC energies [3] and $\sqrt{eB}\sim 0.5\text{ GeV}$ – for typical Pb+Pb collisions – at LHC energies [4].

Besides strong magnetic fields, high densities can also be created in HICs. High-density strongly interacting matter produced in relativistic collisions is relevant not only for ongoing experiments, such as RHIC, but also for future ones, like NICA and FAIR. From a theoretical viewpoint, high densities also provide a path to understanding the different phases of the strong interactions via its underlying microscopic theory: Quantum Chromodynamics (QCD).

At vanishing density, the QCD phase transition is a smooth crossover [5], whereas at high densities it is conjectured to become first order. The presence of a first-order line would imply the existence of a second-order transition point, known as a Critical Endpoint (CEP). This hypothesis has been corroborated by analytical approaches, such as the Functional Renormalization Group [6], Dyson-Schwinger equations [7], holographic models [8], etc.

On the one hand, first-principles numerical simulations of the strong interactions on a lattice – Lattice QCD (LQCD)– successfully account for all non-perturbative features of QCD at vanishing and low densities. On the other hand, the approach becomes unfeasible at moderate and high densities due to the sign (or complex action) problem. Although various techniques have been used to circumvent this issue (see Ref. [9] for a recent review), to this day, the determination of the phase structure of QCD at high $\mu_{B}$ from first principles remains a widely open problem.

Theoretically, one can extend the QCD phase diagram to the $\mu_{B}^{2}<0$ (imaginary) domain, where lattice simulations are sign-problem-free. Interestingly, at imaginary $\mu_{B}$, the phase diagram has a second-order transition point called Roberge-Weiss CEP [10]. Also, the QCD phase diagram in the $T$-$B$ plane is directly accessible within lattice simulations. In the presence of a magnetic field, the QCD crossover becomes sharper and the transition temperature is lowered, suggesting that the crossover would eventually become a first-order transition starting at a CEP. Moreover, lattice simulations of asymptotically strong magnetic fields found further indication for this scenario [11], thus supporting previous expectations based on generic arguments [12] (see also Ref. [13]). Recent LQCD results constrained this CEP to the region: $4<eB<9\text{ GeV}^{2}$, $65<T<95\text{ MeV}$ [14]. The impact of magnetic fields on the Roberge-Weiss CEP has also been studied using LQCD, where the findings indicate that this point approaches the $T$-$B$ plane, thus opening the possibility of a connection between the two critical points [15].

Another way of probing thermodynamic features of QCD is via the Equation of State (EoS). Moreover, knowing the EoS as a function of thermodynamic variables ($T$, $\mu_{B}$, $B$, etc.) is essential to model the hydrodynamic evolution of relativistic collisions across the phase diagram. In this regard, the present knowledge of the interplay between finite densities and magnetic fields is limited. Currently, the EoS at $B=0$ is known from LQCD in a broad range of temperatures at vanishing $\mu_{B}$ with 2+1 flavors [16, 17], as well as at non-vanishing $\mu_{B}$ with 2+1 flavors [18], and 2+1+1 flavors [19]. At $B\neq 0$, the EoS has been calculated with 2+1 flavors at non-zero $\mu_{B}$ but vanishing strangeness chemical potential [20]. For a recent review on lattice investigations of the phase diagram and the EoS at nonzero magnetic fields, see Ref. [21].

Hence, in this proceedings article, we aim to shed light on the EoS of dense and magnetized QCD matter with four dynamical quark flavors. Our method of choice to circumvent the sign problem is analytic continuation from simulations at imaginary $\mu_{B}$, thus extending our previous work [22], where we carried out simulations at vanishing $\mu_{B}$. Our EoS is taken along a trajectory satisfying two constraints motivated by the HIC setup, namely, strangeness neutrality and isospin asymmetry. Our discussions revolve around two main quantities: the chemical potentials $\mu_{Q}$ and $\mu_{S}$ required by our constraints, and the leading-order behavior of the EoS. We show that the magnetic field significantly impacts the latter, in particular, near the crossover temperature. In Fig. 1, we illustrate the conjectured QCD phase diagram in the $T$-$\mu_{B}^{2}$-$B$ space and indicate where we carried out simulations.

This article is organized as follows: In Sec. 2, we describe the simulation setup, including our lattice action and simulation parameters. In Sec. 3, we introduce the definitions of the observables we compute as well as the constraints that we impose on our EoS. This is followed by Sec. 4, in which we explain our method to fulfill the constraints at non-zero chemical potential. Finally, in Secs. 5 and 6, we present our results and conclusions, respectively.

In the fermion sector, we employed $2+1+1$ flavors of rooted-staggered fermions with four stout-smearing steps and masses tuned to the physical point. For details on our Line of Constant Physics (LCP), see Ref. [23]. In the gauge sector, we used the tree-level Symanzik-improved gauge action [24]. In this setup, we generated gauge configurations on a $32^{3}\times 8$ lattice for various temperatures, imaginary chemical potentials, and magnetic field strengths. We varied the temperature via the inverse gauge coupling $\beta$ according to $T=[a(\beta)N_{t}]^{-1}$, where $a$ is the lattice spacing and $N_{t}$ the number of sites in the temporal extent, and the imaginary chemical potential via $\mu_{B}=i\pi j/8$, with $j\in\mathbb{Z}$. We scanned the temperature in the range $T=135$–$200$ MeV in steps of $5$ MeV, and the following chemical potentials: $j=0,3,4,5$. Finally, without loss of generality, we considered a uniform magnetic field pointing in the $z$ direction. To simulate the magnetic field, we multiplied the gluon links by the corresponding U(1) factors at each lattice site. For the exact form of these factors, see [25].

In a finite volume, the magnetic field flux satisfies the following quantization condition

$$eB=\frac{6\pi N_{b}}{L_{x}L_{y}}\,,\hskip 14.22636ptN_{b}\in\mathbb{Z}\,,$$

(1)

where $L_{x}$ and $L_{y}$ are the lattice extents along the $x$ and $y$ directions, respectively. Due to QED-related charge- and magnetic-field renormalizations, the amplitude of the field is typically expressed in terms of the renormalization-group-invariant quantity $eB$. In this fashion, we simulated $eB=0.3\,,0.5\,,0.8$ GeV² magnetic fields, which encompass the values expected in the physical systems described in Sec. 1. To achieve the desired values of $eB$ – given that the change of flux can only be carried out in integer steps $N_{b}$ – for each of the aforementioned field strengths, we simulated two values of $N_{b}$ that surround the desired value. Thus, we obtained results at the desired $B$ using a linear interpolation.

Here, we give the definitions of the observables that we compute. Let us start by introducing the pressure in terms of the grand canonical partition function of QCD in the rooted-staggered formalism

$$\hat{P}\equiv\frac{P}{T^{4}}=\frac{1}{VT^{3}}\ln\mathcal{Z}\,,\hskip 28.45274pt\mathcal{Z}(\{\hat{\mu}_{i}\},T,B)=\int\mathcal{D}Ue^{-S_{g}}\prod_{f}\det[\not{D}(\hat{\mu}_{f},B)+m_{f}]^{1/4}\,,$$

(2)

where $V$ is the spatial volume, $\{\mu_{i}\}=\mu_{u}\,,\mu_{d}\,,\mu_{s}$, and $\hat{\mu}_{f}\equiv\mu_{f}/T$ are the reduced chemical potentials in the quark basis. Notice that we do not include the charm chemical potential although the charm quark is dynamically included in our simulations. Moreover, we define the quark-number susceptibilities as

$\displaystyle\chi^{uds}_{ijk}(\{\hat{\mu}_{i}\},T,B)=\frac{1}{VT^{3}}\frac{\partial^{i+j+k}}{\partial\hat{\mu}_{u}^{i}\partial\hat{\mu}_{d}^{j}\partial\hat{\mu}_{s}^{k}}\ln\mathcal{Z}\,.$

(3)

These susceptibilities can be converted to the physical basis – involving baryon, charge, and strangeness quantum numbers – yielding, for instance, the following number densities

$$n_{B}=\frac{1}{3}\chi^{uds}_{100}+\frac{1}{3}\chi^{uds}_{010}+\frac{1}{3}\chi^{uds}_{001}\,,\hskip 14.22636ptn_{Q}=\frac{2}{3}\chi^{uds}_{100}-\frac{1}{3}\chi^{uds}_{010}-\frac{1}{3}\chi^{uds}_{001}\,,\hskip 14.22636ptn_{S}=-\chi^{uds}_{001}\,.$$

(4)

Using these quantities, we can now introduce the constraints that we impose on the EoS, namely,

$\displaystyle n_{S}=0\,,\hskip 28.45274pt\frac{n_{Q}}{n_{B}}=0.4\,.$

(5)

These are known respectively as strangeness neutrality and isospin asymmetry. The ratio $n_{Q}/n_{B}$ and the value of $n_{S}$ are chosen to match the experimental values expected for the byproducts of typical Pb+Pb collisions, for instance. These constraints imply that there is only one independent chemical potential, since $\mu_{Q}$ and $\mu_{S}$ can be expressed in terms of $\mu_{B}$.

Taking these constraints into account, we can compute the pressure change at finite chemical potential along the line of strangeness neutrality and isospin asymmetry using the integral method

$\displaystyle\Delta\hat{P}(\hat{\mu}_{B})=\hat{P}(\hat{\mu}_{B})-\hat{P}(\hat{\mu}_{B}=0)=\int_{0}^{\hat{\mu}_{B}}\differential\hat{\mu}_{B}^{\prime}\derivative{\hat{P}}{\hat{\mu}_{B}^{\prime}}\,,$

(6)

where $\hat{P}$ was introduced in (2). In terms of the number densities (4), the total derivative of the pressure with respect to $\hat{\mu}_{B}$ can be written as

$$\derivative{\hat{P}}{\hat{\mu}_{B}}=\quantity(n_{B}+n_{Q}\partialderivative{\hat{\mu}_{Q}}{\hat{\mu}_{B}}+n_{S}\partialderivative{\hat{\mu}_{S}}{\hat{\mu}_{B}})\,,$$

(7)

where the partial derivatives of $\hat{\mu}_{Q}$ and $\hat{\mu}_{S}$ with respect to $\hat{\mu}_{B}$ appear due to the constraints. Knowing the right-hand side of Eq. (7) allows us to compute the pressure shift $\Delta\hat{P}(\hat{\mu}_{B})$ via the integral method. Next, we discuss how to tune the charge and strangeness chemical potentials in terms of $\hat{\mu}_{B}$ such that strangeness neutrality and isospin asymmetry are fulfilled.

This section focuses on realizing the experimental conditions of HICs on the lattice. One way is to precisely tune the simulation parameters. Another way is to compute additional derivatives and extrapolate to strangeness neutrality afterwards. To avoid large computational costs during tuning while still simulating close to the strangeness neutral case to avoid a complicated extrapolation, we combine both. To simulate close to strangeness neutrality at $\hat{\mu}_{B}(j)$ we estimate $\hat{\mu}_{Q}(\hat{\mu}_{B}(j))$ and $\hat{\mu}_{S}(\hat{\mu}_{B}(j))$ from all simulations up to $\hat{\mu}_{B}(j-1)$. The shift into exact strangeness neutrality, which is a small difference in $\hat{\mu}_{Q}$ and $\hat{\mu}_{S}$. The determination of the shift is performed after the simulation. We determine the new simulation parameters at $\hat{\mu}_{B}(j)$ by performing extrapolations in $\hat{\mu}_{Q/S}$ using Taylor-expansions up to a certain order, depending on the number of simulation in $\hat{\mu}_{B}$ direction already performed. This means we make an ansatz such as Eq. (8), where $a$ are the parameters we want to determine.

$\displaystyle\hat{\mu}_{Q/S}(\hat{\mu}_{B})=\sum_{i=1}^{N}a_{2i+1}\hat{\mu}_{B}^{2i+1}$

(8)

To estimate systematic effects we calculate the parameters $a$ using different approaches. Combining simulations with estimated parameters with a correction to strangeness neutrality afterwards has the advantage that we keep the computational cost for both the simulation and analysis in balance with low statistical uncertainties. The small deviation of the estimated parameters from the true ones in strangeness neutrality ensures that by performing a linear extrapolation we can calculate the shift $\Delta\hat{\mu}_{Q}$ and $\Delta\hat{\mu}_{S}$ in the chemical potentials to reach the experimental conditions while keeping the impact on the uncertainties in our observables as small as possible. The susceptibilities in strangeness neutrality after the shift, denoted as $\tilde{\chi}^{BQS}_{ijk}$ read as written in Eq. (9).

$\displaystyle\tilde{\chi}^{BQS}_{ijk}=\chi^{BQS}_{ijk}+\Delta\hat{\mu}_{Q}\chi^{BQS}_{i(j+1)k}+\Delta\hat{\mu}_{S}\chi^{BQS}_{ij(k+1)}$

(9)

We will now discuss the different approaches to estimating the chemical potentials for the simulations

In this proceedings article, we investigated the interplay between finite density and nonzero magnetic fields in QCD in the light of observables that can be related to the EoS. We performed simulations at imaginary chemical potential with 2+1+1 dynamical quarks and a physical pion mass. We used different expansion schemes to constrain our simulation parameters to a strangeness-neutral and isospin-asymmetric line, where we defined $n_{Q}/n_{B}=0.4$ as the isospin asymmetry ratio. These constraints are motivated by phenomenological expectations in HICs.

Using a multidimensional spline fit, we took the first steps towards the analytic continuation of our data to the real-chemical-potential axis. At this stage, our results show that the magnetic field has a significant impact on the observables, especially around the crossover temperature. Specifically, we observed that the quantity $\hat{\mu}_{B}^{-1}\differential\hat{P}/\differential\hat{\mu}_{B}$ ceased to be monotonic with temperature in the presence of strong enough magnetic fields, namely, for $eB\gtrsim 0.5$ GeV². This behavior has previously been observed in second-order susceptibilities and in the leading-order coefficient of the EoS [22]. See also Ref. [33] for further discussion about the impact of $B$ on these observables.

In order to make contact with continuum QCD physics, the continuum limit of the results presented here must still be carried out, as well as the full analytic continuation of our data. Nevertheless, our current findings reveal interesting features of the impact of strong magnetic fields on hot and dense QCD matter which might be insightful for the modeling of systems where such medium is created: HICs.

This work is supported by the MKW NRW under the funding code NW21-024-A. The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V. for funding this project by providing computing time on the GCS Supercomputer SuperMUC-NG at Leibniz Supercomputing Centre. The authors are also grateful to Gergely Markó for discussions that shaped our approach to the results presented here. Moreover, we acknowledge Arpith Kumar and Jin-Biao Gu for helpful discussions on dense and magnetized QCD during this conference. Gergely Endrődi acknowledges funding by the Hungarian National Research, Development and Innovation Office (Research Grant Hungary 150241) and the European Research Council (Consolidator Grant 101125637 CoStaMM).