Magnetic fields in heavy ion collisions: flow and charge transport

At the earliest times after a heavy-ion collision, the magnetic field created by the spectator nucleons will generate an extremely strong, albeit rapidly decreasing in time, magnetic field. The impact of this magnetic field may have detectable consequences, and is believed to drive anomalous transport effects like the Chiral Magnetic Effect (CME). We detail an exploratory study on the effects of a dynamical magnetic field on the hydrodynamic medium created in the collisions of two ultrarelativistic heavy-ions, using the framework of numerical ideal MagnetoHydroDynamics (MHD) with the ECHO-QGP code. In this study, we consider a magnetic field captured in a conducting medium, where the conductivity can receive contributions from the electromagnetic conductivity $\sigma$ and the chiral magnetic conductivity $\sigma_{\chi}$. We first study the elliptic flow of pions, which we show is relatively unchanged by the introduction of a magnetic field. However, by increasing the magnitude of the magnetic field, we find evidence for an enhancement of the elliptic flow in peripheral collisions. Next, we explore the impact of the chiral magnetic conductivity on electric charges produced at the edges of the fireball. This initial $\sigma_\chi$ can be understood as a long-wavelength effective description of chiral fermion production. We then demonstrate that this chiral charge, when transported by the MHD medium, produces a charge dipole perpendicular to the reaction plane which extends a few units in rapidity. Assuming charge conservation at the freeze-out surface, we show that the produced charge imbalance can have measurable effects on some experimental observables, like $v_1$ or $\langle \sin \phi \rangle$. This demonstrates the ability of a MHD fluid to transport the signature of the initial chiral magnetic fields to late times.


INTRODUCTION
The collision of two ultrarelativistic nuclei deposits enough energy such that the constituents of the nucleons become liberated, forming a strongly interacting plasma, the Quark Gluon Plasma (QGP); phenomenological studies suggest that this is the most perfect fluid created in nature [1]. In the earliest moments after the collision, the system is subjected to what is expected to be the strongest magnetic field created in nature [2][3][4], of the order of eB ∼ (m π ) 2 ∼ 10 14 T. This magnetic field is generated by the spectator protons in the collision, and may be captured if the medium produced has finite electric conductivity [5]. In order to quantify the effects of the magnetic field, one must study the mutual interactions of this magnetic field and the matter produced in these collisions in a consistent framework. For the QGP, ab initio studies in terms of the fundamental theory, Quantum Chromodynamics (QCD), remain elusive.
Instead, a long-wavelength effective description of the QGP phase of the system in terms of relativistic hydrodynamics has proven very successful [6]. Coupling this mutual interaction of relativistic hydrodynamics and a magnetic field results in the effective description known as MagnetoHydroDynamics (MHD). MHD is an established and very successful tool to describe the complex dynamics of the electromagnetic fields in astrophysics: e.g. the magnetic fields around and within compact objects [7][8][9][10], protostellar accretion disks [11], and jets [12] or in the solar wind [13,14]. In this work we extend the use of the MHD framework also to the realm of heavyion collisions. In this context, the magnetic fields exhibit very short characteristic timescales, of order ∼ 1 fm or even less, a circumstance that might put into question the applicability of MHD. Actually, especially in our case, in which, for simplicity, we embrace the ideal MHD limit, that is the limit of infinite electrical conductivity (compared to the time scale of interest), subtle inconsistencies might appear [15]. The basic assumption of infinite electrical conductivity translates into imposing a null electric field in the fluid comoving frame at any time, implying a zero relaxation time and, therefore, hidden causality problems. This issue calls for the future addition of resistive effects into the code [16], including a proper numerical treatment [17] of the terms associated to the electric field with characteristic evolution time much shorter than the rest of the system. However, we recall that ideal hydrodynamics [18] has been successfully exploited to explain most of the features of the observed anisotropic collective flows [19,20] before being superseded by viscous hydrodynamics [21,22] and before investigating its foundations and the extent of its applicability [23][24][25][26][27]. Similarly, we find legitimate to use the ideal MHD framework for some preliminary analysis, leaving an improvement of the model to follow up projects and keeping a close eye on the progress toward a deeper understanding of its foundations [28,29]. The use of fully fledged MHD in heavy-ion collisions would open a new possibility to explore the electromagnetic properties of the QGP, like the electrical conductivity [30][31][32][33][34], σ, through the comparison with experimental data. However, despite the limits of our model, in this paper we will elucidate a number of novel implications of interacting magnetic fields and a perfect fluid on the observables of heavy-ion collisions.
Another exciting prospect is the ability to develop a greater understanding of anomalous transport phenomena in heavy-ion collisions, like the Chiral Magnetic Effect (CME) [2,35], the interest of the present study, and the Chiral Magnetic Wave [36]. In heavy-ion collisions, topological configurations of the gluon fields generate chiral fermions through the chiral anomaly of QCD [37]. In the presence of a strong external magnetic field, local chirality imbalance will then generate an electric current; this is the CME. While there are a number of potential sources of these topological configurations throughout the pre-equilibrium [38][39][40][41] and QGP stage of the evolution [42], due to the very short lifetime of the magnetic field, it is likely that the best window for the generation of a CME signal is during the pre-equilibrium stage. To this end, recent non-perturbative studies suggest that topological sphaleron configurations may be produced abundantly enough at the earliest time after a heavy-ion to create a viable CME signal [40,43]. In line with the longwavelength description of the QGP in terms of a fluid, we can consider that the chiral charge produced during the pre-equilibrium stage is effectively encoded in σ χ . Its effects are then incorporated into the MHD description, whereby we can study the charge separation in the fluid and the effects on the produced particles. In particular, we will demonstrate that a clear electric charge dipole, which extends a few units in rapidity, is observed.
Previous simplified studies have investigated the effects of the electromagnetic fields on the hydrodynamic description of heavy ion collisions [44][45][46], albeit neglecting the back-reaction of the fluid on the fields. Only recently has the back-reaction from the conducting current been consistently taken into account within the Magne-toHydroDynamics (MHD) framework [47,48], albeit in the ideal limit of infinite electrical conductivity. We here assume the same framework, based of the numericallyvalidated ECHO-QGP code [47]. While the description of the initial magnetic field is still an active topic of investigation (c.f. [49,50]), we consider initial magnetic fields which are the result of an electrically and/or chiralmagnetically conducting medium, with electric conductivity σ and CME conductivity σ χ . As the magnetic field is expected to decrease rapidly as a function of time, it is interesting to understand what effect a strong magnetic field may have on the produced particles and observables like elliptic flow. While only the order of magnitude of the initial field is known, signatures of this strong field may aid in determining the relevant field strength. To this end, in this paper, we study ideal (3+1)D MHD simulations using geometrical Glauber initial conditions [51] and we explore how varying basic parameters like the impact parameter, the conductivity of the medium in the pre-equilibrium phase, the freeze-out temperature and the magnitude of the initial magnetic field affects the elliptic flow of pions at mid-rapidity. The primary interest of this exploratory work is to understand the implications that dynamical magnetic fields may have on QGP observables; we do not wish to obfuscate this by including now-standard phenomenological features, such as viscosity or hadronic cascades. We therefore refrain from making theory-to-data comparison.
The outline of this paper is as follows: first, in Section I, we recall the framework that we are using and we briefly present how we compute the initial conditions for the magnetic field. In Section II, we present and discuss the results of our simulations. In particular, we perform exploratory studies of the magnetic field influence on the elliptic flow of pions. Then in Section III, we consider the electric charge imbalance created by an initially helical magnetic field, and discuss the implications for searches for the CME. In Section IV, we discuss known issues with our approach, and how these may be addressed. Finally, in Section V, we discuss the limits of our current approach and discuss future directions which will enable direct comparison with experimental data through appropriate observables.

I. SETUP OF THE NUMERICAL SIMULATIONS
We perform (3+1) dimensional numerical simulations of magnetic fields coupled to ideal relativistic hydrodynamics in Milne/Bjorken coordinates. The formalism implemented into the code as described in Ref. [47], is summarized in App. A. The initial energy density distribution is computed according to a geometrical Glauber model [52,53], whose details can be found in App. B.

A. Initial magnetic field
We compute the initial electromagnetic field produced by a point charge moving at constant velocity in a medium with electric conductivity σ and chiral magnetic conductivity σ χ as in Ref. [54]: We would like to stress that assuming that the conductivity is just a scalar quantity, constant both in space and time, clearly represents an oversimplification [32][33][34] of the properties of a system out of equilibrium and undergoing an extremely fast dynamical evolution. Moreover, this procedure relies on a semi-classic scheme that lacks some quantum features [55]. However, we recall that the scope of the present works is not to provide a complete description of the QGP formation and evolution, but just to perform an exploratory evaluation of how the magnetic fields might modify the collective flows of charged pions and induce an electric charge separation with respect to the reaction plane. Therefore, despite the strong approximations at the basis of their derivation, we consider Eqs. (1-3) sufficiently good for our purposes. As we mention also in Section II A, we consider the possibility that by using this initialization procedure we might underestimate the magnitudes of the magnetic fields, so, to partially compensate this issue, we perform several additional simulations with the values of the magnetic field components amplified up to four times. Then, we approximate the two nuclei as uniformly charged spheres that are Lorentz-contracted to disks and, by numerical integration over each of them, as explained in detail in Ref. [56], we compute the spatial distribution of the initial magnetic field. We assume that the collision does not affect the motion and the distribution of the electric charges contained in the nuclei. If we also assume that the fluid has initial null velocity in Milne/Bjorken coordinates, the transformation of the magnetic field components B i from Cartesian to Milne/Bjorken coordinates componentsB i is given by: The derivation of Eqs. (4) is provided in App. F. Of course, in the context of a future more refined model, all these approximations introduced in the pre-equilibrium dynamics should be reconsidered [57], possibly allowing for a non null initial velocity field influenced by the electromagnetic fields. Nevertheless, the development of an adequate framework for the initial conditions, combining both QCD [58] and QED, possibly with random topological charge densities [59] and assuming non constant and non uniform conductivity, requires a considerable effort, which goes far beyond the goals of this study and that must be tackled within a separate dedicated project.
We note that the addition of a magnetic field of chiral origin explicitly breaks the rotational symmetry of B x,y in the transverse plane, as the comparison between the top (no chiral B) and the bottom (with chiral B) rows of Fig. (1) shows. However in the regions where the differences between the plots are more evident, the magnitude of the magnetic field is small and it is reasonable to expect only minor effects on bulk observables.

B. Parameter set
The aim of this study is to develop an understanding of the impact of a dynamical magnetic field to the standard hydrodynamics description of heavy-ion collisions. For the present study, we do not fine tune the parameters to exactly reproduce experimental data and instead employ standard values taken from the existing literature, in particular from Ref. [45]. Unless specified otherwise, we adopt the same values for the conductivities from Ref. [54]. At the beginning of the numerical simulations, we assume that the fluid contains the same magnetic field computed as explained in the previous section, but the electrical conductivity of the medium becomes infinite, thus allowing us to apply the ideal MHD formalism. Unfortunately, right now our code is not able to handle a finite conductivity and, rather than unnecessarily extend this approximation also to the initial conditions, we preferred to maintain a more realistic scenario for them, albeit at the cost of introducing an inconsistency in the time evolution of the properties of the medium. By virtue of the ideal MHD formulation, we neglect the contribution of the initial electric fields and we assume that the fluid has vanishing initial velocity. The parameters used in the simulations are summarized in Table I. Our standard grid resolution is 0.2 fm along x and y transverse directions and 0.2 along the longitudinal η direction. The grid resolution is doubled (0.1) when we evaluate also the electric charge density.

II. COMPUTATION OF AZIMUTHAL ANISOTROPY OBSERVABLES
So-called 'flow' observables, characterized by the Fourier harmonics of the particle spectra [60], have be-   come key measurements for characterizing the QGP. It is expected that these quantities characterize both the expansion as well as azimuthal anisotropies developed through the fluid evolution. Therefore, it is a reasonable starting point to understand the influence of dynamical magnetic fields interacting with the fluid by studying these flow harmonics. Our main observables are the elliptic flow v 2 (p T ) at mid-rapidity and the directed flow v 1 (y) of charged pions. We compute the thermal spectra with the Cooper-Frye prescription [61] as described in Ref. [52]. We extend this formulation by preserving local charge conservation across the freezeout hypersurface in the local rest frame, resulting in a charge dependent spectrum; this will be discussed in greater detail later and in App. H.

A. Elliptic flow as a function of centrality
Let us start with the exploration of the centrality dependence of the effect of the magnetic fields. To this aim, we simulate Au+Au collisions at √ s NN = 200 GeV for the impact parameter b = 2, 5, 8, 10, 12 fm. Given the uncertainties in the properties and in the dynamics of the medium in the pre-equilibrium phase, it is possible that the magnitude of the initial magnetic fields may be underestimated by Eqs. (1)(2)(3). To investigate this uncertainty in the strength of the initial magnetic field, we perform systematic studies by increasing the initial field B 0 (obtained from Eqs. (1-3)) by factors of 2 through 4. Fig. (2) explores the elliptic flow in more detail. The first row of Fig. (2) shows that the effect of the magnetic field on v 2 for almost central collisions (b = 2 fm) is very small, with only a slight suppression for p T 3 GeV, with B = 4B 0 . For more peripheral collisions, however, we observe an enhancement of v 2 , but only if we increase our initial magnetic field, otherwise its effects are negligible. This is expected, as similar results have already found in simulations with transport models [62] and in analytic estimates [63]. Although the model, the numbers and many details are quite different, essentially we observe an increasing trend of v 2 by increasing the magnitude of the magnetic field as noticed in Ref. [46]. Most likely, this behavior is due to the spatial distribution of the magnetic field along the transverse plane ( Fig. 1), which is rather similar to the spatial distribution of the pressure (Fig. 13). Since the magnetic field produces a pressure B 2 /2, it gives an additional contribution to the initial pressure anisotropy of the fluid, whose gradient induces a correspondent momentum anisotropy, of which the elliptic flow is just the second momentum of the Fourier decomposition with respect to the azimuthal angle. Going on with the examination of Fig. (2), we note how in the intermediate centrality class (b = 5 fm) for RHIC energies we can already observe an enhancement of the elliptic flow which is almost absent at LHC energies. For peripheral collisions, shown in the last rows of Fig. (2)), we observe an enhancement of v 2 both at RHIC and LHC energies, however the effect is stronger at √ s NN = 200 GeV. serve that at RHIC energies the maximum enhancement is obtained for b = 8 fm and it tends to decrease for more peripheral collisions, while at LHC energies the enhancement steadily increases up to b = 12 fm. It is important to recall that these results are for a specific choice of parameters which we have not tuned to experimental data. In this case the effects of the magnetic field seem to be very modest both at RHIC and LHC energies, albeit in the second case we observe a kind of ripple at |y| ≈ 2. However, if we modify the initial conditions by introducing a tilting in the initial energy density distribution [53,64] as described in App. B 1, so to obtain a directed flow with a negative slope at mid-rapidity dN dy | y=0 < 0 , more similar to the experimental results [65,66], then the situation changes. As shown in Fig figure) is still negligible, but v 2 is suppressed, an effect opposite to magnetic field of classic origin. It is likely that the different behavior is a consequence of the different orientation of the initial fields. It is important to note that this situation is one where only the magnetic field of chiral origin is present; physically this is an unlikely situation, as the magnitude of the standard magnetic field is expected to be an order of magnitude larger. Moreover, we are treating the evolution of a magnetic field of chiral origin as a classic one, completely neglecting the evolution of axial charges. It is reasonable to expect that any axial charge produced will cascade into magnetic fields [67], further increasing the magnetic field of chiral origin. However, in a more realistic scenario, it remains unclear how strong the contribution of magnetic fields of chiral origin on v 2 is.
C. Impact of the electrical conductivity in the pre-equilibrium phase Next, we explore the impact of the electrical and chiral magnetic conductivities of the medium in the preequilibrium phase, which enter in Eqs. (1-3) and determine the magnitude and the spatial distribution of the initial magnetic field. In the context of our model, these are the only parameters which are directly related to the intrinsic properties of the medium and not to the collision energy or geometry. We consider the cases in which the conductivities of the medium in the pre-hydro phase (σ 0 ) are both scaled by a factor 0.1, 0.2, and 5. We consider the case of the void, as well. We find that simple variations of σ 0 are not sufficient to trigger significant modifications in v 2 at mid-rapidity, compared to the case with no magnetic field, both at RHIC and LHC energies, not even for peripheral collisions. However, if we enlarge the initial B fields by a factor 3, we notice that the enhancement of the elliptic flow depends on the chosen values of σ and σ χ to initialize B, with a dependence neither linear nor monotonic, as

III. CHARGE PRODUCTION AND THE CME
The basic mechanism behind the Chiral Magnetic Effect (CME) has long been understood, and has been confirmed in a variety of theory calculations, such as with lattice studies [43,68,69], and observed in condensed matter experiments [70][71][72][73]. In astrophysics, anomalous transport phenomena may explain certain observations, e.g. neutron star kicks [74]. However, a clear signal at current heavy-ion collision experiment results remains elusive [75][76][77][78]. From the chiral anomaly, it has long been known that topological configurations of the gauge fields will produce domains of axial charge (chiral fermions); this axial charge is then the seed for the CME current when an external magnetic field is present. Recent progress has been made however towards understanding sources of axial charge production in heavyion collisions. At the earliest times, the initially very strong gauge fields feature a longitudinal 'flux-tube'-like structures which may source axial charge [38,39,41]. These configurations quickly break up, whereby an overoccupied non-Abelian plasma known as the Glasma is formed [39]. In this regime, real-time topological transition, called sphaleron transitions, can produce axial charge. While this stage is very short, as the system is rapidly thermalizing, it has been determined that a few sphaleron transitions can be expected in the average collision at RHIC or the LHC [40]. As opposed to field strength fluctuations, these domains of topological charge will produce chiral charge coherently, and may allow for the growth of a sizeable CME current. During the QGP stage, it is also possible that chiral fermions are anomalously produced via sphaleron transitions [42].  However, as the magnetic field during the QGP stage of the evolution is expected to be small compared to that of the pre-equilibrium stage, it is likely that CME production after the pre-equilibrium stage is negligible. It was also argued in Ref. [40] that the sphaleron rate out of equilibrium is parametrically larger than in equilibrium, reinforcing the assumption of the dominance of pre-equilibrium CME generation. However, a full chiral-MHD description is needed to study the dynamics of chirality production during the hydrodynamic evolution. The pre-equilibrium axial charge production can be encoded, for simplicity, in terms of a uniform chiral chemical potential, µ 5 . By virtue of MHD being a longwavelength effective description, it is natural to encode this initial axial charge as a chiral magnetic conductivity, σ χ . In essence, the pre-equilibrium axial charge is cast in macroscopic language of MHD in terms of a helical (chiral) magnetic field sourced by σ χ . This however neglects any pre-equilibrium CME current generation; previous real-time lattice studies have demonstrated that these may be sizeable [43]. The motivation for this effective description of microscopic axial charge in terms of macroscopic magnetic helicity can also be understood due to a self-similar cascade [67], and may be driven by a chiral plasma instability [79].
At present, there are no quantitative calculations of the amount of axial charge produced during the preequilibrium stage of a heavy-ion collision. In light of this, we consider a range of reasonable values for the chiral magnetic conductivity, from σ χ = 1.5 MeV (as before; this corresponds to µ 5 322 MeV) to σ χ = 15 MeV. We then can study the local electric charge density generated in the fluid by the magnetic field. In Figures 9, we show this charge density in the transverse plane at midrapidity (top row) and in the y − η plane for x = 0.1 fm (bottom row) at τ = 5.4 fm. In the left column, only a magnetic field of classical origin is considered. We see that the charge densities accumulated at the edges of the fireball along the magnetic field (y) direction are manifestly symmetric. Thus no charge separation is seen. However, in the top right column, we consider only a magnetic field of chiral origin and can clearly see an electric charge dipole along the magnetic field direction. In the bottom right figure, we consider a magnetic field of both classical and chiral origin, and observe that an electric charge dipole still remains manifestly present. Re- sults for the electric charge from magnetic fields of both classical and chiral origin in the transverse plane at midrapidity are shown in Fig. 10. While the dipole structure is less apparent than the top right figure of Fig. 9, it is nevertheless still present.
Furthermore, in Fig. 10, the left and the right panels present results obtained using two different methods to compute the electric charge density in the fluid comoving frame. In the bottom left, the charges are determined directly from the Maxwell equations, while in the bottom right they are determined from the vorticity (Fig. 9 uses this method). These methods, explained in App. G, involve the computation of different partial derivatives. In principle the resulting charge density should be the same, however, in practice they can disagree, especially in regions with strong gradients. It is important to stress two points about the charge densities shown in Figs. 9 and 10. First, they are located in the region (in Bjorken coordinates) close to the expansion wave, under the as-sumption that the whole computation domain is filled by a fluid, but in general the region where the QGP transforms into particles, i.e. the freezeout hypersurface, is different. Second, the detectors measure a particle distribution in the momentum space, therefore it is crucial to understand how to translate this spatial charge distribution into some experimental observables.
In our model we compute the spectra of the emitted particles by using the Cooper-Frye prescription [61], so we should modify this recipe in such a way to take into account a contribution due to the electromagnetic fields. Unfortunately, the underlying assumption of the ideal MHD approach that the electric field in the comoving frame of the fluid is always zero, with an infinitely small relaxation time, makes modifying the distribution function in the denominator of the Cooper-Frye equation in a consistent way very difficult [80,81]. It will be possible to overcome this limitation with resistive MHD. However, in this work we appeal to a simple, physically motivated solution based on the charge density in the comoving frame of the fluid. Essentially, we compute the electric charge density, as measured in the comoving frame of the fluid, for each cell of the freeze-out hypersurface. We then modify the number of the produced positive and negative pions in such a way to obtain the same net charge density. With this procedure, the distribution of the momenta of the emitted hadrons in each cell is still insensitive to the electric charge. Details for this procedure are discussed in App. H. While physically motivated, we acknowledge that this method is imperfect. Nevertheless, our aim is to perform an exploratory study on the effects of the dynamics of magnetic fields on a fluid and not quantitative phenomenology, so such an approximation seems justified.
We evaluate the effects of our charge dependent freezeout procedure for charged pions in Au+Au collisions at √ s NN = 200 GeV with impact parameter b = 8 fm, initializing the energy density distribution with tilting parameter η m = 2. We look at the directed and elliptic flow and at the sin φ observables, where φ is the angle between the direction of transverse momentum of the pion and the reaction plane.
On the left side of Fig. 11 we show the directed flow v 1 of pions versus the rapidity, while on the right side we show the difference ∆v 1 = v 1 (π + ) − v 1 (π − ) between the directed flow of positive and negative pions (multiplied by a factor 10 3 ). The right side of Fig. 11 demonstrates a fair agreement between the results obtained by computing the electric charge density using the Maxwell equations or vorticity method.
We note that the sign of ∆v 1 with respect to rapidity seems to be the opposite of what the preliminary experimental measurements suggest [82], albeit these exper- imental data refer to unidentified charged particles (so, not only pions) with respect to pseudorapidity at a different collision energy ( √ s N N = 5.02 TeV). We will discuss this further in Sec. IV. In the left side of Fig. 12 we show the elliptic flow v 2 of pions versus rapidity, while in the right side we show the difference ∆v 2 between the elliptic flow of positive and negative pions (multiplied by a factor 10 3 ). The right side of Fig. 12 shows some clear discrepancies between the results obtained by computing the electric charge density with the Maxwell equations or with the vorticity method, probably due to the numerical errors in summing over two different sets of gradients and to the limited accuracy (second order) in solving the evolution equations. However, in any case the difference between the v 2 of positive and negative pions seems to be rather small. The sin φ 1 observable (with φ the angle between the direction of transverse momentum of the pion and the reaction plane) aims to isolate particle production in the direction of the magnetic field, perpendicular to the reaction plane. Therefore, it can act as a simple signal of magnetic field driven effect like the CME. If the average of sin φ is different from 0, it means that there is an unbalance between the number of particles emitted up or down along the direction orthogonal to the reaction plane. The results of this calculation for different choices of magnetic field are reported in Table II. We observe that when considering an initial magnetic field which includes a contribution from chiral charges (second row), there is a clearly non-zero sin φ . However, if the magnetic field is generated by classical vector currents only, sin φ 0 (first row). Thus, we can clearly demonstrate the conversion of initial chiral magnetic fields in final state charge asymmetry. The change in the sign of sin φ might be explained by the following considerations. The CME current is proportional to µ 5 and magnetic field. The chiral magnetic field is equivalent locally to some value of µ 5 -the chiral charge of the chiral knot of magnetic field corresponds locally to a certain chiral charge, and thus to a chiral chemical potential µ 5 . So the CME current for the chiral magnetic field only is J chiral ∼ µ chiral 5 B chiral . Now, for chiral + classical, the total magnetic field changes by about an order of magnitude, while µ 5 stays the same, so the CME current is J chiral+classical ∼ µ chiral 5 B chiral+classical . Since B chiral+classical is an order of magnitude larger than B chiral , the CME current and the resulting charge separation should be an order of magnitude larger as well. Moreover, if the direction of B chiral+classical is different from the direction of B chiral , the sign of separation can be different as well. In Ref. 40 it is shown that, in the transverse plane, the components of classic magnetic fields have an azimuthal orientation, while chiral magnetic fields have a radial orientation. Going beyond the aims of the current work, it would be interesting to study more elaborate magnetic field driven charge dependent observables, like those of Refs. [83][84][85], but probably it is wiser to dive into this project only after the development of a more refined version of our model.   sin φ (π + ) sin φ (π − ) ∆ sin φ Classic B −1.49 · 10 −15 −1.67 · 10 −15 1.8 · 10 −16 ≈ 0 Chiral B −2.82 · 10 −6 +2.82 · 10 −6 −5.64 · 10 −6 Classic + chiral B 8.16 · 10 −5 −5.54 · 10 −5 1.37 · 10 −4 Table II: First two columns: sin φ at |y| 0 for positive and negative pions, third column: ∆ sin φ = sin φ (π + ) − sin φ (π − ). The results refer to numerical simulations of Au+Au collisions at √ sNN = 200 GeV with impact parameter b = 8 fm. In these computations the electric charge density was computed by using the vorticity method (App. G 2). When including an initial magnetic field with a contribution from chiral charges, the results turn from a roughly null value (first row) to a non zero signal (second row)

IV. KNOWN ISSUES
In the previous section, we noted that there is tension between our results in right panel of Fig. 11 for the rapidity dependence of ∆v 1 (y) = v 1 (π + ) − v 1 (π − ) compared to recent ALICE results [82]. A similar discrepancy between theoretical and experimental results has been observed in a complementary framework in Refs. [5,81], whereby the electromagnetic fields are added to the fluid evolution perturbatively (see also [86][87][88][89] for related calculations for heavy flavors). There, the authors considered finite conductivity of the medium and therefore have both electric and magnetic fields. They identified two classical magnetic field driven mechanisms which resulted in charge-and rapidity-odd contributions to the Fourier harmonics: I) the Faraday effect, resulting from changing magnetic fields, II) the Lorentz force, resulting from a charged particle moving in a magnetic field. They observed that these effects contribute with opposite signs to charge-dependent v 1 (y). Using this as a basis to speculate on the nature of the apparent disagreement between theory and experiment on the sign of ∆v 1 (y), it is possible that in our formalism, and potentially others, there is a delicate interplay between the properties of the magnetic field in the medium, related to the rate of expansion in comparison to the decrease in the magnetic field with time, which the inclusion of a temperature dependent conductivity and viscosity might drastically alter. There also several other possible explanations for the behavior of ∆v 1 (y) obtained in the present simulations. A possible source of error might be the prescription to determine charge dependent spectra, which assumes a kind of chemical potential to modify the particle species abundances, without taking into account any modification of the momentum distribution due to the electromagnetic field. Moreover, in this study the local charge imbalance is attributed only to pions, neglecting the contributions of other less abundant hadrons. The computation of the electric charge density itself shows some evidence of a limited accuracy, e.g. in Fig. (12), which might be improved either adopting constrained transport schemes [90,91] or other frameworks in which this quantity is explicitly evolved [28,29]. There are also notable improvements possible for the initialization procedure. To begin, by virtue of the ideal MHD framework, the electric field is only a derived quantity and therefore, by setting the initial velocity field to zero, the initial electric field is also set to zero; we thus neglect its contribution. However, in Ref. [81] the electric field has been shown to be of the same order of magnitude of the magnetic field. In addition to that, the magnetic field generated by the protons is assumed to be instantaneously transferred to the QGP fluid. This instant quench should be replaced with a characteristic interaction time profile. Another contribution, albeit probably less important, might come from nucleon and sub-nucleon fluctuations in the initial conditions, which we neglect in the current implementation. This may produce non-zero localized charge densities, inducing non-linear effects, which are not captured by the present approach. Finally, the solution utilized to reproduce the experimental v 1 slope, tilting the initial energy density distribution [53,64], deserves further study and improvement. If we do not utilize this tilting, for LHC energies we obtain the same sign of the slope of ∆v 1 (y) as seen in experiment, while for RHIC energies the slope remains unchanged.
Given the considerable effort required in refining this model and extensively testing a wider set of alternative initial conditions, we postpone this task to a future work. We believe that this current framework is an instructive first step towards the ultimate goal of a dissipative resistive chiral MHD formulation, which will require the development of a more advanced numerical framework.

V. SUMMARY AND OUTLOOK
In this study, we performed (3+1)D one fluid ideal MHD simulations of relativistic collisions of gold nuclei at √ s NN = 200 GeV, for RHIC, and of lead nuclei at √ s NN = 2.76 TeV, for the LHC. We adopted a somewhat simplified approach to make possible a first systematic study of the evolution and influence of backreacting magnetic fields on common heavy-ion observables. Taking into account the possibility of an underestimation of the initial magnetic field, we studied the effects of amplifying the magnitude of the magnetic field components by up to a factor four. Compared to "pure" hydrodynamics, we found an enhancement of v 2 in peripheral collisions, with the greatest change found for the largest values of the magnetic field. Otherwise, this enhancement is almost negligible, in overall general agreement with Refs. [46,48,62,63]. At RHIC energies, this enhancement tends to be stronger, and with a sufficiently large magnetic field, it is evident already in semi-peripheral collisions. In the initialization of the magnetic field B, we took into account also a contribution of chiral origin, although we neglected any fluctuations or dissipation of axial charge. We studied how different values of the electric and chiral magnetic conductivities of the medium in the pre-equilibrium phase can impact the final v n . We found that a simple modification of these parameters is not sufficient to produce significant changes in v 1 and v 2 . However, the initialization of the electromagnetic fields, based on the simple, but unrealistic assumption of constant scalar conductivities for all time [30][31][32][33][34][92][93][94], should be replaced by more sophisticated modeling.
We then computed the electric charge density in the comoving frame of the fluid. When initializing the magnetic field with a contribution originating from a preequilibrium CME, represented in terms of a chiral conductivity σ χ , we detected an electric charge separation in the direction orthogonal to the reaction plane. However, since this is not a quantity that can be directly measured in the experiments, it is necessary to translate this charge imbalance into an observable which can be experimentally evaluated. The ideal MHD assumption prevents a simple modification of the Cooper-Frye prescription to take into account the different electric charges of the hadrons, therefore we have presented an alternative method based on electric charge conservation, whereby an effective electric charge chemical potential is determined and used to ascertain the charge dependent pion spectra. We observed splitting between the directed flow of positive and negative pions that is of the same order of magnitude as recent experimental results [82], however with a different ordering of π + compared to π − . We note however that our results are in qualitative agreement with a number of recent theoretical calculations [5,81,[86][87][88][89]. It is therefore a challenge to all of these descriptions to understand the physical mechanism which may result in this different ordering. Finally, within our simplified setup we were able to compute a non-zero average value of the sine of the emission angle of the particles with respect to the reaction plane, a necessary signature of the CME. This motivates further refinements of our model, in hopes of making more quantitative predictions of the CME.
The most obvious improvement to our approach is the inclusion of resistive effects, already present in many relativistic MagnetoHydroDynamic codes for astrophysics [7,91,[95][96][97][98][99], together with the viscous hydrodynamic corrections (these exist already in ECHO-QGP [52], but are presently kept separate from the MHD module). From the theory side, there are recent advancements towards a deeper understanding of the connections between the kinetic theory and MHD [28,29] that can help in improving the consistency of the formalism. For a proper quantitative investigation of the CME within the chiral-MHD framework, it will be necessary to consider the evolution of axial charges coupled with the fluid [100][101][102], which one may be able to derive from chiral kinetic theory [103][104][105][106][107][108][109][110][111][112]; a possible intermediate step could be the simplified, but still interesting, formalism illustrated in Ref. [113].

Appendix A: Ideal relativistic MagnetoHydroDynamics
We follow the relativistic MagnetoHydroDynamics (MHD) [114,115] approach, based on the conservation laws for a one-fluid current N µ and for the total (matter and electromagnetic fields) energy-momentum tensor of the fluid T µν , namely with d µ being the covariant derivative. In addition, we consider the second law of thermodynamics where s µ is the entropy current, and, for the electromagnetic fields, the Maxwell's equations where F µν is the Faraday tensor and F µν = 1 2 µνλκ F λκ is its Hodge dual. When polarization and magnetization effects are neglected, the electromagnetic contribution to the energy-momentum tensor can be expressed as: In the ideal limit all dissipative fluxes can be neglected, local equilibrium is assumed and we can write a single fluid four-velocity u µ (u µ u µ = −1). After introducing the electric and magnetic fields measured in the comoving frame of the fluid: the Faraday tensor and its dual can be decomposed with respect to u µ as Since we are dealing with electromagnetic fields strongly coupled with the fluid, we need a relation (Ohm's law ) between the electric currents and the fluid, which is adequately modeled, in many cases, by a linear dependence: whereρ e is the electric charge density in the comoving frame, j µ the conduction current (j µ u µ = 0), and σ µν the fluid conductivity tensor. The ideal MHD approximation, adopted in this work, consists in assuming σ → ∞, so that, to avoid the onset of infinite currents, this requires e µ = 0 (A12) must hold, as well. The adoption of Eq. (A12) brings many simplifies the Faraday tensor and its dual, as well as the structure of the evolution equations; at the end are given by the following system (see Ref. [47] for details): in the unknowns n, e, p, u µ , and b µ , all defined in the comoving frame of the fluid. In Eqs. (A13-A15), n = −N µ u µ is the baryon density, e = T µν m u µ u ν the fluid energy density, and p = 1 3 ∆ µν T µν m the kinetic pressure, 1 2 b 2 = 1 2 b µ b µ the energy density of the magnetic field, and g µν is the metric tensor.
However, to solve numerically the system of equations (A13-A15), we need to perform a 3+1 splitting [116,117] of time and spatial components and rewrite the equations in conservative form [118]. For more details on this procedure, see Refs. [52,119,120]; here we summarize the most important steps. We consider two possible metrics: diag(g µν )=(−1, 1, 1, 1) for Minkowski space-time or diag(g µν )=(−1, 1, 1, τ 2 ) for Milne (Bjorken) space-time. Given the fluid velocity v i for an Eulerian observer in the laboratory frame, the fluid four velocity can be expressed as where γ = (1 − v 2 ) −1/2 is the Lorentz factor of the bulk flow and v 2 = v k v k . We now consider the three spatial components of the electric field E i and magnetic field B i as measured in the laboratory frame, which are related to the four vectors e µ and b µ by: where ε ijk is the Levi-Civita pseudo-tensor of the spatial three-metric, i.e. ε ijk = |g| therefore in ideal MHD the electric field is a not an independent quantity, but is derived from v i and B i . We can then rewrite Eqs. (A13-A15) in a form appropriate for numerical integration by clearly separating time and space derivatives and tensor components: where are respectively the set of conservative variables and fluxes, while the source terms are given by The components of energy momentum tensor T µν are: in which u em = 1 2 (E 2 + B 2 ) is the energy density. Finally, the solenoidal condition coming from the time component of Eq. (A15), is enforced by using a modified version of the hyperbolic divergence cleaning method [16,99,[121][122][123][124]. This procedure has been validated through several tests [47]. However, we verified that using or not the hyperbolic divergence cleaning method does not have a noticeable impact in many results presented in this paper, although its use has a positive influence on the stability of the code. Nevertheless, since this method controls the error on the divergence of B by transporting it away with a dumping effect, it might introduce some spurious effects in the computation of the electric charge density and, therefore, we did not use it in the simulations in which this quantity was evaluated. In the future, it may be desirable to use the Constrained Transport approach [90,91].

Appendix B: Initialization of the energy density distribution
We define the thickness function: where n 0 is the nuclear density, δ the width and R the radius of the nuclear Fermi distribution, which depends on the specific nucleus. Then, we define the following functions: where σ is the inelastic NN cross section, A the mass number of the colliding nuclei, and where x T = (x, y) is the vector of the transverse plane coordinates and b is the impact parameter vector, connecting the centers of the two nuclei. The function W N , which gives the contribution of the wounded nucleons, is then simply: The initial proper energy density distribution is given by: where the total weight function W (x, y, η) is defined as: (B7) where n BC (x, y) is the mean number of binary collisions: and α is the collision hardness parameter, which can vary between 0 and 1. The energy density profile in the longitudinal direction is modulated by while η f lat and σ η are parameters.

Tilted initial energy density distribution
In some cases, in order to obtain a directed flow in a better agreement with the experimental results [64], we modify Eq. (B5) and redefine the wounded nucleons weight function W N as: where: while Eq. (B9) becomes: Essentially, this kind of initialization introduces a tilting in the initial energy density distribution, regulated by the parameter η m (see Ref. [64] for details). As the left side of Fig. (14) shows, this tilting changes significantly the directed flow of pions and, with η m = 2, it is able to reproduce quite well the experimental data by the STAR collaboration [125] in the central rapidity region, although, of course, for a proper comparison with experimental data viscous effects, hadronic rescattering and resonance decay feed-down should be included in the model as well. However, as shown in the right side of Fig. (14), the elliptic flow at mid-rapidity seems not to be influenced by the tilting, at least for the values of η m that we tested. Given the weak sensitivity of v 2 to the tilting of the initial energy density distribution, we employed this initialization only in a few cases regarding v 1 , in particular when studying charge dependent spectra, like in Figs. (11,12).
Appendix C: Effects of the energy density pedestal on v1 and v2 At the beginning of the simulation, the energy density created by the collision between two nuclei is very high inside the fireball at the center of the grid, but, at least in principle, it should be zero outside of it (see, e.g., Fig. (13)). However, we cannot run hydro simulations in the void, therefore we must add a minimum energy density pedestal. The magnetic field, instead, has a large magnitude even in the space regions outside the fireball, as can be noticed in Fig. (1). When the magnetic pressure is much larger than the thermal pressure, the code can easily crash. To avoid this situation, the artificial additional energy density pedestal must be considerably larger than in the pure hydro case. Fig. (15) shows the ratio between the initial magnetic pressure and thermal pressure. It is important to verify that this larger minimum energy density layer does not change appreciably the Fourier harmonics v n compared to the case when a much smaller value is used. We make the test in pure ideal hydrodynamical simulations, assuming that the sit-uation does not significantly change in the MHD case, and we look at the p T integrated directed flow vs rapidity and at the elliptic flow at mid-rapidity vs transverse momentum of charged pions. At RHIC energies, whose results are shown in Fig. (16), we can notice a small, but appreciable deviation in v 1 (y) for |y| > 4 (left), while the effects are almost negligible for v 2 (p T , y = 0). At LHC energies, whose results are shown in Fig. (17), the effects are almost negligible both for v 1 and v 2 .
Appendix D: Effects of the reconstruction algorithm on v1 and v2 ECHO-QGP inherits from the ECHO code many high order reconstruction algorithms. In general lower order algorithms tends to be more diffusive, but, when stability is a problem, this weakness becomes an advantage in term of robustness of the code. In this work we use the second order TVD2 algorithm. In this section we check that, for our purposes, the choice of the reconstruction algorithm does affect significantly the final results. This statement is demonstrated by Fig. (18), showing the integrated directed flow vs rapidity (left) and the elliptic flow at midrapidity vs transverse momentum (right) of charged pions produced in Pb+Pb collisions at √ s NN = 2.76 TeV with an impact parameter b = 8 fm, neglecting the effects of magnetic field (pure ideal hydro simulations).
Appendix E: Effects of the freeze-out temperature In all our simulations we choose as freeze-out temperature 154 MeV, the best value given in Ref. [126], with  an uncertainty of ±9 MeV. We want to check that a different choice of the freeze-out temperature within these limits does not affect significantly the final results. To this aim, we run two additional simulations per collision energy with impact parameter b = 12fm, using as values of the freeze-out temperature 145 MeV and 163 MeV. In Fig. (19) we compare the elliptic flow of pions computed in these cases with the standard case. The results show, as expected, that the magnitude of the elliptic flow depends inversely on the value of the freeze-out temperature (the larger the freeze-out temperature, the smaller v 2 ), however the magnetic fields always produce a similar enhancement of the elliptic flow. Temperature dependent effects, which might certainly appear in a more quanti-tative analysis, are out of the scope of the present work.
Appendix F: Transformations from Minkowski to Milne/Bjorken coordinates For the reader's convenience, we derive the transformation laws of the components of the three dimensional electromagnetic field from Minkowski coordinates x µ = (t, x, y, z), with metric diag(g ij )=(−1, 1, 1, 1) to Milne coordinatesx µ = (τ, x , y , η), with metric diag(g ij )=(−1, 1, 1, τ 2 ). The Milne coordinates are re-  lated to the Minkowski coordinates by: while the opposite transformation is given by: From now on a tilde character will be used to distinguish a tensor or a vector in Milne/Bjorken coordinates (F ) from the same object in Minkowski coordinates (F ). Keeping this convention in mind, we express the Faraday tensor in Minkowski and Milne coordinates as: In both cases, we run simulations with a fixed impact parameter b = 12 fm. We assume to have an initial magnetic field 4 times larger than our estimates. We explore the effect of different freeze-out temperatures. We notice how the presence of a magnetic field always enhances the elliptic flow compared to the pure hydro case.
We recall that the formula to change coordinates for a rank 2 contravariant tensor is: The only non trivial derivatives are: We obtain:Ẽ x = cosh ηE x − sinh ηB y (F9) E y = cosh ηE y + sinh ηB x (F10) In the ideal MHD framework: therefore: Now we want to express the four velocity u µ M , in Minkowski coordinates, but using Milne/Bjorken variables, of a fluid with Cartesian velocity components v x = v y = 0 and v z = z/t: which is the initial velocity profile used in this work. In this case the transformation laws (F9-F11) simplify to: while Eqs. (F14-F13) simplify to: B y = cosh ηB y − sinh η(tanh η)B y = cosh 2 η − sinh 2 η cosh η B y = B y cosh η .
where the B k are the components of the three dimensional magnetic field, or, with a small rearrangement of the terms: Assembling the previous equation together, we arrive to the final expression for electric charge density in the fluid comoving frame as: (G9) We recall that, in ideal MHD, the electric field in the lab frame is a quantity given by Eq.(A19).
In cartesian coordinates (but still using natural units), Eq. (G9) has the simplified expression:

Method based on the vorticity
Here we briefly summarize the same method described in the appendix of Ref. [91]. After defining the kinematic vorticity four-vector as [127] ω λ = µνλκ ∇ µ u ν u κ = µνλκ ∂ µ u ν u κ , we split the covariant derivative of the fluid velocity as where a µ = (u ν ∇ ν )u µ is the acceleration. We notice that ω µ u µ = 0 and a µ u µ = 0, i.e. both the kinematic vorticity and the acceleration are orthogonal to fluid velocity. If we use the definitions of e µ and b µ , i.e. the electric and magnetic fields in the fluid comoving frame, given in Eqs. (A7,A8), we find the relation F µν ∇ µ u ν = e µ a µ + b µ ω µ .
We recall that, according to Eq. (G4), −I µ u µ =ρ e , whereρ e is the electric charge density in the fluid comoving frame, therefore: If we replace the last term of the previous equation with Eq. (G13), we end up with the expressioñ ρ e = ∇ µ e µ − e µ a µ − b µ ω µ .
From Eq. (A12) we already know that, in the ideal MHD approximation e µ = 0. Moreover, if we assume that Ohm's law e µ = η r j µ (G17) holds, with the resistivity η r constant, then ∇ µ e µ = η r ∇ µ j µ . Therefore, if, as in the ideal MHD limit, we assume that η r = 0, Eq.(G16) simplifies intõ Appendix H: Charge dependent freeze-out In the charge dependent freeze-out we modify the distribution functions of pions f 0 according to a factor λ, dependent on the local charge density, so that the distribution function of positive and negative pions become f + = (1 + λ)f 0 and f − = (1 − λ)f 0 , obtaining at the end the same total number of pions 2f 0 , but with different abundancies of π + and π − .
From a practical point of view, first we evaluate the expected average π ± density n π in the fluid rest frame according to a Bose-Einstein distribution Assuming a chemical potential µ 0, given our chosen freezeout temperature T = 154 MeV, for the charged pion mass m = 139.6 MeV we get n π 0.04462, which corresponds to a charge density ρ 0 c = ±n π √ 4πα c 0.006 (GeVfm) 1/2 . If the net charge density ρ c in a freezeout hypersurface cell is non zero, we assume that this due to a charged pion production imbalance, i.e. for one half to an enhancement of pion production with the same sign of the charge and for one half to a suppression of pion production with the opposite sign charge by a factor 1 + s · ρ c /(2ρ 0 c ), where s is the sign of the electric charge of the pion. We do not take into account other hadron species. This approach has several shortcomings: 1) as we already mentioned, it deals only with the number of produced pions, but not with their momentum distribution, assumed to be isotropic in the fluid comoving frame, thus neglecting all local electric currents 2) the computation of the electric charge density is imperfect, as we have demonstrated by the differences in the two methods employed here, which are particularly evident in the right side of Fig. (12). Moreover, there is no mechanism to prevent unrealistic charge densities, with enhancing/suppressing factors larger than 1, albeit typically we are far from this situation, 3) although pions are the most abundant species of hadrons produced in high energy heavy ion collisions, many other kinds are also produced (for a more complete hadronization which also respects local charge conservation, see e.g. [128]) 4) since the ions are positively charged, the QGP is not neutral (this, of course, is a general issue of the RMHD approach).