Transverse expansion of hot magnetized Bjorken flow in heavy ion collisions

We argue that the existence of an inhomogeneous external magnetic field can lead to radial flow in transverse plane. Our aim is to show how the introduction of a magnetic field generalizes the Bjorken flow. We investigate the effect of an inhomogeneous weak external magnetic field on the transverse expansion of in-viscid fluid created in high energy nuclear collisions. In order to simplify our calculation and compare with Gubser model, we consider the fluid under investigation to be produced in central collisions, at small impact parameter; azimuthal symmetry has been considered. In our model, we assume an inhomogeneous external magnetic field following the power-law decay in proper time and having radial inhomogeneity perpendicular to the radial velocity of the in-viscid fluid in the transverse plane; then the space time evolution of the transverse expansion of the fluid is obtained. We also show how the existence of an inhomogeneous external magnetic field modifies the energy density. Finally we use the solutions for the transverse velocity and energy density in the presence of a weak magnetic field, to estimate the transverse momentum spectrum of protons and pions emerging from the Magneto-hydrodynamic solutions.


Introduction
Collisions of two heavy nuclei at high energy produce a hot and dense fireball. Quarks and gluons could reach the deconfined state, called quark gluon plasma (QGP), in a very short time (∼1 fm/c) after the initial hard parton collisions of nuclei. A very handy model which describes the typical motion of partons after collision is the Bjorken flow model [1]. Based on some assumptions such as boost invariance along beam line, translation and rotation invariance in the transverse plane, one can show that all quantities of interest only depend on the proper time τ and not on the transa e-mail: afarzaneh@hsu.ac.ir verse (x ⊥ , φ) coordinates, nor on the rapidity η. Using the above assumptions, together with invariance under reflection η → −η, one can determine the four-velocity profile. The four-velocity is u μ = (1, 0, 0, 0) in the (τ, x ⊥ , φ, η) coordinate system. Besides, it is straightforward to show that the energy density decays as τ −4/3 in the local rest frame if the medium is equilibrated and the equation of state of the medium is p = /3.
Based on the size of the colliding nuclei, one realizes that assuming translational invariance in the transverse plane is not realistic [2]. Using the Bjorken model one often assumes that in medium the radial flow (u ⊥ ) is zero. However, this assumption is not correct even for central collisions, and it might mislead the subsequent hydrodynamical flow, on which much of heavy-ions phenomenology depends.
The aim of our work is to generalize the Bjorken model by considering an inhomogeneous external magnetic field acting on the medium. We show that the presence of the magnetic field leads to non-zero radial flow. In order to simplify our calculation, we consider central heavy ions collisions. We still consider rotational symmetry around beam line, as well as boost invariance along the beam line. However, we assume that translational invariance in the transverse plane is broken by the magnetic field. Then we obtain a four-velocity profile which has a non zero radial component. In the present study for central collisions (small impact parameter), we provide an analytical solution for the transverse expansion of a hot magnetized plasma, based on perturbation theory.
We concentrate on the special case of a (1 + 2) dimensional, longitudinally boost-invariant fluid expansion as the Bjorken flow; the fluid also radially expands in the transverse plane, under the influence of an inhomogeneous external magnetic field which is transverse to the radial fluid velocity (this proceeds according to the so called transverse MHD).
We consider an inviscid fluid coupled to an external magnetic field. As one expects in central collisions, we assume that the external magnetic fields is small compared to the fluid energy density [3]. Therefore, we can neglect the coupling to the Maxwell's equations and solve the conservation equations perturbatively and analytically [4].
Moreover the presence of external magnetic field may induce internal electromagnetic fields of the fluid. The internal magnetic fields are dictated by Maxwell's equations and one should solve the conservation equations and Maxwell's equations coupled to each other by numerical methods [25]. In this work we neglect the effects of such internal magnetic field. Hence we will consider the system with an inhomogeneous external magnetic field and will investigate the anisotropic transverse flow and the modified energy density of the fluid induced by the external magnetic fields. As in Gubser flow, the finite size of the colliding nuclei leads to non-zero radial velocity (u ⊥ ); we show that the inhomogeneous weak external magnetic field also leads to nonzero radial velocity and can produce modifications on the radial expansion of the plasma in central collisions.
We remind the reader that recently a wide range of studies has shown that relativistic heavy-ion collisions create also huge magnetic field due to the relativistic motion of the colliding heavy ions carrying large positive electric charge [5][6][7][8][9][10][11][12][13][14][15][16]. The interplay of magnetic field and QGP matter has been predicted to lead to a number of interesting phenomena. One can see recent reviews on this topic in Refs. [17][18][19][20] for more details.
Previous theoretical studies show that the strength of the produced magnetic field depends on the center of mass energy ( √ s N N ) of the colliding nuclei, on the impact parameter (b) of the collision, on the electrical and chiral conductivities (σ el , σ χ ) of the medium [6,10,11,15]. Moreover, the magnetic field in central collisions becomes non-zero due to the fluctuating proton position from event to event [3,13]. It has been found that the ratio of magnetic field energy to the fluid energy density (σ = eB 2 /2 ) in central collisions is much smaller than in peripheral collisions [3]. The authors of Ref. [3] computed the fluid energy density and electromagnetic field by using the Monte Carlo Glauber model. The initial energy density for the fluid at proper time τ i = 0.5 fm was fixed to ∼40 GeV/fm 3 . They found σ 1 for most of the events, at the center of the collision zone and for impact parameter b = 0, while for large b as compared to central collisions, σ becomes larger as a result of the increase in magnetic field and decrease in fluid energy density. In a plasma σ = 1 indicates that the effect of magnetic field in the plasma evolution can not be neglected, but it is worth observing that in some situations, even σ ∼ 0.01 may affect the hydrodynamical evolution [3].
We found that this coupling can indeed affect the solutions with respect to the ones of Ref. [4]. In the present work, we show that the perturbative approach of Ref. [4] can be applied to the case of central collision, in order to find analytical solution for the transverse expansion of QGP matter in the presence of an external magnetic field.
The paper is organized as follows. In Sect. 2, we introduce the ideal relativistic magnetohydrodynamic equations in their most general form, considering them in the case of a plasma with infinite electrical conductivity. In Sect. 3 we present our perturbative approach and the analytical solutions we found. Section 4 illustrates and discusses the general results obtained. Section 5 contains a calculation of the transverse momentum spectrum together with a comparison of this quantity with experimental results obtained at RHIC. Conclusions and subsequent outlook can be found in the last section.

Ideal relativistic magneto-hydrodynamic
We deal with the case of an ideal non-resistive plasma, with vanishing electric field in the local rest-frame (e μ = 0), which is embedded in an external magnetic field (b μ ) [31,32]. The energy momentum conservation equations read: where In the above g μν is the metric tensor, and P are the energy density and pressure, respectively. Moreover d μ is the covariant derivative, defined later in Eq. (9). The four velocity is defined as satisfying the condition u μ u μ = −1.
Canonically one takes projections of the equation  For the transverse projection we use the definition μν = g μν + u μ u ν ; then α ν d μ (T μν pl + T μν em ) = 0 gives: Notice that α should be a spacelike index. Moreover

Ideal transverse MHD setup in the transverse expansion
We assume that the medium has a finite transverse size and expands both radially and along the beam axis, the only nonzero components of u μ = (u τ , u ⊥ , 0, 0) being u τ , which describes the boost invariant longitudinal expansion, and u ⊥ , which describes the transverse expansion. For the sake of simplicity we suppose that u φ = 0 because we claimed that the system is rotationally symmetric. It is more convenient to work in Milne coordinates, x m = (τ, x ⊥ , φ, η), such that: Moreover we suppose that the external magnetic field is located in the transverse plane as Fig. 1.
The metric for the coordinates (τ, x ⊥ , φ, η) is parameterized as follows: Correspondingly In this configuration it is found that u τ = −u τ = −u 0 and ∂ τ = −∂ τ . We have to take care of the following covariant derivative (instead of the usual one): where the Cristoffel symbols are defined as follows: Here we frequently take advantage of the following formula: Hence the only non-zero Christoffel symbols, here, The constraint u 2 = u τ u τ + u ⊥ u ⊥ = −u 2 0 + u 2 ⊥ = −1 must be satisfied as well.
We now look for the perturbative solution of the conservation equations in the presence of a weak external inhomogeneous magnetic field pointing along the φ direction in an inviscid fluid with infinite electrical conductivity and obeying Bjorken flow in z -direction. Our setup is given by: where c is the energy density at proper time τ 0 . Then the energy conservation and Euler equations [Eqs. (4), (5)] reduce to two coupled differential equations. Up to O(λ 2 ), they are: The combination of the two above equations yields a partial differential equation depending on u ⊥ and b φ : For b φ = 0, Eq. (19) is a homogeneous partial differential equation, which can be solved by separation of variables. The general solution is where k can be real or imaginary numbers, c k 1,2 and c k 1,2 are integration constants.
For non-vanishing b φ we assume a space-time profile of the magnetic field in central collisions in the form: We see that the magnitude of b φ is zero at x ⊥ = 0. In order to find solutions for transverse velocity u ⊥ and energy density consistently with the assumed magnetic field, we found it convenient to first expand the magnetic field, Eq. (21) into a series of x ⊥ -dependent functions: where k ≥ 1 are now real integers and B 2 k are constants. For simplicity, we have assumed the time dependence of the magnetic field square as τ n with n < 0, which approximately characterizes the decay of the magnetic field in heavy ion collisions. This is our key to convert the solution of the partial differential Eq. (19) into a summation of solutions of ordinary differential equations.
Moreover we replace the solution (20) for the radial velocity u ⊥ (τ, x ⊥ ), which is valid for b φ = 0, with the following Ansatz [4]: It mantains the x ⊥ dependence of Eq. (20), but embodies the τ dependence in the coefficients of the Bessel functions. Note that from the initial condition u ⊥ (τ, Now we can substitute the Eqs. (22) and (23) into Eq. (19) and end up with the equation (at fixed k): Here we can apply separation of variables, thus obtaining the following ordinary differential equation for the function f (kx ⊥ ): Its general solution is given by where 1 F 2 is the hypergeometric function. The first term is a well-defined function, but the second one diverges in x ⊥ = 0 for any n except n = −4/3 which will be considered in details; hence, d 1 must be zero. For two values of the parameter n, n = −1 and n = −4/3, the solution of Eq. (25) takes a simple form: and In order to implement the orthogonal properties of the Bessel functions, for the case n = −4/3 we set d 2 = 0 in Eq. (28). Then we can easily describe the external magnetic field as a series of Bessel functions, by restricting ourselves to the cases n = −1 and n = −4/3.
We write the solution for n = −1 as where the coefficients B 2 k are given by β 1k being the kth zero of J 1 . For n = −4/3 the solution for the magnetic field can be written as where the coefficients B 2 k are given by β 0k being the kth zero of J 0 ; in the above k = β ik /a (i = 0, 1). Finally the coefficients a k (τ ) in Eq. (23) can be obtained by solving the following ordinary differential equation: The analytical solution for n = −1 is while for n = −4/3 the solution is In the above G pq mn is the Meijer function. The transverse velocity then takes the form In order to completely determine the function u ⊥ (τ, x ⊥ ) we must fix the integration constants c k 1 and c k 2 . It is conve-nient to consider the boundary conditions at τ → ∞.
By making late-time expansion of u ⊥ , one finds that u ⊥ takes the asymptotic form f (τ )τ 1/6 where f (τ ) is an oscillatory function. In order to prevent divergencies of the transverse velocity one has to impose that the coefficient of τ 1/6 is equal to zero. The solutions satisfying these boundary condition at τ → ∞ are shown in the following.

Results and discussion
In this section we will present the transverse velocity and energy density numerically obtained from our perturbation approach: this two quantities will help in understanding the space time evolution of the quark-gluon plasma in heavy ion collisions. The typical magnetic field produced in Au-Au peripheral collisions at √ s N N = 200 GeV reaches |eB| ∼ 10m 2 π . The estimate ∼ 5.4 GeV/fm 3 at a proper time of about τ = 1 fm is taken from [2]. By taking m π ≈ 150 MeV and e 2 = 4π/137, one finds B 2 / c ∼ 0.6. This value in central collisions is much smaller than in peripheral collisions; therefore, in our calculations we assumed the even smaller value B 2 c / c = 0.1, which correspond to σ ∼ 0.015. Note that in our calculations any change in the ratio B 2 c / c will only scale the solutions. We will use cylindrical coordinates whose longitudinal component is chosen to be the third component of Cartesian coordinate, e.g., x = (x ⊥ , φ, z).

Numerical solution for the case n = −1
The external magnetic field profile Eq.  Figure 2 shows a comparison between the approximated magnetic field in Bessel series and the assumed magnetic profile Eq. (21). Note that the Fourier expansion matches the assumed magnetic profile in the whole region of x ⊥ , hence the solutions for the radial velocity and the energy density are valid in the entire region x ⊥ ∈ (0, ∞).
is displayed, at either fixed τ or fixed x ⊥ , respectively. From Fig.3, one finds that v x ⊥ (τ, 0) = 0 and the radial velocity v x ⊥ first increases from x ⊥ = 0, has a maximum at intermediate x ⊥ and then gradually decreases with x ⊥ . As shown in Fig.4, v x ⊥ at fixed x ⊥ becomes smaller at late times, due to the decay of the magnetic field, in agreement with the curves displayed in Fig.3. Figure 5 shows the correction energy density 1 (τ, x ⊥ ) as a function of x ⊥ for different values of τ ; we remind the reader that the total energy density is = 0 (τ ) + 1 (τ, x ⊥ ) and the latter is the component which is truly affected by the magnetic field. Figure 6 shows the correction energy density as a function of τ for different values of x ⊥ . Here we find that for x ⊥ = 0 the correction energy density is positive, starting from zero at proper time τ = 1 fm and showing a shallow maximum; for x ⊥ = 0.5 fm the correction energy density is negative at τ = 1 fm and increases with τ reaching zero at approximately τ = 4.5 fm, becoming then slightly positive. For the other values of x ⊥ the correction energy density is negative at any time and monotonically increases toward zero.  [2] has nearly the same behavior as in Fig. 8, stemming from a similar trend of the correction energy density as a function of τ , like the one illustrated in our Fig. 6. In the Gubser work for τ < 4.6 fm the energy density is positive for x ⊥ ≤ 3 fm and negative for x ⊥ ≥ 4 fm and it is negative for any x ⊥ for τ > 4.6 fm.
It is interesting to investigate variations of the spatial width of the external magnetic field: this affects the Fourier series which reproduces the assumed distribution for the magnetic field; moreover we find that v x ⊥ (τ, x ⊥ ) and 1 (τ, x ⊥ ) have an important dependence on the parameter α (with dimension square of inverse length), which characterizes the spatial width of the magnetic field. In Fig. 9, we plot the external magnetic profile at τ = 1 fm for several different values of α. In Figs. 10 and 11, we plot v x ⊥ and 1 at τ = 1 fm for references. The v x ⊥ gets smaller when α is increased. It seems that the parameter α plays the role of  the parameter 1/q 2 in Ref. [2]: indeed the radial flow velocity (versus x ⊥ at τ = 0.6 fm) becomes smaller when 1/q increases.   21), one may take the first 100 terms of the series in the calculation. Fig. 12 shows a comparison between the approximated magnetic field by the Bessel series and the assumed magnetic profile Eq. (21). Figures 13 and 14 show v x ⊥ (τ, x ⊥ ) at either fixed τ or fixed x ⊥ , respectively. The qualitative behaviors of v x ⊥ (τ, x ⊥ ) in both figures are different from the case n = −1 and the amplitude is smaller. While for n = −1, the direction of the fluid velocity is always positive, for n = −4/3 the direction of fluid velocity changes during the expansion of the fluid. Figure 15 shows the correction energy density as a function of x ⊥ for different values of τ . Figure 16 shows the correction energy density as a function of τ for different val-   Figures 17 and 18 show the normal and Log-Log plots of the total (τ, x ⊥ ) as a function of τ for several values of x ⊥ , respectively. In Fig. 15, we find that for x ⊥ = 0 the correction energy density is always positive and then it decreases from the value at x ⊥ = 0 with increasing x ⊥ . From Fig. 16 one also finds that for the case n = −4/3 the correction energy density is always positive, at variance with the    In the previous sections we have obtained as analytical solution the transverse velocity and energy density in the presence of a weak magnetic field. Now we can use these results to From the local equilibrium hadron distribution the transverse spectrum is calculated at the freeze out surface via the Cooper-Frye (CF) formula: We note that T f is the temperature at the freeze out surface. The latter is the isothermal surface in space-time at which the temperature of inviscid fluid is related to the energy density as T ∝ 1/4 . It must satisfy T (τ, x ⊥ ) = T f . In our convention, where is the longitudinal proper time, x ⊥ the transverse (cylindrical) radius, η = 1 2 log t+z t−z the longitudinal rapidity (hyperbolic arc angle), the azimuthal angle ϕ p belonging to the spacetime point x μ . Similarly u ⊥ is the transverse flow velocity and ϕ is its asimuthal angle. Finally p T is the detected transverse momentum, m T = m 2 + p 2 T the corresponding transverse mass, while Y is the observed longitudinal rapidity, which gives our final expression for the CF formula is the solution of the T (τ f , x ⊥ ) = T f and the degeneracy is g i = 2 for both the pions and the protons. The above integral over x ⊥ on the freeze-out surface is evaluated numerically. The spectrum Eq. (47) is illustrated in Figs. 22 and 23 for three different values of the freeze out temperature (140, 150 and 160 MeV) and compared with experimental results obtained at PHENIX [33] in central collisions. Our proton spectrum appear to underestimate the experimental data, except at low p T , but their behavior with p T has the correct trend of a monotonical decrease. The pion spectrum, instead, appears in fair agreement with the experimental results, which are very close to the theoretical curves. This is an indication that hadrons with different masses have dif-  black, purple and red lines correspond to a freeze out temperature of 140, 150 and 160 MeV, respectively. Circles: PHENIX data [33] ferent sensitivities to the underlying hydrodynamic flow and to the electromagnetic fields. Indeed, the difference between the charge-dependent flow of light pions and heavy protons might arise because the former are more affected by the weak magnetic field than the heavy protons [34].
For comparison, we also show the results obtained by Gubser, which appear to be more flat and typically overestimate the experiment. We also notice that, for the proton case, the highest value of the freeze out temperature we employed (as suggested, e.g. in Ref. [35]) slightly brings (for protons) the calculation closer to the experimental data; however it also shows a kind of saturation phenomenon and points to the need of including other effects not considered in the present work.

Conclusions
In the present work, we investigated central heavy ion collisions in the presence of a transverse external magnetic field. Making use of Milne coordinates, in our setup the medium is boost-invariant along the z direction and the magnetic field, which is a function of τ and x ⊥ , points along the φ direction. The energy conservation and Euler equations reduced to two coupled differential equations, which we solved analytically in the weak-field approximation. We showed in detail how the fluid velocity and energy density are modified by the magnetic field. The solutions obtained by our numerical calculations assume an initial energy density of the fluid at time τ = 1 fm fixed to ∼ 5.4 GeV/fm 3 and a ratio of the magnetic field energy to the fluid energy density, σ , fixed to ∼ 0.015. We consider two different decays with time of the magnetic field: τ n , with n = −1 or n = −4/3. We remark that in Ref. [4] the external magnetic field was approximated by a Fourier cosine series and, due to the oscillatory behavior of the cosine function, the magnetic field reduces to zero in the fringes for |x| = π . Consequently, these authors had to focus on the valid region −π < x < π and the behavior of the transverse velocity and of the correction energy density was difficult to analyze near the fringes. In the present work the magnetic field is approximated by a series of Bessel functions and the solutions are valid for the entire region of x ⊥ , i.e., (0, ∞).
Another point concerning the choice of the τ dependence of the magnetic field is related to the ratio σ between magnetic and fluid energy densities: in Ref. [3], it was found that in central collisions, at the center of the collision region, σ 1 for most of the events; nevertheless, large values of σ were observed in the outer regions of the collision zone. Therefore, our assumption for the spatial distribution of the external magnetic field for the case n = −1 may be more realistic at face of the physical conditions.
In general, our study in a simple setup, which includes an azimuthal magnetic field in the matter distribution, is worthwhile to check the possible effect of this change on the transverse expansion of the fluid. We showed that by combining the azimuthal magnetic field with the boost symmetry along the beam direction, a radial flow perpendicular to the beam axis is created and the energy density of the fluid is altered. We stress that the present work presents an approximated calculation which can be useful for cross checking current and future numerical calculations in some limiting region. Indeed, the effect of such a scenario on hadronic flow in heavy ion collisions requires more pragmatic debates.
Our study can be generalized in many directions: breaking of rotational symmetry can be introduced, the conservation equation can be coupled to Maxwell's equations and solved consistently. Of course in this case only numerical solutions can be found, while in the present paper we were able to obtain analytical solutions.

Data Availability Statement
This manuscript has no associated data or the data will not be deposited. [Authors' comment: Data will be sent by the authors upon request to anybody interested.]