Non-relativistic approximate numerical ideal-magneto-hydrodynamics of (1+1D) transverse flow in Bjorken scenario

In this study, we investigate the impact of the magnetic field on the evolution of the transverse flow of QGP matter in the magneto-hydrodynamic (MHD) framework. We assume that the magnetic field is perpendicular to the reaction plane and then we solve the coupled Maxwell and conservation equations in (1+1D) transverse flow, within the Bjorken scenario. We consider a QGP with infinite electrical conductivity. First, the magnetic effects on the QGP medium at mid-rapidity are investigated at leading order; then the time and space dependence of the energy density, velocity and magnetic field in the transverse plane of the ideal magnetized hot plasma are obtained.


Introduction
It is commonly accepted nowadays that collisions of relativistic heavy-ions create a hot and dense fireball matter. Quarks and gluons are in a deconfined state, called a quark-gluon plasma (QGP), for a very short time (∼ 1 fm/c) after the initial hard parton collisions of nuclei. The hydrodynamic approach has given one of the best descriptions for the QGP matter: especially for estimating the lowest shear viscosity over the entropy ratio, this theoretical framework has shown acceptable consistencies with many experimental results [1][2][3][4][5][6].
Recently it has been shown that in the peripheral AAcollisions such as Pb-Pb at the center of mass energy √ s = 2.76 TeV and Au-Au at the center of mass energy √ s = 200 GeV a huge magnetic field is created, of the order of eB∼ 10 18 − 10 19 G, which is 10 13 times larger than the strongest steady magnetic field ever realized in the laboratory. It has been claimed that the existence of such strong a e-mail: afarzaneh@hsu.ac.ir fields may be important for a variety of new phenomena like the Chiral Magnetic Effect (CME), Chiral Magnetic Wave (CMW), Chiral Electric Separation Effect (CESE), Chiral Hall Separation Effect (CHSE), pressure anisotropy in QGP, influence on the direct and elliptic flow, shift of the critical temperature. A series of reviews and more references can be found in Refs. [7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25][26]. Hence, it will be worth to further investigate the properties of the QGP in the presence of EM fields.
There have been several works which have explored the behavior of the space-time evolution of electromagnetic fields created by the colliding charged beams moving at relativistic speed in z-direction, as a solution of the Maxwell equations [27][28][29][30][31][32][33]. In Refs. [27][28][29][30][31][32], for the sake of simplicity, the classical Ohm law for induced currents in QGP has been suggested. In many studies, it has been assumed that there are no couplings between the electrodynamic and hydrodynamic equations in a QGP medium. Based on this assumption, it has been shown that the electromagnetic field depends only on the impact parameter of the colliding nucleons b, on the center of mass energy √ s and on the electric and chiral magnetic conductivities of the QGP; besides, its decrease with time is much slower than in vacuum. In addition, in several works (see, e.g., Refs. [27][28][29][30][31][32][33]) the electromagnetic field is derived from the Maxwell equations without coupling to the velocity of the fluid, assuming that the latter has negligible influence on the field itself.
According to this point, considering the Bjorken flow four velocity, the electromagnetic response of QGP in a quantum regime has been investigated in Ref. [30], and it has been concluded that the induced electric current in the plasma fireball cannot generate a classical electromagnetic field. In [35] charged dependence of flow coefficients has been discussed, and the effects of the EM field on directed flow has been studied, showing that these effects are negligible. However, we claim that this controversial result has been obtained by imposing the requirement that the velocity of charged particles v is smaller than the velocity of the expanding plasma u.
Other work obtained a series of preliminary results, by estimating the significance of strong EM fields on the QGP medium [36][37][38][39][40]. In most cases, it has been assumed that the Maxwell equations decouple from the time evolution of the QGP, and then the evolution of the EM fields and their influence on the flow coefficients have been studied. Results have revealed that after collision, the strength of the EM fields decrease. In addition, it was found that the ratio of magnetic pressure over the thermal pressure b 2 /P of the hot in-viscid fluid is negligible. However, the presence of a medium with finite electrical conductivity can substantially delay the decay of the magnetic field [40].
It is obvious that the resulting EM field is a solution of a complicated magneto-hydrodynamic problem [1][2][3][4][5][6]. In fact, the relativistic magneto-hydrodynamic (RMHD) setup is one of the necessary tools in order to describe the hot plasma in the presence of EM fields [41,42]. For this purpose, one needs a numerical code that solves the equations of (1+3) dimensional relativistic magneto-hydrodynamics (RMHD).
Recently, in Refs. [43][44][45][46], some efforts have been made toward both numerical and analytical approaches aimed to solve the RMHD setup, by considering some constraints, specific of high energy heavy ion collisions. In Ref. [46] the main goal was to obtain an analytical solution for a (1+1) dimensional Bjorken flow within ideal transverse RMHD; these authors have neglected (consistently with their hypothesis) the coupling to Maxwell's equations and have analytically solved the energy-momentum conservation equations in a perturbation framework.
Another recent work employs a (1+3) dimensional RMHD code [47]: these authors have used the initial conditions according to the solutions obtained from the Maxwell equations in the early time of the collision: there are, however, many uncertainties in the conditions of the pre-equilibrium phase.
In this paper we improve previous research by removing some of the above-mentioned restrictions: in particular we simulate (1 + 1) dimensional ideal magneto-hydrodynamics in the Bjorken scenario to determine the effect of the magnetic field on the behavior of an inviscid fluid. Here, we consider the combination of non relativistic hydrodynamic equations with Maxwell equations and solve numerically in (1+1) dimensions a set of coupled MHD equations. This improves some previous, analogous, work, where the coupling between Maxwell equations and conservation equations has been neglected or treated perturbatively. For the purpose of numerical calculations, we have supplemented a relatively simple code which incorporates the contribution of a coupled electromagnetic field in (1+1) dimensions. One important novelty is that we use the boundary conditions at late time (τ → ∞): indeed the late-time dynamics has been governed by ideal hydrodynamics and is known, while the early-time conditions are unknown. In order to check our code, we compare our results with the analytical solutions of Ref. [46]. We find, indeed, that their results can be recovered by the numerical solutions, at least in (1+1)-D transverse evolution.
The paper is organized as follows: in Sect. 2, we introduce the ideal relativistic magneto-hydrodynamic equations in their most general form, considering them in the case of a plasma with infinite electrical conductivity. In the end of this section we restrict the formalism to the case of nonrelativistic transverse flow, which will be employed in the subsequent calculations. In Sect. 3 we present our numerical procedure with details in the setup; results obtained with the spatial initial condition are shown in Sect. 4. Finally, we summarize our conclusions and present a possible subsequent outlook in the last section.

Ideal transverse MHD setup in (1+1D) expansion
We consider the relativistic magneto-hydrodynamic (RMHD) framework, in order to describe the interaction of matter and electromagnetic fields in quark-gluon plasmas [41,42]. For the sake of simplicity, we assume an ideal relativistic plasma with massless particles and infinite electrical conductivity. In addition, the fluid is considered to be ultra relativistic, thus implying that the rest mass contributions to the equation of state (EOS) can be neglected, and the pressure is simply proportional to the energy density: P = c 2 s = 1 3 where c s = 1 3 is the speed of sound. For an ideal fluid with infinite electrical conductivity, the equations of RMHD can be written in the form of the covariant conservation laws, where and Here F μν is the dual tensor of electromagnetic field. and P are energy density and pressure, respectively. b μ is the magnetic field four vector in the local rest-frame of the fluid, which is related in the standard way to the one measured in the lab-frame. In the present paper we assume a fluid with infinite electrical conductivity, so the electric field four vector in the local rest-frame equals zero (e μ = 0). Besides, the single fluid four velocity u μ (u μ u μ = −1) is defined as follows: In Eqs. (1) and (2) the covariant derivative is given by where i jk are the Christoffel symbols and g i j is the metric tensor. It is more convenient to work with Milne coordinates rather than the standard Cartesian coordinates for a longitudinally boost-invariant flow: Here, the metric is given by Working in Milne coordinates, one can easily obtain the Christoffel symbols: the only non-zero ones being τ ηη = τ and η τ η = 1/τ . Then four distinct conservation equations can easily be derived from d μ T μν = 0 in the Milne coordinate system. They are given by In contrast with the energy-momentum tensor T μν , the dual electromagnetic tensor F * μν is antisymmetric; hence the homogeneous Maxwell equation, d μ F * μν = 0, leads to the following equations: In order to simplify the problem, we assume that the magnetic field is perpendicular to the reaction plane, pointing along the y direction in an inviscid fluid with infinite electrical conductivity, following the Bjorken expansion along the z direction and moving, in the transverse plane, only in the x direction. The boost invariance of the Bjorken expansion allows us to restrict the discussion to the z = 0 plane, where symmetry reasons impose u z = 0. Then In our setup, the energy-momentum and dual electromagnetic tensors are given by When Eqs. (22)-(23) are plugged into Eqs. (13)- (20), one obtains the following coupled equations: In the following, we will assume that the transverse velocity u x is non-relativistic, so we will keep only first-order terms in u x . Since we consider (1+1D) flow, all thermodynamical variables depend only on τ, x coordinates. Applying the definition of u 2 = −1 one finds that u μ = (1, u x , 0, 0) andγ → 1. Using all the above assumptions the set of equations (24)-(30) reduce to As one expects in ideal MHD, the energy conservation equation (31 does not include the B field. In the next section we present a numerical method to solve the above coupled equations (31)-(33) simultaneously.

Numerical calculation
In this section we will solve the coupled non-relativistic hydrodynamic and Maxwell equations, which are summarized in Eqs. (31)- (33). The solutions of the three coupled differential equations will be obtained by using the numerical method of lines (MOL). This method is a technique for solving partial differential equations (PDE) by discretizing one variable in one of the two dimensions and then by integrating the semi-discrete problem as a system of ordinary differential equations (ODE). Here we discretize the partial derivatives with respect to the space variables and obtain a system of ODEs in the time variable: then the initial value software Mathematica has been used to solve this ODE system. It is necessary that the partial differential equation problem be well posed as an initial value problem in at least one dimension, since these are the conditions for an appropriate use of the employed ODEs integrators. Hence we discretize the coordinate x with N (N even) uniformly spaced grid points x i = (i − 1)h, x N +1 = π, i = 1, 2, . . . , N and h = π/N . We use a second-order finite difference formula for the first derivative in x. In this configuration v i (τ ) indicates v(τ, x i ). In Fig. 1, the lines along which the discrete quantities v i (τ ) are defined, are shown. Using the second-order difference approximation for the first derivative in x results in After substituting the first derivatives with respect to x for In Ref. [46] the profile of the magnetic field has been defined by where n is a negative value which governs the decay of magnetic field with increasing time. The Fourier expansion of the above square magnetic field is approximated as b y (τ, x) 2 = B 2 c τ n (0.28 + 0.44 cos x + 0.21 cos 2x + 0.06 cos 3x + 0.01 cos 4x) (36) and Fig. 2 shows a comparison between the Fourier cosine series and the Gaussian distribution at the late time τ 1 = 20 fm. 1 Due to the oscillatory property of the cosine function, the solutions are valid only in the region −π < x < π. The spatial width of the magnetic field depends on the impact parameter of the considered peripheral collision. Following the method which has been presented in Ref. [46] one can obtain analytical solutions for transverse velocity v(τ, x) and energy density ratio These are 3N +3 first-order, coupled ODEs with boundary conditions:

Numerical results of MHD
We shall now show the results obtained by numerically solving the above outlined system of equations. To resume the procedure, we recall that with the method of lines we have fixed discrete values for the variable x (3N + 3 values) and defined derivatives with respect to x via the second-order difference method; then the original set of equations reduces to a system of 3N + 3 coupled ordinary differential equations for the quantities i (τ ), b i (τ ), v i (τ ).
The above-mentioned ODEs have been solved by using an ODE solver of Mathematica, with respect to the time variable. The initial boundary conditions for the ODE integrator have been fixed at the late time τ 1 = 20 fm and derived from the analytical solutions of Ref. [46].
This procedure allowed us to solve the (R)MHD equations in 1+1 dimensions (Eqs. (31)- (33)) and to obtain the space-time evolution of the magnetic field, the energy density and the velocity of plasma, v(τ, x), (τ, x) and b(τ, x). Our results for these functions are presented in following figures, where they are plotted versus x at fixed τ or versus τ at fixed x, for three different values of n ( n < −1). Figures 5 and 6 show the variation of the fluid velocity in terms of either τ or x with different values of n. Figure 5 shows that |v| at fixed x is large at early times end becomes small in late times. This behavior is strikingly evident for the smallest |n| employed here. The transverse velocity v in terms of x at the fixed time τ = 1 fm has been plotted in Fig. 6. v increases from x = 0 to a maximum at intermediate x and gradually decrease with x. Besides, one can see that when the |n| increases, the v at fixed τ becomes smaller due to the faster decay of the magnetic field. Figures 7 and 8 show the energy density in terms of x at fixed τ or in terms of τ at fixed x, for three different values of n ( n < −1). Figure 7 indicates that grows from x = 0 up to some intermediate value of x, where it seems to saturate: the increase is more rapid when n > −2. The behavior of the energy density as a function of τ is monotonically decreasing; notice the different scale between Figs. 7 and 8.  Finally, in Figs. 9 and 10 the magnetic field b y (τ, x) is plotted versus x at fixed time (τ = 1) or in terms of τ at fixed x. Figure 9 shows that b y (τ, x) becomes small for large x. Besides, when |n| < 2 the magnetic field is quite larger, in the central area, than for the other values of n.
Next we wish to validate our numerical work by comparing it with the approximate analytical calculation of Ref. [46]: as already stated these authors assume an external, timedecreasing magnetic field B with a Gaussian distribution in x.
In Fig. 11. we show the profile of the magnetic field for different values of n, for the fixed time τ = 1: the black curves correspond to our solutions and the blue curves correspond to the approximate analytical solution of Ref. [46], where the magnetic field is considered as the one given in Eq. (35). At τ = 1 obviously the profiles of magnetic field Eq. (35) is independent of the value of n, while in the present work we obtain different results, as already seen in Fig. 9. The case n = −2 corresponds to ideal magneto-hydrodynamics and is referred to as the "frozen-flux condition", which stems from the Maxwell equations with conservation of the entropydensity current: in this case the analytical solution coincides with the one obtained in the present calculation. Instead for n = −7/3 the magnetic field decreases by nearly a factor 0.5 and for n = −5/3 it increases by a factor 1.5 relative to case n = −2. Figure 12 shows the ratio b(τ, 0) 2 / (τ, 0) as a function of τ at x = 0, for different values of n: it is seen that all the analytical solutions of Ref. [46] (blue lines) converge to the value 0.1 at τ = 1 for any value of n, as expected. Our results, instead, reach the same limit only for the value n = 2, while for n = −7/3 (n = −5/3) the ratio is typically smaller (larger). This shows the effect of the coupling between Maxwell's and conservation equations. Moreover, for n < −1, the ratio decreases with increasing time: this implies that the energy density of the magnetic field decays much faster than the fluid energy density in relativistic heavy ion collisions. Finally, in order to show the influence of the MHD equations on the modification of velocity, magnetic field and energy density we consider the specific case n = −4/5, which corresponds to a weaker dependence on time of the magnetic field, Eq. (35). This is illustrated in Figs. 12 and 13, where again we compare our solutions with the approximated analytical ones of Ref. [46]. The latter implies that v(τ, x) = 0 for n = −1; hence one may expect a change in the direction of the transverse velocity. Figure 13a, b show the transverse velocity results from the numerical solutions of the present work and from the analytical solutions, respectively, at different times, for n = −4/5. From Fig. 13b one finds that the magnitude of the transverse velocity decreases with respect to the proper time τ 0 = 1 fm, as expected, and the velocity profile has a similar shape compared with the case n < −1, but the direction becomes negative. In addition (notice the small numbers in the vertical scale) it is nearly zero, since in the approximate analytical solution the fluid velocity is only modified by the spatial gradient of the external magnetic field. On the contrary, Fig. 13a shows that the direction of the fluid velocity is positive and decreases with time, until τ = 10 fm, where the sign changes; moreover, its magnitude is much larger than the analytical solution at early times. As a conclusion we can state that, in the analytical solution, the transverse flow led by a Gaussian magnetic field points outward for n < −1 and inward for n > −1, while the results of the present work, where the MHD equations are solved numerically, for both cases (n < −1 and n > −1) the transverse flow points outward in early time, though they have opposite direction at late time. Figure 14a shows the behavior of b(τ, 0) 2 / (τ, 0) versus x, for different proper times: it decreases with time from the value 2.5 at proper time τ 0 = 1 fm to the value 0.6 at the late time τ 1 = 10 fm, similarly to the case n < −1 (see Fig. 12). We also plot the analytical solution for b(τ, 0) 2 / (τ, 0) in Fig. 14b: in this case the considered ratio increases with time from the value τ 0 = 0.1 at proper time τ 0 = 1 fm to the value 0.6 at the late time τ 1 = 10 fm. This again stresses the importance of the coupling between Maxwell's and the conservation equations.

Conclusion and outlook
In this work, we present a numerical method for the solution of coupled non-relativistic hydrodynamic equations and Maxwell's equations, i.e., non-relativistic magnetohydrodynamics (MHD), which has recently become of growing interest for the study of relativistic heavy ion collisions. By solving the coupled conservation and Maxwell's equations, we obtain numerical results for the fluid velocity, the energy density and the magnetic field.
We work in the (1+1D) dimensional MHD model where the transverse magnetic field, fluid velocity, and energy density are considered as a function of one spatial (x) and one temporal (τ ) variable; besides, the magnetic field points along the orthogonal y direction. In our setup, the medium is boostinvariant along the z direction. It turns out that the transverse velocity is rather small at all times and for different parameterizations of the initial magnetic field: hence we treated the transverse flow in the non-relativistic approximation.
The core of our method is twofold: (1) the adoption of a discretized spatial variable, in terms of which the derivatives are expressed with the method of the second-order finite difference formula, (2) the adoption of suitable boundary conditions for the numerical solution of the resulting system of ordinary differential equations in the time variable.
Although it is well known that during relativistic heavy ion collisions intense magnetic fields are developed, their knowledge in the initial times of the collision, where the QGP has been formed, could not be used here as the desired boundary conditions. Indeed for the purpose of comparing with the results of Ref. [46], it turned out more convenient to introduce initial conditions at a late time and solve numerically the coupled equations inversely in time We found it appropriate to assume for the late-time quantities (fluid velocity, energy density and magnetic field) the approximate analytical solutions found by the authors of Ref. [46]: in contrast with our approach, these solutions where found by neglecting the back reaction from the internal fields, i.e., the combination of Maxwell equations and conservation equations was discarded. A weak external, magnetic field is adopted, with Gaussian distribution dependence in space and power-law decay dependence in time. For our method these analytical solutions are appropriate for the late-time boundary conditions, since in the final stage of the plasma evolution, the magnetic field is indeed quite small and the back reaction from the internal fields can be safely neglected.
After presenting our results for the fluid velocity, the energy density and the magnetic field as functions of both space and time, we have validated our numerical calculation by making a comparison with the approximate analytical solutions of Ref. [46]. As expected the two approaches give similar results at late times as well as for specific choices of the time evolution of their magnetic field (which is governed by the parameter n), but in other conditions the coupling between Maxwell and conservation (RMHD) equations, here taken into account, appears to be quite relevant.
In particular we notice that the transverse velocities have the same direction in early and late time for the case n < −1 (faster decay of the magnetic field), while for the case n > −1 the transverse velocities appear to change sign in late time. Hence, the transverse flow propagates on the same direction, for any value of n only in the early stages of the collision. It should be noticed, however, that a weak decay of the magnetic field is probably less realistic than a strong one.
According to the estimated conditions of heavy ion collision experiment at RHIC, one find b 2 / = 0.17 − 0.68 at τ = 0.6 fm. As a result, for the validity of the weak-field expansion, in Ref. [46] b 2 / = 0.1 was chosen at the proper time τ 0 = 1. We find that only for the case n = −2, b 2 / converges to this value at the proper time τ 0 = 1. For the casen < −2, b 2 / is smaller than 0.1 and decreases when n becomes smaller. For the case n > −2, b 2 / is bigger than 0.1 and increases when n becomes larger. For both cases n < −1 and n > −1, b(τ, x) 2 / (τ, x) decreases with time at fixed x. This is a preliminary result of our approach, which is potentially interesting and deserves further investigation also from experimental point of view.
As a final remark and outlook for future developments we wish to consider the possibility of extending the present calculation to the case of (1+2D) dimensions, taking into account both transverse dimensions. This will also allow one to investigate differences in the azimuthal distribution of the velocities and hence the so-called elliptic flow. It is well known that this is one of the crucial characteristics of the deconfined plasma. At present there is an interesting debate about the influence of the magnetic field on the v 2 coefficient of the elliptic flow (see for example Ref. [48,49]) and the role of the magnetic field in this connection still deserves additional efforts.