Directed flow of charged particles within idealized viscous hydrodynamics at energies available at the BNL Relativistic Heavy Ion Collider and at the CERN Large Hadron Collider

Following the Boz$\dot{\textrm{e}}$k-Wyskiel parametrization tilted initial condition, an alternative way to construct a longitudinal tilted fireball based on the Glauber collision geometry is presented. This longitudinal tilted initial condition combined with the Ideal-CLVisc (3 + 1)D hydrodynamic model, a nonvanishing directed flow coefficient $v_{1}$ in a wide rapidity range is observed. After comparing the model's results with experimentally observed data of directed flow coefficient $v_{1}(\eta)$ from $\sqrt{s_{NN}}$ = 200 GeV Cu + Cu, Au + Au collisions at RHIC energy to $\sqrt{s_{NN}}$ = 2.76 TeV and $\sqrt{s_{NN}}$ = 5.02 TeV Pb + Pb collisions at the LHC energy. We find that the directed flow measurements in heavy-ion collisions can set strong constraints on the imbalance of forward and backward incoming nuclei and on the magnitude asymmetry of pressure gradients along the $x$ direction.


I. INTRODUCTION
In the past decades, relativistic heavy-ion collisions provided methods to explore and understand the deconfined quark-gluon plasma (QGP). The strongly collective flow of such hot dense medium is one of the key observations in the physics of high-energy heavy-ion collisions [1][2][3][4].
The rapidity-odd directed flow (shorted as directed flow or v 1 ) refers to the collective sideward deflection of particles and is the first-order harmonic of the Fourier expansion of the particle azimuthal distribution with respect to the reaction plane [5,6]. The directed flow is believed to be created at very early stage (during the nuclear passage time: 2R/γ ≈ 0.1 fm/c) at large rapidity η (in the fragmentation region) [7]. Therefore, it may keep trace of the bulk collective dynamics and the subsequent evolution into a thermalized hot QCD matter, which provides a unique insight to investigate the initial condition.
At current stage, the directed flow v 1 and the splitting ∆v 1 between particles and anti-particles are measured for both light charged hadrons and heavy quark productions (e.g., J/ψ, D 0 ,D 0 ) at the BNL Relativistic Heavy Ion Collider (RHIC) and at the Large Hadron Collider (LHC) [8][9][10][11][12]. A striking characteristic of the measured directed flow for charged particle is the large value of v 1 at RHIC energy [8] and large value of splitting ∆v 1 of charm hadrons produced at LHC energy [11]. The directed flow has been investigated by different models and mechanisms, such as the transport model [13], three di-the CLVisc(3+1)D ideal hydrodynamic simulation are presented. In Sec. III, our numerical results on the directed flow of charged particles are presented. Finally, in Sec. IV, we summarize the results and present a short outlook.

II. SETUP OF THE NUMERICAL SIMULATIONS
In this work, following the Refs. [15,21,27], the viscous corrections are not included in the present study 1 and the viscosity effect will be included in future studies. The (3 + 1)-dimensional numerical simulations are performed with in Milne/Bjorken coordinates.

A. Initial condition
The initial energy density distribution is computed according to a modified optical Glauber model [30,40].
In the optical Glauber model, the nucleus thickness function T (x, y) from the Woods-Saxon distribution is where n 0 is the average nuclear density, d is the diffusiveness, x, y, z is the space coordinates and R is the radius of the nuclear Fermi distribution, which depends on the specific nucleus. The parameters used for nucleus Cu, Au and Pb in current study are listed in Table. I. T 1 (x T ) and T 2 (x T ) are the densities of participants from the two nuclei, 1 The directed flow is generated at very early stage in the evolution and it is given by the initial tilted source. The ideal fluid approximation could give a sizable directed flow coefficient v 1 and directly reflect the evolution of pressure [15].
where A is the mass number of the colliding nuclei, σ N N is the inelastic-scattering cross section, σ N N are set to 40 mb for √ s N N = 200 GeV Cu + Cu, Au + Au, 64mb for √ s N N = 2.76 TeV Pb + Pb and 67mb for Pb + Pb √ s N N = 5.02 TeV collisions [41], and where x T = (x, y) is the vector of the transverse plane coordinates and b is the impact-parameter vector, connecting the centers of the two nuclei.
The function W N gives the contribution of the wounded nucleons. To generate a tilted fireball along the longitudinal direction, the wounded nucleons [40] weight function W N is modified as follow, here the longitudinal tilted parameter H t is a free parameter that reflects the imbalance in the emitting contributions from forward-going and backward-going participant nucleons. As we will see in the following section, varying the parameter H t results in strong dependencies in the magnitude of pion's and charged particle's directed flow. The parameter η t is a constant and chosen to be 8.0 for all the collision systems. The energy density distribution at the hydrodynamic starting time τ 0 is given by where ε 0 the maximum energy density given in Table. II. For most-central collisions at RHIC and the LHC energy, the total weight function W (x, y, η s ) is defined as [27,33] W (x, y, η s ) = (1 − α)W N (x, y, η s ) + αn BC (x, y) with α = 0.05 being the collision hardness parameter 2 and n BC (x, y) being the mean number of binary collision: The energy density profile in the longitudinal direction is modulated by [30], pact parameter b, which can be obtained by interpolation [41]. The impact parameter using in current study is presented in Table. III.  With the above parametrizations, we illustrate the profile of energy density and magnitude of pressured gradient distribution for 0-5% Au + Au collisions at  show the magnitude of initial pressure gradients in the transverse plane at τ = 0.2 fm and forward-rapidity η s = 2.1. The magnitude of pressure gradient (−∂ x P and −∂ y P ) is calculated from the initial energy density and equation of state (EoS). We find that the magnitude of the pressure gradient shows an asymmetry along the x > 0 and x < 0 direction in the transverse plane. The value of −∂ x P shows a maximal value of around 16 GeV/fm 4 (14 GeV/fm 4 ) at the x > 0 (x < 0) panel. With the hydrodynamic expansion, considering the contribution from the transverse and longitudinal pressures, the final directed flow coefficient becomes negative for positively rapidity [15,39]. TeV with the parameter H t = 0.70. Figure 2(b) and 2(c) show the magnitude of initial pressure gradients in the transverse plane at τ = 0.2 fm and forward-rapidity η s = 2.1. The magnitude of pressure gradient (−∂ x P and −∂ y P ) shows an asymmetry along the x > 0 and x < 0 direction. Furthermore, the value of −∂ x P shows a maximal value of around 48 GeV/fm 4 (40 GeV/fm 4 ) along the x > 0 (x < 0) direc-tion in the transverse plane.
The information on how the above initial spatial anisotropy is transferred to the momentum space [15,39] is encoded on the directed flow coefficient, which will be presented in Sec. III.

B. Hydrodynamic equations/simulation
The hydrodynamic equations from the literature read [42][43][44][45][46] where T µν = εu µ u ν − P ∆ µν is the energy-momentum tensor for ideal hydrodynamics, u µ = γ(1, v) denotes the fluid velocity four-vector, ε is the energy density and P is the pressure. The pressure P is given as a function of energy density ε by the equation of state (EoS). The lattice QCD equation of state from the Wuppertal-Budapest group (2014) [47] is used in the current study. Notice that the viscosity and net baryon density are set to zero in current study. The energy-momentum conservation equations [Eq. (10)] are solved numerically by using Kurganov-Tadmor (KT) algorithm [30], which was introduced to the field of high-energy physics by the McGill group [48]. The 3D partial differential equations are solved by updating the values of fluid cells at each time step. For each collision system, we run the Ideal-CLVisc (3 + 1)D hydrodynamic simulation with number of cells N cells = N x ×N y ×N ηs = 201 × 201 × 105 for 1600 time steps on GPU NVIDIA GeForce RTX 2080TI (Turing Features) and server CPU Intel Xeon E5-1620v2. For more details about this GPU simulation part, see Refs. [30,49]. The freeze-out condition used in current calculation is the isothermal freezeout condition [30], where one assumes the hypersurface is determined by a constant temperature T f rz = 137 MeV.
Based on the Cooper-Frye formula [50], the Ideal-CLVisc provided a "smooth method" to compute the different particle spectra on the freeze-out hypersurface, where the numerical integration is performed over the freeze-out hyper-surface and smooth particle spectra are [30]. p T and φ are chosen to be Gaussian quadrature nodes to simplify the calculation of p T or φ integrated spectra. Hadron spectra from resonance decays are also computed via integration and parallelized on the GPU.
where Ψ 1 is the first-order event plane of the collision [15].

III. NUMERICAL RESULTS
In this section, numerical results from ideal CLVisc hydrodynamic simulation with experimental data for RHIC and the LHC energies are presented.

A. Cu+Cu and Au+Au
√ sNN = 200 GeV collisions Figure 3 shows a comparison of charged hadron pseudo rapidity distributions dN/dη between our model and the PHOBOS measurement [51] in the most central centrality 0-6% 3 Au+Au and Cu+Cu collisions at GeV. We find that the dN/dη obtained from the Ideal-CLVisc with modified initial conditions under three different settings of H t are almost indistinguishable on the plot.  GeV 4 . Our hydrodynamic simulations with modified initial conditions give a same distribution of v 2 (p T ) under different H t compared with the STAR measurement [54].
Our theoretical calculations of pseudorapidity distribution dN/dη and elliptic flow v 2 (p T ) for different settings of tilted parameter H t show that different initial longitudinal tilted fireball almost does not affect the pseudo rapidity distribution and the elliptic flow coefficient v 2 (p T ). . 4 Please notice that our model following those assumptions: (1) modified optical Glauber model, which miss event-by-event eccentricity fluctuations; (2) dissipative effect is not included during the hydro expansion; (3) the freeze-out temperature T f rz = 137 MeV. Above assumptions makes our model sightly overestimate the elliptic flow coefficient v 2 (p T ) and v 2 (η) in several centrality classes at current stage [17,52,53]. We will focus on the event-by-event dissipative hydrodynamic simulations in future work.
However, a non-zero directed flow v 1 can be generated in the hydrodynamic simulation with the modified initial conditions. Figure 5 shows the result for the directed flow coefficient v 1 (η) of charged particles emitted after a hydrodynamic evolution. The dashed (-dotted) curves are the results for Au + Au and Cu + Cu √ s N N = 200 GeV collisions. One finds the experimental data in the centrality classes 0%-5% and 5%-40% are reproduced well in a large pseudorapidity region. For the peripheral collisions (in the centrality bin 30%-60%), our model overestimates the directed flow coefficient v 1 at large rapidity (fragmentation region). The directed flow in such a region maybe has different origins, such as the baryon stopping effect and fluctuations [13,15] Figure 6 shows only the most central pseudorapidity distribution (0-5%) for charged particles from the Ideal-CLVisc is in comparison with experimental data from the ALICE Collaboration [55,56]. The hydrodynamic simulation with different longitudinal tilted initial condition (H t = 0.0, 0.5, 2.0) gives almost the same charged multiplicity distribution dN/dη for the most central collisions. We find that the elliptic flow v 2 (p T ) is indistinguishable for different H t (dashed and dashed-dotted curves). The experimental data are from the CMS Collaboration [57] and ATLAS Collaboration [58].
In Figure 8 we plot the pseudorapidity dependence of the directed flow v 1 (η) for pions (π + ). Figure 8  TeV in the centrality classes 10%-20% and 30%-40%. We see that hydrodynamic simulation with modified initial conditions reproduce the experimental data from ALICE Collaboration [59] in the centrality bin 10%-20% with tilted parameter H t = 0.7 and 30%-40% with H t = 0.8. Figure 8 (bottom panel) shows the result for Pb + Pb √ s N N = 5.02 TeV in the centrality class 5%-40%.
We see that the ALICE Collaboration [11] points for the centrality class 5%-40% are for positive charged particles, and results for the pion (π + ) directed flow from CLVisc simulation agree with the experimental data to a reasonable level. This is because (a) a distinction for the whole positive and negative particle yields needs a proper theory, which has not been developed totally in the current stage; (b) π + plays a vital role for the yield of the total positive charge particles and mainly comes from the hydrodynamic evolution [55,56].
To get information that how the longitudinal tilted structure introduced in Eq. (5) is related to the final azimuthal asymmetry measured by the directed flow, we plot the H t dependence of the directed flow of charged particles in Figure 9 for Pb + Pb  ficient at large rapidity. The value of H t extracted from the STAR and ALICE data show that (a) the initial spatial pressure gradient asymmetry in the transverse plane of bulk medium at RHIC energy is larger than at the LHC energy at the initial proper time; (b) the larger the impact parameter b, the larger the magnitude asymmetry of initial pressure gradient in the transverse plane at the initial stage.

IV. SUMMARY
In this work, following the previous work Refs. [15,21,27,39], an alternative parametrization to construct longitudinal tilted initial condition based on the Glauber model is presented. A data-driven phenomenological rapidity-dependent weight function [Eq. (5)] is obtained to generate the initial tilting longitudinal expansion of the strongly-coupled QCD matter, and the results show that such a tilting source leads to a magnitude asymmetry of the pressure gradient along the x-direction and longitudinal η s direction [15]. The existence of such asymmetries push the light quarks and gluons toward the x > 0 (and x < 0) direction at different rapidities, which gives a nonzero directed flow coefficient.
The pseudorapidity distribution dN/dη and the elliptic flow coefficient v 2 (p T ) of charged particles are presented for Cu + Cu, Au + Au and Pb + Pb collisions for serval longitudinal tilted parameters H t . The longitudinal tilted structure in hot-QCD matter affect both the multiplicity density distribution and elliptic flow coefficient weakly and almost negligible, consistent with previous studies [15,39].
Our modified initial conditions with the ideal CLVisc hydrodynamic simulation is used to fit the directed flow coefficient v 1 (η) of charged particle and pions measured by the STAR Collaboration and the ALICE Collaboration. We find that the directed flow coefficient is generated at a very early stage in the evolution and is given by the initial tilted source. The directed flow coefficient is decreased with the center-of-mass energy √ s N N , the reasons are (a) with increasing range of longitudinal rapidity the imbalance between the thickness function T 1 and T 2 goes down; (b) the contribution of binary collisions to the Glauber model is believed to increase with the √ s N N [27,39].
We also remark that, besides the initial fireball spatial asymmetry contributing to the directed flow, the following aspects are also important and can be extended to future work.
(1) In real heavy-ion collisions, the fluctuations of the initial energy density also contribute to the local pressure gradient asymmetry, e.g., T R ENTo-3D initial condition model [60]. It is possible that the dynamical fluctuations or the initial flow in the z direction can generate a large pressure gradient asymmetry for both angle and magnitude [15-17, 28, 60].
(2) During the hydrodynamic expansion, the viscosity corrections reduce the longitudinal pressure and increase the transverse pressure. As a result, a smaller directed flow coefficient is observed after a viscous hydrodynamic simulation. This means that the ideal hydrodynamics is a suitable tool for the study of the presence of directed flow in heavy ion collisions [39]. The expansion of the strongly coupled matter is affected by the shear viscosity and bulk viscosity in both transverse plane and longitudinal direction. The observables such as the elliptic flow and p T spectra are more sensitive to the viscosity effect according to viscous hydrodynamics calculations [30,52,53,61,62]. Moreover, the v 2 (η) coefficient is sightly overestimated in our model as the lack of the shear viscosity effect [30,39].
Because v 2 (η) places severe tight constraints to the initial shape of the fireball away from midrapidity, and thus to the parametrization of the initial state [30,52,53]. In the future study, viscous correction will be included and the v 2 (η) coefficient will be used to constrain our initial condition.
(3) For non-central heavy ion collisions, an extremely strong magnetic-field is created by the colliding charged beams moving at relativistic speed (almost 10 16 − 10 20 Gauss) [63][64][65][66]. Due to the expansion along the beam axis, the Lorentz force is directed along the negativex direction in the forward rapidity region for positively charged quarks, which generates a directed transverse flow. In addition to the above Hall effect, the time dependence of the magnetic-field generates an electric field due to the Faraday effect. The induced Faraday current provides a no-zero finite drift velocity in the transverse plane due to the magnetic field. Lots of work has investigated the contribution of the combination of the above two effects on the directed flow coefficient, such as the MHD [27] and decoupled hydro-magnetic frame [22,24]. And one has found that the contribution of the magnetic field effect on the soft hadron directed flow coefficient is less than 5 × 10 −4 at large rapidity, which is almost 80 -100 times smaller than the contribution from the tilted initial conditions [27].
(4) Taking into account the asymmetry between forward and backward moving participants, the noncentral heavy ion collisions produce not only strong angular momentum, strong magnetic-field but also global and local vorticity and hyperon polarization [28,31].
(5) The directed flow of heavy hadrons could be a great probe to investigate the initial pressure gradient asymmetry of bulk matter [16,23]. Vice versa, such kind of tilted medium could be the background of heavy quark propagation, and its effect on open heavy flavor production v 1 and R AA might be interesting problems [21,[67][68][69][70]. Furthermore, D 0 directed flow found at both STAR and ALICE still contains large statistical uncertainly and systemic uncertainly, one suggests that hydrodynamic + transport model together may put more constraints on the directed flow for better understanding of the initial stage in heavy ion collisions [26]. (6) We also need to take into account the hadronic cascade. For a proper comparison with more experimental data, one should include such interactions with the addition of hadronic transport models such as UrQMD [71,72] and SMASH [73,74]. (7) Recently, a large number of studies used the parametrizations for the longitudinal structure of the fireball in Ref. [15] to investigate the directed flow coefficient of heavy meson [19][20][21]28]. Collision geometry-based 3D initial conditions with hydrodynamic simulations from the group of C. Shen et al. [17,18] also described the directed flow coefficient of π + at RHIC energies well. It will be interesting to see which one of the models can lead to the reproduction of the experimental data with the least amount of tweaking of the parameters.
These important aspects will be studied in the future.

Acknowledgments
The authors thank Xin-Nian Wang for helpful comments and providing the GPU computing platform at the initial stage of this study. Z-F. Jiang would like