The fundamental structures and interactions of all visible matter in our universe are well described by the Standard Model of particle physics together with gravitation. In the Standard Model, gauge fields arising from underlying gauge symmetry principles provide the strong, weak and electromagnetic forces that hold the physical world together. A fascinating aspect of gauge fields is related to gauge field configurations with nontrivial topological windings. They emerge in various physical systems of different dimensions and play crucial roles in many important phenomena. Well-known examples include magnetic vortices in superconductors, monopoles in the electroweak theory, as well as instantons and sphalerons in non-Abelian gauge theory. See reviews in e.g. [1][2][3][4].
Let us focus on the instantons and sphalerons of non-Abelian gauge theories in four-dimensional spacetime [5,6]. These topological configurations constitute crucial ingredients for our very existence. On the very small scale of subatomic dynamics, they lead to the spontaneous chiral symmetry breaking in the strong interaction vacuum [2,3,7,8] and thus provide the dominant source of mass for visible matter. On the very large scale of cosmic evolution, they allow the violation of baryon number conservation and thus possible emergence of the large matter-anti-matter asymmetry in today's universe [9,10]. The essence of such topological configurations is the tunneling transitions across energy * liaoji@indiana.edu barriers between the topologically distinct vacuum sectors of a non-Abelian gauge theory characterized by different Chern-Simons numbers. In doing so, they themselves "twist" topologically around spacetime boundary and can be characterized by their topological winding numbers, defined as:
where the integrand q(x) is the local topological charge density of a given gauge field configuration described by gauge field strength tensor G µν (x) with g being the corresponding coupling constant.
Despite their significance, the topological configurations are elusive experimentally. A direct detection of their presence and consequences in laboratories could substantially advance our understanding of the underlying tunneling mechanism that is at the heart of quark confinement and baryon asymmetry, but has not been achieved so far. A concrete proposal [11][12][13][14] toward this goal is to look for the parity-odd "bubbles" (i.e. local domains) arising from the topological transitions of gluon fields in Quantum Chromodynamics (QCD). Specifically, these bubbles could occur in the hot quark-gluon plasma (QGP) created by relativistic heavy ion collisions. The parity-odd nature of such a bubble can be quantified by the macroscopic chirality N 5 generated for the light quarks in the plasma. Indeed, this is enforced by the famous chiral anomaly relation for massless fermions (i.e. the light flavor quarks in the case of QCD):
In the above J µ 5 is the local chiral or axial current for each quark flavor while the Eq.( 3) is the spacetime-integrated version of Eq.( 2), with N R and N L being the number of right-handed (RH) and left-handed (LH) quarks. This latter equation has its deep mathematical root in the celebrated Atiyah-Singer index theorem and physically means that each topological winding generates two units of net chirality per flavor of light quarks. Therefore, measuring the net chirality of QGP provides a unique way of directly accessing the fluctuations of gluon field topological windings in heavy ion collision experiments. Note though the winding number Q w in principle can not be an accurately measurable observable according to the Wigner-Araki theorem, for which the practical constraint could nevertheless be made arbitrarily small [15]. What one can hope for is to quantify the average winding number that occurs in the fireball spacetime volume from event-by-event topological fluctuations.
The quark net chirality, however, is also challenging to detect due to the fact that the QGP born from collisions would expand, cool down and eventually transition into a low temperature hadron phase where the spontaneous breaking of chiral symmetry makes the net chirality unobservable. Fortunately, there is a way out by virtue of the so-called chiral magnetic effect (CME) [16,17]. The CME is an anomalous transport phenomenon where an electric current J is induced along an external magnetic field B under the presence of net chirality in a system with massless fermions of charge Q e :
The µ 5 is a chiral chemical potential that quantifies the net chirality N 5 . The study of CME has attracted significant interests and activities from a broad range of physics disciplines such as high energy physics, condensed matter physics, astrophysics, cold atomic gases, etc. See recent reviews in e.g. [18][19][20][21][22]. In the context of heavy ion collisions, the CME current (4) leads to a charge separation in the quark-gluon plasma that results in a specific hadron emission pattern and can be measured via chargedependent azimuthal correlations [14]. Extensive experimental efforts have been carried out over the past decade to look for its signatures at the Relativistic Heavy Ion Collider (RHIC) and the Large Hadron Collider (LHC): see more details in recent reviews [23][24][25].
In short, there is a promising pathway for experimental probe of gauge field topology in heavy ion collisions: the winding number Q w of gluon fields ⇒ net chirality of quarks ⇒ CME current ⇒ correlation observables. Here we demonstrate, for the first time, how this strategy actually works quantitatively for counting the topological winding numbers of gauge field with measurements of the CME. To be precise, the Q w fluctuates from event to event with equal chance of being positive or negative and what one can hope for is to determine its average variance σ w ≡ Q 2 w event associated with such fluctuations over the fireball spacetime domain. This quantity σ w is ultimately related to the two-point functions of the topological charge density and is what one aims to extract in this study. This has become possible due to two developments. Experimentally, great progress has been made in data analysis methods to separate background contamination and extract the CME signal with much reduced uncertainty [24,25]. Theoretically, we have developed a new framework for accurately computing the CME transport while taking into account the stochastic dynamics of gluon field topological fluctuations during the evolution in these collisions.
The various ingredients and flow chart of our framework is illustrated in Fig. 1. This is built upon a stateof-the-art modeling tool, the anomalous-viscous fluid dynamics (AVFD) [26][27][28], which provides the quantitative link by simulating the CME transport during the dynamical evolution of a heavy ion collision from initial topological winding Q w to final experimental signal. However, random topological fluctuations occur both at the very beginning of a collision and during the course of its evolution. They result in flipping of chirality and impose important impact on the CME transport, which was previously neglected. To precisely count the initial topological windings, one must account for such missing key ingredient. In this work we've successfully implemented both the event-by-event initial topological fluctuations and the stochastic dynamics of gluon field topological fluctuations into the AVFD simulations (as indicated by the green blocks in Fig. 1), thus paving the way for counting Q w with CME measurements. See Method section for a more detailed discussion about this framework.
Here we present details of the various components (as illustrated in Fig. 1) for the computational framework developed for this study.
Heavy ion collisions at RHIC and LHC energies are well described by relativistic viscous hydrodynamic simulations, which have been thoroughly vetted with extensive experimental data. The "backbone" of our framework for describing such bulk evolution is based on the MUSIC (2+1)D code package [29,30] with initial time τ 0 = 0.4fm/c and shear viscosity parameter η/s = 0.08. The event-wise initial conditions (e.g. for entropy density profiles) of the bulk hydro are generated with Monte-Carlo Glauber simulations.
The anomalous-viscous fluid dynamics (AVFD) [26][27][28] is the key component for implementing the dynamical CME transport in the realistic environment of a relativistically expanding viscous fluid. This state-of-theart tool numerically solves the anomalous hydrodynamic equations for the coupled evolution of quarks' vector and axial currents together on top of the bulk collective flow. Starting from a given axial charge initial condition, the magnetic-field-induced CME currents eventually lead to a charge separation effect in the fireball, which turns into a dipole term ±a 1 sin(φ -Ψ RP ) in the azimuthal angle distribution of positive/negative charged hadrons with respect to the reaction plane orientation Ψ RP . Such a dipole signal is experimentally measurable through a difference between same-sign (SS) and opposite-sign (OS) charged hadron pair correlations, H SS-OS = 2a 2 1 . In short, the AVFD is a hydrodynamic realization of CME for quantifying its signal in heavy ion collisions. For more details, see refs. [26][27][28]. It may be noted that the lifetime of the magnetic field is a major source of theoretical uncertainty. To nail this down would require quantitative understanding of dynamically evolving magnetic field in the medium along the line of studies in e.g. [31,32]. The influence of this factor on CME signal is thoroughly investigated in [27]. We follow [27] to use a lifetime of 0.6 fm which is reasonable and supported by phenomenological analysis [33].
For the goal of this work, one crucial new ingredient has been introduced into our framework. During the hydrodynamic evolution, there exist random topological fluctuations of the gluon fields that would necessarily influence the axial current evolution. These fluctuations eventually amount to a relaxation effect toward equilibrium with vanishing topological charge on long time scale. To account for such effect, one needs to introduce the resulting relaxation term into the anomalous hydrodynamic equation for the axial current [34][35][36][37]:
In the above, J µ f,5 is the axial current of each quark flavor with electric charge Q f and color number N c = 3, while n f,5 is the corresponding axial charge density. The E • B in Eq.( 5) comes from the Abelian anomaly due to electromagnetic fields. The stochastic dynamics of gluon topological charge density q(x) gives rise to a new contribution i.e. the second term in Eq.( 5), where the τ cs is an important relaxation time for the topological charge fluctuations. Its physical meaning is simple: over this time scale, the q approaches equilibrium value. The τ cs is controlled by the Chern-Simons diffusion rate Γ cs , i.e.
In the above χ is the total quark number susceptibility, T is the temperature and N f = 3 is the light flavor number. Some discussions are now in order concerning the two key quantities here, i.e. τ cs and Γ cs . First, it is useful to compare τ cs with different time scales involved in the problem. The microscopic scale relevant for bulk equilibrium in heavy ion collisions is presumably on the order of τ micro ∼ 0.1fm/c while the magnetic field lifetime scale is on the order of τ B ∼ 1fm/c. In the scenario of τ cs τ micro , the topological fluctuations behave as random noise that would wash out any meaningful CME transport on macroscopic scale. On the other hand, if τ cs τ B then the relaxation of initial topological fluctuations would be too slow to have an impact on the CME transport. It is only when the τ cs is on a par with τ B that the stochastic dynamics of topological charges will play an important role. Therefore, a good estimate of τ cs is crucial, which however is highly nontrivial and strongly dependent on the estimate of Γ cs . For the latter, the best theoretical guidance is the perturbative result from [38]: Γ cs ≈ 30(α s T ) 4 . One has to recognize that the applicability of this result, when extrapolated to the regime relevant for heavy ion collisions, is just anyone's guess. It is therefore to clearly demonstrate the significant theoretical uncertainty present in this parameter and its consequences for the outcome of our calculations. In Table . I, we list the values of Γ cs and τ cs for a wide range of choices for α s based on the aforementioned formula. Needless to say, this formula may become totally inapplicable when α s increases beyond certain range. However, we notice that the resulting τ cs value (which is the actual input into our simulation framework) already spans a rather generous range of possibilities. Comparing simulation results from these different choices will help calibrate the theoretical uncertainties associated with the relaxation time scale of topological charges, as we shall discuss later in the Results section. While such a widespread range of τ cs appears superficially discouraging, it shall nevertheless be noted that many studies of the initial stages in heavy ion collisions have found that the α s = 0.3 is quite reasonable. Since the CME transport occurs mainly during the early stage, it is therefore a plausible choice to use this value for estimating τ cs in the present work. Furthermore, the large uncertainty helps remind us that the very initial motivation of studying CME is to utilize this phenomenon for providing unique experimental constraints on such a highly uncertain property of gauge fields. It may also be noted that the small but finite masses for the light flavor quarks (i.e. u, d, and s quarks) would also contribute to the diffusion rate. However recent analysis [39] has shown convincingly that such mass contribution, even for the strange quarks with ∼ 100MeV mass, is smaller than the above Chern-Simons diffusion rate by a few orders of magnitude and thus negligible. That means the strange quarks will also contribute to the anomalous transport just like the u and d flavors.
Last but not least, it is important to properly generate the initial conditions for the axial charge density based on the initial topological charge density of the gluon fields. At the very early stage of high energy heavy ion collisions, many flux tubes of strong chromo electric and magnetic fields in parallel or anti-parallel configurations are formed. In such a glasma picture [40][41][42][43], these flux tubes extend along the collision beam axis and are localized on the transverse plane. The chromo field strength inside the tubes is on the order of Q 2 s while their transverse size is on the order of 1/Q s , with Q s being the saturation scale. Depending on whether the chromo electric and magnetic fields are in parallel or anti-parallel configurations, each flux tube possesses randomly positive or negative topological charge density. They seed the generation of initial axial charge density n 5,i in the collisions [12,44,45]. We develop the following procedure to sample the initial axial charge density. For each collision event, we randomly sample a total of N coll glasma flux tubes on the transverse plane, where N coll is the binary collision number for each event generated from Monte Carlo Glauber simulations. For the i-th tube located at position x i , the chromo fields generate the following local axial charge upon integrating the Eq.( 2) up to the hydro initial proper time τ 0 :
where x, y are transverse spatial coordinates. In the above, the plus or minus sign is randomly sampled for each tube. For our computations we use an average value of Q 2 s = 1.5 GeV 2 in the above for collisions at RHIC energy [40][41][42][43]. The Gaussian factor in the bracket represents a smearing of the generated axial charge over a spread size of σ. In our calculations we vary σ in a reasonable range between (0.5 ∼ 1)fm as an estimator of systematic uncertainty due to the initial condition. The λ is a dimensionless strength parameter reflecting the fact that we only know these chromo fields up to the order of magnitude. Clearly λ controls the amount of initial topological windings and consequentially the axial charges.
As a last step, we superpose all the flux tubes in a given event together to obtain the overall axial charge initial profile: n 5 (τ 0 ) = N coll i=1 n 5,i . This can then be used as the initial condition for solving previously discussed evolution equation for the axial current in our framework. Finally one can integrate the n 5 over the fireball to obtain the total axial charge N 5 in each event.
Procedurally, one varies the value of control parameter λ and perform event-by-event simulations to generate the CME signal. This establishes a mapping between the initial topological fluctuations σ w and the final CME signal H SS-OS . Then by comparison with experimentally extracted value for H SS-OS , one can determines the corresponding σ w . It may be worth mentioning that the correlation measurements typically contain both CME signal and background contamination. Methods have been developed and demonstrated to be able to extract the signal part [46][47][48][49][50][51][52][53].
Using the above framework, we first show results for the AuAu collisions at 200GeV beam energy and (20 ∼ 50)% centrality in Fig. 2, where the dependence of CME-induced correlation signal H SS-OS on the initial topological fluctuations σ w is shown. (For this calculation the τ cs is calculated with α s = 0.3.) The blue circles (along with statistical uncertainty bars and systematic uncertainty boxes) are obtained from numerical simulations by varying the control parameter λ for initial chirality generation (-see Method section for details). A theoretically-expected quadratic dependence is clearly observed, with the dashed blue curve showing an excellent fitting: H SS-OS = 1.038 × 10 -9 σ w 2 . The star symbol and grey band show the current experimental constraints from analysis of the STAR collaboration data [24,25,54,55]. Comparison between calculations and measurements suggest an optimal value of σ w = 119 (as indicated by the red dot). This demonstrates how the initial topological fluctuations can be quantitative extracted from any definitive CME signal from future experimental measurements.
To demonstrate the impact of the stochastic dynamics, we have compared the current result with a reference calculation in which the stochastic dynamics is turned off and the initial axial charges do not dissipate away. We find that the correlation signal H SS-OS is reduced by a factor of 3.6 as compared with the reference calculation. Clearly the impact of stochastic topological fluctuations is significant and must be accounted for. To calibrated the uncertainty associated with the τ cs estimate, we've also done calculations with different input α values. We find that H SS-OS = 3.578 × 10 -9 σ w 2 for α s = 0.1 and H SS-OS = 7.4 × 10 -11 σ w 2 for α s = 0.5. The comparison confirms the sensitivity of the calculated signal to the choice of τ cs , which also implies that future measurement of CME signal could help constrain this key property of gauge field topological fluctuations.
As noted before, the finite quark masses have only negligible contributions to the axial charge dynamics and the strange quarks are expected to participate in CME transport just as the u and d flavors. A nontrivial consequence of this fact is a large charge correlation signal for strange mesons such as charged kaons, predicting a ratio H SS-OS K ± /H SS-OS π ± ≈ 1.2 which can be tested with CME measurements for identified hadrons. The star symbol and grey band show the current experimental constraints from analysis of the STAR collaboration data.
The red dot at σ w = 119 indicates the value of σ w that gives a correlation signal in consistence with experimental data.
Using the above result for AuAu collisions at (20 ∼ 50)% centrality as the benchmark, we next extend the calculations to all centrality classes as well as to other colliding systems at RHIC 200GeV beam energy, including the CuCu collisions as well as the isobaric RuRu and ZrZr collisions [56,57]. The ratio of σ w to initial total entropy S is found to increase from central to peripheral collisions and also increases from larger to smaller colliding systems. This behavior implies that when the fireball created in the collisions becomes smaller (due to either centrality or nucleus size), the σ w decreases less sensitively than the entropy S. Physically, the σ w is correlated with initial seeds of local gluon field flux tubes, the number of which roughly scales with the binary collision number N binary . The N binary is not directly measurable but is closely correlated with measurable charged hadron multiplicity N ch . To gain further insight, we plot in Fig. 3 the dependence of σ w on N ch with all centrality and all colliding systems from our computations. Remarkably the results demonstrate a universal scaling behavior, as is visible from the black straight line on the log-log scale. A fitting analysis reveals the following scaling relation:
It would be exciting to test this novel finding with future experimental analyses from all these colliding systems.
Finally, we utilize our finding above to analyze the isobar collisions and offer insights into the interpretation of the latest isobar data from STAR Collaboration [58], specifically focusing on the (20 ∼ 50)% centrality class. Based on the measured multiplicities and our finding in Eq. ( 8), one obtains σ Ru w = 36.32 and σ Zr w = 35.35. Using these inputs for AVFD simulations, one arrives at the following predictions for CME signals: H Ru = 1.41 × 10 -5 and H Zr = 1.17 × 10 -5 , with a signal ratio H Ru /H Zr = 1.2. For a reasonable estimate of the background correlations into the experimental observable known as the ∆γ-correlator, one could utilize an approximate scaling relation ∆γ bkg gv2 N ch [52,59,60] where the coefficient g is estimated to be 0.41 by a comparison with experimental data for ∆γ = 2H + ∆γ bkg consistently in both isobar systems. One can therefore further obtain the CME signal fraction f CM E in the overall ∆γcorrelator of both isobar systems, with f Ru CM E ≈ 0.040 and f Zr CM E ≈ 0.032, which appear to be consistent with an independent analysis in [60]. It however should be noted that a definitive answer for any potential CME signal in isobar data is still lacking, due to a range of uncertainties involved in both theoretical baseline estimates and experimental analyses.
In summary, we report a newly developed framework that has implemented a key ingredient, the stochastic dynamics of gauge field topological fluctuations, into the EBE-AVFD simulations for CME transport in these collisions.Such a framework provides the necessary tool to quantify initial topological fluctuations of QCD gluon fields from any definitive CME signal to be extracted from experimental measurements of the chiral magnetic effect in relativistic heavy ion collisions. By further applying this tool toward a variety of centrality class and colliding systems, a universal scaling relation between the initial gluon field topology fluctuations and the final charged hadron multiplicity has been identified in Eq.( 8), which can be tested by experimental data. With combined theoretical and experimental efforts in the near future, one can look forward to opportunities for extracting the relevant topological transition rate and understanding the mechanism underlying nontrivial phenomena such as cosmic baryon asymmetry as well as sponta-neous chiral symmetry breaking in QCD vacuum.