Investigating the effects of QCD matter’s electrical conductivity
on charge dependent directed flow

At high temperatures, hadronic matter smoothly transitions into a strongly interacting quark-gluon plasma (QGP). Relativistic heavy-ion collision experiments at the Relativistic Heavy Ion Collider (RHIC) and the Large Hadron Collider (LHC) have produced and continue to investigate that deconfined hadronic state of QGP. Experimental evidence favors a strongly coupled fluid behavior for QGP, which leads to phenomenological models of the QGP expansion using relativistic viscous hydrodynamics [1].

Simultaneously, the relativistic charges of the collision spectators create intense electromagnetic (EM) fields in and around the collision region. Estimates using Liénard-Wiechert potentials place the initial magnetic field for RHIC energies of $\sqrt{s}=200$ GeV at $e|\vec{B}|\sim\mathcal{O}(m^{2}_{\pi})$ and LHC energies of $\sqrt{s}=2760$ GeV at $e|\vec{B}|\sim\mathcal{O}(10m^{2}_{\pi})$ [2, 3]. Because of the collision geometry, those intense EM fields penetrate the created QGP, which modifies the life-time of the EM fields. Lattice QCD, perturbative QCD (pQCD), and AdS/CFT calculations predict QGP to be a resistive medium [4, 5, 6, 7, 8, 9]. Using that information, the interplay between currents in QGP and the EM fields can be parametrized by the electric conductivity $\sigma$ transport coefficient.

A few collision observables have also been suggested to investigate the electric conductivity of QGP, including using electroweak probes [10, 11]. An example relevant to this paper is the charge dependent directed flow in asymmetric collisions [12]. Theoretically, in asymmetric collisions, an electric field is directed from the larger ion to the smaller ion because of the difference in the number of protons. Then the electric field induces a charged current inside QGP via the Coulomb effect that modifies the directed flow of final state electrically charged hadrons. The argument is that the difference between the directed flow of electrically charged hadrons $\Delta v_{1}\equiv v_{1}^{+}-v_{1}^{-}$ is sensitive to the electric conductivity of QGP. At RHIC, the STAR detector has reported a difference in the directed flow between negatively and positively charged hadrons for Cu-Au collisions at $\sqrt{s_{\text{NN}}}=200$ GeV [13]. However, the uncertainty of the result is too large and the understanding of other non-EM effects is not complete enough for any definite statements to be made about QGP’s conductivity.

Instead, we will consider the charge dependent directed flow $\Delta v_{1}$ of symmetric collisions. Charge dependent directed flow was suggested to be a signal of electromagnetic fields in symmetric collisions [14, 15]. Recent STAR data for the directed flow of protons and kaons supports that theory [16], but results for pions and D mesons indicate non-EM effects are important [17, 18]. Directed flow $v_{1}$ in symmetric collisions has also been studied using relativistic magneto-hydrodynamics [19, 20] and a resistive model was used to study the $\Delta v_{1}$ for pions [21]. The result of the resistive model illustrated the $\Delta v_{1}$ dependence on the electric conductivity of QGP.

For this work, we have improved our relativistic resistive magneto-hydrodynamic (RRMHD) model [22] to study different final state hadrons across collision centralities¹¹1Previously only peripheral centralities and final state pions were possible.. We aim to show the qualitative role of QGP’s electric conductivity in $\Delta v_{1}$. This contrasts with other works that have focused on baryon charges and currents to explain the STAR data [16]. Those works have suggested the positive slopes of $\Delta v_{1}$ can be solely explained by the transported quark effect. But, as was pointed out in Ref. [15], the initial protons in the collision mean the plasma has a net positive charge. We will demonstrate with our RRMHD model that the initial positive charge can also explain positive slopes for $\Delta v_{1}$. This demonstration should be thought as a first step because the initial non-equilibrium dynamics of QGP have not been settled, hence an exact quantitative role of the initial positive charges is difficult to establish.

Physically, an imbalance in $\Delta v_{1}\neq 0$ is thought to be a caused by the Faraday induction + Coulomb effects and the Hall effect acting on charged particles inside QGP. If the combined Faraday induction + Coulomb effects are larger than the Hall effect, $\Delta v_{1}$ will feature a negative slope over the mid-rapidity region [15].

In the next section, Sec. II, we introduce our relativistic resistive magneto-hydrodynamic (RRMHD) model. Then in Sec. III we will show the time evolution of the magnetic field in our model. Next, how we fix the initial conditions of the EM fields as solutions of Maxwell’s equations in a medium and the initial conditions of QGP are discussed in Sec. IV. Finally, in Sec. VI we present the results for charge dependent directed flow using our RRMHD model. Our conclusions are summarized in Sec. VII and we discuss the necessary future improvements to move this type of study from qualitative to quantitative.

For this work, we have built upon the relativistic resistive magneto-hydrodynamic (RRMHD) model of [21, 23, 22] such that we can study central collisions. The overall algorithms of the model remain the same, but some implementations have been changed to avoid sources of numerical instability; unrealistic negative pressures and imaginary Lorentz factors. We have also reduced the computation time by changing the numerical implementation to use the parallel computation architecture called Message Passing Interface (MPI). Because of the implementation changes, we have rebenchmarked the model using the same tests as described in the original paper and found only small differences. That model simultaneously solves the relativistic ideal-hydrodynamic equations and Maxwell’s equations with Ohm’s law as a current source. Reviewing the relevant details, we start with the conservation laws of ideal-hydrodynamics:

Here $\nabla_{\mu}$ is the covariant derivative. Then those ideal-hydrodynamic equations are augmented with Maxwell’s equations,

where $\epsilon^{\mu\nu\rho\sigma}=\tilde{\epsilon}^{\mu\nu\rho\sigma}/\sqrt{-g}$ is the Levi-Civita tensor and tilde epsilon is the antisymmetric symbol. Then to close the system of equations we use the first order covariant form of Ohm’s law,

We have assumed the electric conductivity $\sigma$ is a constant and $q=-J^{\mu}u_{\mu}$ is the electric charge density in the fluid’s comoving frame.

The feature of RRMHD compared to ideal relativistic magneto-hydrodynamics (MHD) is a finite conductivity $\sigma=1/\eta$, which is the inverse of $\eta$ the resistivity. In ideal relativistic MHD models of QGP, the conductivity is taken to be infinite, in other words, QGP is treated as a perfect conductor. However, Lattice QCD and pQCD calculations suggest the conductivity should be $\sigma=\mathcal{O}(0.1\sim 10)$ for the highest energy regions produced at RHIC and the LHC [7]. Additionally, a finite conductivity reduces the lifetime of the EM fields compared to idea MHD and can introduce new phenomena like magnetic reconnection. Modeling QGP with a finite conductivity leads to interesting physics beyond the ideal MHD scenario.

The magnetic field plays a major role in the charge-dependent directed flow $\Delta v_{1}$ at symmetric collisions. One role the magnetic field plays is to induce Faraday and Hall effect currents in the QGP. Then the interplay of Lorentz + Coulomb + Faraday currents within the plasma leads to an electric chemical potential. During the time evolution, the chemical potential is frozen into the freezeout hypersurface. Thus, an understanding of the time evolution of the magnetic field is important when calculating an observable like $\Delta v_{1}$ which depends on freezeout hypersurface.

Numerous calculations have been done to estimate the magnitude and lifetime of the magnetic field in relativistic heavy-ion collisions under simplifying assumptions [15, 14, 20, 19, 24, 25]. In our model, we have three assumptions:

• •

  We assume the electric conductivity of QGP is a constant value over all spacetime. Strictly speaking, dimensional analysis alone suggests the electric conductivity is temperature-dependent, which means it should change as the plasma expands and cools.

• •

  We ignore the pre-equilibrium dynamics that can create an initial electric charge distribution and change the configuration of the initial EM fields [26, 27].

• •

  Last, the initial EM fields are created by the spectators of the collision, but during the time evolution we only include the plasma charges as a source of the EM fields. In other words, our model focuses on the behavior of charges inside the plasma and how they alone affect the time evolution of the EM fields.

For a more complete description, we plan to resolve the first and third items in a future version of our model. The second item depends on an understanding of the pre-equilibrium dynamics, which has recently been an active area of research (e.g.,  [28, 29]). Some possible directions include Color Class Condensate (CGC), IP-Glasma, flux tube, and kinetic models.

With these assumptions in mind, our model captures how a finite electric conductivity modifies the time evolution of the magnetic field. For an intuition of the modifications we derive the magnetic wave equation in Milne coordinates, starting from Faraday’s and Ampere’s laws:

In the usual way, we take the $\partial_{0}$ of Faraday’s Law and substitute it into Ampere’s law to find a wave equation for the magnetic field components,

where $\Box_{M}$ is the d’Alembert operator in Milne coordinates defined to be,

Those are the same equations that have been discussed in the literature [30]. As introduced in Sec. II, we augment Faraday’s and Ampere’s laws Eq.(6) and Eq.(7) with Eq.(5) Ohm’s law. Because our goal is a qualitative feel for the magnetic field’s time evolution in Milne coordinates, we simplify Ohm’s law by neglecting any particle-density gradients and considering the fluid’s rest frame. Then the diffusion current becomes $j^{i}=\sigma E^{i}$ and the magnetic field wave equations are,

Since the sigma dependent terms on the right-hand side of those equations are Faraday’s law, we rewrite them to arrive at:

Equation (17) is for a non-linearly damped magnetic wave in the $\eta$ direction. Equations (15) and (16) are similar to Eq.(17) except for a driving term proportional to the electric field. Even though those three equations are non-linearly damped, the damping is simple enough to gain a qualitative feeling of the magnetic field time evolution.

• •

  For, $\tau^{-1}+\sigma>1$ the magnetic field is over damped and exponentially decays.

• •

  At, $\tau^{-1}+\sigma=1$ the magnetic field is critically damped and will reach an asymptotic value as quickly as possible.

• •

  Then for, $\tau^{-1}+\sigma<1$ the magnetic field is under damped and will oscillate about an asymptotic value.

Because of the non-linearity of the damping, it is possible for the magnetic field to transition from under damped to over damped if the conductivity is between $0<\sigma<1$. As an example, for, $\sigma=0.0294\text{ fm}^{-1}$ the transition occurs close to $\tau=1\text{ fm/c}$. For conductivities larger than 1 the transition would occur at negative $\tau$, which cannot be realized. Thus, the magnetic field is always over damped when $\sigma>1\text{ fm}^{-1}$.

For relativistic heavy-ion collisions, the symmetry of the system is best captured with Milne coordinates $(\tau,x,y,\eta)$, which have the metric $g_{\mu\nu}=\text{diag}(-1,1,1,\tau^{2})$. Accordingly, our RRMHD model solves Eqs.(1)-(5) using those coordinates with the volume factor $\sqrt{-g}=\tau$. We have adopted the optical Glauber model for the initial energy density profile and included an initial longitudinal tilt [31]. Even though we initialize the model after QGP is created at $\tau_{0}=0.4$ fm/$c$, the source of the initial EM fields is dominated by the collisions spectators [24]. Therefore, we calculate the initial EM fields using the analytic solution from [25], which is a solution of Maxwell’s equations in a medium with the ion charges as the source. These are the same initial conditions as the previous works [23, 21, 20] and are summarized in table 1.

The inelastic cross-section $\sigma_{\text{NN}}$ is calculated with the parametrization $\sigma_{\text{NN}}(s)=28.84+0.0458\ln^{2.374}[s]$ from Ref. [32].

Numerically, we use a constant grid size of $168\times 168\times 72$ with $-16.8\text{ fm}<(x,y)<16.8\text{ fm}$ and $-7.2<\eta<7.2$, except for calculating the initial EM fields. Instead, for the initial EM fields, we adopt the following approach based on the analytic solution [25].

1. 1.

   Assume the nuclei are very thin in the beam direction $\eta$ because of Lorentz contraction, such that their charge is distributed only in the transverse plane $x,y$. Then we can fix the $\eta$ position of the charges based on the collision kinematics.

2. 2.

   At each grid point $\begin{pmatrix}x,&y,&\eta\end{pmatrix}$, calculate the EM fields by integrating over the contribution from all possible sources. Practically, this means we create a new grid for the EM sources in the transverse plane of size $200\times 200$ for $-15.0\text{ fm}<(x,y)<15.0\text{ fm}$ and integrate over it with the trapezoidal method. This size ensures the area of both nuclei are integrated over.

The second step is repeated for all the original grid points $\begin{pmatrix}x,&y,&\eta\end{pmatrix}$. The numerical Courant-Friedrichs-Lewy (CFL) value is fixed to be $C_{\text{CFL}}=0.05$.

Lastly, we calculate the final state hadron distribution using the Cooper-Frye formula [33]. We are ignoring the final state interactions that would be included in an afterburner, because our current goal is to focus on the effects of the QGP conductivity. Our implementation of the Cooper-Frye formula starts by constructing a 4-D hypersurface at the energy density boundary $\epsilon_{f}=0.15\text{ GeV/fm}^{3}$ [34]. We then construct a distribution using the fluid velocity $v_{i}$, and charge density $q$ that is perpendicular to the hypersurface $d\Sigma(t,x,y,z)$. Next, the accumulated charge density $q$ is converted to the electric chemical potential using the linear approximation $\mu_{Q}\approx q(\hbar c)^{2.5}/T_{f}^{2}$. Where the final temperature is calculated using the ideal gas equation of state,

Then we integrate over the entire hypersurface to arrive at a distribution. The details of the implementation are given in appendix B.

As discussed at the beginning of Sec. III, the time evolution of the magnetic field induces Lorentz + Faraday currents in QGP. Coulomb currents caused by the electric field of the positive spectators also appear, but for symmetric collisions largely cancel each other. Then the EM field component in the charge dependent directed flow $\Delta v_{1}$ for symmetric collisions is dominated by the Lorentz + Faraday currents, specifically the magnetic field. In particular, the $B_{y}$ component of the field has the largest contribution to $\Delta v_{1}$, which approximately appears as $\gamma B_{y}$ in the Faraday current and $-\gamma u_{z}B_{y}$ in the Lorentz current [15]. Where $\gamma$ is the Lorentz factor of QGP in the local fluid rest frame and $u_{z}$ is the $z$ component of the fluid velocity.

In our RRMHD model, the time evolution of $B_{y}$ depends on the initial conditions and the electric conductivity $\sigma$ of QGP. As a starting point, we have used smooth initial solutions to Maxwell’s equations in a medium as discussed in Sec. IV that predicts a large initial field intensity for $B_{y}$. That being said, solving the time evolution of the magnetohydrodynamics Eqs. (1)-(5) with different $\sigma$’s has an impact on $\Delta v_{1}$ through changes of $B_{y}$, which will be shown in Sec. VI.

Figure 1 illustrates the impact of $\sigma$ on the time evolution of $B_{y}$ for our RRMHD model of QGP + EM fields. Following the setup in Sec. IV we prepare a peripheral Au-Au collision, with an impact parameter of $b=10$ fm, at the collision energy of $\sqrt{s_{\text{NN}}}=200$ GeV. The RRMHD simulation was run from $\tau_{0}=0.4$ fm/$c$ until the energy density for all grid cells below the freezeout value $\epsilon_{f}=0.15\text{ GeV/fm}^{3}$. Focusing on the center of the grid $(x=0\text{ fm},y=0\text{ fm},\eta=0)$, in Fig. 1 we find a different time evolution for $\tau B_{y}$ that depends on the electric conductivity. Larger values of the electric conductivity result in a longer lifetime of the magnetic field. Because we have assumed the initial conditions are the same, the point at the earliest time in Fig. 1 is equal. Additionally, from equations Eq.(15), Eq.(16), and Eq.(17) we expect the magnetic field to behave as a non-linear oscillator. In particular, values of $\sigma<1\text{ fm}^{-1}$ result in oscillations around zero following under-damped behavior. Whereas, those of $\sigma>1\text{ fm}^{-1}$ never cross zero with $\tau B_{y}$ acting over-damped.

The reduction in differences between the results for larger sigma values, like between $\sigma=10\text{ fm}^{-1}$ and $\sigma=100\text{ fm}^{-1}$, is partly due to numerical diffusion overwhelming the physical diffusion. This is a long-standing issue for resistive (relativistic) MHD models in astrophysics, where the conductivities can become relativity large, see Ref. [35] as an example of a detailed analysis. Here for each simulation, we keep the typical length-scale $L$ of the problem fixed and the order of the fluid velocities $u$ remains largely unchanged, but the conductivity increases by an order of magnitude. This means the magnetic Reynolds number $S\equiv\sigma L|u|$ also increases by an order of magnitude for each conductivity. When studying QGP, we limit the conductivity to be $\sigma<5\text{ fm}^{-1}$. This limit is justified by Lattice, pQCD, and AdS/CFT calculations and should keep the numerical diffusion smaller than the physical [4, 5, 6, 7, 8, 9].

A useful feature of the analytic solutions used for the initial EM fields, see Sec. IV and Ref. [25], is that it can be applied to all times. As a comparison in Fig. 2 that analytic solution for $-B_{y}$ follows a similar time evolution as our RRMHD solution. The deviations after $t\sim 1.25$ fm/$c$ between the analytic and RRMHD solutions can be attributed to different EM sources used in each calculation. At early times, both our RRMHD model and the analytic calculation only use the spectators as sources of the EM fields. Then, during the time evolution, our RRMHD model only considers QGP as a possible EM source; whereas, the analytic calculation only includes the spectators.

Overall, we can find general agreement about the nature of the magnetic field between our RRMHD model solutions and expectations. The magnetic field rapidly decreases from its initial peak within a time frame of $\sim 2$ fm/$c$. How the magnetic field decreases depends on the support of the plasma, with larger conductivity values increasing the life-time of the magnetic field, as seen in Fig. 1. We expect this dependence to be reflected in the charge dependent directed flow $\Delta v_{1}$, because of changes to the plasma’s electric charge density brought on by the Lorentz plus Faraday currents.

Suggested as an observable for the EM fields, charge dependent directed flow $\Delta v_{1}$ can be a measure of the electric currents inside QGP [15, 14]. It is defined as the difference between the observed directed flows of positive electrically charged particles and negative electrically charged particles,

Each charge dependent directed flow is calculated with the method [36],

which is dependent on $Y$ the momentum-space rapidity. As explained in Sec. II and appendix B we have adopted the Cooper-Frye formula [33, 34] for hadronization at the linear order approximation for the electric chemical potential $\mu_{Q}\simeq q/(g_{h}T_{f}^{2})$. Where $q$ is the electric charge density and $T_{f}=140$ MeV is the freezeout temperature. Finally, we use the Gauss-Legendre method to integrate over the transverse momentum $p_{T}$ and angle $\phi$ in Eq.(20).

For our model, the source of non-vanishing $\Delta v_{1}$ is a net electric charge $q$ on the freezeout hypersurface, which depends on the conductivity $\sigma$ through Eq.(5). As we have already discussed, the net electric charge is produced by the Lorentz and Faraday currents, and the initial positive charge of the plasma. Figures 3 and 4 illustrate the net electric chemical potential on the freezeout hypersurface. For a larger conductivity, the magnitudes of $\mu_{Q}$ increase because the lifetime of the EM fields have also increased, see Fig. 1 and arguments about the Lorentz and Faraday currents.

Additionally, an asymmetry in the chemical potential $\mu_{Q}$ is present in the longitudinal plane projected on Fig. 4. This asymmetry is created in peripheral collisions, i.e., $b\neq 0$ fm, because the initial magnetic fields of the spectators break the reflection symmetry across $\eta=0$. Consequently, the charged directed flow of Eq.(20) is different for positively vs. negatively charged final state hadrons.

An observable difference between the charge of protons and anti-protons is illustrated in Fig. 5, where the conductivity of QGP is varied across the panels. The rapidity $Y$ dependence of (anti-)proton charge dependent directed flow Eq.(20) separates in panels b and c. That separation is captured by a nonzero difference $\Delta v_{1}\neq 0$ from Eq.(19). This is a similar result as the observations at the STAR experiment [16].

In panel a of Fig. 5 for zero conductivity, $\sigma=0\text{ fm}^{-1}$, the difference between flows is zero because we have assumed the initial charge density is in equilibrium, i.e., $q_{0}(x,y,\eta)=0$. This means the solution for Ohm’s law in Eq.(5) is always zero. By increasing the conductivity in panels b and c the proton and anti-proton directed flow separate. This is expected from ohm’s law Eq.(5) and our discussions about Figs. 3 and 4 that larger conductivities increase $\mu_{Q}$ on the freezeout hypersurface, which results in more significant differences between the charged directed flow.

The negative linear slope of $\Delta v_{1}$ in panels b and c of Fig. 5 is a sign that the Faraday current has a greater net effect than the Lorentz current, which is in agreement with Ref. [15]. This shows the approximation of the Faraday current as $\gamma B_{y}$ and the Lorentz current as $-\gamma u_{z}B_{y}$ at the beginning of Sec. III holds for our RRMHD model. Moreover, an order increase in the electric conductivity $\sigma$ brought a similar increase in the slope, $\Delta(dv_{1}/dy)$ suggesting this observable could be used to understand $\sigma$, the electric conductivity of QGP. Although, a more refined model is required to draw any quantitative conclusions.

We have demonstrated how our relativistic resistive magneto-hydrodynamic (RRMHD) model can be used to study the charge dependent directed flow $\Delta(dv_{1}/dy)$ across collision centralities. The negative slope of $\Delta(dv_{1}/dy)$ was found to strongly depend on QGP’s electric conductivity, $\sigma$, which is a unique feature of our model. The slope of $\Delta(dv_{1}/dy)$ is caused by an electric charge imbalance on the freezeout hypersurface as illustrated in Figs. 3 and 4 in combination with the Cooper-Frye formula. This is the first work to establish such behavior and test is across centralities, which is shown in Fig. 6.

In future work, we aim to improve our model in several ways. Our results suggest that a complete modeling of $\Delta(dv_{1}/dy)$ requires a combination of baryon and electric charges. This means incorporating the Transported quark effects of Refs. [37, 38] with the electric conductivity effects of our RRMHD model. It is possible to accomplish that incorporation using a phenomenological approach as described in appendix A, or a more sophisticated model of the initial collision stages. Moreover, a dedicated analysis of how different models of the initial stages change the initial EM fields + charge density would be interesting on its own.

Another important improvement is the inclusion of vortical effects on the magnetic field. It has been established by experimental data with phenomenological models that collision participants and the resulting QGP have a large angular momentum. Theoretical models suggest that angular momentum could positively modify the time evolution of the magnetic fields. Incorporating that effect in our model should result in a better description of the electromagnetic fields.

The previous improvements focus on improving the accuracy of calculating the charge density $q$. For a consistent calculation, we should use an equation of state (EoS) that includes the electric chemical potential, for example [39]. Although, considering in this work the charge density was found to be small on the freezeout hypersurface, the impact of changing the EoS to include the electric chemical potential may also be small.

In this work, we have assumed the electric conductivity of QGP is a scalar. But, dimensional analysis shows the conductivity should be temperature-dependent and dedicated calculations of the conductivity agree [7]. Additionally, for intense magnetic fields, the conductivity splits into a component parallel with the magnetic fields and a component perpendicular. Considering both of though effects would change the conductivity from a scalar to a tensor. A future investigation about how a tensor conductivity could affect experimental observables maybe interesting.

We have chosen a single freezeout temperature and then used the Cooper-Frye formula to calculate the final particle distributions. This is quite an oversimplification of the final hadron dynamics that would be better described with an afterburner + EM field calculation. Possible candidates are the kinetic solvers UrQMD, SMASH, or JAM, and ideally, we would like our model to be able to interface with all three.

Lastly, it has been well established that viscous hydrodynamics is required to explain the elliptical flow $v_{2}$ of final state hadrons. Although, it may seem natural to include viscosity in our model, because of how it could modify the Lorentz current, we do not believe it to be as important as other improvements. For example, our results agree with Ref. [15] even though they use iEBE-VISHNU + a perturbative EM field calculation.

Even without those changes, we have demonstrated a relationship between the negative slope of $\Delta(dv_{1}/dy)$ and QGP’s electric conductivity, $\sigma$ in Figs. 5 and 6. This contrasts with other works [37, 38] that have focused on initial baryon charges and currents to explain the experimental data. In addition, we also demonstrated how initial positive charge can also cause positive slopes for $\Delta v_{1}$, as was originally discussed in Ref. [15]. All of this work suggests an understanding of charge dependent directed flow at STAR [16] and ALICE [18] calls for a consideration of QGP’s electric conductivity.