Numerical magneto-hydrodynamics for relativistic nuclear collisions

We present an improved version of the ECHO-QGP numerical code, which self-consistently includes for the first time the effects of electromagnetic fields within the framework of relativistic magnetohydrodynamics (RMHD). We discuss results of its application in relativistic heavy-ion collisions in the limit of infinite electrical conductivity of the plasma. After reviewing the relevant covariant $3\!+\!1$ formalisms, we illustrate the implementation of the evolution equations in the code and show the results of several tests aimed at assessing the accuracy and robustness of the implementation. After providing some estimates of the magnetic fields arising in non-central high-energy nuclear collisions, we perform full RMHD simulations of the evolution of the Quark-Gluon Plasma in the presence of electromagnetic fields and discuss the results. In our ideal RMHD setup we find that the magnetic field developing in non-central collisions does not significantly modify the elliptic-flow of the final hadrons. However, since there are uncertainties in the description of the pre-equilibrium phase and also in the properties of the medium, a more extensive survey of the possible initial conditions as well as the inclusion of dissipative effects are indeed necessary to validate this preliminary result.


I. Introduction
High-energy nuclear collisions, studied by several experimental collaborations at RHIC and at the LHC, allow one to explore the QCD phase-diagram in the high-temperature region, from high to almost vanishing baryonic density. Strong evidence, coming both from soft and hard observables, was obtained for the onset of a deconfined phase in the RHIC and LHC energy regime. Furthermore, at the experimentally accessible conditions (i.e. slightly above the deconfinement phase-transition), the produced system, with a lifetime ∼ 10 fm/c, was found to behave like a collective, stronglyinteracting medium, rather opaque to penetrating probes, in contrast to the expected gas of weakly-interacting quarks and gluons. Relativistic hydrodynamic models (nowadays including also dissipative effects) were developed to describe the evolution -driven by pressure gradients -of the produced matter and turned out to reproduce the data quite well [1][2][3][4][5][6][7][8], in particular the various flow-harmonics arising from the collective response of the system to the anisotropies and fluctuations in the initial conditions.
While the main purpose of relativistic heavy-ion experiments is the study of strong interactions at extreme energy densities similar to the early universe, it was recently realized that during the collisions of high-Z nuclei (Z = 82 for Pb) at ultra-relativistic energies, one can also produce the strongest magnetic fields reached in our universe, with initial values of B ∼ 10 15 T and oriented mainly in the direc- * inghirami@fias.uni-frankfurt.de tion perpendicular to the reaction-plane [9]. In the last years it was suggested [9,10] that, besides leading to the production of a strongly-interacting deconfined system, the presence of these strong magnetic fields in relativistic heavy-ion collisions opens also the possibility of exploring peculiar nonperturbative features of QCD, such as the appearance of nontrivial topological configurations of the color-field. Once coupled to quarks, these configurations characterized by a nonvanishing winding number lead to an excess of quarks of a given chirality (chiral anomaly), depending on the value of the topological charge, and hence, on an event-by-event basis, to a violation of parity (clearly preserved after an event-average).
In the presence of strong magnetic fields this can give rise to observable effects, with a separation of oppositely-charged particles with respect to the reaction-plane. Since for massless particles with a fixed handedness (e.g. right handed quarks) the chirality coincides with the helicity (i.e. the projection of the spin along the particle momentum) and since particles tend to align their magnetic moments along the B-field, one would have an excess of positively-charged u-quarks moving in the direction of the magnetic field and an excess of negative d-quarks moving in the opposite direction. Clearly, averaging over a large sample of events, each one with a different excess of right or left-handed quarks, the effect should cancel at the level of single-particle distributions; however, it should leave its fingerprints in multi-particle correlations, as suggested in [11]. Due to the interplay between a nonperturbative feature of strong interactions (the chiral anomaly) and the role of the magnetic field, such a phenomenon was called Chiral Magnetic Effect (CME) and is currently studied by different experimental collaborations at RHIC and at arXiv:1609.03042v3 [hep-ph] 11 Jan 2017 the LHC [12][13][14]. Analogous effects have been recently observed also in astrophysics (as an explanation of Neutron Stars kicks) [15] and in solid-state physics, placing Dirac semimetals in parallel magnetic and electric fields [16][17][18][19]. Other related phenomena (Chiral Magnetic Wave [20], Chiral Separation Effect [21], Chiral Vortical Effect [22]), all arising from an unbalance among right and left-handed particles and from the presence of a strong magnetic field or angular momentum, were suggested to occur in non-central heavy-ion collisions: for an overview we refer the reader to [22].
An unambiguous observation of the CME in heavy-ion collisions would be clearly a result of deep theoretical interest, since it would represent a manifestation of the nontrivial topological structure of a Yang-Mills theory. However, in order to separate opposite-sign charges with respect to the reaction-plane, the initial magnetic field generated by the colliding nuclei must be sufficiently long lived. The lifetime of the magnetic field depends strongly on the nature of the produced medium. In the vacuum the initial magnetic field decays rather rapidly. On the contrary in the opposite limit, in the presence of an ideal plasma with infinite electric conductivity, the freezing of the magnetic-flux makes the field survive much longer and may allow for the manifestation of signatures of the possible chiral unbalance in the final charged-hadron spectra, even though, at the same time, a large conductivity would also tend to compensate any local charge excess. Unfortunately, so far in the literature one can find only semi-analytic estimates of the time-evolution of the magnetic field in heavy-ion collisions, based on simplifying assumptions [23][24][25][26][27][28]. A fully realistic calculation would require to solve the the Maxwell equations together with the continuity equations for the energy-momentum tensor (closed by some form of Ohm's law), i.e. it calls for a full Relativistic Magneto-HydroDynamic (RMHD) description of the medium, in which the evolution of the electromagnetic field is consistently coupled with the evolution of the plasma: this is the challenge we address with the present paper.
For this first study we consider the case of an ideal plasma, with no dissipative effects and, in particular, an infinite electric conductivity, which makes the electric field in the local restframe of the medium vanish. We also neglected any anomalous term in the currents, although previous studies [29,30] in simplified models showed that they would not to contribute to entropy production, being in this sense "ideal": the inclusion of dissipative and anomalous terms (necessary for the description of the CME) in our setup is left for future work. In light of the small experimental uncertainties reached at the LHC and RHIC on flow measurements the development of a code able to consistently treat the coupled evolution of the plasma and Z-enhanced electromagnetic fields represents in any case a necessary baseline for any claim that CME (and other related phenomena that we will be able to address after including anomalous currents) can be disentangled from possible other confounding electromagnetic effects that could lead to charge separation.
Our paper is organized as follows. In Sec. II we present the RMHD equations in their most general form, focusing then on their ideal limit, i.e. on the case of a plasma with infinite electrical conductivity (and neglecting other dissipative effects such as viscosity and thermal conduction). Only the ideal case is considered for the present paper. In Sec. III we discuss the numerical implementation of the ideal-RMHD equations, written in a conservative form, within our improved ECHO-QGP code. In Sec. IV we discuss the results of a large variety of numerical tests to prove the accuracy and the robustness of the implementation: the shock-tube problem, the description of Alfvén waves, the rotor test, the reproduction of the one-dimensional Bjorken expansion in a magnetic field and the accurate treatment of the in-vacuum self-similar expansion in transverse-MHD. In Sec. V we show the results obtained from the code with simplified (but reasonable) initial conditions for non-central nucleus-nucleus collisions. At least in the context of this simplified approach, the magnetic field is not able to modify the elliptic flow of the final hadrons substantially. Nevertheless, further and more realistic investigations are needed before solid conclusions can be drawn. Finally, in Sec. VI we discuss our findings and the future perspectives of our work, with the idea of performing 3D+1 simulations based on a much broader pool of different initial conditions, possibly including dissipative effects. The appendix is devoted to a discussion of the propagation of linear perturbations in RMHD, focusing on the case of fast-magnetosonic and Alfvén waves, which are the ones relevant for the analysis carried out in this paper.

II. Ideal relativistic magnetohydrodynamics
Relativistic MHD (RMHD hereafter) is a one-fluid description of the interaction of matter and electromagnetic fields in plasmas [31,32]. In general, as in the Newtonian limit of classical MHD, one assumes that there is a dominant species determining a main fluid current, while a secondary species must be responsible for the conduction current, namely the source for the electromagnetic field. The RMHD evolution equations describing the dynamics of the overall system are the conservation laws for this fluid current N µ (associated to the net-baryon current or to any other conserved charge, if any) and for the total (matter and fields) energy-momentum tensor of the plasma T µν , namely with d µ being the covariant derivative, thus to be supplemented by the second law of thermodynamics where s µ is the entropy current. On the other hand, the electromagnetic field obeys Maxwell's equations where F µν is the Faraday tensor and F µν = 1 2 µνλκ F λκ is its dual. Notice that here we have neglected possible polarization and magnetization effects of the plasma, therefore we do not make a distinction between microscopic and macroscopic fields [33]. Under this assumption, the electromagnetic contribution to the energy-momentum tensor is known to be for which d µ T µν f = J µ F µν , from Maxwell equations. Introducing the matter contribution to the energy-momentum tensor T µν m and letting T µν = T µν m + T µν f , Eq. (2) gives where the right-hand-side is the Lorentz force acting on the plasma.
In the ideal limit all dissipative fluxes can be neglected and local equilibrium is assumed. A single fluid four-velocity u µ (u µ u µ = −1) can be thus defined and we write where we have introduced the projector ∆ µν = g µν + u µ u ν (∆ µν u ν = 0). In the above zeroth-order relations n = −N µ u µ is the main charge density, e = T µν m u µ u ν the fluid energy density, and p = 1 3 ∆ µν T µν m the kinetic pressure, all quantities are defined in the comoving frame. The Faraday tensor and its dual can also be split with respect to u µ as where e µ = F µν u ν , (e µ u µ = 0), are the electric and magnetic fields measured in the comoving frame of the fluid. Since the electromagnetic fields do not evolve in vacuum, but are strongly coupled with the fluid, we must now provide an appropriate Ohm law relating the current with the fields. In the simplest case one usually assumes the linear form where ρ e is the electric charge density in the comoving frame, j µ the conduction current ( j µ u µ = 0), and σ µν the plasma conductivity tensor. The presence of a finite conductivity in the plasma gives rise to (anisotropic) magnetic dissipation and Joule heating, as well as to topological field line changes known as magnetic reconnection. Recent theoretical and numerical results may be found in [34] and references therein. In the ideal MHD approximation considered in the present paper we assume a conductivity high enough to avoid the onset of huge currents in the plasma. We can then replace the Ohm law with its limiting case When the above condition holds, the expressions for the Faraday tensor and for its dual are simplified, and the number of unknowns is reduced. In particular, Eq. (4) will be used to derive the current, if needed, while Eq. (5) will become the evolution equation for b µ . Moreover, the electromagnetic energymomentum tensor becomes where b 2 = b µ b µ , which can be plugged into Eq. (2) together with the corresponding matter contribution in Eq. (9). Summarizing, the system of ideal RMHD equations is in the unknowns n, e, p, u µ , and b µ . Non-conservative versions of the above equations can also be found. It is useful to decompose the covariant derivative as where D ≡ u µ d µ indicates derivation along u µ (reducing to the Eulerian time derivative in the nonrelativistic limit), and ∇ µ = ∆ ν µ d ν is the derivative transverse to the flow (reducing to the spatial gradient in the nonrelativistic limit). The charge conservation (baryon-number in the case of heavy-ion collisions) becomes where θ ≡ d µ u µ = ∇ µ u µ is the expansion factor. The energy equation is derived by projecting the d µ T µν = 0 conservationlaw along the flow u ν , where, we remember, the total energymomentum tensor is given by the sum of the matter and field components: T µν = T µν m + T µν f . From Eq. (7) we get which leads to Written in the above form, the energy equation is rather general, the right-hand side representing the Joule heating of the fluid. However, as previously discussed, in ideal MHD the electric field in the local rest-frame vanishes, e µ = 0, thus one simply has independent of b µ , as in ordinary relativistic hydrodynamics. This form of the energy equation will be exploited in discussing the Bjorken-flow of a magnetized plasma in Sec. IV D. However, if the two contributions are kept together, we may also write The relativistic extension of the MHD Euler equation is retrieved by projecting the total energy-momentum conservation law transverse to the flow, that is Several expressions may be derived from the last RMHD equation for the evolution of b µ , here we choose to rewrite it as where we have used the relation d µ b µ = b µ Du µ . Finally, the system of ideal RMHD equations must be closed by choosing an equation of state (EoS), for instance of the form p = P(e, n), under the assumption that in the ideal case each local equilibrium state can be completely determined by u µ and two thermodynamical variables (e and n in this case). The Euler and Gibbs-Duhem relations read e + p = T s + µn, de = T ds + µdn, where we have defined the local temperature T = (∂e/∂s) n and the chemical potential µ = (∂e/∂n) s . Eqs (29), (25), and (22) allow us to write We then retrieve the expected result that in the ideal case, when all dissipative terms are neglected, there is no entropy production and Eq. (3) holds as an equality. Notice that the entropy current is conserved even in the case of vanishing charge (baryon-number) density and chemical potential n = µ = 0, as appropriate for high-energy heavy-ion collisions and an ultrarelativistic EoS with p = P(e).

III. The RMHD module in ECHO-QGP
We now rewrite the evolution equations for ideal RMHD in a form suitable for numerical integration, for which we need a clear separation between time and space components (the so-called 3 + 1 split) and the preservation of the original conservative character of the equations, since shock-capturing numerical codes such as ECHO-QGP require to solve a series of balance laws. Here we will provide the basic expressions, for further formal and technical details details see [6,35,36] and references therein.
Neglecting curvature effects due to gravitational fields, we consider here a metric in special relativity (though not necessarily Minkowskian) of the form where the three-metric coefficients g i j may depend both on space x i and time x 0 , in general. It is first useful to introduce the fluid velocity v i and electric and magnetic fields E i and B i as measured in the laboratory frame, which are spatial vectors (vanishing time component). 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 , whereas the fields are, respectively where ε i jk is the Levi-Civita pseudo-tensor of the spatial three-metric, namely ε i jk = |g| which is known once v i and B i have been determined. In this case the b µ field is Let us now rewrite Eqs. (18)(19)(20) in a form appropriate for numerical integration, by clearly separating time and space derivatives and tensor components. We find the system where are respectively the set of conservative variables and fluxes, while the source terms are given by where the symmetric and antisymmetric properties of T µν and F µν , respectively, have been exploited in deriving the above balance laws. The components of T µν appearing in the expressions for the conserved variables and fluxes are where we have defined the electromagnetic energy density u em = 1 2 (E 2 + B 2 ). We recall that while B i is a dynamical variable, E i is a derived quantity, obtained from Eq. (35).
One final constraint comes from the time component of Eq. (20), that is the solenoidal condition which, if valid at the initial time of the evolution, should be preserved analytically by the last equation of the above RMHD system. From a numerical point of view, however, this constraint needs some specific techniques to be actually enforced. In fact, the accumulation of the numerical errors associated to the computation of the derivatives of the magnetic field may lead to the violation of the solenoidal (i.e. "null-B divergence") condition (44), implying the formation of unphysical magnetic monopoles and fictitious forces. There are several methods to avoid, or at least to limit, this issue [37][38][39][40]. We adopted the method proposed by Dedner for MHD and later extended to the cases of special and general relativity [41][42][43][44][45][46].

A. Numerical procedures
ECHO-QGP is based on finite difference schemes. At the beginning of the simulation, the initial values of the primitive variables n (the baryon density), v i (the contravariant components of the velocity of the fluid in the lab frame), p (the pressure of the fluid in the comoving frame) and B i (the contravariant components of the magnetic induction field in the lab frame) are discretized on the computational grid by evaluating them at the center of each cell. Time integration of conservative variables is performed using a second or third order Runge-Kutta algorithm, then, at each sub-timestep: • the values of the primitive variables are reconstructed at cell borders, for each direction (several algorithms are implemented and can be selected [36]: TVD2, CENO3, WENO3, WENO5, PPM4, MPE3, MPE5, MPE7), • fluxes in Eq. (39) are computed, • the Riemann problem for fluxes at cell interfaces is solved using the HLL (Harten-Lax-Van Leer) [47] approximate method, • the divergence of these numerical fluxes and source terms in Eq. (40) are computed at cell centers, allowing to integrate the discretized evolution equations for the conservative variables, • the new primitive variables are retrieved from the evolved conservative ones.
This last step above implies to solve a system of non-linear equations and currently there is no known algorithm which guarantees a global convergence to the solutions. The system is more easily solved by providing an initial guess for the solution, usually chosen as the values of the primitive variable at the previous timestep. However, in a rapidly evolving system as in the case of heavy ion collision, this guess may not be close enough to the real solution and the algorithm may fail or converge to other (unphysical) solutions. Nevertheless, if we restrict to the use of a specific analytic Equation of State (EoS), then the system of non linear equations may be considerably simplified and it is possible to develop very robust inversion routines [36,48].
For the present study we focus for sake of simplicity on the ultra-relativistic gas EoS p = e/3, using an "ad hoc" version of the method described in [36], hereafter shortly summarized. We exploit Eq. (35) to rewrite equations (41) and (43), then we compute S 2 = S i S i and S i B i , which are known since B i is both a conservative and primitive variable (the difference is only in the factor |g| 1 2 ). After introducing the new variables x = v 2 = v i v i and y = 4pγ 2 , with some algebraic manipulations we can formulate the following system of equations: These coupled non-linear equations are solved through a nested procedure: Eq. (45) is solved for x with a one dimensional iterative hybrid Newton-Raphson/bisection method [49] with bracketing between 0 and 1; at each iteration of this routine, the y variable is obtained by finding the (unique) positive root of the third order polynomial of Eq. (46) multiplied by y 2 with x = x(y). The solution of the system allows then to compute the primitive variables through the relations: For EoS where the pressure p depends also on the baryon density n, like the ideal gas EoS used in [36] and in the shock tube test presented here, note that the latter quantity can be easily obtained by dividing the corresponding conserved variable by the Lorentz factor γ. However, for a comparison to high energy HIC data, a lattice QCD based equation of state should be employed [50] (in contrast to the simplified EoS used for the present study), which unfortunately does not allow to simplify the system of non linear equations on which the inversion routine is based and needs a more careful (and slower) numerical treatment as discussed above.

IV. Tests
In this section we present some numerical test problems selected in order to validate the code. We avoid to repeat tests aimed at simply measuring the accuracy of the "core" algorithms, since ECHO-QGP for relativistic hydrodynamics [6,7] has been already validated against basic benchmarks, and many additional tests have been performed on the original ECHO code [36], from which ECHO-QGP has been derived sharing the same base structure. Instead, here we focus on checking the correctness of its results in the ideal RMHD context. We use the ultrarelativistic EoS p = e/3, if not mentioned otherwise.
We will use either Minkowski (t, x, y, z) or Milne [τ, x, y, η s ] coordinates, where τ ≡ √ t 2 − z 2 is the longitudinal propertime and η s ≡ 1 2 ln t+z t−z the space-time rapidity. In the following, in writing four-vector components in Milne coordinates, we will employ square brackets. Notice that in the first case the three-metric is g i j = diag{1, 1, 1}, with |g| 1 2 = 1, whereas for Milne coordinates g i j = diag{1, 1, τ 2 }, with |g| 1 2 = τ. In both cases ∂ j g ik = 0 and the source terms in the evolution equations simplify considerably. Notice that in Milne coordinates, where g 33 = τ 2 , the source term for the energy equation contains a non-vanishing term proportional to 1 2 ∂ 0 g 33 = τ.

A. Magnetized shock tube
In order to test the shock-capturing properties of ECHO-QGP for relativistic MHD, we run a 1D shock-tube test in Minkowski coordinates comparing the numerical results against the solutions of the same problem computed by the exact Riemann solver developed by Giacomazzo and Rezzolla [51]. Since the cited solver works for an ideal-gas EoS, for the present test we impose with an adiabatic index Γ = 4/3, where ρ = nm stands for the mass density in the comoving frame (m is the rest mass and n is the number density of the conserved species), in a situation in which particle creation/annihilation is negligible, so that (e − ρ) is the thermal energy density. To employ Eq. (48) in this test, when retrieving the primitive variables we used the same method described in Ref. [36]. The initial conditions for the non-vanishing quantities are provided in [52] and listed in Table (I) using proper dimensionless units.
The test runs from an initial time t = 0 to a final time t = 4, with a grid resolution of 0.0025 (400 cells per unit of length). Results are displayed in Fig. (1). The comparison shows excellent agreement between the RMHD implementation in ECHO-QGP and the exact result.
B. Large-amplitude CP Alfvén-wave A multi-dimensional relativistic MHD test with an exact 1 solution is provided by the propagation along the diagonal of a square numerical domain of a large-amplitude Circularly Polarized (CP) Alfvén-wave [36].
We consider a Cartesian X − Y − Z frame, rotated along Z ≡ z in the x − y plane in such a way that X coincides with the diagonal y = x of the numerical domain. A relativistic MHD CP Alfvén wave is defined by the magnetic field and velocity components where B 0 is the uniform background field, the dimensionless parameter η = B 2 X + B 2 Y /B 0 sets the scale of the perturbation, and φ is the phase. For propagation along X we have φ = k(X − v A t), where k = 2π/λ is the wave-number and the relativistic Alfvén velocity for arbitrary large amplitudes η is given by [36]: (50) We remind that in our ideal MHD approach the electric field is given by Eq. (35) and we notice that the quantities Here we use the ultrarelativistic EoS p = e/3, where p and e remain constant to their initial uniform values p 0 and e 0 . Notice that, as expected, for small amplitudes the Alfvén speed in Eq. (50) correctly reduces to the expression derived in Appendix A 2 for the linearized case. With the above assumptions the CP Alfvén wave has a period T = λ/v A , so that at time t = nT , with n any integer number, the numerical solution is expected to assume the same configuration as at t = 0. In the following we will consider the case of a perturbation with wavelength λ = L/2, where L is the length of the diagonal of the x − y domain. We perform the test in a square numerical domain [0, 2π , so that L = 4π, discretized with a grid of 512 x 512 cells, choosing p 0 = e 0 /3 = B 2 0 = 1 and also a large amplitude of the wave η = 1 and a unit wave number k = 1 (so that λ = 2π = L/2).
In Fig. (2) we compare the z components of the velocity and of the magnetic field for t = 5T , that is after n = 5 periods, along the diagonal of the grid y = x. Neither deformations nor phase lags are observed for the depicted components as well for the other quantities not shown here. The accuracy obviously depends on an adequate numerical resolution and on the order of time and spatial integration. Further details can be found in Ref. [36]. Finally, note that large-amplitude Alfvén    waves, even if exact solutions of MHD equations, may be unstable on long timescales due to coupling with compressive modes [53,54].

C. Rotor test
We now describe a modified version of the 2D "rotor" test [35,43], here both in Minkowski and in Milne coordinates, using the the ultrarelativistic EoS p = e/3.
An initially rigidly rotating disk of radius r 0 is threaded by a constant magnetic field, causing a rapid slow down of the motion. In the previous examples found in the literature the disk is denser than the surrounding medium, but, since in our case the density does not have any influence on the evolution of the system, because the EOS does not depend on it, we assume that the region inside the disk has an initial thermal pressure larger than the region outside. After this modification, the new test proposed here becomes a sort of mixture between the "rotor" and the "magnetized cylindrical blast wave" tests [35].
The initial velocity of the fluid is null outside of the disk, while inside the disk its components are: in the case of Minkowski coordinates, while, in Milne coordinates, v z is substituted by v η = 0, which amounts to assume a longitudinal Bjorken expansion v z = z/t.
The values of the parameters chosen for the test are listed in Table (II). The major difference between the results in the two coordinate systems is the decay of the thermal and magnetic pressures in the case of Milne coordinates, which occurs in every region of the grid, due the longitudinal expansion of the system. Then, in both cases we observe a compression wave, due to the larger initial inner pressure and due to the motion of rotation of the disk, forged into an asymmetric shape by the effects of the magnetic field.

D. Bjorken flow
This test consists in a comparison with the analytical solution for the temporal evolution of a one-dimensional boost-invariant flow, obtained extending the model by J.D. Bjorken [55] to the case of transverse MHD [56]. We consider the relativistic flow along the z-direction of an ideal magnetized fluid, with pressure p and energy density e, related by the ultrarelativistic EoS p = e/3, both constant in the transverse x − y plane and independent from the space-time rapidity η s (one employs Milne coordinates). For the flow profile one considers a longitudinal boost-invariant Hubble-law expansion v z = z/t, leading to a four velocity u µ = (cosh η s , 0, 0, sinh η s ). In Milne coordinates the fluid velocity reads simply u µ = [1, 0, 0, 0], so that for the comoving derivative and the expansion rate one has D = ∂ τ and θ = 1/τ. The transverse MHD hypothesis, i.e. the assumption of having a magnetic field b µ = (0, b x , b y , 0) orthogonal to the fluid velocity u µ , so that u µ b µ = 0, allows one to derive from Eq. (26) the energy-conservation equation [56] ∂ However, under the hypothesis of infinite conductivity, one has also from Eq.(25) This allows one to obtain the evolution equation for the magnetic field: Considering the case of an ultrarelativistic p = e/3 EoS, it is possible to derive from the above the time evolution of the energy density and of the magnetic field: and Notice that, in an ideal plasma, due to the flux-freezing condition, the magnetic field decreases according to the same law as the conserved charges or of the entropy. We perform the test for three different values of the initial magnetization σ 0 = b 2 0 /e 0 : 0, 1 and 10 (in adimensional units). The comparison with the analytic results (shown in Figures  (5) and (6)) shows perfect agreement between the simulation and the exact solution.

E. Self-similar expansion in vacuum
With the purpose of performing a non-trivial validation of our numerical code, here we consider an exact solution of  Within ideal RMHD, we recall that the relation between the magnetic field in the comoving (b µ ) and laboratory (B) frames is given by Eq. (36):   and if all quantities are constant in the transverse plane, the set of equations reduces to The above equations have to be solved together with the one providing the evolution of the magnetic field in the plasma which, in ideal MHD where E = −v × B, leads to Writing explicitly the derivatives one gets: We now wish to address the case of a plasma, initially at rest, with magnetic field, pressure, energy and entropy density b 0 , p 0 , e 0 and s 0 for z < 0 and vanishing on the right. We want to study how the system evolves in time, extending the study performed in [57] to the case of an ultra-relativistic plasma of massless particles. For this purpose, it is useful to introduce the self-similar variable ξ ≡ z/t, which allows one to rewrite the equations as: In Ref. [57] the system was closed by combining the induction equation for the magnetic field with the one for mass conservation. Actually, in the case of heavy-ion collisions, such a choice would not be meaningful, since one deals with an ultrarelativistic plasma of massless particles, in which particleantiparticle pairs are continuously created/annihilated. However, in the absence of dissipative effects, one can replace the conservation equation for the mass with the one for the entropy. One can write the conservation law d µ s µ = ∂ µ s µ = 0 for the entropy current Entropy conservation can be expressed by Eq. (30) or, here more conveniently, in terms of its density in the laboratory frame: Introducing the Lagrangian derivative d/dt ≡ (∂ t + v ·∇) and combining Eqs. (60) and (65) one gets: In the transverse one-dimensional MHD case we are addressing one has then: This allows one to rewrite the set of RMHD equations as (the prime index denotes the derivative with respect to the selfsimilar variable ξ) The equation for the entropy, together with the rather general EoS p = c 2 s e (here c s is the sound speed), leads to The system has a non-trivial solution only if the determinant vanishes, i.e. if A rarefaction wave propagates from the outside inside the plasma. The position of the rarefaction front, characterized by a vanishing value of the fluid velocity v = 0 and with all the other quantities equal to their initial unperturbed values is given by One gets then which, in the case of and ideal ultrarelativistic gas EoS, reduces to in agreement with what obtained for the fast magnetosonic speed in Eq. (A8) of A 1. Hence, with the initial condition we chose, the position of the rarefaction front propagates backwards with a velocity equal to the fast magnetosonic speed: z rf (t) = −c f t. We now look for an explicit solution written in terms of the ratio B between the initial thermal and magnetic pressure. We will try to follow an approach as close as possible to the one employed by Lyutikov and Hadden [57]. In the case of an ideal ultra-relativistic plasma one has p ∼ T 4 and s ∼ T 3 , so that One gets then Exploiting Eq. (68c) one obtains In the approach by Lyutikov (generalized to our ultrarelativistic case) one writes the above equation in terms of the parameter and variable Hence, we get which we can recast as The latter is equivalent to Eq. (6) in the paper by Lyutikov, except that now it depends only on the parameter B (thermal pressure and particle/entropy density are not independent variables in an ultra-relativistic plasma) and it is does not include the term arising from the mass density. The above equations can be equivalently written in terms of the variables One obtains From the first equation we define From the second equation one gets then The latter can be easily integrated, obtaining ln δ ξ (s 1 ) where δ ξ 0 can be fixed through the initial condition, namely the development of a left-propagating rarefaction-wave, with velocity equal to the fast magnetosonic speed: parameter. We used a grid of 801 cells, the reconstruction algorithm MPE5, the approximate Riemann solver HLL and the time integration algorithm was a second order Runge-Kutta. The initial pressure was: left side (z ≤ 0) p 0 = 1000, right side (z > 0) p 0 = 5 · 10 −5 ≈ 0 (due to numerical reasons).
In Fig. (7) we display a comparison between the above semi-analytic solution and the numerical result provided by our code. The graph shows s 1 = s/s 0 vs ξ = z/t at t = 20 for three different values ( 10, 1 and 0.1 ) of the B = 2p 0 /B 2 0 parameter. We used a grid of 801 cells, reconstruction algorithm: MPE5, approximate Riemann solver: HLL, time integration algorithm: second order Runge-Kutta. Initial pressure was: left side (z ≤ 0) p 0 = 1000, right side (z > 0) p 0 = 5 · 10 −5 ≈ 0 (due to numerical reasons, since ECHO-QGP cannot run with true null pressure). Again we observe excellent agreement between the numerical implementation and the analytical results for a large variety of parameters.

V. Results of RMHD simulations for HIC
We plan to present a more extensive study of the QGP evolution in a subsequent article, nevertheless here we present some preliminary results to evaluate the impact that the interplay between magnetic field and hydro evolution may have on some experimental observables. Although the whole 3D+1 formalism has been already implemented into the code, for simplicity here we will show a basic 2D+1 application to Heavy Ion Collisions.  where x ± = t ± v/z and e = √ 4πα, α being the fine structure constant. We mention that the B field has dimensions [GeV 1/2 fm −3/2 ], so that B 2 has the same dimensions as the pressure, i.e. [GeV/fm 3 ].
Then, we approximate the electric charge distribution inside the two colliding nuclei as being uniform and spherical and we perform an integration over it to get the total magnetic field in each point of our computational grid. We assume that the motion and the distribution of the electric charges are unaffected by the collision between the nuclei. A detailed description of the whole procedure can be found in Ref. [23]. Since at the moment our code is not able to handle configurations where the magnetic pressure is much larger than the thermal pressure, which is the case in regions outside the fireball, where the initial energy density is less than 30 MeV/fm 3 we rescale the magnetic field so that the ratio between the magnetic and the thermal pressure does not exceed 0.1. This procedure does not affect the final results because at such low temperature there is no participating QCD matter and the hydrodynamic description of the medium would cease to be valid anyway.
Our choices of the parameters for the initial conditions are summarized in Table (III). The initial distribution of the thermal pressure, the magnetic field and the ratio of thermal to magnetic pressure in the transverse plane are shown in Figs. (8), (9) and (10) for Au+Au, b=10 fm reactions at √ s NN =200 GeV. We always use the same initial conditions for the initial energy density distribution, but for the initial magnetic field we consider two cases: 1. B = 0 (no magnetic field) 2. B 0 and σ = 5.8 MeV In the first case we consider a pure hydrodynamical simulation, without magnetic field. In the second case we assume that in the pre-equilibrium phase there is a medium with finite constant electrical conductivity σ = 5.8 MeV, which allows to compute an initial magnetic field distribution as in Ref. [23], shown in Fig. (9).
We assume that at the time τ 0 the fluid is in local thermal equilibrium, its electrical conductivity σ becomes infinite and that the magnetic field generated by the fast moving electric charges contained in the protons of the nuclei is converted into  Table (III). The parameters are for the reaction Au+Au, b=10 fm at √ s NN =200 GeV.
the magnetic field of the fluid, while, consistently with the hypothesis that initially the fluid is at rest and it has infinite electrical conductivity, we assume that there is no initial electric field in the fluid frame (otherwise, for Eq. (35), we should have also initial non null fluid velocity). We neglect dissipative effects and we assume that the fluid obeys the e = p/3 EoS. We run the simulation until thermal freeze-out, when the energy density is below 150 MeV/fm 3 . Then we compute the spectra and the elliptic flow of the pions produced. Here we adopt the Cooper-Frye prescription [6,60], without any modification to the distribution function due to the electromagnetic interaction.

B. Results
In Fig. (11) and (12) we compare the decay of the magnetic field in the ideal 2D+1 RMHD simulation in the center of the of overlap region of the two nuclei (i.e. in the center of the grid: x = y = z = η = 0) with some common analytical models. Fig. (11) shows the comparison between the decay of the magnitude of B in the center of the grid during the 2D+1 RMHD evolution and the decay expected for a Bjorken flow, following the analytic law τ 0 /τ. Fig. (12) show the comparison of the time evolution of the magnitude of the magnetic field (in neutral pion mass units squared) at the center of the grid in five different cases: (a) ECHO-QGP 2D+1 RMHD evolution starting from initial conditions as described in this section, with σ = 5.8 Figure 9: (color online) The initial spatial distribution of the components of the magnetic field B, computed using the method described in Ref. [23]. The parameters are for the reaction Au+Au, b=10 fm at √ s NN =200 GeV. We notice that the expansion of the fluid in the transverse plane leads to a faster decrease compared to the case of a pure longitudinal Bjorken-flow [56] and tends to become roughly exponential. However, the decay of the magnetic field of the fluid is still slower than in the case that the fields are generated by two electric charges moving in opposite directions in a uniform medium with constant finite electrical conductivity, as in Ref. [23], especially if there is no medium at all and the electric charge propagates in empty space. We stress that this comparison between different decay rates is based on a simplified model of HIC. In a 3D+1 simulation, adopting a more realistic EoS and including dissipative effects, the decay rate of the B-field might be considerably quantitatively different. In Fig. (13) and (14) we compare the elliptic flow and the transverse momentum distribution of pions, computed with the Cooper-Frye prescription [6,60], with and without the presence of an initial magnetic field, computed as described in the previous section of this article. According to our current results, the presence of a magnetic field with a magnitude and spatial distribution evaluated according to Ref. [23] seems to have a negligible impact both on the pion spectra and on the elliptic flow. This is in contrast to Ref. [62] where was suggested that the magnetic field might substantially influence the anisotropic flow. In Ref. [63] it was indeed found that a significant enhancement of the elliptic flow might be possible. A direct comparison with our results is however not possible because of the many differences compared to our approach. However, in contrast to Ref. [62,63] and the present study, Ref. [61] reported the opposite result, namely a reduction of the anysotropic flow. This was attributed to the effects of the magnetic squeezing. However the model at Ref. [61] does not satisfy the divergence-free condition for the magnetic field. There the magnetic field has a rather large magnitude and it is not completely coupled with the fluid.

VI. Conclusions, discussion and outlook
We presented the extension of the ECHO-QGP code to the relativistic magnetohydrodynamic regime, in the limit of infinite electrical conductivity, i.e. without taking into account any resistive effect. In the present version, the code has been tested with an ideal-gas EoS, either in the presence of a finite mass-density or in the ultrarelativistic regime (p = e/3). After introducing the physics equations on which the code is based, we gave an overview of their numerical implementation. Then, we illustrated the results of several tests to validate the implementation. Since our final aim is to exploit the code to study the evolution of the Quark-Gluon Plasma formed in Heavy-Ion collisions, we showed first applications in this context, adopting simplified initial conditions. Due to the (on average) small ratio of the magnetic to thermal pressure, the magnetic field does not seem to significantly affect the fluid evolution and we observed only a tiny effect on inclusive hadronic observables such as the elliptic flow and transverse momentum spectra of pions. However, in our approach the magnitude of the initial magnetic field could have been underestimated, possibly because in the pre-equilibrium phase we considered the electrical conductivity σ as constant, while there are some evidences that it increases with the temperature [64][65][66][67]. Other authors, employing different initial conditions for the magnetic field, found a non-negligible effect of the latter on the hadron elliptic-flow [61][62][63]. Clearly this would affect the estimate of the viscosity-to-entropy η/s ratio obtained by comparison of hydrodynamic results with experimental data: if, for example, part of the hadron v 2 in noncentral collisions arose from the magnetic field, one should reduce the contribution from the hydrodynamic expansion, via for instance a larger value of η/s.
Our preliminary results suggest also that the formation of a deconfined conductive plasma, compared to the case of the vacuum, might slow down the decay of the initial magnetic field generated by the colliding nuclei, possibly affecting nonperturbative phenomena relying on the presence of huge magnetic fields to show up. Since our study refers to the case of an ideal plasma, with infinite electrical conductivity, our results have to be considered as an upper limit on the lifetime of the magnetic field produced in heavy-ion collisions.
However, the recent estimates both from lattice QCD computations [64][65][66] and fitting of experimental data [67] point toward high, but finite value for the electrical conductivity of the QGP. For a quantitative comparison with experimental data this has to be taken into account including the effects of the electrical resistivity. We expect a considerably acceleration of the decay of the magnitude of the magnetic field compared to our studies.
As a next step, we plan to evaluate better the role of the initial magnitude and spatial distribution of the magnetic fields, performing full 3D+1 simulations, already possible with the present setup. This will allow one to explore a broader range of possible initial conditions under different models [68], using a more realistic EoS.
The next development of the code will involve the inclusion of dissipative effects (shear and bulk viscosity and a finite electric conductivity), using the numerical techniques presented in [6] and [34] and already implemented in previous versions of the ECHO code, going beyond the approximation of an ideal plasma. A major conceptual achievement would be represented by the inclusion in our setup of anomalous currents, allowing one to provide a consistent description of the CME and to estimate the possibility of disentangling it from other charge-separating effects related to the presence of strong electromagnetic fields.
Then, indeed, it would be necessary to modify the Cooper-Frye formula by taking into account the presence of an electromagnetic field and of a non uniform spatial distribution of electric charges. After that, for a proper comparison with experimental data, one should compute the effects on the final particle spectra and on collective flows, in the post-freeze-out phase, of decays, elastic collisions and of magnetic deflections by the Lorentz force.
Finally, we deem that applications of numerical calculations performed with the present relativistic MHD version of the ECHO-QGP code could be also relevant for cosmological (generation of the primordial magnetic fields [69]) or astrophysical studies. For instance, the sudden transition from an hadronic to a QGP-like equation of state in a protomagnetar (phase transition to a quark star) has been recently suggested as a possible explanation for the observed cases of (long) Gamma-Ray Burst events with double prompt emission peaks [70]. In this appendix we want to present a study of the propagation of small perturbations in a relativistic plasma embedded in a constant magnetic field. Although this represents a standard MHD subject, we think it is useful for the reader to explicitly re-derive the main results for the case of an ultrarelativistic plasma addressed in this paper, with no conservation equation for the mass density, at variance with usual astrophysical studies. We then perform small fluctuations around a homogeneous background, keeping in the equations only terms linear in the fluctuations. Taking into account that γ ∼ O(δ 2 ) one has (we consider the case of a one-dimensional flow along the z-axis) u µ = [1, 0, 0, δv], p = p 0 + δp, e = e 0 + δe, b µ = b µ 0 + δb µ . (A1) Notice that the index 0 in the magnetic field is used to denote its unperturbed background value and not as a covariant index.
Clearly, fluctuations in the pressure and energy density are related by the Equation of State.

Magnetosonic waves
We firs want to evaluate the velocity of propagation of magnetosonic disturbances. This will be relevant for the study of the self-similar one-dimensional flow described by the Lyutikov solution given in Sec. IV E. We focus then on the propagation along the z-axis (i.e. δ = δ(t, z)) of the following perturbations u µ = [1, 0, 0, δv]+O(δ 2 ), b µ = [0, b 0 +δb, 0, 0]+O(δ 2 ) (A2) where, to linear order in the fluctuations, B 0 = b 0 and δB ≈ δb, so that one can identify the magnetic field in the laboratory and in the comoving frame. The system of RMHD equations reduces to Let us now perform a Fourier analysis of the fluctuations, inserting in the above the ansatz δ = δ ω,k e −iωt+ikz . From the last equation, one gets for the magnetic field (turning out to fluctuate in phase with the velocity) δb ω,k = b 0 (k/ω)δv ω,k , , which corresponds to the zero mass-density limit of Eq. (3) of the paper by Lyutikov and Hadden. In terms of the thermal to magnetic-pressure ratio one gets . (A10)

Alfvén waves
Alfvén waves are MHD excitations which propagates along the lines of the unperturbed magnetic field. In full generality we will consider the evolution of the following perturbations (still neglecting O(δ 2 ) terms in the fluctuations) where we take δ = δ(t, x ⊥ ): we will see that only the dependence on x, i.e. the direction of the unperturbed magnetic field, matters. We start considering the equations for the evolution of the components of the magnetic field. To linear order in the fluctuations we have: If initially absent, no field perturbation develops along the x and y directions, perpendicular to the velocity fluctuation. Hence, in the following we set δb x = δb y = 0. On the other hand, from Faraday's law one has so that, employing the Fourier ansatz δ = δ ω,k e −iωt+ik x x+ik y y (k ⊥ = (k x , k y ) = (k cos θ, k sin θ)), one gets δb ω,k = −b 0 (k cos θ/ω)δv ω,k .
The magnetic field develops a z-component, fluctuating in opposition of phase with respect to the velocity. Let us now move to the equation for the energy and the fluid velocity.
For the Euler equation one gets instead: In Fourier space one has then (e 0 + p 0 + b 2 0 ) ω δv ω,k + b 0 k x δb ω,k = 0, which, employing Eq. (A14), leads to The perturbation propagates then along the x-axis (the direction of the unperturbed magnetic field) with group velocity equal to the Alfvén speed v x g = (dω/dk which corresponds to the weak-fluctuation (η → 0) limit of the exact result quoted in Eq. (50). Assuming an ideal-gas EoS e 0 = 3p 0 , the latter can be expressed in terms of B as in agreement with Eq. (3) of Lyutikov paper [57], once setting to zero the contribution from the mass-density.