Electromagnetic fields from the extended Kharzeev-McLerran-Warringa model in relativistic heavy-ion collisions
EElectromagnetic fields from the extended Kharzeev-McLerran-Warringa modelin relativistic heavy-ion collisions
Yi Chen,
1, 2
Xin-Li Sheng, and Guo-Liang Ma
4, 1, ∗ Shanghai Institute of Applied Physics, Chinese Academy of Sciences, Shanghai 201800, China University of Chinese Academy of Sciences, Beijing 100049, China Key Laboratory of Quark and Lepton Physics (MOE) and Institute of Particle Physics,Central China Normal University, Wuhan 430079, China Key Laboratory of Nuclear Physics and Ion-beam Application (MOE),Institute of Modern Physics, Fudan University, Shanghai 200433, China
Based on the Kharzeev-McLerran-Warringa (KMW) model that estimates the strong electromag-netic (EM) fields generated in relativistic heavy-ion collisions, we generalize the formulas of EMfields in the vacuum by incorporating the longitudinal position dependence, the generalized chargedensity distribution and retardation correction. We further generalize the formulas of EM fieldsin the QGP medium by incorporating a constant Ohm electric conductivity. Using the extendedKMW model, we observe a slower time evolution and a more reasonable impact parameter b de-pendence of the magnetic field strength than those from the original KMW model in the vacuum.The constant electric conductivity from lattice QCD data helps further prolong the time evolutionof magnetic field, so that the magnetic field strength at the thermal freeze-out time matches therequired magnetic field strength for the explanation of the observed difference in global polarizationsof Λ and ¯Λ hyperons in Au+Au collisions at RHIC. These generalized formulations in the extendedKMW model could be potentially useful for many EM fields relevant studies in relativistic heavy-ioncollisions at lower colliding energies and for various species of colliding nuclei. I. INTRODUCTION
The study of strongly interacting matter and its properties in the presence of strong electromagnetic (EM) fieldshas been a hot topic for more than a decade [1–34]. The ultra-relativistic heavy-ion collisions provide a unique way forcreating and exploring the strongly interacting matter at extremely high temperature and high energy density, wherethe matter is expected to be deconfined and reach a new state of matter, which is so-called the “quark-gluon plasma”(QGP) [35–38]. Properties of strongly interacting matter are governed by quantum chromodynamics (QCD), whichhave been widely and experimentally studied both on the Relativistic Heavy Ion Collider (RHIC) at BrookhavenNational Laboratory (BNL) and on the Large Hadron Collider (LHC) at CERN. In heavy-ion collisions, e.g., at topRHIC energy of √ s = 200 GeV or top LHC energy of √ s = 5 .
02 TeV, the two oppositely fast moving (almost at thespeed of light) colliding nuclei in non-central nucleus-nucleus (A-A) collisions can generate hitherto the strongest EMfields [1–34], which are usually estimated to be at the order of magnitude of eB ∼ eE ∼ m π ∼ G at top RHICenergy, or ∼ G at LHC energies, where m π is the pion mass. Here we note that the event-by-event fluctuationsof generated EM fields have been widely studied in Refs. [9–11, 15, 20–22, 28–33] in recent years, which can give riseto non-vanishing components of EM fields such as | B x | and | E x | due to the fluctuations of proton positions in thetwo colliding nuclei. In the QCD vacuum, topologically non-trivial gauge field configuration with non-zero windingnumber Q w [39–41] of deconfined QGP in the presence of such a strong magnetic field B will induce a non-conservedaxial current j µ = (cid:80) f q f (cid:104) ¯ ψ f γ µ γ ψ f (cid:105) A as well as a vector current j µ = (cid:80) f q f (cid:104) ¯ ψ f γ µ ψ f (cid:105) V along the direction ofmagnetic field, which are respectively called the “charge separation effect” (CSE) and the “chiral magnetic effect”(CME) [2, 6, 42–50].Since the axial current j µ requires a charge imbalanced C -odd environment while the vector current j µ requires achirality imbalanced P -odd environment, an asymmetry between the amount of positive and negative charges alongthe direction of magnetic field B in heavy-ion collisions is expected. Experimental observations of CME can beregarded as direct evidences of topologically non-trivial gluon field configurations, and furthermore can be interpretedas indications of event-by-event local P and CP violation of QCD at quantum level [2, 42]. Besides the CME and CSE,it is well known that the strong magnetic field could also influence many QCD processes [48], e.g. the induction ofchiral symmetry breaking [51], the influences on chiral condensation [52], and the modification of in-medium particlemass [53–59]. As an important consequence, the QCD phase diagram may be dynamically modified by such a strong ∗ Electronic address: [email protected] a r X i v : . [ nu c l - t h ] J a n magnetic field [60–64], e.g. color-superconducting phases at very high baryon densities could also be strongly affectedby the strong magnetic field [65]. Coupling with some anomaly processes, many interesting effects [48], e.g. theformation of π -domain walls [66], could also be generated by such a strong magnetic field.Due to the fact that the generated magnetic field B in heavy-ion collisions can not be directly event-by-eventmeasured, it is an enormous challenge to measure the magnetic field B induced chiral anomalous effects in currentexperimental measurements. In heavy-ion collisions, the magnetic field B is generated along the direction preferentiallyperpendicular to the reaction plane, hence experimental measurements are usually conducted using the two-particlecorrelators γ αβ = (cid:104)(cid:104) cos( φ α + φ β − RP (cid:105)(cid:105) firstly proposed by Voloshin [67], where α and β denote the electric chargesign of particles α and β , φ α and φ β are their azimuthal angles respectively, and Ψ RP is the azimuthal angle ofthe constructed reaction plane for a given event, and (cid:104)(cid:104)· · · (cid:105)(cid:105) denotes the average over all particle pairs and thenover all events. Therefore, the same-sign (SS) and opposite-sign (OS) correlators can be respectively written as γ ss ≡ ( γ ++ + γ −− ) / γ os ≡ ( γ + − + γ − + ) /
2. Based on [5, 67], the magnetic field driven CME is expected tocontribute to a negative γ ss but a positive γ os . The STAR Collaboration [68–73] and ALICE Collaboration [74] haveindependently measured the γ ss and γ os correlators, which indeed show the expected features of CME. However,there exit some ambiguities [75–83] in the interpretation of the experimental data due to the large backgroundcontaminations, potentially arising from the elliptic-flow ( v ) driven background contributions, such as the transversemomentum conservation (TMC) [78, 82, 83], local charge conservation (LCC) [79–81]. Hence, a dedicated run ofthe Zr + Zr and Ru + Ru isobar colliding systems at RHIC has been proposed [5], which is expected toyield unambiguous evidence for the CME signal by varying the signal but with the v -driven backgrounds roughlyfixed [5, 84–89].In this paper, we will review some generic features of an analytical model for the estimation of EM fields generatedin relativistic heavy-ion collisions from the original work initiated by Kharzeev, McLerran and Warringa [2], whichwe refer to as the original KMW model. On the ground of it, we formulate our generalizations of the estimated EMfields. We first start from the generalization of charge distributions used in heavy-ion collisions, based on which wepoint out that the formulas of EM fields in the original KMW model can be properly extended by incorporating thelongitudinal position dependence through the generalized relativistic three-dimensional charge distribution models,e.g., three-parameter Fermi model, for both spherically symmetric colliding nuclei and axially deformed ones. Also, wemake the retardation correction (RC) to the estimated EM fields contributed by participants, such that our formulationof the estimated EM fields embedded with both generalized charge distributions and retardation correction (RC) canbe potentially applied to lower energy region, such as the current beam-energy-scan (BES) program at RHIC, theunder planning FAIR, NICA and J-PARC programs. It is because the Lorentz contraction is not so large thatthe “pancake-shaped disk” approximation used in the original KMW model [2] is no longer valid in these energyregions. Besides, we further extend the formulas of EM fields to incorporate the medium feedback effects accordingto the Faraday’s induction law, where a constant Ohm electric conductivity σ is properly and analytically embedded.Finally, we make some numerical evaluations of the generalized EM fields formulas in the extended KMW modelfor the detailed comparisons of time evolution, impact parameter b dependence as well as the prediction of possiblelifetime t B of the estimated magnetic field B ( t, r ) in comparison with those from the original KMW model.This paper is organized as follows. We present detailed formulations of the estimated EM fields in the extendedKMW model for heavy-ion collisions in Sec. II, which consists of three subparts. We first present in Sec. II A aformal generalization of the charge distributions from the widely used three-parameter Fermi model in which therelativistic Lorentz contraction effects on the geometries of both spherically formed and axially deformed collidingnuclei are taken into account. We then in Sec. II B make a generalization of EM fields in the vacuum from the widelyused Li´enard-Wiechert equations, and naturally generalize the EM fields in the vacuum from the original KMWmodel with the generalized charge distributions and retardation correction. We further extend the EM fields withmedium feedback effects by incorporating a constant Ohm electric conductivity σ into the extend KMW model forhigh-energy heavy-ion collisions in Sec. II C. Some evaluations and comparisons about the time evolution, centrality(impact parameter b ) dependence, as well as the prediction of possible lifetime of the estimated magnetic field fromthe extended KMW model in comparison with those from the original KMW model are presented in Sec. III. Wefinally summarize the main processes of such generalization along with the conclusions and discussions in Sec. IV.The notation we use in this paper is the rationalized Lorentz-Heaviside units within the natural units, with (cid:126) = c = 1and µ = (cid:15) = 1. II. ELECTROMAGNETIC FIELDS FROM EXTENDED KMW MODEL
Before the detailed discussion, let us first emphasize that the physical situations of heavy-ion collisions consist ofevent-by-event P and CP violation processes due to the effects of topological charge fluctuation with non-trivial QCDgauge field configurations (characterized by the topological invariant, the winding number Q w ) in the vicinity of thedeconfined QGP phase [2]. This kind of P and CP violation processes can only locally happen under some specialand even extreme conditions in QCD, such as in the instantaneous deconfined QGP phase with chirality imbalance atextremely high temperature T ∼ Λ QCD in the presence of a strong magnetic field eB ∼ m π , which will equivalentlygenerate event-by-event locally non-vanishing θ angle for the effects of so called “ θ -vacuum” [90–92]. In other words,the fluctuations of local P and CP violated metastable domains of the QCD θ vacuum on the event-by-event basis areintrinsically connected to the amount of charges Q separated by the simultaneous magnetic field B for the non-trivialgluon field configurations with Q w (cid:54) = 0. It means that the expectation value of the amount of separated charges Q is proportional to the strength of magnetic field | B | , therefore the charge asymmetry fluctuations ∆ ± will directlyproportional to | B | [2]. It is therefore of crucial importance to quanitfy the strength of EM fields, especially the timeevolution and the possible lifetime t B of the magnetic field B ( t, r ) in relativistic heavy-ion collisions [93].Similar to the conventional setup of colliding system for heavy-ion collisions, we choose the x axis along the impactparameter b , and z axis along the beam direction of projectile such that the x − z plane is exactly the reaction plane(RP). The y axis is then chosen to be perpendicular to the reaction plane, as illustrated in Fig. 1. PT FIG. 1: (Color online) Illustration of the initial geometry of the colliding system projected on the transverse plane ( z = 0) atan impact parameter b for non-central high-energy nucleus-nucleus collisions. The centers of the projectile (denoted as P) andtarget (denoted as T) nucleus are respectively located at ( b/ , ,
0) and ( − b/ , ,
0) at t = 0 when the two colliding nuclei arecompletely overlapping with each other, moving parallel or anti-parallel to the beam z direction. A. Generalization of Charge Distributions for Heavy-ion Collisions
Since in symmetric heavy-ion collisions such as gold-gold collisions at top RHIC energy with the center of massenergy √ s = 200 GeV per nucleon pair along the z direction, the Lorentz contraction factor γ = √ s/ (2 m N ) (cid:39) . m N is the average nucleon mass), which corresponds to the beam rapidity Y = cosh − ( γ ) (cid:39) .
36. In thiscase, the two gold nuclei will be Lorentz contracted to be less than 1% of their original size in the z direction, becauseof which the original KMW model [2] approximates the two colliding nuclei as “pancake-shape disk”. It is furtherassumed in this model that the charge distribution can be treated as uniformly distributed within the two collidingnuclei and further be limited only onto the transverse plane with a two-dimensional surface number density. Thereforein this paper, we will first formulate our generalization of the charge distributions, and hence formulate the estimationof generated EM fields in the extended KMW model.It has been widely acknowledged that final state anisotropies of emitted charged particles in heavy-ion collisionsare very sensitive to the initial conditions, such as the experimentally measured final state charge particle elliptic flow v is actually very sensitive to the initial eccentricity ε [94]. Hence, the initial conditions may dominantly cause thecontamination of the experimentally measured CME signal. Therefore, a better quantification of the initial geometrywill be of crucial importance, especially for the study of CME in heavy-ion collisions since the initial geometry ofcharge distributions may significantly modify not only the centrality dependence but also the time evolution of thegenerated magnetic field B ( t, r ).Besides, due to relativistic motion of colliding nuclei, the overall size of each nucleus will be Lorentz contractedby a factor of γ along the beam z direction, such that a spherically symmetric nucleus will become an ellipsoidalnucleus with the volume shrunken by a factor of γ while the charge density magnified by a factor of γ simultaneously.Hence, we first start from generalizing the charge distributions into a relativistic three-dimensional form with thecorresponding standard parameters obtained from the high-energy electron-scattering measurements [95]. Moreover,what the original high-energy electron-scattering experiments measure is actually the charge spatial distribution inthe so-called “Breit Frame” [96] rather that the nuclear density profile with relativistic motion. Hence, we notethat the originally inferred charge distribution parameters from electron-scattering experiments should therefore beproperly modified [94] so as to be unambiguously implanted into the widely used Monte-Carlo simulations of theinitial conditions for the finite-number nucleons sampling of each nucleus [97].For non-central high-energy nucleus-nucleus (A-A) collisions at center of mass energy √ s = 2 γm N per nucleon pairwith the same setup of colliding system as illustrated in Fig. 1, the charge distribution used in the original KMWmodel can be generalized from the widely used charge distribution models listed in [95], i.e. the widely accepted andused three-parameter Fermi model (3pF) for heavy ions, with the relativistic Lorentz contraction effect taken intoaccount. For a source point r (cid:48) = ( x (cid:48) , y (cid:48) , z (cid:48) ) = ( r (cid:48)⊥ , z (cid:48) ) within the colliding nucleus (with charge Ze , and e = | e | ) withits center located at r (cid:48) c = (0 , , ρ ( r (cid:48) ) = γρ ω [ r (cid:48) /R (cid:48) ( θ, φ )] (cid:8) [ r (cid:48) − R (cid:48) ( θ, φ )] /d (cid:9) Θ( r (cid:48) ) , (1)where r (cid:48) = r (cid:48) ( γ ) ≡ (cid:113) r (cid:48)⊥ + ( γz (cid:48) ) . (2)In the expression of ρ ( r (cid:48) ) in Eq. (1), the polar and azimuthal angle variations of axially deformed R (cid:48) ( θ, φ ) is definedas [21, 84, 86] R (cid:48) ( θ, φ ) = R (cid:2) β Y ( θ ) + β Y ( θ ) (cid:3) , (3)where R is the spherical monopole charge radius, Y ml ( θ, φ ) denotes the spherical harmonics, with θ being the polarangle with respect to the symmetry axis of each nucleus, and φ being the corresponding azimuthal angle in eachnucleus, and they are defined as φ = sin − ( y (cid:48) / | r (cid:48)⊥ | ) and θ = cos − ( z (cid:48) /r (cid:48) ). The diffusion depth of the nuclear surface d can roughly be related to the skin thickness t by t = 4 ln(3) d (cid:39) . d , where t equals the distance in which the chargedensity falls from 90% to 10% of its central value ρ .Note that ω in Eq. (1) is an additional parameter [95], introduced to characterize the charge central depression orelevation phenomena of some nuclide such as Ca, relative to the traditional two-parameter Fermi model (2pF) withparameter c (which is in analogue to R (cid:48) in 3pF model) and d . Only in 3pF model with ω = 0 that the half densityradius R / = R (cid:48) ( θ, φ ), which is the distance from the center of the nucleus to the point at which the charge densityequals half its central value ρ . ρ is usually regarded as the charge number normalization constant, which can beobtained from the following charge number normalization condition Z = (cid:90) + R A − R A d x (cid:48) (cid:90) + R A − R A d y (cid:48) (cid:90) + R A /γ − R A /γ d z (cid:48) ρ ( r (cid:48) ) . (4)where R A is the radius of each nucleus. One should note that the γ factor in front of ρ in Eq. (1) is due to thefact that since the volume of each colliding nucleus will be Lorentz contracted by a factor of γ in the z direction, inaccordance to which the charge number density will be simultaneously enlarged by a factor of γ at the same time.The function Θ( r (cid:48) ) is introduced so as to constrain that the charge density is just within the ellipsoidal collidingnucleus, and it can further be used to separate the spectators from the participants when restricting it just withinthe transverse plane, namely Θ( r (cid:48) ) → ˜Θ( r (cid:48)⊥ ), which are separately defined asΘ( r (cid:48) ) = θ (cid:104) R − r (cid:48)⊥ − ( γz (cid:48) ) (cid:105) , ˜Θ( r (cid:48)⊥ ) = θ (cid:104) R − r (cid:48)⊥ (cid:105) . (5)Here we note that Θ( r (cid:48) ) and ˜Θ( r (cid:48)⊥ ) are similar to the function θ ( r (cid:48)⊥ ) introduced in the original KMW model [2]. Thefunction θ ( x ) in Eq. (5) is the standard unit step function, which reads θ ( x ) = (cid:40) x < x > . (6)For a spherically symmetric nucleus ( β = β = 0), i.e. the Au nucleus, R (cid:48) is actually independent of polar andazimuthal angle distributions and is just the monopole charge radius, namely R (cid:48) = R (similar to the parameter c in 2pF model), as a consequence of which charge number density ρ ( r (cid:48) ) = ρ ( r (cid:48) ) is only a function of r (cid:48) ( γ ). In therelativistic situation under this condition, the relativistic three-dimensional charge number density ρ ( r (cid:48) ) according toEq. (1) reads as ρ ( r (cid:48) ) = γρ ω ( r (cid:48) /R ) (cid:2) ( r (cid:48) − R ) /d (cid:3) Θ( r (cid:48) ) . (7)However, for a axially deformed nucleus, e.g. the Zr nucleus, its shape may deviate from the ideal sphere andthe charge distribution depends on the polar and azimuthal angle distributions, for which Eq. (3) is to some extentintroduced so as to account for such kind of charge distribution inside the 3pF model in Eq. (1). In the charac-terization of deformations in Eq. (3), β Y and β Y characterize the angular variations of R (cid:48) ( θ, φ ) of quadrupoleand hexadecapole deformations respectively, where β and β are the corresponding quadrupole and hexadecapoledeformation parameters, and Y ( θ ) and Y ( θ ) are the corresponding spherical harmonics. B. Electromagnetic Fields from Extended KMW Model in the Vacuum
Since we have generalized the charge number density in a relativistic three-dimensional form in Eq. (1) and Eq.(7), let us now first quote a simple situation where EM fields are generated by a pair of two constantly but oppositelymoving point-like charged particles of the same electric charge Ze (where e is the elementary electric charge, e = | e | )with rapidity Y n = ± Y ( Y > v n = (0 , , tanh Y n ), parallel or anti-parallel tothe z direction. We now set that at t = 0 the two point-like charged particles are located at r (cid:48) n ( t = 0) = ( x (cid:48) n , y (cid:48) n , z (cid:48) n ).Here the index n = ± respectively represents the charged particles moving in the positive and negative z -directions.Also, it is worth to note that since in the three-dimensional situation in heavy-ion collisions for a charge point r (cid:48) n within the two colliding nuclei at time t = 0 when the two colliding nuclei are usually assumed to be completelyoverlapping with each other, the z coordinate z (cid:48) n for this charge point in general will be non-zero. One can start fromusing the Li´enard-Wiechert potentials generated by this two point-like charge particles to evaluate the EM fields fora field point r = ( x, y, z ) at observation time t , which gives rise to the following widely used L-W equations of EMfields [3, 5, 8–11, 15, 16, 20–22, 28–33] e E ( t, r ) = α EM (cid:88) n = ± Z n R n (cid:0) − v n (cid:1)(cid:16) R n − [ v n × R n ] (cid:17) / ,e B ( t, r ) = α EM (cid:88) n = ± Z n v n × R n (cid:0) − v n (cid:1)(cid:16) R n − [ v n × R n ] (cid:17) / , (8)where R n = r − r (cid:48) n ( t ) are the relative positions of the field point r to the source points r (cid:48) n at the same observationtime t , r (cid:48) n ( t ) = ( x (cid:48) n , y (cid:48) n , z (cid:48) n ) + v n t are the source positions of the two point-like charge particles at time t , and α EM isthe EM fine-structure constant, defined as α EM = e / π ≈ / Y , in the following form, e E ( t, r ) = Zα EM (cid:88) n = ± cosh( Y n ) · R n (cid:110) ( x − x (cid:48) n ) + ( y − y (cid:48) n ) + [ t sinh Y n − ( z − z (cid:48) n ) cosh Y n ] (cid:111) / ,e B ( t, r ) = Zα EM (cid:88) n = ± sinh( Y n ) · e z × R n (cid:110) ( x − x (cid:48) n ) + ( y − y (cid:48) n ) + [ t sinh Y n − ( z − z (cid:48) n ) cosh Y n ] (cid:111) / , (9)which in nature are well consistent with the derivation in [2] for the magnetic field generated by a point-like charge whenreducing z (cid:48) n → z (cid:48) n (cid:54) = 0 is intended for the three-dimensional form of charge distributionsuch as Eq. (1) or (7) where the z coordinate of a point-like charge inside the colliding nuclei is not necessarily limitedonly within the transverse plane, hence we can generalize the “pancake shaped disk” approximation as well as thetwo-dimensional surface number densities used in the original KMW model [2].Now one can estimate the strength of EM fields generated in heavy-ion collisions. For a field point r = ( x, y, z ) =( r ⊥ , z ) at observation time t , the EM fields generated by two identical colliding nuclei with charge Ze and beamrapidity ± Y ( Y >
0) parallel or anti-parallel to the z direction can be estimated by combining the generalizedrelativistic three-dimensional charge number density ρ ( r (cid:48) ) in Eq. (1) or Eq. (7) with the EM fields generated by apair of two point-like charges in Eq. (8) or Eq. (9), and then performing the integration of charge number densitiesover the corresponding coordinate space within the two colliding nuclei. Nevertheless, the situation will be a littledifferent from the point-like charge case since colliding two heavy ions will separate the charged bulk matter intospectators (denote as S) and participants (denote as P). The participants in general will slow-down to a certain extentdue to the baryon-junction stopping effect [100], while the spectators in general are assumed to be moving with thesame beam rapidity Y b = ± Y and do not scatter at all. Therefore, we follow the method proposed in [2] by likewisesplitting the contributions to the EM fields generated by the two colliding nuclei into four parts, which simply readas E = E +S + E − S + E +P + E − P , B = B +S + B − S + B +P + B − P . (10)For a charge position vector r (cid:48) = ( x (cid:48) , y (cid:48) , z (cid:48) ) = ( r (cid:48)⊥ , z (cid:48) ) with the same setup of colliding system as in Fig. 1, thecharge number density ρ ( r (cid:48) ) in Eq. (1) or Eq. (7) as well as the corresponding Θ( r (cid:48) ) and ˜Θ( r (cid:48)⊥ ) functions should beaccordingly shifted by a half impact parameter vector b , namely ρ ± ( r (cid:48) ) = ρ ( r (cid:48) ∓ b / ± ( r (cid:48) ) = Θ( r (cid:48) ∓ b /
2) and˜Θ ± ( r (cid:48)⊥ ) = ˜Θ( r (cid:48)⊥ ∓ b / e E ± S ( t, r ) = α EM cosh( Y ) (cid:90) V ± d r (cid:48) ρ ± ( r (cid:48) ) (cid:104) − ˜Θ ∓ ( r (cid:48)⊥ ) (cid:105) R ± (cid:110) ( r ⊥ − r (cid:48)⊥ ) + [ t sinh( ± Y ) − ( z − z (cid:48) ) cosh Y ] (cid:111) / ,e B ± S ( t, r ) = α EM sinh( ± Y ) (cid:90) V ± d r (cid:48) ρ ± ( r (cid:48) ) (cid:104) − ˜Θ ∓ ( r (cid:48)⊥ ) (cid:105) e z × R ± (cid:110) ( r ⊥ − r (cid:48)⊥ ) + [ t sinh( ± Y ) − ( z − z (cid:48) ) cosh Y ] (cid:111) / , (11)where V ± donate the ellipsoidal volumes occupied by the two colliding nuclei, which can be related to the integrationdomains in Eq. (4) for charge number normalization condition.For EM fields contributed by participants, we at present do not consider contributions by newly created particleseither as [2] because the numbers of newly produced positively and negatively charged particles are nearly equaland their expansion is estimated to be almost spherical. Such an assumption should be reasonable, especially forperipheral collisions. Therefore one only needs to take into account the contributions of baryon-junction stoppingbulk matter which are initially there. The normalized rapidity Y distribution of the baryon-junction stopping bulkmatter of projectile (+) and target ( − ) can be empirically estimated [2, 17, 27, 100] as f ± ( Y ) = α y α y Y ) e ± α y Y , − Y ≤ Y ≤ Y . (12)Here we note that the parameter α y can be estimated as α y (cid:39) α , where α is the Regge trajectory intercept, definedas α (0) (cid:39) α B (0) − α R (0)) (cid:39) / α B (0) (cid:39) α R (0) (cid:39) . α B (0) should be used if nucleon exchange effectively dominates, which will somehow give rise to a relativelysmaller α y . Recently, the ALICE Collaboration has reported the data at different center of mass energies √ s forthe anti-baryons to baryons ratio on baryon stopping [101], which seems to support that α y (cid:39) . α y (cid:39) . α y between 0.48-0.50 does not largely affect the final results for the time evolution of the magnetic field strength.The contributions to EM fields by participants can therefore be estimated by inserting Eq. (12) into Eq. (9), andthen combining with Eq. (1), which follow as e E ± P ( t, r ) = α EM (cid:90) V ± d r (cid:48) (cid:90) Y − Y d Y Ψ ± ( Y ) cosh Y · ρ ± ( r (cid:48) ) ˜Θ ∓ ( r (cid:48)⊥ ) R ± (cid:110) ( r ⊥ − r (cid:48)⊥ ) + [ t sinh Y − ( z − z (cid:48) ) cosh Y ] (cid:111) / ,e B ± P ( t, r ) = α EM (cid:90) V ± d r (cid:48) (cid:90) Y − Y d Y Ψ ± ( Y ) sinh Y · ρ ± ( r (cid:48) ) ˜Θ ∓ ( r (cid:48)⊥ ) e z × R ± (cid:110) ( r ⊥ − r (cid:48)⊥ ) + [ t sinh Y − ( z − z (cid:48) ) cosh Y ] (cid:111) / , (13)where the refined functions Ψ ± ( Y ) for baryon-junction stopping effect are introduced to account for whether theretardation time t ret of a participating charge at observation time t is ahead of the collision time t c or simply afterthe collision time t c , which readsΨ ± ( Y ) = θ ( t ret − t c ) f ± ( Y ) + θ ( t c − t ret ) δ ( Y ∓ Y ) . (14)Here t ret is obtained by solving the retardation relation t ret = t − (cid:12)(cid:12) r − r (cid:48) − t ret · e z tanh Y (cid:12)(cid:12) for each rapidity Y withinthe normalized rapidity integral. We refer to this kind of treatment of participants in Eqs. (13-14) as retardationcorrection (RC) in this paper. For the case without retardation correction, the functions Ψ ± ( Y ) are simply previouslyused normalized rapidity distribution functions f ± ( Y ) in [2, 17, 27], namely Ψ ± ( Y ) = f ± ( Y ).The expressions for the estimated EM fields in Eqs. (11) and (13) are in general consistent with the derivation inthe original KMW model [2]. One distinct difference is that the depth in the z direction is now taken into accountin the extended KMW model rather than the infinitely thin depth (zero z -depth) in the original KMW model,which actually extends the dimensions of the charge density integration. The other distinct difference is that theretardation correction (RC) from the perspective of field propagation for the treatment of participating charges asin Eq. (14) through Ψ ± ( Y ) rather that f ± ( Y ) in [2, 17, 27] is now explicitly considered. Hence, we think that ourgeneralization in extended KMW model may give rise to some observable differences for the total estimated EM fieldssince the retardation due to field propagation is most relevant to the z coordinate. Again, we note that the chargenumber densities ρ ± ( r (cid:48) ) are not longer the uniformly distributed surface densities, which however are generalized intorelativistic 3pF model forms for the incorporation of Lorentz contraction effect and the effect due to deformationswith the standard parameters obtained from high-energy electron-scattering measurements [95]. Therefore, Eqs. (11)and (13) will be very easily and appropriately applied to lower energy regions and various colliding systems, wherethe “pancake-shaped disk” approximation used in the original KMW model [2] is no longer valid. C. Electromagnetic Fields from Extended KMW Model with Medium Feedback Effects
Up to now, we limit our discussions on the EM fields without considering any medium feedback effects, namely theOhm electric conductivity σ = σ Ohm = 0 in above discussion for the extended KMW model. For high-temperature( ∼ Λ QCD ) hot QCD matter, the real QGP is a conducting medium and the electric conductivity σ is generallyacknowledged to be proportional to the plasma temperature T [102–108], namely σ ∝ T . Since it is known accordingto the Faraday’s induction law that the electric conductivity σ may (to some extent) substantially slow down thedecrease of the generated EM fields, and therefore largely prolong the lifetime of EM fields, which is verified tobe crucial for the estimation of CME in heavy-ion collisions [93, 109]. We therefore implant a constant electricconductivity σ into the following Maxwell’s equations so as to incorporate the medium feedback effects, ∇ · E = ρ ext ε , ∇ · B = 0 , ∇ × E = − ∂ t B , ∇ × B = ∂ t E + σ E + j ext . (15)We then obtain the following partial differential wave equations for EM fields, (cid:0) ∇ − ∂ t − σ∂ t (cid:1) B = −∇ × j ext , (cid:0) ∇ − ∂ t − σ∂ t (cid:1) E = ∂ t j ext + ∇ (cid:16) ρ ext ε (cid:17) . (16)The solutions for above EM wave equations have been analytically solved in [23], in which the authors have embeddedboth electric conductivity σ and chiral magnetic conductivity σ χ for the EM fields generated by a relativistic point-likecharge. In the absence of chiral magnetic conductivity σ χ (since in general, σ χ (cid:28) σ ), we obtain the following solutionin the cylindrical coordinates for the magnetic field generated at a field point r = ( r ⊥ , z ) = ( x, y, z ) as B r B φ B z ( t, r ) = Q π γvR ⊥ ∆ / (cid:20) γ | v | σ √ ∆ (cid:21) e A , (17)with R ⊥ = | r ⊥ − r (cid:48)⊥ | = (cid:112) ( x − x (cid:48) ) + ( y − y (cid:48) ) , ∆ = γ ( vt + z (cid:48) − z ) + R ⊥ ,A = γvσ γ ( vt + z (cid:48) − z )] − γ | v | σ √ ∆ , (18)for a point-like charged particle with electric charge Q = + Ze moving with velocity v = (0 , , + v ) along the positive z direction, which is located at r (cid:48) = ( r (cid:48)⊥ , z (cid:48) ) = ( x (cid:48) , y (cid:48) , z (cid:48) + vt ) at observation time t . The radial and longitudinalcomponents of magnetic field B r and B z are both vanishing when only embedded with the electric conductivity σ ,namely B r = B z = 0.The transformations of EM fields F (where F = B , E ) in the Cartesian coordinates ( F x , F y , F z ) with that in thecylindrical coordinates ( F r , F φ , F z ) simply read F = F x F y F z = F r cos φ − F φ sin φF r sin φ + F φ cos φF z , (19)where φ is the azimuthal angle (with respect to x -axis) of the transverse relative position vector, namely R ⊥ = r ⊥ − r (cid:48)⊥ .Note that since we have e z × R = ( y (cid:48) − y, x − x (cid:48) , R = r − r (cid:48) is the relative position vector, the magneticfield B in Eq. (17) in the Cartesian coordinates therefore can be rewritten in a more compact form, which reads B ( t, r ) = Q π γv e z × R ∆ / (cid:20) γ | v | σ √ ∆ (cid:21) e A , (20)which shares almost the same form as the L-W form of magnetic field in Eq. (9) except for the factor (1+ γ | v | σ √ ∆ / e A due to the inclusion of electric conductivity σ for the QGP medium response. Therefore, Eq. (20) can naturallyreproduce the well-known L-W form of magnetic field in the vacuum, i.e., Eq. (8) or Eq. (9), when the electricconductivity σ is vanishing. For a special situation in Eq. (17) when the azimuthal angle φ = 0, the x component ofmagnetic field B x is consequently vanishing and the tangential φ component of magnetic field B φ is actually B y andalong the y direction, as is usually mentioned in literature [14]. Also, we note that the expression of B in Eq. (17)was firstly derived in Ref. [17].For the electric field E , the tangential component of electric field E φ can be easily obtained by applying the relation E φ = − vB r = 0. Meanwhile, the authors in [23] have obtained the following explicit expressions for radial as well aslongitudinal components E r and E z in the leading linear order of σ for high-energy heavy-ion collisions as E r ( t, r ) = Q π (cid:26) γR ⊥ ∆ / (cid:18) γ | v | σ √ ∆ (cid:19) − σ | v | R ⊥ (cid:20) γ | v |√ ∆ (cid:18) t − z − z (cid:48) v (cid:19)(cid:21) e − σ [ t − ( z − z (cid:48) ) /v ] (cid:27) e A ,E φ ( t, r ) = 0 ,E z ( t, r ) = Q π (cid:26) σ v | v | e − σ [ t − ( z − z (cid:48) ) /v ] Γ(0 , − A ) + e A ∆ / (cid:20) γ ( z − z (cid:48) − vt ) − v | v | A √ ∆ − γvσv ∆ (cid:21)(cid:27) , (21)where Γ(0 , − A ) = (cid:82) ∞− A d t exp( − t ) /t is the incomplete gamma function. According to the transformations in Eq. (19),the radial and longitudinal components of electric field E r and E z are respectively along the x and z directions whenthe azimuthal angle φ = 0. Here we note that Eqs. (17) and (21) can naturally reduce to the well-known L-Wequations for EM fields in the vacuum when the electric conductivity σ is vanishing ( σ = 0). It has been furtherproved in [23] that Eqs. (17) and (21) are self-consistent and well satisfy the Maxwell’s equations in Eq. (15). Hencein the following, we will generalize above EM fields for point-like charge in Eqs. (17) and (21), and make them moreapplicable for high-energy nucleus-nucleus collisions.We combine the electric conductivity σ embedded EM fields in Eqs. (17) and (21) with the relativistic three-dimensional charge number density in Eq. (1) or Eq. (7), and then apply the decomposition of EM fields as inEq. (10). In analogue to Eqs. (11) and (13), we finally obtain the following explicit expressions for magnetic fieldembedded with a constant electric conductivity σ . For the contributions of spectators, we obtain e B ± S ( t, r ) = lim Y →± Y α EM sinh Y (cid:90) V ± d r (cid:48) ρ ± ( r (cid:48) ) (cid:104) − ˜Θ ∓ ( r (cid:48)⊥ ) (cid:105) e z × R ± ∆ / (cid:20) σ sinh | Y | √ ∆ (cid:21) e A , (22)with the definitions in Eq. (18) accordingly rewritten as R ⊥ = | r ⊥ − r (cid:48)⊥ | = (cid:112) ( x − x (cid:48) ) + ( y − y (cid:48) ) , ∆ = [ t sinh Y − ( z − z (cid:48) ) cosh Y ] + R ⊥ ,A = σ sinh Y t sinh Y − ( z − z (cid:48) ) cosh Y ] − σ sinh | Y | √ ∆ , (23)where R ± = r − r (cid:48) ( t ) are the relative positions of the field point r to the source point r (cid:48) at the observation time t .Again, the positive and negative signs ± correspond to the nuclei moving in the positive and negative z -directions.For the contributions of participants, we have e B ± P ( t, r ) = α EM (cid:90) V ± d r (cid:48) (cid:90) Y − Y d Y Ψ ± ( Y ) sinh Y · ρ ± ( r (cid:48) ) ˜Θ ∓ ( r (cid:48)⊥ ) e z × R ± ∆ / (cid:20) σ sinh | Y | √ ∆ (cid:21) e A . (24)For the electric field E generated by two colliding nuclei, we first notice that E φ is vanishing in the point-like chargecase, therefore there is no contribution from both spectators and participants for the tangential component electricfield E φ . The cylindrical coordinates are usually not extensively used in heavy-ion collisions, we therefore directlytransform the expressions for electric field into the Cartesian coordinates after using Eq. (19) as follows, (cid:18) eE ± x, S eE ± y, S (cid:19) ( t, r ) = lim Y →± Y α EM (cid:90) V ± d r (cid:48) ρ ± ( r (cid:48) ) (cid:104) − ˜Θ ∓ ( r (cid:48)⊥ ) (cid:105) (cid:26) R ⊥ cosh Y ∆ / (cid:18) σ sinh | Y | √ ∆ (cid:19) − σR ⊥ tanh | Y | (cid:20) | Y |√ ∆ (cid:18) t − z − z (cid:48) tanh Y (cid:19)(cid:21) exp (cid:18) − σ (cid:20) t − z − z (cid:48) tanh Y (cid:21)(cid:19) (cid:27) e A R ⊥ (cid:18) x − x (cid:48) y − y (cid:48) (cid:19) , (cid:18) eE ± x, P eE ± y, P (cid:19) ( t, r ) = α EM (cid:90) V ± d r (cid:48) (cid:90) Y − Y d Y Ψ ± ( Y ) ρ ± ( r (cid:48) ) ˜Θ ∓ ( r (cid:48)⊥ ) (cid:26) R ⊥ cosh Y ∆ / (cid:18) σ sinh | Y | √ ∆ (cid:19) − σR ⊥ tanh | Y | (cid:20) | Y |√ ∆ (cid:18) t − z − z (cid:48) tanh Y (cid:19)(cid:21) exp (cid:18) − σ (cid:20) t − z − z (cid:48) tanh Y (cid:21)(cid:19) (cid:27) e A R ⊥ (cid:18) x − x (cid:48) y − y (cid:48) (cid:19) , (25)which are respectively contributions of spectators and participants for the x and y components of electric field E x and E y . For the z component of electric field E z contributed by spectators and participants, we have eE ± z, S ( t, r ) = lim Y →± Y α EM (cid:90) V ± d r (cid:48) ρ ± ( r (cid:48) ) (cid:104) − ˜Θ ∓ ( r (cid:48)⊥ ) (cid:105) (cid:26) sgn( Y ) σ tanh Y exp (cid:18) − σ (cid:20) t − z − z (cid:48) tanh Y (cid:21)(cid:19) Γ(0 , − A )+ e A ∆ / (cid:20) ( z − z (cid:48) ) cosh Y − t sinh Y − sgn( Y ) A √ ∆ − σ sinh Y tanh Y ∆ (cid:21) (cid:27) ,eE ± z, P ( t, r ) = α EM (cid:90) V ± d r (cid:48) (cid:90) Y − Y d Y Ψ ± ( Y ) ρ ± ( r (cid:48) ) ˜Θ ∓ ( r (cid:48)⊥ ) (cid:26) sgn( Y ) σ tanh Y exp (cid:18) − σ (cid:20) t − z − z (cid:48) tanh Y (cid:21)(cid:19) Γ(0 , − A )+ e A ∆ / (cid:20) ( z − z (cid:48) ) cosh Y − t sinh Y − sgn( Y ) A √ ∆ − σ sinh Y tanh Y ∆ (cid:21) (cid:27) . (26)Here the sign function sgn( Y ) = Y / | Y | denotes the sign of rapidity Y . Note that Eqs. (22-26) quantify the x , y and z components of EM fields contributed by spectators and participants in heavy-ion collisions with a constant electricconductivity σ implemented, where we also incorporate the generalized relativistic charge distributions, e.g., the 3pFmodel in Eqs. (1) or (7), and the retardation correction (RC) through Ψ ± ( Y ). On the one hand, we currently in Eqs.(22-26) do not include the chiral magnetic conductivity σ χ into the Maxwell’s equations in Eq. (15) since in general σ χ (cid:28) σ . Also, the value of electric conductivity σ from recent lattice calculations [102–108] still shows very largeuncertainties, which may easily submerge the contributions from chiral magnetic conductivity σ χ . On the other hand,as has shown in [23] that even when the chiral magnetic conductivity σ χ is included, the dominative components suchas B φ , E r and E z will not be changed at all. Therefore, we hold the point of view that Eqs. (22-26) can well quantifythe dominant contributions of medium feedback effects for the time evolution and centrality (impact parameter b )dependence of EM fields generated in heavy-ion collisions.Here we note that the chiral magnetic conductivity σ χ in general will be complex due to the spatially anti-symmetricpart of the off diagonal photon polarization tensor, as has been formally discussed and evaluated in [44] by using thelinear response theory. It has also been shown in [44] that the chiral magnetic conductivity σ χ has very strongfrequency ω and temperature T dependence, and the induced current j ( t ) can even change from positive to negativein the time evolution. We therefore currently postpone the inclusion of the chiral magnetic conductivity σ χ into theextended KMW model, but leave it to our future study.Before we move forward, let us make some short remarks here. The merits of above generalization for the chargenumber densities as well as the formulations of estimated EM fields both in the vacuum in Eqs. (10-14) and in themedium in Eqs. (22-26) in the extended KMW model can be at least boiled down to two points: one is that it is morerealistically akin to the physical situations for charge distributions with relativistic motion, for which both sphericalformed and axially deformed charge number densities are explicitly elaborated, e.g., Eqs. (1) or (7). Therefore, the0formulations of estimated EM fields based on the generalized charge number densities can be easily applied to variouscolliding systems, such as the asymmetric Cu + Au colliding system or the ongoing STAR experiments for isobarcolliding systems in which the Zr and Ru nuclei are widely regarded as an axially deformed nucleus [84–87]. Herewe note that a lot of discussions and the related estimations for the EM fields generated in isobar collisions have beenpresented [30, 32, 84, 86, 88, 89, 110].Meanwhile, since the charge number density ρ ( r (cid:48) ) in Eqs. (1) or (7) is formally embedded with relativistic effectdue to Lorentz contraction, it will be much more appropriate and adaptable for the situation at very low energyregion when the Lorentz contraction effect is not that significantly large, such as at the lower RHIC energy with √ s = 7 . γ ∼
4. As a comparison, the “pancake shape disk” approximation used in [2]will only be a good approximation in high energy region where the Lorentz factor γ is substantially large. Thereforethe model proposed in this paper will be much more appropriate for applications to the STAR-BES program, andsome even lower energy region like the under planning FAIR, NICA and J-PARC programs.The other point is that since the estimation of generated EM fields are explicitly embedded with medium feedbackeffects through a constant electric conductivity σ and refined baryon-junction stopping effect with proper retardationcorrection in Eq. (14), the generalized formulas of EM fields from Eqs. (10-14) in the vacuum or Eqs. (22-26) in themedium hence can further be used for many EM fields related studies, such as the CME related charge asymmetryfluctuations like a ++ or a + − correlators [2, 68, 87, 110], CME current J = σ χ B and its related studies for chiralmagnetic conductivity σ χ [44], in-medium particle mass [53–59] as well as the QCD phase diagram under strongmagnetic field [60–65], and so on.Due to the ellipsoidal-type geometry and also the peculiar charge domains occupied by spectators and participants,the numerical integrations of the charge distribution for the estimated EM fields from Eqs. (10-14) or Eqs. (22-26)are actually not very easily and quickly performed. Meanwhile, we notice that the L-W equations of EM fields in thevacuum, i.e. Eq. (8), are widely accepted and used. We therefore propose one possible simplification of the formulasof EM fields in the conducting QGP medium from Eqs. (22-26) by replacing the continuous integration of the chargedistribution with the discrete summation of all point-like charges, which turns out to be an alternative solution to doMonte-Carlo numerical simulations of generated EM fields from the L-W equations [3, 5, 8–11, 15, 20–22, 28–33] asfollows e B ± S ( t, r ) = lim Y →± Y α EM (cid:88) N S Z n sinh( Y ) · e z × R ± ∆ / (cid:20) σ sinh | Y | √ ∆ (cid:21) e A ,e B ± P ( t, r ) = α EM (cid:88) N P Z n (cid:90) Y − Y d Y Ψ ± ( Y ) sinh Y · e z × R ± ∆ / (cid:20) σ sinh | Y | √ ∆ (cid:21) e A , (27)for the decomposed magnetic field B contributed by spectators and participants. For that of the electric field E , weobtain the following expressions for x and y components of the electric field E x and E y , (cid:18) eE ± x, S eE ± y, S (cid:19) ( t, r ) = lim Y →± Y α EM (cid:88) N S Z n (cid:26) R ⊥ cosh Y ∆ / (cid:18) σ sinh | Y | √ ∆ (cid:19) − σR ⊥ tanh | Y | (cid:20) | Y |√ ∆ (cid:18) t − z − z (cid:48) tanh Y (cid:19)(cid:21) exp (cid:18) − σ (cid:20) t − z − z (cid:48) tanh Y (cid:21)(cid:19) (cid:27) e A R ⊥ (cid:18) x − x (cid:48) y − y (cid:48) (cid:19) , (cid:18) eE ± x, P eE ± y, P (cid:19) ( t, r ) = α EM (cid:88) N P Z n (cid:90) Y − Y d Y Ψ ± ( Y ) (cid:26) R ⊥ cosh Y ∆ / (cid:18) σ sinh | Y | √ ∆ (cid:19) − σR ⊥ tanh | Y | (cid:20) | Y |√ ∆ (cid:18) t − z − z (cid:48) tanh Y (cid:19)(cid:21) exp (cid:18) − σ (cid:20) t − z − z (cid:48) tanh Y (cid:21)(cid:19) (cid:27) e A R ⊥ (cid:18) x − x (cid:48) y − y (cid:48) (cid:19) , (28)1and also that for the z component of the electric field E z , eE ± z, S ( t, r ) = lim Y →± Y α EM (cid:88) N S Z n (cid:26) sgn( Y ) σ tanh Y exp (cid:18) − σ (cid:20) t − z − z (cid:48) tanh Y (cid:21)(cid:19) Γ(0 , − A )+ e A ∆ / (cid:20) ( z − z (cid:48) ) cosh Y − t sinh Y − sgn( Y ) A √ ∆ − σ sinh Y tanh Y ∆ (cid:21) (cid:27) ,eE ± z, P ( t, r ) = α EM (cid:88) N P Z n (cid:90) Y − Y d Y Ψ ± ( Y ) (cid:26) sgn( Y ) σ tanh Y exp (cid:18) − σ (cid:20) t − z − z (cid:48) tanh Y (cid:21)(cid:19) Γ(0 , − A )+ e A ∆ / (cid:20) ( z − z (cid:48) ) cosh Y − t sinh Y − sgn( Y ) A √ ∆ − σ sinh Y tanh Y ∆ (cid:21) (cid:27) . (29)Here the expressions for R ⊥ , ∆ and A are the same as the definitions in Eq. (23), and Ψ ± ( Y ) in Eq. (14). N S and N P denote the total numbers of spectators and participants of point-like charges, respectively. Also, we note that asimilar Monte-Carlo simulation using above formulas of EM fields for point-like charges in Eqs. (27-29) has firstly beenperformed in [23], since such an alternative formulation can be more directly obtained from the EM fields of point-likecharge in [23] when further embedded with the refined baryon-junction stopping effect with retardation correction asin Eq. (14). If one assumes that each proton at nucleon level in the two colliding nuclei can be treated as point-likecharge, as has been extensively assumed in literatures [3, 5, 8–11, 15, 20, 22, 28–33], then above generalization of EMfields in Eqs. (27-29) for Monte-Carlo simulations in heavy-ion collisions can have at least two merits: one is that theEM fields from Eqs. (27-29) are analytically embedded with a constant electric conductivity σ for the incorporation ofmedium effects rather than the original L-W equations in the vacuum; the other is that the EM fields in Eqs. (27-29)are also implanted with the refined baryon-junction stopping effect with retardation correction as in Eq. (14). III. RESULTS AND DISCUSSIONS
Since we have explicitly formulated the estimations of generated EM fields for heavy-ion collisions in the vacuum inEqs. (10-14) or in the medium in Eqs. (22-26), as well as an alternative solution to the currently used Monte-Carlosimulations of point-like charges in the medium in Eqs. (27-29), let us first give some pre-analysis before performingthe ab initio integration of charge distribution or event-by-event simulations of generated EM fields at specific space-time points. Due to the mirror and centrosymmetric symmetries of the colliding system (as is illustrated in Fig. 1),the total electric field E will vanish while the total magnetic field B will only remain its y component pointing in thenegative y direction at the center point r = of the overlapping region, which further results in a charge separationdue to the CME. Therefore, we only focus on the y component of magnetic field B y in this paper.Let us first declare some abbreviations that we will use in the following. The original KMW model is incorporatedwith a two-dimensional (2D) surface charge number density, which is abbreviated as 2D KMW model. The extendedKMW model is generalized by incorporating the longitudinal position dependence and generalized three-dimensional(3D) volume charge number density, which is likewise abbreviated as 3D KMW model but including two cases: with(w/ ) and without (w/o) retardation correction (RC). The collision time t c for simplicity is chosen at fully overlappingtime, namely t c = 0. Besides, we do check that even when one considers the collision time t c totally from the geometryconsideration regardless of the causality reason according to the theory of special relativity, namely the collision time t c can be estimated as t c = − t d , where t d is a departure time defined in Eq. (30), our results basically will not bechanged due to a much earlier retardation time t ret .It should also be noted that the formulation of estimated magnetic field in the original KMW model is actually only inthe vacuum without retardation correction. In order to make a better comparison of the two models in the conductingmedium, we also show the generalized 2D KMW model embedded with a constant Ohm electric conductivity σ , whichcan be easily obtained from the reduction of extended KMW model by removing the longitudinal position dependenceand replacing the corresponding charge number density ρ . Such a generalization of the formulas of the magnetic fieldin the conducting medium actually is consistent with that in Refs. [17, 27].Moreover, since it is mainly the magnetic field strength that has been numerically evaluated in the original KMWmodel [2], and also the following papers [17, 27]. For a comparison consideration, we only show the results of estimatedmagnetic field in this paper. The results of estimated electric field can be similarly evaluated from the extended KMWmodel in the same way as we present here.2 A. Time Evolution of the Magnetic Field Strength
Let us now make a sample comparison of the time evolution of the total magnetic field strength estimated fromEq. (10) in the vacuum. In Fig. 2, we show the estimation of time evolution of the total magnetic field strength e B ( t, r ) in the vacuum at the center point r = of the overlapping region in Au+Au collisions at √ s = 200 GeVwith different impact parameters b = 4 , ,
12 fm. From the inserted plot in the upper panel, one can clearly see thatthe magnetic field from 3D KMW model with (w/ ) and without (w/o) retardation correction (RC) at t = 0 willsystematically yield a relatively smaller magnetic field strength (especially for larger impact parameter b ), comparedwith that too sharp peak of magnetic field strength around t ∼ t = 0 actually are indeed consistent with the numerical simulations from HIJING model in [10],which will be clearly demonstrated in the next subsection. -4 -3 -2 -1
3D Vacuum 2D 3D -4 -3 FIG. 2: (Color online) Comparison of the time evolution of the total magnetic field strength in the vacuum from extendedKMW model with (w/ ) and without (w/o) retardation correction (RC) and that from original KMW model at different impactparameters b = 4 , ,
12 fm for the center point r = of the overlapping region in Au+Au collisions at √ s = 200 GeV. The twodashed horizontal lines indicate the later-time constraints of the magnetic field strength estimated in [93] from the differencebetween global polarizations of Λ and ¯Λ hyperons in Au+Au collisions at √ s = 200 GeV in [112] at one standard deviation(1 σ ) and three standard deviations (3 σ ), respectively. However, compared with that “seemly enhanced” magnetic field strength from the extended (3D) KMW modelwithout retardation correction, it is also noticeable that the original KMW model will generally give rise to a weakermagnetic field strength just shortly ( ∼ .
03 fm / c ) after the collision, and the differences in between become rathernoticeable at later times for more central (or smaller impact parameter b ) collisions. Thus the extended KMWmodel without retardation correction will yield a longer lifetime t B of the magnetic field than the original KMWmodel. When incorporated with retardation correction, however, the magnetic field strength from the extended (3D)KMW model will firstly decrease very fast, and then decrease slowly and even grows up to a certain extent forsome time range at some impact parameters, e.g., b = 4 , f ± ( Y ) in Eq. (14) are basically ruled out by the θ ( t ret − t c )function during the comparison of retardation time t ret with the collision time t c . The contributions from chargedparticipants are mainly dominated by those with beam rapidity ± Y , hence the magnetic field strength at early stagedecreases much faster than those without retardation correction; At middle or later stage time evolution, however,the contribution from baryon-junction stopping bulk matter grows dominant since the spectators have fly away. Theresulting magnetic field strength contributed by those participants will increase first and then decrease. Hence theresulting total magnetic field strength will show non-monotonic behavior, especially for smaller impact parameter b -3 -2 -1
3D Medium 2D 3D
FIG. 3: (Color online) Comparison of the time evolution of the total magnetic field strength in the medium with constant Ohmelectric conductivity σ = 5 . b = 4 , ,
12 fm for the center point r = of theoverlapping region in Au+Au collisions at √ s = 200 GeV. The two dashed horizontal lines indicate the later-time constraintsof the magnetic field strength estimated in [93] from the difference between global polarizations of Λ and ¯Λ hyperons in Au+Aucollisions at √ s = 200 GeV in [112] at one standard deviation (1 σ ) and three standard deviations (3 σ ), respectively. The electric conductivity σ has been calculated by using a first principal lattice QCD approach in the quenchedapproximation in Ref. [105], according to which we set the electric conductivity σ = 5 . e B ( t, r ) in the conducting medium with σ = 5 . r = of the overlapping region in Au+Au collisions at √ s = 200 GeV with different impactparameters b = 4 , ,
12 fm. One can clearly see that the magnetic field strength in the conducting QGP medium casein Fig. 3 decays much slower than that in Fig. 2 in the vacuum case. There clearly show non-monotonic behaviors ofthe magnetic field strength at early time stage for different impact parameters in the conducting QGP medium ratherthan the monotonically decreasing behaviors of magnetic field strength in the vacuum. Besides, compared with theoriginal KMW model (2D), we find the extended KMW model (3D) will in general systematically yield a relativelylarger magnetic field strength in the conducting QGP medium except at very early stage around t ∼
0. Also, we noticethat the differences between any two cases become more significant around the non-monotonic peaks ( t ∼ . / c ),especially for smaller impact parameter b .Moreover, we notice that the time evolution of the estimated magnetic field strength B ( t, r ) in the conductingmedium at the later-time stage evolution in Fig. 3 is quiet different from that in the vacuum case in Fig. 2 for theimpact parameter b dependence. This can be intuitively understood as follows: since in the vacuum case there is nomedium feedback effect due to the vanishing electric conductivity σ , the magnetic field decays quickly and there isonly the baryon-junction stopping effect that may to some extent hinder the moving charged bulk matter from flyingaway from the field point r = . Thus the baryon-junction stopping effect will dominate in the later stage of timeevolution especially for smaller impact parameter b due to larger bulk matter of the participants, the resulting totalmagnetic field strength therefore becomes larger for smaller impact parameter b at later time stage. In the conductingmedium, however, since the electric conductivity σ is non-negligible, it is the Faraday’s induction effect rather thanthe baryon-junction stopping effect that dominates for the later time evolution. The resulting total magnetic fieldstrength in the conducting medium therefore shows rather distinct impact parameter b dependence at later time stagefrom that in the vacuum. We also notice that the original (2D) KMW model and the extended (3D) KMW model can4give almost identical results after t (cid:38) . / c , and the difference in between grows little during the time evolutionwhich is quite different from that in the vacuum case.The “kinks” marked as the points at different impact parameters b for the extended (3D) KMW model in Fig. 3are due to the fact that the two colliding nuclei start to separate from each other. The corresponding vertical linesindicate the departure time t d when the two colliding nuclei just depart from each other, which can be estimated fromthe following geometry relation in the x − z plane (RP) at an impact parameter b as follows t d = t d ( b ) ≡ R A γ (cid:115) − (cid:18) b R A (cid:19) . (30)The original (2D) KMW model due to the same reason actually also have “kinks”, which however are all located at t = 0 (we do not show in Fig. 3), since the charges in the original KMW model are restricted within the transverseplane with zero depth in the z direction.By using the interesting idea proposed in the recent Ref. [93], one can use the experimentally measured polarizationdata [111, 112] on Λ and ¯Λ hyperons to make some constraints on the allowed later-time magnetic field strength, whichhas been also recently explored in Ref. [113]. It has been estimated in [93] that in the one standard deviation (1 σ ) limitof the recent measurements of the global polarization of Λ and ¯Λ hyperons in Au+Au collisions at √ s = 200 GeV [112], the magnetic field strength at thermal freeze-out can be estimated as e | B | = eT s | ∆ P H | | µ Λ | < . × − m π . (31)where T s ≈
150 MeV is the temperature of the emitting source, ∆ P H = P Λ −P ¯Λ is the difference in global polarizationsof Λ and ¯Λ hyperons, and µ Λ = − µ ¯Λ = − . µ N . In the three standard deviations (3 σ ) limit, the magnetic fieldstrength is estimated as e | B | < . × − m π [93]. Here we choose m π as the mass of π rather than that of π ± ,which is consistent with the pion mass that we use throughout this paper. Therefore, we can immediately use thesepossible constraints to compare with the magnetic field strength estimated from the extended KMW model and thatfrom the original KMW model both in the vacuum in Fig. 2 and in the conducting medium in Fig. 3, which couldenable us to make some estimations of the possible lifetime t B of the generated magnetic field in Au+Au collisions at √ s = 200 GeV.By comparing with the possible constraints shown by two dashed horizontal lines at 1 σ and 3 σ limits in Fig. 2 andFig. 3, we notice that the magnetic field strength in the vacuum in Fig. 2 at a possible lifetime t s (cid:39) /c of the QGPin Au+Au collisions at √ s = 200 GeV [114], is extremely smaller than the 1 σ bound of the possible constraint fromthe polarization data. This means that although we have improved the estimation of time evolution of the magneticfield strength in the vacuum in the extended KMW model, which results in a prolonged lifetime of magnetic fieldcompared with that from the original KMW model, but the strength of the magnetic field at thermal freeze-out time t = t s still appears quite inadequate (about two orders of magnitude smaller) for explaining the observed differencebetween global polarizations of Λ and ¯Λ hyperons. Hence, this indicate that the generated QGP in heavy-ion collisionsneed to maintain the magnetic field stronger for a longer lifetime than that in the vacuum, which further requiresthat the generated QGP is more like a conducting medium.In Fig. 3, as our expectation above, we surprisingly notice that the magnetic field strength in the conducting mediumwith a constant electric conductivity σ = 5 . < t < /c almost lies between the 1 σ and 3 σ bounds. Wefind that the magnetic field strength at impact parameter b = 4 −
12 fm lie between eB (cid:39) (2 . − . × − m π at t = t s (cid:39) /c , which meets the required condition of magnetic field strength for the explanation of observeddifference between global polarizations of Λ and ¯Λ hyperons. Our results support that if the generated QGP is aconducting medium, the enough strength of magnetic field can be reached at the thermal freeze-out time. B. Impact Parameter Dependence of the Magnetic Field Strength
Let us now compare the impact parameter b dependences of magnetic field strength in the vacuum and in theconducting QGP medium.In Fig. 4, we show the impact parameter b dependences of the estimated magnetic field strength in the vacuum fromthe extended (3D) KMW model and from the original (2D) KMW model with (w/ ) and without (w/o) retardationcorrection (RC) at the center point r = of the overlapping region in Au+Au collisions at √ s = 200 GeV at t = 0.Meanwhile, we also show the results from numerical simulations with HIJING model in [10, 50] and that from asample analytical approach (denote as LGK model) initially proposed in [16], both of which have been used for thecomparison of impact parameter dependence of magnetic field in Ref. [50].5 R a ti o HIJING3D LGK2D KMW w/o RC2D KMW w/ RC3D KMW w/o RC3D KMW w/ RC
Vacuum (b)(a)
FIG. 4: (Color online) (a) Comparison of impact parameter b dependence of the estimated total magnetic field strength in thevacuum at the center point r = of the overlapping region in Au+Au collisions at √ s = 200 GeV at t = 0 fm/ c from sixdifferent approaches. (b) The corresponding ratio (2D/3D) of the total magnetic field strength from 2D and 3D KMW modelswith (w/ ) and without (w/o) retardation correction (RC). The dashed vertical line indicates the boundary at b = 2 R A , abovewhich the two colliding nuclei will miss each other and there is no contribution from participants, as indicated by the coloredbands. The corresponding “kinks” are also marked by the points along the dashed vertical line. At first glance, one may notice that the magnetic field strength in the vacuum at smaller impact parameter b from the 3D KMW model seems to deviate a lot from the result obtained with HIJING model in [10]. Meanwhile,the LGK model seems to work well and be consistent with the HIJING result. But we should point out that inthe evaluation of generated magnetic field with HIJING model in [10], at t = 0 when the two colliding nuclei arecompletely overlapping with each other, all nucleons inside the two colliding nuclei are assigned with the same beamvelocity v z = 1 − (2 m N / √ s ) , and in the LGK model presented in [50] all the moving charges inside the two collidingnuclei are regarded as spectators and assigned with the beam velocity ± v z , but they both neglect the baryon-junctionstopping effect during the overlapping process. Unlike these two approaches, in the extended KMW model, we haveincluded the baryon-junction stopping effect through the experimentally supported f ± ( Y ) in Eq. (12) along with thegeneralized charge distribution. Therefore, it is not surprising that the simulation with HIJING model and the LGKmodel in [50] will roughly yield similar and consistent results, and both of them will give slightly larger magnetic fieldstrength, compared with the extended (3D) KMW model without (w/o) retardation correction (RC), especially atsmaller impact parameter b due to a larger charged bulk matter of participants which is affected by f ± ( Y ). When theretardation correction (RC) is included, however, the extended (3D) KMW model with (w/ ) retardation correction(RC) give a slightly larger magnetic field strength at larger impact parameter b than the HIJING and LGK models,which can be mainly attributed to the generalized charge distribution for the incorporation of Lorentz contractioneffect on the geometries of the two colliding nuclei used in the extended KMW model.In Fig. 4, one may notice the “kinks” at the boundary b = 2 R A from both original and extended KMW modelswhere the two colliding nuclei begins to miss each other and there is no contribution from participants (which willbe demonstrated in the following Fig. 5), compared with the other two models that seem to smoothly across theboundary. The “kinks” are due to the fact that the charge number densities in both 2D and 3D KMW models areaccompanied with the step functions Θ ± so as to constrain that the charges are limited within the two colliding nuclei(with radius R A before the Lorentz contraction), i.e. Eq. (1). Here we note that the HIJING model uses an relativelylarge upper limit of the nucleus radius for the corresponding sampled protons inside, which actually largely exceedsthe generally supposed nucleus charge radius R c estimated from the nucleus root mean square (RMS) charge radius R rms measured in experiments [115] such as R A (cid:39) R c = (cid:112) / R rms , while the LGK model used in [50] also employsan relatively large nuclear charge radius R c in the treatment of effective charge Z eff .Besides, one can clearly notice that, compared with other three 3D models used in Fig. 4, the original (2D) KMW6 Spec.Part. w/o RCPart. w/ RC
Vacuum (2D)
Spec.Part. w/o RCPart. w/ RC
Vacuum (3D)(a)(c) (b)(d)
FIG. 5: (Color online) Comparison of impact parameter b dependence of the magnetic field strength in the vacuum contributedby spectators (S) in the upper panels (a) and (b), and that contributed by participants (P) in the lower panels (c) and (d)at the center point r = of the overlapping region in Au+Au collisions at √ s = 200 GeV at t = 0 fm / c from original (2D)KMW model (left panels) with that from extended (3D) KMW model (right panels). The dashed vertical lines indicate theboundary at b = 2 R A , above which the two colliding nuclei will miss each other and there is no contribution from participants,as indicated by the colored bands. The corresponding “kinks” are also marked by the points along the dashed vertical lines. model indeed largely overestimates the magnetic field strength at t = 0, as shown in Fig. 2 that the original (2D)KMW model will give too sharp peaks of magnetic field strength around t ∼ ∼ . .
5) at relatively smaller impact parameter b with (without) retardationcorrection (RC), especially at the boundary b = 2 R A by a factor of ∼ . . b dependence of the estimated magnetic fieldstrength in the vacuum contributed by spectators (S) and participants (P) from the original (2D) KMW model (leftpanels) with that from the extended (3D) KMW model (right panels) with (w/ ) and without (w/o) retardationcorrection (RC). From the upper panels of Fig. 5, one can clearly see that it is mainly the contributions of spectatorsin the original (2D) KMW model that cause the too shape peak of the estimated magnetic field strength (left),compared with that result from the extended (3D) KMW model (right), which is also consistent with observed tooshape peaks in the time evolution in Fig. 4 at t ∼ f ± ( Y ) in Eq. (14) is basically ruled out by the θ ( t ret − t c ) function during the comparison between retardation time t ret and collision time t c . The contributions of charged participants actually are all calculated with beam rapidity ± Y instead of the normalized rapidity distribution f ± ( Y ). Hence this kind of enhancement is quite significant at t = 0, especially for the 2D KMW model due to the more localized two-dimensional surface charge number density.It is quite noticeable that the magnetic field strength contributed by participants in the 3D KMW model for a middleimpact parameter b , e.g., b = 8 fm, is almost one third of that in the 2D KMW model for the case with retardationcorrection (RC).Different from the situation in the vacuum in Fig. 4, the centrality dependence of the generated magnetic fieldstrength in the conducting medium should be realistically evaluated at some time after the collision, since the formationof QGP medium certainly needs some time after the collision time t c for a more realistic consideration. At present,it is still not very clear how long it takes for the formation of QGP medium in heavy-ion collisions, so we estimate it7
2D w/o RC 2D w/ RC 3D w/o RC 3D w/ RC
Medium (KMW) R a ti o FIG. 6: (Color online) (a) Comparison of impact parameter b dependence of the total magnetic field strength in the conductingQGP medium at hydrodynamical (pre-) equilibration time t = τ = 0 . /c [114, 116] with a constant Ohm electric conductivity σ = 5 . r = of theoverlapping region in Au+Au collisions at √ s = 200 GeV. (b) The corresponding ratio (2D/3D) of the total magnetic fieldstrength from 2D and 3D KMW models with (w/ ) and without (w/o) retardation correction (RC). The dashed vertical lineindicates the boundary at b = 2 R A in Fig. 4, compared to which there is no “kink” at all since the two colliding nuclei havealready missed each other before t = τ and the corresponding shaded color bands are also not shown. from the widely used (pre-) equilibration time in hydrodynamical simulations, namely t = τ = 0 . /c [114, 116].In Fig. 6, we show the impact parameter b dependence of the estimated total magnetic field strength in theconducting QGP medium at t = τ = 0 . /c [114, 116] with a constant electric conductivity σ = 5 . r = in Au+Au collisions at √ s = 200 GeV. One can see that the difference between the original (2D) KMWmodel and the extended (3D) KMW model at t = τ in the conducting medium is not as significant as that at t = 0 inthe vacuum as shown in Fig. 4. From the lower panel of Fig. 6, we can also see that the extended (3D) KMW modelsystematically yields a relatively stronger total magnetic field strength than the original (2D) KMW model for theconducting medium at t = τ , which is distinctly different from the situation in the vacuum at t = 0 shown in Fig. 4.This is resulted from the delay of the generated magnetic field in the conducting medium according to the Faraday’sinduction law since a constant Ohm electric conductivity σ has been explicitly embedded into the Maxwell’s equationsin Eq. (15). Moreover, there is no longer any “kink” at t = τ in Fig. 6, since the two colliding nuclei have alreadydeparted from each other at the estimated t d in Eq. (30).Likewise, in Fig. 7 we show a detailed comparison of the impact parameter b dependence of the estimated magneticfield strength in the conducting QGP medium contributed by spectators (S) and participants (P) at t = τ =0 . /c [114, 116] with a constant electric conductivity σ = 5 . b dependent magnetic field strength contributed byspectators, although the result in the extended (3D) KMW model is slightly stronger. In the lower panels of Fig. 7, wefind that the retardation correction (RC) can only mildly modulate the impact parameter b dependent magnetic fieldstrength contributed by participants in the conducting medium, which is more noticeably stronger in the extended(3D) KMW model than that in the original (2D) KMW model, similar to the situation of spectators in the upperpanels of Fig. 7. These behaviors in the conducting medium, however, are distinctly different from those in thevacuum at t = 0 fm / c in Fig. 5 where the original (2D) KMW model will systematically yield a stronger magneticfield strength than the extended (3D) KMW model.8 Spec. Part. w/o RC Part. w/ RC
Medium (2D)
Spec. Part. w/o RC Part. w/ RC
Medium (3D)(a) (b)(c) (d)
FIG. 7: (Color online) Comparison of impact parameter b dependence of the magnetic field strength in the conducting QGPmedium contributed by spectators (S) in the upper panels (a) and (b), and that from participants (P) in the lower panels (c)and (d) at hydrodynamical (pre-) equilibration time t = τ = 0 . /c [114, 116] with a constant Ohm electric conductivity σ = 5 . r = fm of the overlapping region in Au+Au collisions at √ s = 200 GeV from original (2D)KMW model (left panels) with that from extended (3D) KMW model (right panels). The dashed vertical lines indicate theboundary at b = 2 R A in Fig. 4, compared to which there is no “kink” at all since the two colliding nuclei have already missedeach other before t = τ and the corresponding shaded color bands are also not shown. IV. SUMMARY
In summary, based on the Kharzeev-McLerran-Warringa (KMW) model [2] for the estimation of strong EM fieldsgenerated in heavy-ion collisions, we make an attempt to generalize the formulas of estimated EM fields in the originalKMW model, which eventually turns out to be the above formulations in the extended KMW model.We first start from generalizing the widely used charge distribution (or charge number density) models by incorpo-rating the Lorentz contraction effect for the ellipsoidal geometry of each colliding nucleus due to relativistic motion,e.g., the three-parameter Fermi (3pF) model, for both spherically symmetric colliding nuclei as well as axially deformedcolliding nuclei used in heavy-ion collisions. We then combine the widely used L-W equations of EM fields generatedby point-like charges with the generalized relativistic three-dimensional charge number densities. By employing thedecomposition method for the contributions of spectators and participants of the two colliding nuclei as that in theoriginal KMW model [2], we then naturally generalize the formulas of EM fields in the vacuum by incorporating thelongitudinal position dependence and retardation correction (RC) for heavy-ion collisions.Since the medium feedback effects in the high-temperature QCD matter or the conducting QGP medium maysubstantially modify the time evolution of the generated EM fields according to the Faraday’s induction law, weexplicitly and analytically embed a constant electric conductivity σ into the Maxwell’s equations for the incorporationof medium effects, and eventually formulate the estimation of generated EM fields with medium feedback effects in theextended KMW model for heavy-ion collisions. For simplicity, we also propose one possible kind of simplification of theformulas of EM fields in the conducting QGP medium, which turns out to be an alternative way to the currently usedMonte-Carlo simulations of generated EM fields from widely used L-W equations in heavy-ion collisions, where weincorporate both conducting medium effects through a constant Ohm conductivity σ and the refined baryon-junctionstopping effect with retardation correction (RC).Our formulations of the estimated EM fields in the extended KMW model will result in a slower time evolution(or a longer lifetime t B ) and a more reasonable impact parameter b dependence of the magnetic field strength inthe vacuum. The formulations of estimated EM fields embedded with medium feedback effects through a constantelectric conductivity σ = 5 . t = t s can meet the9required magnetic field strength for the explanation of the observed difference in global polarizations of Λ and ¯Λhyperons in Au+Au collisions at √ s = 200 GeV. This supports that if the generated QGP is a conducting QCDmedium, an enough strength of magnetic field can be maintained until the thermal freeze-out time.It should be emphasized that we only applied the generalization of charge distributions with the three-parameterFermi model (3pF) for Au+Au collisions at √ s = 200 GeV in this work, which is just a simple case of spheric nucleus,various colliding nuclei such as axially deformed isobaric nuclei (zirconium vs ruthenium experiment at RHIC) can besimilarly generalized in the same way as we present here. Also, the extended KMW model can be applied to variouscolliding energies, especially for the lower energy regions like the ongoing STAR-BES program and the under planningFAIR, NICA and J-PARC programs. Hopefully, these generalized formulations also could be employed for many EMfields relevant studies, such as the charge asymmetry correlators or fluctuations related to the experimental searchingfor the CME, in-medium particle mass, as well as the QCD phase diagram under strong magnetic field. Acknowledgments
The authors would like to gratefully thank Yu-Gang Ma, Jin-Hui Chen, Song Zhang and Qi-Ye Shou for theirstimulating and helpful discussions and comments, and Chen Zhong for maintaining high-quality performance ofcomputational facilities. Y. C. thankfully acknowledges the helpful and fruitful discussions with Wei-Tian Deng,Huan-Xiong Yang, Heng-Tong Ding, Xu-Guang Huang, and Qun Wang, as well as the timely discussions and helpsfrom Xin-Li Zhao, Xin Ai, Bang-Xiang Chen, Yu Guo, Xiao-Liang Xia, Hui Li and Xian-Gai Deng. Y. C. and G.-L.M. are supported by National Natural Science Foundation of China under Grants Nos. 11890714, 11835002, and11961131011, and the Key Research Program of the Chinese Academy of Sciences under Grant No. XDPB09. X.-L.S. is supported by National Natural Science Foundation of China (NSFC) under Grant Nos. 11890714 and 12047528. [1] J. Rafelski and B. Muller, Phys. Rev. Lett. , 517 (1976).[2] D. E. Kharzeev, L. D. McLerran, and H. J. Warringa, Nucl. Phys. A , 227 (2008).[3] V. Skokov, A. Y. Illarionov, and V. Toneev, Int. J. Mod. Phys. A bf 24, 5925 (2009).[4] Kirill Tuchin, Phys. Rev. C , 034904 (2010); Phys. Rev. C , 039903 (2011).[5] Sergei A. Voloshin, Phys. Rev. Lett. , 172301 (2010).[6] M. Asakawa, A. Majumder, and B. M¨uller, Phys. Rev. C , 064912 (2010).[7] V. Voronyuk, V. D. Toneev, W. Cassing, E. L. Bratkovskaya, V. P.Konchakovski, and S. A.Voloshin, Phys. Rev. C ,054911 (2011).[8] L. Ou and B. A. Li, Phys. Rev. C 84, 064605 (2011).[9] A. Bzdak and V. Skokov, Phys. Lett. B , 171 (2012).[10] Wei-Tian Deng and Xu-Guang Huang, Phys. Rev. C , 044907 (2012).[11] V. D. Toneev, V. P. Konchakovski, V. Voronyuk, E. L. Bratkovskaya, and W. Cassing, Phys. Rev. C , 064907 (2012).[12] V. D. Toneev, V. Voronyuk, E. L. Bratkovskaya, W. Cassing, V. P. Konchakovski, and S. A. Voloshin, Phys. Rev. C ,034910 (2012).[13] Kirill Tuchin, Phys. Rev. C , 024911 (2013).[14] Kirill Tuchin, Adv. High Energy Phys. , 490495 (2013).[15] John Bloczynski, Xu-Guang Huang, Xilin Zhang, and Jinfeng Liao, Phys. Lett. B , 1529 (2013).[16] Yunpeng Liu, Carsten Greiner, and Che Ming Ko, arXiv: 1403.4317 [nucl-th].[17] Umut G¨ursoy, Dmitri Kharzeev, and Krishna Rajagopal, Phys. Rev. C , 054905 (2014).[18] L. McLerran and V. Skokov, Nucl. Phys. A , 184 (2014).[19] B. G. Zakharov, Phys. Lett. B , 266 (2014).[20] Wei-Tian Deng and Xu-Guang Huang, Phys. Lett. B (2014) 296-302.[21] John Bloczynski, Xu-Guang Huang, Xilin Zhang, and Jinfeng Liao, Nucl. Phys. A , 85 (2015).[22] V. Roy and S. Pu, Phys. Rev. C , 064902 (2015).[23] Hui Li, Xin-li Sheng, and Qun Wang, Phys. Rev. C , 044903 (2016).[24] Robert Holliday, Ryan McCarty, Balthazar Peroutka, and Kirill Tuchin, Nucl. Phys. A , 406 (2017).[25] Balthazar Peroutka and Kirill Tuchin, Nucl. Phys. A , 64 (2017).[26] Duan She, Sheng-Qin Feng, Yang Zhong, and Zhong-Bao Yin, Eur. Phys. J. A , 48 (2018).[27] Umut G¨ursoy, Dmitri Kharzeev, Eric Marcus, Krishna Rajagopal, and Chun Shen, Phys. Rev. C , 055201 (2018).[28] Xin-Li Zhao, Yu-Gang Ma, and Guo-Liang Ma, Phys. Rev. C , 024910 (2018).[29] Xin-Li Zhao, Guo-Liang Ma, and Yu-Gang Ma, Phys. Lett. B , 413 (2019).[30] Xin-Li Zhao, Guo-Liang Ma, and Yu-Gang Ma, Phys. Rev. C , 034903 (2019). [31] Yi-Lin Cheng, Song Zhang, Yu-Gang Ma, Jin-Hui Chen, and Chen Zhong, Phys. Rev. C , 054906 (2019).[32] Jan Hammelmann, Alba Soto-Ontoso, Massimiliano Alvioli, Hannah Elfner, and Mark Strikman, Phys. Rev. C ,061901(R) (2020).[33] X. G. Deng and Y. G. Ma, Phys. Rev. C , 054610 (2020).[34] Bang-Xiang Chen, and Sheng-Qin Feng, Chin. Phys. C , 024104 (2020).[35] E. V. Shuryak, Phys. Rept. , 71 (1980).[36] P. Braun-Munzinger, K. Redlich, J. Stachel, in Quark Gluon Plasma 3 , ed. by R. C. Hwa, X.-N. Wang, World Scientific,Singapore (2004).[37] J. Adams et al . (STAR Collboration), Nucl. Phys. A , 102 (2005).[38] Adam Bzdak, ShinIchi Esumi, Volker Koch, Jin-feng Liao, Mikhail Stephanov, and Nu Xu, Phys. Rept. , 1 (2020).[39] A. A. Belavin, A. M. Polyakov, A. S. Shvarts and Yu. S. Tyupkin, Phys. Lett. B , 85 (1975).[40] G. ’t Hooft, Phys. Rev. Lett. , 8 (1976); Phys. Rev. D , 3432 (1976); Phys. Rev. D , 2199 (1978).[41] R. Jackiw and C. Rebbi, Phys. Rev. Lett. , 172 (1976).[42] Kenji Fukushima, Dmitri E. Kharzeev, and Harmen J. Warringa, Phys. Rev. D , 074033 (2008).[43] Harmen J. Warringa, J. Phys. G: Nucl. Part. Phys. (2008) 104012.[44] Dmitri E. Kharzeev and Harmen J. Warringa, Phys. Rev. D , 034028 (2009).[45] Kenji Fukushima, Dmitri E. Kharzeev, and Harmen J. Warringa, Nucl. Phys. A , 311 (2010).[46] Kenji Fukushima, Dmitri E. Kharzeev, and Harmen J. Warringa, Phys. Rev. Lett. , 212001 (2010).[47] Dmitri E. Kharzeev, Ann. Phys. , 205 (2010).[48] Dmitri Kharzeev, Karl Landstriner, Andreas Schmitt, and Ho-Ung Yee (ed.), Lect. Notes Phys., vol. , (2013).[49] Dmitri E. Kharzeev, Prog. Part. Nucl. Phys. , 133 (2014).[50] Koichi Hattori and Xu-Guang Huang, Nucl. Sci. Tech. , 26 (2017).[51] V. P. Gusynin, V. A. Miransky and I. A. Shovkovy, Phys. Rev. Lett. , 3499 (1994); Phys. Rev. Lett. , 1005 (1996).[52] A. J. Mizher, M. N. Chernodub, and E. S. Fraga, Phys. Rev. D , 105016 (2010).[53] Hao-Lei Chen, Kenji Fukushima, Xu-Guang Huang, and Kazuya Mameda, Phys. Rev. D , 104052 (2016).[54] C. Bonati, M. D’Elia, M. Mariti, M. Mesiti, F. Negro, A. Rucci and F. Sanfilippo, Phys. Rev. D , 074515 (2017).[55] G. S. Bali, B. B. Brandt, G. Endrodi and B. GlaBle, Phys. Rev. D , 034505 (2018).[56] M. Coppola, D. GAmez Dummc, and N. N. Scoccola, Phys. Lett. B ,155 (2018).[57] Shijun Mao, Phys. Rev. D , 056005 (2019).[58] E. J. Ferrer and A. Sanchez, Phys. Rev. D , 096006 (2019).[59] H.-T. Ding, S.-T. Li, A. Tomiya, X.-D. Wang,and Y. Zhang, arXiv: 2008.00493 [hep-lat].[60] K. Fukushima and H. J. Warringa, Rev. Mod. Phys. , 025001 (2016).[61] Stefan Rechenberger, Phys. Rev. D , 054013 (2017).[62] Massimo D’Elia, Floriano Manigrasso, Francesco Negro, and Francesco Sanfilippo, Phys. Rev. D , 054509 (2018).[63] Guo-yun Shao, Wei-bo He, and Xue-yan Gao, Phys. Rev. D , 014020 (2019).[64] Shile Chen, Jiaxing Zhao, and Pengfei Zhuang, arXiv: 2005.08473 [nucl-th].[65] J. O. Andersen, W. R. Naylor and A. Tranberg, Phys. Rev. Lett. , 032007 (2008).[66] D. T. Son and M. A. Stephanov, Phys. Rev. D , 014021 (2008).[67] S. A. Voloshin, Phys. Rev. C , 057901 (2004).[68] I. V. Selyuzhenkov, (STAR Collaboration), Rom. Rep. Phys. , 049 (2006).[69] B. I. Abelev et al . (STAR Collaboration), Phys. Rev. Lett. , 251601 (2009).[70] B. I. Abelev et al . (STAR Collaboration), Phys. Rev. C , 054908 (2010).[71] S. A. Voloshin, (STAR Collaboration), Indian J. Phys. , 1103 (2011).[72] G. Wang, (STAR Collaboration), Nucl. Phys. A , 904 (2013).[73] L. Adamczyk et al ., (STAR Collaboration), Phys. Rev. Lett. , 052302 (2014).[74] B. I. Abelev et al ., (ALICE Collaboration), Phys. Rev. Lett. , 012301 (2013).[75] F. Wang, Phys. Rev. C , 064902 (2010).[76] A. Bzdak, V. Koch, and J. Liao, Phys. Rev. C , 031901 (2010).[77] J. Liao, V. Koch, and A. Bzdak, Phys. Rev. C , 054902 (2010).[78] A. Bzdak, V. Koch, and J. Liao, Phys. Rev. C , 014905 (2011).[79] S. Pratt, arXiv: 1002.1758 [nucl-th].[80] S. Schlichting and S. Pratt, arXiv: 1005.5341.[81] S. Schlichting, S. Pratt, Phys. Rev. C , 014913 (2011).[82] S. Pratt, S. Schlichting, and S. Gavin, Phys. Rev. C , 024909 (2011).[83] A. Bzdak, V. Koch, and J. Liao, Lect. Notes Phys. , 503 (2013).[84] Wei-Tian Deng, Xu-Guang Huang, Guo-Liang Ma, and Gang Wang, Phys. Rev. C , 041901(R) (2016).[85] V. Koch, S. Schlichting, V. Skokov, P. Sorensen, J. Thomas, S. Voloshin, G. Wang, and H.-U. Yee, Chin. Phys. C ,072001 (2017).[86] Xu-Guang Huang, Wei-Tian Deng, Guo-Liang Ma, and Gang Wang, Nucl. Phys. A , 736 (2017).[87] Shuzhe Shi, Yin Jiang, Elias Lilleskov, and Jinfeng Liao, Ann. Phys. , 50 (2018).[88] Hao-jie Xu, Xiaobao Wang, Hanlin Li, Jie Zhao, Zi-Wei Lin, Caiwan Shen, and Fuqiang Wang, Phys. Rev. Lett. ,022301 (2018).[89] Shu-zhe Shi, Hui Zhang, De-fu Hou, and Jin-feng Liao, Phys. Rev. Lett , 242301 (2020).[90] Dmitri Kharzeev, Phys. Lett. B , 260 (2006). [91] Dmitri Kharzeev and Ariel Zhitnitsky, Phys. Lett. B , 67 (2007).[92] C. G. Callan, R. F. Dashen, D. J. Gross, Phys. Lett. B , 334 (1976).[93] Berndt M¨uller and Andreas Sch¨afer, Phys. Rev. D , 071902(R) (2018).[94] Q. Y. Shou, Y. G. Ma, P. Sorensenc, A. H. Tang, F. Videbæk, and H. Wang, Phys. Lett. B , 215 (2015).[95] H. De Vries, C. W. De Jager, and C. De Vries, At. Data Nucl. Data Tables , 495 (1987).[96] F. J. Ernst, R. G. Sachs and K. C. Wali, Phys. Rev. , 1105 (1960).[97] M. L. Miller, K. Reygers, S. J. Sanders, and P. Steinberg, Annu. Rev. Nucl. Part. Sci. , 205 (2007).[98] G. C. Rossi and G. Veneziano, Nucl. Phys. B , 507 (1977);[99] G. C. Rossi and G. Veneziano, Phys. Rep. 63, (1980).[100] D. Kharzeev, Phys. Lett. B , 238 (1996).[101] E. Abbas et al . (ALICE Collaboration), Eur. Phys. J. C , 2496 (2013).[102] P. B. Arnold, G. D. Moore, and L. G. Yaffe, J. High Energy Phys. , 100 (2003).[103] S. Gupta, Phys. Lett. B , 57 (2004).[104] G. Aarts, C. Allton, J. Foley, S. Hands, and S. Kim, Phys. Rev. Lett. , 022002 (2007).[105] H. T. Ding, A. Francis, O. Kaczmarek, F. Karsch, E. Laermann, and W. Soeldner, Phys. Rev. D , 034504 (2011).[106] A. Francis and O. Kaczmarek, Prog. Part. Nucl. Phys. , 212 (2012).[107] Alessandro Amato, Gert Aarts, Chris Allton, Pietro Giudice, Simon Hands, and Jon-Ivar Skullerud, Phys. Rev. Lett. , 172001 (2013).[108] B. B. Brandt, A. Francis, H. B. Meyer, and H. Wittig, J. High Energy Phys. , 34 (2013).[109] E. Stewart and K. Tuchin, Phys. Rev. C , 044906 (2018).[110] Yin Jiang, Shuzhe Shi, Yi Yin, and Jinfeng Liao, Chin. Phys. C , 011001 (2018).[111] L. Adamczyk et al . (STAR Collaboration), Nature (London) , 62 (2017).[112] J. Adam et al . (STAR Collaboration), Phys. Rev. C , 014910 (2018).[113] Yu Guo, Shu-Zhe Shi, Sheng-Qin Feng, and Jin-Feng Liao, Phys. Lett. B , 134929 (2019).[114] Hui-Chao Song and Ulrich Heinz, Phys. Rev. C , 064901 (2008).[115] I. Angeli and K. P. Marinova, At. Data Nucl. Data Tables , 69 (2013).[116] Peter F. Kolb, and Ralf Rapp, Phys. Rev. C67