The signature of charge dependent directed flow observables by electromagnetic fields in heavy ion collisions

We discuss the generation of the directed flow $v_1(p_T,y_z)$ induced by the electromagnetic field as a function of $p_T$ and $y_z$. Despite the complex dynamics of charged particles due to strong interactions generating several anisotropies in the azimuthal angle, it is possible at $p_T>m$ to directly correlate the splitting in $v_1$ of heavy quarks with different charges to some main features of the magnetic field, and in particular its values at formation and freeze-out time. We further found that the slope of the splitting $d\Delta v_1/dy_z|_{y_z=0}$ of positively and negatively charged particles at high $p_T$ can be formulated as $d\Delta v_1/dy_z|_{y_z=0}=-\alpha \frac{\partial \ln f}{\partial p_T}+\frac{2\alpha-\beta}{p_T}$, where the constants $\alpha$ and $\beta$ are constrained by the $y$ component of magnetic fields and the sign of $\alpha$ is simply determined by the difference $\Delta[tB_y(t)]$ in the center of colliding systems at the formation time of particles and at the time when particles leave the effective range of electromagnetic fields or freeze out. The formula is derived from general considerations and is confirmed by several related numerical simulations; it supplies a useful guide to quantify the effect of different magnetic field configurations and provides an evidence of why the measurement of $\Delta v_1$ of charm, bottom and leptons from $Z^0$ decay and their correlations are a powerful probe of the initial e.m. fields in ultra-relativistic collisions.


I. INTRODUCTION
The ultrarelativistic heavy ion collisions (uRHICs) at both the BNL Relativistic Heavy Ion Collider (RHIC) [1,2] and the CERN Large Hadron Collider (LHC) [3] have created a new state of matter, the quark-gluon plasma (QGP), during their early stage and showed that such a matter is the most perfect fluid in nature [4][5][6]. In the last decades, there are numerous studies in the search for the parity (P ) and charge conjugate parity (CP ) symmetry breaking processes in quantum chromodynamics (QCD) happened in QGP, mainly via the chiral magnetic effect (CME) [7][8][9][10][11][12] and the chiral vortical effect (CVE) [13]. The strongest ever electromagnetic (e.m.) field and the largest relativistic vorticity [14,15] are created in noncentral heavy ion collisions. The strong e.m. field can lead to many other interesting phenomena such as the chiral magnetic wave (CMW) [16][17][18][19] and the splitting in the spin polarization of hyperons [20][21][22]. However, there are a lot of uncertainties in the calculation of magnetic field in heavy ion collisions, and this inspired the search for a direct probe to the strong e.m. fields by measuring the charge dependent flows (v n ) of charged mesons and baryons as well as neutral charged charmed mesons [23,24]. Though these numerical studies are very meaningful for the understanding of e.m. fields, a general description of charge dependent flows going beyond the details of e.m. fields is also important. In Refs [25], we have found several important features of the directed * Electronic address: sunyfphy@lns.infn.it † Electronic address: greco@lns.infn.it ‡ Electronic address: salvatore.plumari@ct.infn.it flow splitting, ∆v 1 , of positively and negatively charged heavy quarks as well as leptons induced by e.m. fields: (i) It is not very sensitive to the details of the spatial and time configurations of e.m. fields; (ii) d∆v 1 /dη at midpseudorapidity (η) depends on the slope of the transverse momenta spectra, especially at high p T ; (iii) The modification by the interaction with QGP for heavy quarks is negligible at p T higher than 2-3 GeV/c. This persuades us to find the general physics behind these features and extend the findings to other charge dependent flow observables. Since the effect due to the interaction with QGP is small for charged quarks at high p T , the findings should have a general application for high p T heavy quarks and energetic jets as well as leptons of arbitrary p T . The purpose of these studies is not only trying to build the bridge between the spatial and temporal configurations of e.m. fields and final observables in the theoretical side, but also can be used to determine whether the ∆v n observed experimentally has an electromagnetic origin.
The configuration of non-central heavy ion collisions is fixed throughout the paper, where the center of the nuclei moving in positive z direction localizes in positive x axis, which produces a strong magnetic field and large vorticity in negative y direction. When mentioning the formula in this paper, it should be noted that it is deduced by assuming a pure interaction with e.m. fields, which means that it should be applied to high p T heavy quarks and to leptons at arbitrary p T where the strong interaction with QGP leads to negligible modifications. The numerical results from transport simulation for heavy quarks include, however, the strong interaction with QGP according to the current standard approach to their dynamics.
The paper is organized as follows: In Sec. II, we describe the general formula for the charge dependent flow observables induced by e.m. fields, and give the physical meaning to the coefficients in the general formula. In Sec. III, we study specifically the ∆v 1 of positively and negatively charged particles, and present the direct connection between the coefficients in ∆v 1 and the y component of the magnetic field. Sec. IV presents several numerical results that can be understood by the general formula probing the robustness of the formula. We also discuss the importance of measuring the ∆v 1 and propose also a new measurement of the spectra ratio of positively and negatively charged leptons from Z 0 decay. Summary and conclusions are discussed in Sec. V.

II. CHARGE DEPENDENT FLOW HARMONICS BY ELECTROMAGNETIC FIELDS
The effects of the e.m. fields on the phase-space distribution function can be expressed in terms of a transition function T (∆p x , ∆p y , ∆y z , p x , p y , y z ). The function T represents the distribution of the shifts in the transverse momenta p T = (p x , p y ) and rapidity y z due to the electromagnetic field. In order to guarantee the particle number conservation the distribution function T satisfy the following normalization condition d 2 ∆p T d∆y z T (∆p x , ∆p y , ∆y z , p x , p y , y z ) = 1. Starting from a boost-invariant spectra of charged particles f (p T ), after the modification by e.m. fields, which is considered as a small perturbation, the distribution f ′ (p T , y z ) shall be: where f T = f (p T , y z )T (∆p T , ∆y z , p T , y z ) and the derivatives are evaluated at (p T , y z ), and the average shifts ∆p a with a = x, y, z are defined as: By the definition of rapidity, one can further express ∆y z in terms of ∆p x , ∆p y , ∆p z as: where φ = tan −1 (p x /p y ) is the azimuthal angle relative to the reaction plane in momentum space. Since the colliding systems are symmetric with y ↔ −y, in momentum space one should have ∆p x (p T , φ, y z ) = ∆p x (p T , 2π − φ, y z ), −∆p y (p T , φ, y z ) = ∆p y (p T , 2π − φ, y z ) and ∆p z (p T , φ, y z ) = ∆p z (p T , 2π − φ, y z ). Therefore, after a Fourier decomposition with respect to the angle φ, the average shift can be expressed as: In heavy ion collisions, if the chiral magnetic conductivity is zero, there is no B z [26,27], and according to Lorentz force the momentum shifts are: The above set of equations suggests that the coefficients a n , b n and c n have direct relations with the e.m. fields as a function of φ x in coordinate space. Since then the distribution function f ′ (p T , y z ) after the modification of e.m. fields relates to the initial f (p T , y z ) as: One can read from Eq. (7) that the Lorentz force in the longitudinal direction (c n ) leads also to non-zero charge dependent v n that measures the anisotropy in transverse momenta. This is an effect that was not brought to light in the previous studies focused on the numerical simulations [23][24][25][28][29][30], even if they naturally include it.
If p T is larger than the mass of charged particles, Lorentz force is not sensitive to p T any more, because p T /m T ≈ 1 and furthermore the equations of motion do not depend on the p T of the particle implying similar trajectories. All of this together lead to: while in general a n , b n , c n are p T dependent and the specific dependence is determined by the way the strong interaction with the QGP medium acts on the specific particle. Given that such momentum dependence can be discarded in the high p T limit, one can rearrange the terms in Eq.(7) that can be rewritten as: where d n and e n are mixed combinations of a n , b n and c n which can be read from Eq. (7). Moreover, if two types of charged particles have similar formation time, such as charm quarks and leptons from Z 0 decay, then at p T ≫ m all their coefficients a n , b n and c n are differed only by their charges regardless of the complex spatial and temporal configurations of e.m. fields, which provides a strong correlation between their flow observables supplying strong probes of the e.m. fields. One immediate consequence of Eq. (9) is that for charged particles with a peculiar spectra, for example a sudden change in p T like the leptons from Z 0 decay, as long as these coefficients are non-zero, there will be a sudden change in the spectra ratio and the ∆v n of positively and negatively charged particles just at exactly the same p T .
Eq. (7) provides a general p T scaling for all flow observables induced by e.m. fields, and it reduces to Eq. (9) at p T ≫ m, for any configuration of e.m. fields. This scaling is certainly different from the collective flows generated by the strong interaction with QGP for light and heavy quarks, since a n , b n and c n induced by e.m. fields are charge dependent, and become constant at p T ≫ m, while they are not charge dependent but flavor dependent hence with a p T dependence determined by the flavor dependence of the strong interactions. In a preliminary study presented in a recent Proceedings [28], we showed that the interaction with QGP is negligible for the charge dependent flow observables induced by e.m. fields at high p T for heavy quarks and leptons of arbitrary p T , hence the ∆v n should be a general signature of effects induced by e.m. fields. In Section IV the comparison to realistic simulations including the effective hot QCD matter interaction will allow to assess the range of validity of the approximations done to deduce Eq.(9).

III. CHARGE DEPENDENT DIRECTED FLOW INDUCED BY ELECTROMAGNETIC FIELDS
Reading from Eq. (7), under the conditions in Eq.(8), the directed flow v 1 becomes, up to quadrupole moments: In AA collisions in the overlapping region, one finds the quadrupole moment of e.m. fields is smaller than their 0-th order moment, which can be seen by looking at the e.m. fields at the very initial stages of AA collisions.
With this assumption, one has the simplified: It is noted that v 1 has a contribution also from the Lorentz force in longitudinal direction. In AA collisions, since the colliding systems are also symmetric with x, y, z ↔ −x, −y, −z, which leads to a n (y z ) = (−1) n+1 a n (−y z ), b n (y z ) = (−1) n+1 b n (−y z ) and c n (y z ) = (−1) n+1 c n (−y z ), so all terms in Eq. (11) are non-zero and v 1 is odd in rapidity. As shown in Eq.(11) at leading order the v 1 depends only on the two coefficients a 0 and c 1 . In order to understand the role of the e.m. field in the building up of the v 1 , in this section, we study the relation between these coefficients and the e.m. field. Starting from Eq. (5) relating the ∆p x to the time integral of the Lorentz force, we project into the 0-th order in the azimuthal angle both sides, which gives for the coefficient a 0 : where we have used transverse coordinates ρ = x 2 + y 2 and φ x = tan −1 (x/y) while η s is the space-time rapidity η s = 1 2 ln t+z t−z . In the second line, we make use of the approximation η s ≈ y z valid in a boost invariant geometry. Notice that in Eq. (12) E x0 and B y0 are the 0-th order in a Fourier decomposition of the e.m. field . Though in Eq. (12) one has to take into account the coordinate and momentum distributions of charged particles initially produced in the overlap region as well as their corresponding trajectories, it is possible to reduce their complexity. As shown in Ref [25], if the e.m. fields do not change too abruptly with respect to ρ in the overlapping region, which is the case in AA collisions, it is possible to show that the overall effect can be approximated as: where K is a positive constant that depends on the spatial distribution of e.m. fields and the spatial distribution of charged particles when they are initially formed, Θ is the step function, R is a radius about the average of the effective ranges of e.m. fields and overlapping region of colliding systems, and γ is pT mT ( t cosh yz − τ 0 ) with t 0 = τ 0 cosh y z the formation time of charged particles. Because of the Faraday's Law, E x and B y at ρ = 0 are related; we express them as: (the negative sign in B y because its direction stays along the negative y direction). The understanding of the key features of the strength and time dependence of the electromagnetic field that determines the magnitude and the sign of ∆v 1 , is a main aim of the present work. Our strategy has been to reduce the Maxwell equations to a 1D integral equation that can give a quite good approximation of the relations between B y (t, η s ) and E x (t, η s ); this is obtained from the Faraday's Law ∇ × E = −∂B/∂t under the assumption of small space gradients ∂Ez ∂x ∼ 0, hence in the inner part of the QGP fireball this allows to write at ρ = 0 and small η s : with h(t, 0) = 0, i.e. at the collision center E x = 0.
To check the reliability of Eq. (15), we compare E x at ρ = 0 given by h(t, η s ) through Eq. (15), once B y (t, η s ) = −g(t, η s ) has been calculated by solving Maxwell equations, with the E x solution of the full Maxwell equations that includes space gradients. To get an analytical solution of Maxwell equations coupled with conducting medium is quite though in heavy ion collisions, and it is still not yet possible to have reliable solutions in magnetohydrodynamics unless in infinite conductivity case [27]. However, assuming a constant and uniform conductivity, one can obtain the analytical solutions of e.m. fields [23,26,31], which are adopted by several studies recently [23,24,30,32]. The analytical results in this case of E x at ρ = 0 in 5.02 TeV Pb+Pb collisions at impact parameter b = 7.5 fm are shown by the solid lines in Fig. 1, where the different choices of electrical conductivity σ el = 0.0115, 0.023 and 0.046 fm −1 are included. It is seen that with larger conductivity, the maximum magnitude of E x decreases, while the lifetime of it increases. The dashed lines are E x calculated by Eq. (15), i.e. h(t, η s ) , from the time evolution of B y at ρ = 0, that is given by the analytical results in the same conditions. It is seen that Eq. (15) coincides with analytical results at η s = 0.5, and it is a good approximation even at η s = 1. By combining Eqs. (13) and (15), one can relate a 0 directly to B y and its time derivative implicitly including the effect of E x : with t 0 = τ 0 cosh y z and t f (p T ) = (τ 0 +Rm T /p T ) cosh y z . The slope d∆a0 dyz | yz=0 of positively and negatively charged particles and anti-particles is an important quantity of v 1 and it becomes the following simple form: with τ 1 (p T ) = τ 0 + Rm T /p T . Eq. (17) shows that the sign and magnitude of d∆a0 dyz | yz=0 is just determined by the difference of tB y in the center of fireball at the formation time of particles and the time when particles escape the control of e.m. fields, or when particles freeze out. The detailed information of e.m. fields is irrelevant, of course this comes out under the approximation that the gradients of the fields are not too large within the inner part of the fireball and the balance between electric and magnetic field scales in a self-similar way with r−space coordinates. The factorization found in terms of only the t dependence of B y allows a new insight into the generation of the splitting in v 1 between particles with different charges reducing the delicate balance between the magnetic Lorentz force and the Faraday's effect in terms of the time evolution of the magnetic field only.
From Eq. (5), since the system is symmetric with y ↔ −y, which leads to a B x with no 0-th order expansion in φ x , we understand the leading term to c 1 from B x comes from its 2nd order expansion. However, this is expected to be small in the overlapping region of colliding systems. The same happens to the 1st order expansion in φ x of E z , see analytical solutions using constant and uniform conductivity in Ref [26]. The dominant contribution to c 1 is thus from B y ≃ −g(t, y z ) in its 0th order expansion, which can be approximated simply as: dτ g(τ cosh y z , y z ). (18) From Eqs. (16) and (18) one finds a 0 and c 1 are p T independent at p T ≫ m, and this holds for any configurations of e.m. fields as we discuss before. Eqs. (11), (16) and (18) capture essential ingredients for the charge dependent v 1 induced by e.m. fields and tell what information we can extract from the experimental measurement.

IV. NUMERICAL RESULTS
In this section we will use some specific cases to see how Eqs. (11), (16) and (18) can help us better understand numerical results from realistic simulations. We still choose the 5.02 TeV Pb+Pb collision systems at b = 7.5 fm and the e.m. fields are given by the analytical solutions of Maxwell equations assuming a constant and uniform conductivity. It should be noted that there is a discontinuity of e.m. fields in the initial stages of heavy ion collisions, since non-zero conductivity has to appear after collision rather than even before. So the numerical study is just a platform to test our analytical formula rather than making predictions to be tested in experiments. The general analytical formula however can be applied to all charge dependent flow observables induced by e.m. fields in AA systems in relativistic heavy ion collisions.
We include realistic initialization of heavy quarks in transverse momentum space. Charm quarks are formed at same τ 0 = 0.1 fm/c and we initialize the momentum distribution of charm quarks spectra f c in 5.02 TeV Pb+Pb collisions with the prompt distribution obtained within the Fixed Order+Next-to-Leading Log (FONLL) QCD [33,34], that reproduces the D-mesons spectra in pp collisions after fragmentation. We parametrize it as: where the parameters are A = 20.28, n = 1.951, α = 3.137 and B = 0.0752 respectively. The solid red line in Fig. 2 shows the p T dependence of − ∂ ln fc ∂pT , which is seen to approach the maximum value of 0.8 GeV −1 at 3.5 GeV/c for charm and 0.4 GeV −1 at about 6 GeV/c for bottom.
To study quantitatively the dynamics of HQs we solve the relativistic Langevin equation in an expanding QGP background. The background medium is described by the relativistic transport code with fixed shear viscosity to entropy density ratio close to the lower bound 1/4π which was constrained by the experimental data on the collective flows of charged particles [35][36][37][38]. The dynamics of heavy quarks is studied by standard Langevin equations [39][40][41][42][43][44][45][46][47] with the inclusion of Lorentz force [24,32]: where the momentum diffusion coefficient D p is related to the drag coefficient Γ by D p = ΓET , and ξ i is a real number randomly sampled from a normal distribution with ξ i = 0 and ξ i ξ j = δ ij . Before the formation of QGP at about 0.3 fm/c in 5.02 TeV Pb+Pb collisions, Γ and D p are set to zero and heavy quarks interact with only e.m. fields. Γ and D p are derived from a quasi-particle model (QPM) [48][49][50].
The QPM approach accounts for the non-perturbative dynamics by T −dependent quasi-particle masses, with m 2 q = 1/3g 2 (T )T 2 and m 2 g = 3/4g 2 (T )T 2 , plus a T −dependent background field known as a bag constant, with g(T ) tuned to fit the thermodynamics of the lattice QCD [51,52]. This approach has been shown to lead to a good description of the experimental data for both R AA (p T ) and v 2 of charmed and bottomed mesons, both at RHIC and LHC employing an enhancement factor K ∼ 2 of the drag and diffusion coefficients [53][54][55].
A. Transverse momentum dependence of d∆v1/dyz In this section we present our results on the transverse momentum dependence of the directed flow. We first studied charm quarks under e.m. fields generated by conducting medium with an electrical conductivity σ el = 0.023 fm −1 , which is within the bound of LQCD results [56][57][58]. In Fig. 3 by the solid purple line, we show the slope d∆v c 1 /dy z | yz=0 of charm quarks and anti-quarks after the evolution in e.m. fields as well as in QGP. The slope is seen to be negative though very small, and its magnitude increases with p T initially but then decreases at p T >3 GeV/c. The results can be understood exploiting Eqs. (11), (17) and (18), which lead to the following scaling assuming ∂a0 ∂pT close to zero: with where we recall that the function g(t, y z ) = −B y (t, y z ). We see that if B y does not strongly change with η s , β is positive.
We used the scaling above to fit the numerical results at p T > 3 GeV and found that the scaling with α = −3.6 MeV and β = 2.5 MeV agrees quite well with the numerical results. Moreover, with τ 1 g(τ 1 , 0) ≈ 2.5 MeV, τ 0 g(τ 0 , 0) = 11.5 MeV and |q c | = 2/3, one can find K ≈0.6. At low p T , there appears a deviation from the scaling, which is mainly due to the suppression by the strong interactions with the QGP [28]. Still it is remarkable it works also for bottom at p T > 6 Gev/c.

B. The effect of conductivity on d∆v1/dyz
One of key properties of QGP is its electrical conductivity, and increasing the conductivity can increase the lifetime of e.m. fields, which affects the charge dependent flow observables. In Fig. 4, we show the variation of d∆v D 1 /dy z | yz=0 of D 0 (cu) and D 0 (cu) with the variation of σ el , where charmed meson are formed by the Peterson fragmentation as done in [44,53]. It is seen in Fig. 4 that with larger conductivity, the magnitude of the slope becomes smaller. Reading from Eq. (23), which relates α directly to τ 1 g(τ 1 ) − τ 0 g(τ 0 ) = −(τ 1 B y (τ 1 ) − τ 0 B y (τ 0 )), one can see why it is so by looking at the evolution of −tB y at the center of the colliding systems, as shown in Fig. 5. It is seen that increasing σ el will decrease the magnitude of magnetic field initially while increase it at latter time. This leads to the decrease of the magnitude of α due to the decrease in the difference of tB y at τ 1 and τ 0 . Reading from Fig. 5, we can obtain α 1 , α 2 and α 3 from 0.0115 to 0.46 fm −1 taking the ratio (α 1 − α 2 ) : (α 2 − α 3 ) ≈ 0.7, which agrees quite well with the results in Fig. 4. We have to note that such a result is non trivial because one would expect that a larger conductivity inducing a B y with a longer lifetime and a larger strength for nearly all the time evolution, see Fig. 5, would generate a stronger charge/anti-charge splitting of the v 1 . This important aspect is caught by the formula we have derived in the previous Section. The physical reason for this behavior can be understood considering that for small conductivity the quick variation of the magnetic field generates a strong electric field by the Faraday's law that wins over the Lorentz magnetic force that acts in the opposite direction. At increasing conductivity the magnetic field has a slower evolution thus inducing a smaller electric field and this reduces v 1 because there is a nearly exact cancellation between the magnetic field and the electric field. Under the approximation done we have been able to trace back such a delicate dynamics in terms of the variation of the tB y (t) and the slope of the particel spectrum according to Eqs. (22) and (23). It can be expected that with larger conductivity, τ 1 g(τ 1 ) − τ 0 g(τ 0 ) becomes positive, and it should lead to a positive d∆v 1 /dy z | yz=0 that would agree with the experimental measurement in ALICE [59]; this indeed has been seen in Ref. [25].

C. The ∆v1 from charm to bottom quarks
Switching from charm to bottom quarks in the study of ∆v n , one encounters the differences in quark's charge, in the interaction strength with QGP, in the initial spectra, in the mass and in the formation time. As the effect by the difference in charge is trivial, and the interaction strength difference plays a negligible role at high p T , we thus focus on the variations of ∆v n by the other three differences whose impact can be envisaged by our factorized formula in Eq. (22). To isolate the effects induced by these three differences separately, we make only one change each time. The colliding system is still 5.02 TeV Pb+Pb collisions at b = 7.5 fm with a fixed electrical conductivity σ el = 0.023 fm −1 . We first study the effect of particle's spectra on ∆v 1 by using the initial spectra to bottom quarks, which is obtained by FONLL as well. The parameters are found to be A = 0.468, n = 1.838, α = 3.076 and B = 0.0302 respectively, and the − ∂ ln f b ∂pT is shown by the blue dashdotted line in Fig. 2.
As shown in Fig. 6, where the slopes d∆v 1 /dy z | yz=0 at p T = 3 GeV/c are found to be -5.8×10 −3 with charm spectra and -3.1×10 −3 with bottom spectra and -2.9×10 −3 with charm spectra and -2.0×10 −3 with bottom spectra at p T = 10 GeV/c, the spectra taken from bottom quarks decreases the magnitude of the slope of ∆v 1 vs y z , and the suppression is smaller at p T = 10 GeV/c compared to that at p T = 3 GeV/c. The results can be understood by the general scaling in Eq. (22), which is generated by their difference in − ∂ ln f ∂pT of charm and bottom quarks. Specifically, as shown in Fig. 2, − ∂ ln f ∂pT of bottom quarks is always smaller than charm quarks, and their ratio approaches maximum of about a factor of two at p T =3-4 GeV/c, while decreasing at high p T . In Fig. 7, we study the effect of mass of quarks on ∆v 1 vs y z of mesons fragmented by heavy quarks, where we solely change the mass of quarks from m c = 1.3 GeV to m b = 4 GeV, and the mesons from m D = 1.87 GeV to m B = 5.27 GeV. As shown in Fig. 7, the slopes d∆v 1 /dy z | yz=0 at both p T = 5 and 6 GeV/c are found to be about -4×10 −3 for both charm and bottom quarks within the statistical uncertainty of the calculation. This can be expected because, as shown before, the Lorentz force effect is similar at p T > m, and the mass effect will thus be negligible at high p T .
Finally we study the effect of the formation time on ∆v 1 . The formation time of heavy quarks is given by the pair production process that is approximated as 1/2m, and we thus vary τ 0 from 0.1 fm/c to 0.033 fm/c to see how it affects ∆v 1 . The numerical results at p T = 3 GeV/c and 10 GeV/c are shown in Fig. 8, where the slopes d∆v 1 /dy z | yz=0 at p T = 3 GeV/c are found to be -5.7×10 −3 with charm formation time and -3.1×10 −3 with bottom formation time and -2.9×10 −3 with charm formation time and -1.2×10 −3 with bottom formation time at p T = 10 GeV/c. The change is seen to be surpris- ingly large with such a small change in the formation time. The results can be understood again by Eq. (17), or Eq. (22) where tB y at τ 1 does not change, but it changes significantly at τ 0 = 0.1 and 0.033 fm/c, see the blue dash-dotted line in Fig. 5. Our Eq. (22) shows that more generally ∆v 1 is sensitive to τ 0 only when the difference of tB y at τ 1 and τ 0 is dominated by tB y at τ 0 ; if the time dependence of the magnetic field is such that tB y at τ 1 dominates over its value at τ 0 , ∆v 1 should not change significantly by varying τ 0 . In Fig. 9 we show the results for d∆v B 1 /dy z | yz=0 for the case of B mesons at σ el = 0.023f m −1 . We can see the slope of the splitting is about a factor of six smaller than the charm case, see Fig. 4. Such a factor arises from a factor of 2 from the charge, about a factor of 2 from the smaller formation time, and roughly about a factor 1.5 from the smaller slope of the bottom spectrum at p T ≃ 5 − 10 GeV/c D. The correlation between charmed mesons and leptons from the decay of Z 0 As shown in the section above, the charm and bottoms are so different that we may not be able to determine safely whether the experimental measurement of ∆v 1 for B and D have sole e.m. fields origin. Moreover, to extract both α and β from experimental data one needs to use the scaling α −∂ ln f ∂pT + 2α−β pT to fit the data at high p T where the interaction with QGP plays a negligible role on the charge dependent flow observables; However, due to −∂ ln f ∂pT ∝ 1 pT at high p T because of the power law decay of the spectra of quarks at high p T , it is hard to identify α and β separately, because the last is also ∝ pT pT . On the other hand, in Ref [25], leptons from Z 0 decay are found to be an excellent probe of e.m. fields based on the following reasons: (i) Leptons weakly interact with QGP and do not have complex hadronization mechanisms as heavy and light quarks so as to be a cleaner probe; (ii) Leptons from the decay of Z 0 share a similar formation time τ 0 = 1/2.5 GeV −1 as charm quarks, so that all the coefficients a n , b n and c n of leptons can be approximated as 1.5 times those of charm quarks, because they experience the same ∆(tB y (t)) as charm. This makes an important difference with respect to bottom quarks because it resets the uncertainty coming from the difference in the formation time, discussed above. A test of this should be a strong probe of e.m. fields; (iii) Leptons from the decay of Z 0 have a peculiar spectra so that it is easier to identify α, β coefficients; (iv) Since the Lorentz force becomes same at p T ≫ m, a measurement of constant α at such high p T by leptons from Z 0 decay should also be a strong probe of e.m. field.
In this section, we will study how ∆v 1 of leptons correlates with charm quarks, and see if it fulfills our expectation, through numerical simulations for 5.02 TeV Pb+Pb collisions at b = 7.5 fm. We will also show how the general formula (7) instructs us to make general predictions for the leptons.
The spectra of leptons is generated by decaying Z 0 into lepton pairs, where we first obtain the momentum distribution of Z 0 by fitting the experimental measurements [60,61]: The parameters a = 0.6896, n = 0.4283 and ∆ l = 3.034 are found to give quite a good description of p T and y z dependence of Z 0 in 5.02 TeV Pb+Pb collisions [25]. From the spectra of Z 0 , the spectra of lepton can be obtained: where y 2 and y 3 are given by the energy conservation and have two sets of solution: In the above, Γ l /Γ tot is the branching ratio, p f is the magnitude of momentum of each lepton in COM frame of Z 0 (2 p 2 f + m 2 1 = m Z 0 ), m 1 = m 2 are the mass of lepton pairs, and φ is the angle between the transverse momentum of lepton pairs.
The transverse coordinate of Z 0 is given by the binary collisions of colliding nuclei, and the formation time t and longitudinal coordinate z are given by t = τ Z 0 cosh y z and z = τ Z 0 sinh y z with τ Z 0 = 1/m Z 0 = 0.0022 fm/c. Finally the spacetime coordinate of produced leptons is given by their mother Z 0 that moves in a straight line with a decay time having a distribution ρ(∆t) ∝ e − Γ tot ∆t γv with Γ tot = 2.495 GeV and γ v being the Lorentz contraction factor. In Fig. 10, we show − ∂ ln f ∂pT of leptons from the decay of Z 0 deduced by Eq. (25). It is seen that − ∂ ln f ∂pT is negative at p T < m Z 0 /2 = 45 GeV/c, and it jumps to a large and positive value above 45 GeV/c due to the kinematic effect. This peculiar spectra should imprint a signature in the spectra ratio and ∆v n of positively and negatively leptons inspired by the general formula in Eq.
The results of d∆v l 1 /dy z | yz=0 of lepton pairs in 5.02 TeV Pb+Pb collisions at b = 7.5 fm, which is generated by e.m. fields with σ el = 0.0115, 0.023 and 0.46 fm −1 , are shown in Fig. 11. In general, there is a sudden drop of the slope at p T = m Z 0 /2 with a peak structure that is driven by ∂ ln f l ∂pT , according to the Eq.
. However for the conductivity σ el = 0.046 fm −1 such a peak structure disappears, again Eq. (22) allows to understand it; in fact in this case τ 0 B y (τ 0 ) ≃ τ 1 B y (τ 1 ) and hence the α coefficient, multiplying ∂ ln f l ∂pT , becomes quite small. The fittings with our formula are shown by the dashed lines in Fig. 11 agree with the numerical results quite well. The α factor for σ el = 0.046 fm −1 becomes about 50 times than the case for σ el = 0.0115 fm −1 , essentially because the difference in tB y (t) at formation time and escape time are nearly equal. The α ratio of leptons from Z 0 decay and charm quarks, which is about 4.7/3.6 = 1.3 and 8.7/6.3 = 1.38 for electrical conductivity 0.023 fm −1 and 0.0115 fm −1 separately, is close to their charge ratio (implying a quite similar value of the K factor), which confirms the strong correlation between the ∆v 1 of charm and leptons, once subtracting the impact of the very different p T slope of the spectrum.
We close our analysis with a final consideration. As shown in Eq. (7), the peculiar spectra of leptons from Z 0 decay should leave fingerprints in both the spectra f (p T , φ, y z ) and ∆v n of positively and negatively charged particles, as long as a n , b n and c n are non-zero. For example, with the help of the first two lines of Eq. (7) and by knowing that a n , b n and c n do not depend on p T when p T ≫ m, their spectra after the effects of e.m. fields become according to Eq. (7): where a 1 and b 1 become non-zero when E x and E y at η s = 0 have non-zero cos φ x and sin φ x terms respectively. We thus also studied the ratio of the spectra of positively and negatively leptons from Z 0 in 5.02 TeV collisions at b = 7.5 fm using σ el = 0.046 fm −1 , and the results are shown in Fig. 12. It is seen that the ratio, f + /f − − 1, is driven by the term − ∂ ln f ∂pT of leptons shown in Fig. 10, though the ratio is very close to 1 (deviated by 10 −3 ), which means that a 1 and b 1 are non-zero. On the other hand, if one looks at the spatial configurations of E x and E y [26,30,62], one should immediately identify large dipole moments and conclude that they should be nonzero.

V. CONCLUSIONS AND DISCUSSIONS
In this study we have obtained the general formula of the charge dependent flow observables generated by e.m. fields, which has a simple form ∆v n (p T , y z ) = − ∆dn(yz) 2 ∂ ln f ∂pT − ∆en(yz) 2pT at high p T according to Eq. (9) where the specific impact of strong interactions is subdominant. An experimental check of the p T pattern it predicts for the splitting of matter/anti-matter ∆v 1 and the correlations between the charm meson directed flow and the one of leptons from Z 0 decay (yet to be measured) would provide a strong probe of e.m. fields. The coefficients in the formula have a direct relation to the expansion of e.m. fields in φ x , and the Lorentz force in the longitudinal direction contributes also to the charge dependent flow observables that measure the anisotropy in transverse momenta. Moreover the strength of our formula is to trace back the sign and the strength of the splitting ∆v 1 to time dependence of the magnetic field B y (t), in the center of the colliding system, including also the effect of the electric field E x generated by the Faraday's law. This has been obtained under the approximation that the space gradients of the electromagnetic field can be discarded within the core of the QGP fireball created in AA collisions. The formula derived allows to understand that the intial strength of the magnetic field does not determine the sign and/or the strength of ∆v 1 . Furthermore clarify also the relation between the ∆v 1 of charm and bottom, and moreover the one with the leptons from Z 0 decay that appears quite different due to the very different p T dependence of the spectra. To confirm the validity of the formula, we have compared it to the realistic numerical simulation in a relativistic Langevin approach, finding a very good agreement between them. By comparing the numerical results between charm and bottom quarks, we show how the formula serves as a powerful tool to understand the results. Finally, we pointed out the strong correlation between the coefficients a n , b n and c n of charm quarks and leptons from Z 0 decay, where a n , b n and c n are the respective harmonic expansions of the mean variations of the three momenta due to e.m. fields as seen by Eq. (4), which is believed to hold for any configuration of e.m. fields.
The present study utilizes the Peterson fragmentation converting heavy quarks into heavy mesons, and so the shape does not modify much from quarks, that can not be directly probed, to mesons. However, the hadronization mechanism by coalescence plus fragmentation [43,53,55] may modify it quantitatively since the coalescence model combines one heavy quark with light quarks of different p T into mesons by a non random selection of the bulk matter along the hypersurface of hadronization. However, if one looks at the high p T behavior discussed along the present work, the modification should be small since the coalescence contribution significantly decreases with p T .
The present paper is presenting a test of a pocket formula for ∆v 1 by comparison to realistic simulations of AA collisions but under e.m. space-time profile coming from the assumption of the existence of an equilibrated QGP matter at constant σ el . However the analysis could be extended also to other e.m. profiles that currently under consideration for example in the study of the chiral magnetic effect [63]. Our study shows that the value of the magnetic field at the very early time, t < 0.1 fm/c, can significantly modify the relation between the ∆v 1 of D and B mesons due to the relevance of the value of tB y (t) at the the formation time of charm and bottom quarks.