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

TL;DR

This paper models the evolution of magnetic fields and transverse flow in quark-gluon plasma using a non-relativistic, ideal magneto-hydrodynamic approach within the Bjorken scenario, focusing on magnetic effects at mid-rapidity.

Contribution

It introduces a non-relativistic numerical method for ideal MHD in (1+1)D transverse flow in the Bjorken scenario, analyzing magnetic effects on QGP evolution.

Findings

Magnetic fields influence the energy density and velocity profiles of QGP.

Magnetic effects are significant at mid-rapidity in ideal conditions.

The study provides analytical solutions for magnetic field and flow evolution.

Abstract

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.

Equations

$$d_{\mu }{T}^{\mu \nu }$$

$$d_{\mu }{F}^{\star \mu \nu }$$

$$T^{\mu \nu }$$

$$T_{matter}^{\mu \nu }$$

$$T_{EM}^{\mu \nu }$$

$$F^{\star \mu \nu }$$

$$b^{\mu }=F^{\star \mu \nu }{u}_{\nu },(b^{\mu }{u}_{\mu }=0),b^2=b^{\mu }{b}_{\mu }$$

$$u_{\mu }=\gamma (1,\vec{v}),\gamma =\frac{1}{\sqrt{1-v^2}}$$

$$d_{\mu }{A}^{\nu }$$

$$d_p{A}^{\mu \nu }$$

$${\Gamma }_{jk}^i=\frac{1}{2}{g}^{im}\left(\frac{\partial g_{mj}}{\partial x^k}+\frac{\partial g_{mk}}{\partial x^j}-\frac{\partial g_{jk}}{\partial x^m}\right)$$

$$(\tau ,x,y,\eta )$$

$$g^{\mu \nu }=diag(-1,1,1,1/{\tau }^2),g_{\mu \nu }=diag(-1,1,1,{\tau }^2)$$

$${\partial }_{\tau }{T}^{\tau \tau }+{\partial }_x{T}^{x\tau }+{\partial }_y{T}^{y\tau }+{\partial }_{\eta }{T}^{\eta \tau }+\tau T^{\eta \eta }+\frac{1}{\tau }{T}^{\tau \tau }=0,$$

$${\partial }_{\tau }{T}^{\tau x}+{\partial }_x{T}^{xx}+{\partial }_y{T}^{yx}+{\partial }_{\eta }{T}^{\eta x}+\frac{1}{\tau }{T}^{\tau x}=0,$$

$${\partial }_{\tau }{T}^{\tau y}+{\partial }_x{T}^{xy}+{\partial }_y{T}^{yy}+{\partial }_{\eta }{T}^{\eta y}+\frac{1}{\tau }{T}^{\tau y}=0,$$

$${\partial }_{\tau }{T}^{\tau \eta }+{\partial }_x{T}^{x\eta }+{\partial }_y{T}^{y\eta }+{\partial }_{\eta }{T}^{\eta \eta }+\frac{3}{\tau }{T}^{\tau \eta }=0$$

$${\partial }_x{F}^{\ast x\tau }+{\partial }_y{F}^{\ast y\tau }+{\partial }_{\eta }{F}^{\ast \eta \tau }$$

$${\partial }_{\tau }{F}^{\ast \tau x}+{\partial }_y{F}^{\ast yx}+{\partial }_{\eta }{F}^{\ast \eta x}+\frac{1}{\tau }{F}^{\ast \tau x}$$

$${\partial }_{\tau }{F}^{\ast \tau y}+{\partial }_x{F}^{\ast xy}+{\partial }_{\eta }{F}^{\ast \eta y}+\frac{1}{\tau }{F}^{\ast \tau y}$$

$${\partial }_{\tau }{F}^{\ast \tau \eta }+{\partial }_x{F}^{\ast x\eta }+{\partial }_y{F}^{\ast y\eta }+\frac{1}{\tau }{F}^{\ast \tau \eta }$$

$$u_{\mu }=\tilde{\gamma }(1,u_x,0,0),b_{\mu }=(0,0,b_y,0),e_{\mu }=(0,0,0,0).$$

$$T^{\mu \nu }=\left(\begin{array}{cccc}(ϵ+P+b^2){\tilde{\gamma }}^2-P-\frac{b^2}{2}& -(ϵ+P+b^2){\tilde{\gamma }}^2{u}_x& 0& 0\\ -(ϵ+P+b^2){\tilde{\gamma }}^2{u}_x& (ϵ+P+b^2){\tilde{\gamma }}^2{u}_x^2+P+\frac{b^2}{2}& 0& 0\\ 0& 0& P-\frac{b^2}{2}& 0\\ 0& 0& 0& (P+\frac{b^2}{2})(\frac{1}{{\tau }^2})\end{array}\right)$$

$$F^{\ast \mu \nu }$$

$${\partial }_{\tau }\left[(ϵ+P+b^2){\tilde{\gamma }}^2-P-\frac{b^2}{2}\right]+{\partial }_x\left[-(ϵ+P+b^2){\tilde{\gamma }}^2{u}_x\right]$$

$$+\frac{(ϵ+P+b^2){\tilde{\gamma }}^2}{\tau }=0,$$

$${\partial }_{\tau }\left[-(ϵ+P+b^2){\tilde{\gamma }}^2{u}_x\right]+{\partial }_x\left[(ϵ+P+b^2){\tilde{\gamma }}^2{u}_x^2+P+\frac{b^2}{2}\right]$$

$$-\frac{(ϵ+P+b^2){\tilde{\gamma }}^2{u}_x}{\tau }=0,$$

$${\partial }_y\left(P-\frac{b^2}{2}\right)=0,$$

$${\partial }_{\eta }\left((P+\frac{b^2}{2})(\frac{1}{{\tau }^2})\right)=0,$$

$${\partial }_y(\tilde{\gamma }{b}_y)=0,$$

$${\partial }_y(\tilde{\gamma }{u}_x{b}_y)=0.$$

$$\left[-\tilde{\gamma }{\partial }_{\tau }+\tilde{\gamma }{u}_x{\partial }_x\right]b_y+b_y\left({\partial }_x(\tilde{\gamma }{u}_x)-{\partial }_{\tau }(\tilde{\gamma })-\frac{\tilde{\gamma }}{\tau }\right)=0$$

$$(-{\partial }_{\tau }+u_x{\partial }_x)b_y+(\frac{\partial u_x}{\partial x}-\frac{1}{\tau })b_y=0,$$

$$-{\partial }_{\tau }ϵ+u_x{\partial }_xϵ+(1+c_s^2)ϵ(\frac{\partial u_x}{\partial x}-\frac{1}{\tau })=0,$$

$$\displaystyle-\Big{(}(1+c_{s}^{2})\epsilon+b^{2}\Big{)}\partial_{\tau}u_{x}+c_{s}^{2}\partial_{x}\epsilon+b_{y}\partial_{x}b_{y}+\left(\frac{c_{s}^{2}(1+c_{s}^{2})\epsilon+b^{2}}{\tau}\right)u_{x}=0.$$

$$\frac{d{v}_i(\tau )}{dx}=\frac{-3v_i(\tau )+4v_{i+1}(\tau )-v_{i+2}(\tau )}{2h},i=1,....,N+1.$$

$$b_y(\tau ,x)\hat{y}=B_c{\tau }^{n/2}{e}^{-x^2/2}\hat{y}$$

$$B_y(\tau ,x{)}^2=B_c^2{\tau }^n(0.28+0.44\cos x+0.21\cos 2x+0.06\cos 3x+0.01\cos 4x)$$

$$\frac{\partial {ϵ}_i(\tau )}{\partial \tau }$$

$$\frac{\partial b_i(\tau )}{\partial \tau }$$

Full text

1 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 quark gluon plasma(QGP), for a very short time ($\sim$ 1fm/c) after the initial hard parton collisions of nuclei. The hydrodynamic approach has given one of the best description 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]-[6].

Recently, it has been shown that in the peripheral AA-collisions such as Pb-Pb at the center of mass energy $\sqrt{s}=2.76$ TeV and Au-Au at the center of mass energy $\sqrt{s}=200$ GeV a huge magnetic field is created, of the order of eB$\sim 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 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]-[25]. 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 [26]-[32]. In ref [26]-[31], 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 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 $\sqrt{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. [26]-[32]) the electromagnetic field is derived from 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. [29], and it has been concluded that the induced electric current in the plasma fireball cannot generate a classical electromagnetic field. In [34] 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 that the velocity of charged particles $\vec{v}$ is smaller than the velocity of the expanding plasma $\vec{u}$.

Other works obtained a series of preliminary results, by estimating the significance of strong EM fields on the QGP medium [35]-[39]. In most of them, it has been assumed the Maxwell equations decouple from the time evolution of the QGP, and then the evolution of the EM fields as well as 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 [39].

It is obvious that the resulting EM field is a solution of a complicated magneto-hydrodynamic problem [1]-[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 [40]-[41]. For this purpose, one needs a numerical code that solves the equations of (1+3) dimensional relativistic magneto-hydrodynamics (RMHD).

Recently, in refs. [42]-[45] 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. [45] 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 [46]: these authors have used the initial conditions according to the solutions obtained from 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 researches 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 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 ($\tau\rightarrow\infty$): 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. [45]. 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 Section $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 Section $3$ we present our numerical procedure with details in the setup; results obtained with the spatial initial condition are shown in Section 4. Finally, we summarize our conclusions and possible subsequent outlook in the last Section.

2 Ideal transverse MHD setup in (1+1)D expansion

We consider the relativistic magnetohydrodynamic (RMHD) framework, in order to describe the interaction of matter and electromagnetic fields in quark-gluon plasmas [40]-[41]. 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) have been neglected, and the pressure is simply proportional to the energy density: $P=c_{s}^{2}\epsilon=\frac{1}{3}\epsilon$ where $c_{s}=\sqrt{\frac{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

[TABLE]

where

[TABLE]

and

[TABLE]

Here $F^{\star\mu\nu}$ is the dual tensor of electromagnetic field. $\epsilon$ and $P$ are energy density and pressure respectively. $b^{\mu}$ 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 to zero ($e^{\mu}=0$). Besides, the single fluid four velocity $u_{\mu}$ ($u_{\mu}u^{\mu}=-1$) is defined as follow:

[TABLE]

In eqs.(1) and (2) the covariant derivative is given by:

[TABLE]

where $\Gamma^{i}_{jk}$ are the Christoffel symbols

[TABLE]

and $g_{ij}$ the metric tensor.

It is more convenient to work with Milne coordinates rather than the standard Cartesian coordinates for a longitudinally boost-invariant flow:

[TABLE]

Here, the metric is given by:

[TABLE]

Working in Milne coordinates, one can easily obtain the Christoffel symbols: the only non-zero ones being: $\Gamma^{\tau}_{\eta\eta}=\tau$ and $\Gamma^{\eta}_{\tau\eta}=1/\tau$. Then, four distinct conservation equations can be easily derived from $d_{\mu}T^{\mu\nu}=0$ in the Milne coordinate system. They are given by:

[TABLE]

In contrast with the energy momentum tensor $T^{\mu\nu}$, the dual electromagnetic tensor $F^{*\mu\nu}$ is antisymmetric; hence the homogeneous Maxwell equation, $d_{\mu}F^{*\mu\nu}=0$, leads to the following equations:

[TABLE]

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 $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:

[TABLE]

Where $\tilde{\gamma}=\frac{1}{\sqrt{1-u_{x}^{2}}}$, and $u^{\mu}b_{\mu}=0,\ u^{\mu}u_{\mu}=-1$ are satisfied.

In our setup, the energy momentum and dual electromagnetic tensors are given by:

[TABLE]

[TABLE]

When the Eqs.(22)–(23) are plugged into Eqs.(13)–(20), one obtains the following coupled equations:

[TABLE]

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+1)D flow, all thermodynamical variables depend only on $\tau,x$ coordinates. Applying the definition of $u^{2}=-1$ one finds that $u_{\mu}=(1,u_{x},0,0)$ and $\tilde{\gamma}\rightarrow 1$. Using all the above assumptions the set of equations (24)–(30) reduce to:

[TABLE]

As one expects in the 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.

3 Numerical calculation

In this section we will solve the coupled 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}=\pi,\ i=1,2,…,N$ and $h=\pi/N$. We use a second-order finite difference formula for the first derivative in $x$. In this configuration $v_{i}(\tau)$ indicates $v(\tau,x_{i})$. In Fig. 1, the lines along which the discrete quantities $v_{i}(\tau)$ are defined, are shown. Using the second-order difference approximation for the first derivative in $x$ results in

[TABLE]

After substituting the first derivatives with respect to $x$ for the $v_{i}(\tau),\ \epsilon_{i}(\tau),\ b_{i}(\tau)$, in equations (31)–(33) one is left with a set of coupled ODEs: in order to numerically solve these equations the crucial point remains the definition of the boundary conditions.

3.1 Boundary conditions

In our work we consider the boundary conditions at late-time. Indeed the issue under investigation does not allow to have precise information on the early-times conditions. Also the determination of the proper boundary conditions from analytical solutions of the equations(31)–(33) at late-time ($\tau\to\infty$) is still very difficult, due to the coupling between conservation and Maxwell equations. Hence we derived such boundary conditions at late time from the analytical solutions of ref[45].These authors have investigated the magneto hydrodynamics in the presence of an external magnetic field, which follows the power-law decay in proper time and has spatial inhomogeneity characterized by a Gaussian distribution in transverse coordinates. For simplicity they have neglected the coupling of Maxwell’s equations and conservation equations and solved the conservation equations perturbatively and analytically.

In ref. [45] the profile of the magnetic field has been defined by:

[TABLE]

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

[TABLE]

and Fig. 2 shows a comparison between the Fourier cosine series and the Gaussian distribution at the late time $\tau_{1}=20$ fm.111In the ratio illustrated in Fig. 2 the $\tau$ dependence (and hence the precise value of $\tau$)n is irrelevant, since it cancels in the ratio; however in the following this value will be chosen to fix the initial late conditions. Due to the oscillatory property of the cosine function, the solutions are valid only in the region $-\pi<x<\pi$. 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. [45] one can obtain analytical solutions for transverse velocity $v(\tau,x)$ and energy density ratio $\epsilon(\tau,x)/\epsilon_{c}$ at the assumed late time [$\epsilon_{c}$ is the initial energy density of the medium at time $\tau_{0}$]. They are shown in Figs. 3 and 4, for $\tau_{1}=20$ fm and 3 different values of $n$; the value $B_{c}^{2}/\epsilon_{c}=0.1$ has been assumed.

Then we consider the following family of equations:

[TABLE]

These are $3N+3$ first order, coupled ODEs with boundary conditions: $\epsilon_{i}(\tau_{1})\equiv\epsilon_{i}(\tau_{1},x_{i}),\ b_{i}(\tau_{1})\equiv b_{i}(\tau_{1},x_{i}),\ v_{i}(\tau_{1})\equiv v_{i}(\tau_{1},x_{i}),i=1,...,N+1$ and $v_{1}(\tau)=\ v_{N+1}(\tau)=0$. These functions are obtained from analytical solutions of ref. [45] and shown in Figs. 2-4.

4 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 remind 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 $\epsilon_{i}(\tau),\ b_{i}(\tau),\ v_{i}(\tau)$.

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 $\tau_{1}=20$ fm and derived from the analytical solutions of ref. [45].

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(\tau,x),\epsilon(\tau,x)$ and $b(\tau,x)$. Our results for these functions are presented in following figures, where they are plotted versus $x$ at fixed $\tau$ or versus $\tau$ at fixed $x$, for three different values of $n$ ( $n<-1$).

Figs. 5 and 6 show the variation of the fluid velocity in terms of either $\tau$ or $x$ with different values of $n$. Fig. 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 $\tau=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 $\tau$ becomes smaller due to the faster decay of the magnetic field.

Figs. 7 and 8 show the energy density in terms of $x$ at fixed $\tau$ or in terms of $\tau$ at fixed $x$, for three different values of $n$ ( $n<-1$). Fig. 7 indicates that $\epsilon$ 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 $\tau$ is monotonically decreasing; notice the different scale between Fig. 7 and Fig. 8.

Finally in Figs. 9 and 10 the magnetic field $b_{y}(\tau,x)$ is plotted versus $x$ at fixed time ($\tau=1$) or in terms of $\tau$ at fixed $x$. Fig. 9 shows that $b_{y}(\tau,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. [45]: as already stated these authors assume an external, time-decreasing 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 $\tau=1$: the black curves correspond to our solutions and the blue curves correspond to the approximate analytical solution of Ref. [45], where the magnetic field is considered as the one given in Eq. (35). At $\tau=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 magnetohydrodynamics and is referred to as the ”frozen-flux condition”, which stems from Maxwell’s equations with conservation of the entropy-density 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$.

Fig.12 shows the ratio $b(\tau,0)^{2}/\epsilon(\tau,0)$ as a function of $\tau$ at $x=0$, for different values of $n$: it is seen that all the analytical solutions of Ref. [45] (blue lines) converge to the value 0.1 at $\tau=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 RMHD 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. [45]. The latter implies that $v(\tau,x)=0$ for $n=-1$ , hence one may expect a change in the direction of the transverse velocity.

Fig. 13 (a) and (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. 13 (b) one finds that the magnitude of the transverse velocity decreases with respect to the proper time $\tau_{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. 13 (a) shows that the direction of fluid velocity is positive and decreases with time, until $\tau=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 RMHD 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.

Fig. 14 (a) shows the behavior of $b(\tau,0)^{2}/\epsilon(\tau,0)$ versus $x$, for different proper times: it decreases with time from the value 2.5 at proper time $\tau_{0}=1$ fm to the value 0.6 at the late time $\tau_{1}=10$ fm, similarly to the case $n<-1$ (see Fig. 12). We also plot the analytical solution for $b(\tau,0)^{2}/\epsilon(\tau,0)$ in Fig. 14 (b): in this case the considered ratio increases with time from the value $\tau_{0}=0.1$ at proper time $\tau_{0}=1$ fm to the value 0.6 at the late time $\tau_{1}=10$ fm. This again stresses the importance of the coupling between Maxwell’s and conservation equations.

5 Conclusion and outlook

In this work, we present a new numerical method for the solution of coupled relativistic hydrodynamic equations and Maxwell’s equations, i.e., relativistic magnetohydrodynamics (RMHD), 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+1 dimensional RMHD model where the transverse magnetic field, fluid velocity, and energy density are considered as a function of one spatial ($x$) and one temporal ($\tau$) variable; besides, the magnetic field points along the orthogonal $y$ direction. In our setup, the medium is boost-invariant 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: i) 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, ii) the adoption of precise and realistic boundary conditions for the numerical solution of the resulting system of ordinary differential equations in the time variable.

Although it is 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, is very poor and could not be used as the desired boundary conditions: hence we introduce initial conditions at 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. [45]: in contrast with our approach, these solutions where found by neglecting the coupling between the dynamical evolution of the magnetic field and the one of the energy density and velocity of the fluid. A weak external, uncoupled 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 its coupling 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. [45]. 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}/\epsilon=0.17-0.68$ at $\tau=0.6$ fm. As a result, for the validity of the weak-field expansion, in Ref. [45] $b^{2}/\epsilon=0.1$ was chosen at the proper time $\tau_{0}=1$. We find that only for the case $n=-2$, $b^{2}/\epsilon$ converges to this value at the proper time $\tau_{0}=1$. For the case$n<-2$, $b^{2}/\epsilon$ is smaller than 0.1 and decreases when $n$ becomes smaller. For the case $n>-2$, $b^{2}/\epsilon$ is bigger than 0.1 and increases when $n$ becomes larger. For both cases $n<-1$ and $n>-1$, $b(\tau,x)^{2}/\epsilon(\tau,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 2+1 dimensions, taking into account both transverse dimensions. This will also allow 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. [47, Roy]) and the role of the magnetic field in this connection still deserves additional efforts.