Investigation of Experimental Observables in Search of the Chiral Magnetic Effect in Heavy-ion Collisions in the STAR experiment

Figures(8) / Tables(2)

Get Citation
Subikash Choudhury, Xin Dong, Jim Drachenberg, James Dunlop, ShinIchi Esumi, Yicheng Feng, Evan Finch, Yu Hu, Jiangyong Jia, Jerome Lauret, Wei Li, Jinfeng Liao, Yufu Lin, Mike Lisa, Takafumi Niida, Robert Lanny Ray, Masha Sergeeva, Diyu Shen, Shuzhe Shi, Paul Sorensen, Aihong Tang, Prithwish Tribedy, Gene Van Buren, Sergei Voloshin, Fuqiang Wang, Gang Wang, Haojie Xu, Zhiwan Xu, Nanxi Yao and Jie Zhao. Investigation of Experimental Observables in Search of the Chiral Magnetic Effect in Heavy-ion Collisions in the STAR experiment[J]. Chinese Physics C.
Subikash Choudhury, Xin Dong, Jim Drachenberg, James Dunlop, ShinIchi Esumi, Yicheng Feng, Evan Finch, Yu Hu, Jiangyong Jia, Jerome Lauret, Wei Li, Jinfeng Liao, Yufu Lin, Mike Lisa, Takafumi Niida, Robert Lanny Ray, Masha Sergeeva, Diyu Shen, Shuzhe Shi, Paul Sorensen, Aihong Tang, Prithwish Tribedy, Gene Van Buren, Sergei Voloshin, Fuqiang Wang, Gang Wang, Haojie Xu, Zhiwan Xu, Nanxi Yao and Jie Zhao. Investigation of Experimental Observables in Search of the Chiral Magnetic Effect in Heavy-ion Collisions in the STAR experiment[J]. Chinese Physics C. shu
Received: 2020-02-02
Article Metric

Article Views(44)
PDF Downloads(5)
Cited by(0)
Policy on re-use
To reuse of subscription content published by CPC, the users need to request permission from CPC, unless the content was published under an Open Access license which automatically permits that type of reuse.
通讯作者: 陈斌,
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Email This Article


Investigation of Experimental Observables in Search of the Chiral Magnetic Effect in Heavy-ion Collisions in the STAR experiment

    Corresponding author: Yufu Lin,
    Corresponding author: Diyu Shen,
    Corresponding author: Nanxi Yao,
  • 1. Fudan University, Shanghai, 200433
  • 2. Lawrence Berkeley National Laboratory, Berkeley, California 94720
  • 3. Abilene Christian University, Abilene, Texas 79699
  • 4. Brookhaven National Laboratory, Upton, New York 11973
  • 5. University of Tsukuba, Tsukuba, Ibaraki 305-8571, Japan
  • 6. Purdue University, West Lafayette, Indiana 47907
  • 7. Southern Connecticut State University, New Haven, Connecticut 06515
  • 8. State University of New York, Stony Brook, New York 11794
  • 9. Rice University, Houston, Texas 77251
  • 10. Physics Department and Center for Exploration of Energy and Matter, Indiana University, 2401 N Milo B. Sampson Lane, Bloomington, IN 47408, USA
  • 11. College of Physics and Technology, Guangxi Normal University, Guilin, 541004, China
  • 12. Central China Normal University, Wuhan, Hubei 430079, People’s Republic of China
  • 13. Ohio State University, Columbus, Ohio 43210
  • 14. University of Texas, Austin, Texas 78712
  • 15. University of California, Los Angeles, California 90095
  • 16. Department of Physics, McGill University, 3600 University Street, Montreal, QC, H3A 2T8, Canada
  • 17. Wayne State University, Detroit, Michigan 48201
  • 18. Huzhou University, Huzhou, Zhejiang 313000

Abstract: The chiral magnetic effect (CME) is a novel transport phenomenon, arising from the interplay between quantum anomalies and strong magnetic fields in chiral systems. In high-energy nuclear collisions, the CME may survive the expansion of the quark-gluon plasma fireball and be detected in experiments. Over the past two decades, the experimental searches for the CME have aroused extensive interest at the Relativistic Heavy Ion Collider (RHIC) and the Large Hadron Collider (LHC). The main goal of this article is to investigate three pertinent experimental approaches: the $\gamma$ correlator, the R correlator and the signed balance functions. We will exploit both simple Monte Carlo simulations and a realistic event generator (EBE-AVFD) to verify the equivalence in the core components among these methods and to ascertain their sensitivities to the CME signal and the background contributions for the isobar collisions at RHIC.


    • A system is called chiral if it is not invariant under mirror reflection. The imbalance of right- and left-handed particles in a chiral system can be quantified by chiral chemical potential ($ \mu_5 $). In a system of charged fermions with finite $ \mu_5 $, an electric current could be generated in the presence of a strong magnetic field ($ \overrightarrow{B} $),

      $ \overrightarrow{J_e} \propto \mu_5\overrightarrow{B}, $


      which is theorized as the chiral magnetic effect (CME) [1, 2]. The CME physics encompasses a wide range of systems and has generated significant interdisciplinary interests (see some recent reviews in Refs. [3-6]). For example, it has been observed in condensed matter systems using Dirac and Weyl semimetals with emergent chiral quasiparticles (e.g., ZrTe5 [7], Na3Bi [8], TaAs [9] and TaP [10]). In this article, we present method studies in search of the CME in high-energy nuclear collisions.

      Ultra-relativistic heavy-ion collisions have been performed in the experimental facilities, such as the Relativistic Heavy Ion Collider (RHIC) and the Large Hadron Collider (LHC). These experiments aim to create a new phase of hot and dense nuclear matter with a temperature above several trillion Kelvin, consisting of deconfined quarks and gluons, called the quark-gluon plasma (QGP) [11-14]. Inside the QGP, the chiral symmetry is approximately restored, and light-flavor quarks $ u $ and $ d $ become nearly massless and hence chiral. The two preconditions (finite $ \mu_5 $ and $ \overrightarrow{B} $ field) for the CME could be realized in heavy-ion collisions as follows.

      Topological configurations of the color-SU(3) gluon fields in quantum chromodynamics (QCD) can generate chirality imbalance for light quarks through the chiral anomaly [15, 16], thus forming local chiral domains with finite $ \mu_5 $ in a QGP [1, 2, 17-20]. Note that the global chirality imbalance still vanishes when averaged over an infinite number of domains. In non-central heavy-ion collisions, extremely strong magnetic fields (B ~ 1014 T) can be formed [1, 17], by energetic protons. (By convention, the participating nucleons in the overlap region are called participants, and the rest, spectators.)

      Therefore, the CME could happen in heavy-ion collisions, and cause an electric current along the $ \overrightarrow{B} $ direction. Since $ \overrightarrow{B} $ is approximately perpendicular to the reaction plane ($ \Psi_{\rm{RP}} $) that contains the impact parameter and the beam momenta of a collision, the CME will manifest as an electric charge transport phenomenon across the reaction plane.

      In view of the CME-induced charge transport and other modes of collective motion of the QGP, the azimuthal distribution of particles is often Fourier-decomposed for given transverse momentum ($ p_T $) and pseudorapidity ($ \eta $) in an event:

      $\begin{aligned}[b] \frac{dN_{\alpha}}{d\phi^*} \approx& \frac{N_\alpha}{2\pi} [1 + 2v_{1,\alpha}\cos(\phi^*) + 2v_{2,\alpha}\cos(2\phi^*) \\&+ 2v_{3,\alpha}\cos(3\phi^*) + ... + 2a_{1,\alpha}\sin(\phi^*) + ...], \end{aligned} $


      where $ \phi $ is the azimuthal angle of a particle, and $ \phi^* = \phi - { \Psi_{{\rm{RP}}}} $. The subscript $ \alpha $ ($ + $ or $ - $) denotes the charge sign of a particle. Conventionally, the coefficients $ v_1 $, $ v_2 $ and $ v_3 $ are called "directed flow", "elliptic flow", and "triangular flow", respectively. They reflect the hydrodynamic response of the QGP medium to the initial collision geometry and to its fluctuations [21]. Here "RP" does not necessarily mean the reaction plane, but rather a flow symmetry plane obtained from the collective motion of final particles. For simplicity, we still use RP in the following discussions, and RP denotes a specific flow plane. The coefficient $ a_1 $ (with $ a_{1,-} = -a_{1,+} $ in a charge-symmetric system) characterizes the electric charge separation with respect to the flow plane, e.g., due to the CME. It is tempting to average the $ a_{1,\pm} $ coefficient in Eq. (2) over events to measure the CME signal. However, since $ \mu_5 $ flips sign on a domain-by-domain basis with equal probability (the global chirality should be balanced), the event-averaged $ a_{1,\pm} $ is zero by construction. Therefore, most experimental observables have been designed to measure the $ a_{1,\pm} $ fluctuations with respect to the flow plane.

      The experimental confirmation of the CME in heavy-ion collisions will open new windows for accessing fundamental aspects of strong interaction physics such as the QCD chiral symmetry restoration and the topological configurations of non-Abelian gauge fields. New information could also be extracted on the dynamical evolution of the extremely strong magnetic fields in these collisions, which provides crucial inputs for studying many other nontrivial effects from the magnetic field. A community-wide interest has been drawn to search for the CME at RHIC [22-32] and the LHC [33-37] over the past two decades. However, there is no definite conclusion so far for the existence of the CME in heavy-ion collisions. The main challenge in the CME search lies in the background contributions to the experimental observables. For example, mechanisms such as resonances with finite elliptic flow can also enhance charge fluctuations across flow planes [38-47]. To better gauge the background, the STAR experiment at RHIC has collected a large data sample of isobar collisions, and performed the corresponding CME-related analyses.

      The two isobaric systems, namely $ ^{96}_{44} {\rm{Ru}} + ^{96}_{44} {\rm{Ru}}$ and $ ^{96}_{40} {\rm{Zr}} + ^{96}_{40} {\rm{Zr}}$, have the same number of nucleons and hence similar amounts of elliptic flow, but different numbers of protons, which causes a difference in the magnetic field strength and in turn, a difference in the CME signal [49-51]. By keeping the background unchanged and varying the signal level, the two isobaric systems provide an ideal test ground for the CME study. The STAR Collaboration has implemented a blind-analysis recipe [52] to eliminate unintentional biases in data analyses, and all the analysis codes were frozen as part of the blinding procedure before unblinding of the isobar species. Since the isobar collisions will be examined with multiple observables, it is desirable to learn the connection and the difference between them, as well as their sensitivities to the CME signals. The objective of this paper is to do apples-to-apples comparisons between method kernels. For the sensitivity test of the final observables being used in the STAR blind analyses, we shall apply the STAR frozen codes to simulated events with various CME inputs. This work will provide an important reference point for the interpretation of the isobar-collision data that recently came out of the blind analyses [53].

      The paper is structured as follows. In Sec. II, we focus on three observables in search of the CME: the $ \gamma $ correlator [54], the $ R $ correlator [32, 55] and the signed balance functions [56, 57]. We shall uncover the relations between them. Section III describes the simple Monte Carlo calculations, as well as a more realistic event generator: Event-By-Event Anomalous-Viscous Fluid Dynamics (EBE-AVFD) [58-60]. These simulations are deployed in Sec. IV to compare the core components of the experimental observables. In Sec. V, EBE-AVFD is used to estimate the sensitivities of the final observables to the CME signal for the isobar collisions under study. Section VI gives the summary and outlook.

    • Several observables have been proposed to search for the CME in heavy-ion collisions, including the $ \gamma $ correlator [54], the $ R(\Delta S_m) $ correlator [32, 55], and the signed balance functions [56, 57]. It is not a surprise that these methods provide largely overlapping information, since they all make use of similar inputs of particle azimuthal correlations. We will review these approaches, and reveal the relations among them.

    • A.   $ \gamma $ Correlator

    • The three-point $ \gamma $ correlator (later specified as $ \gamma_{112} $) measures the fluctuations of charge separations or $ a_{1,\pm} $ coefficients with respect to a flow plane [54],

      $ \begin{aligned}[b] \gamma_{112} \equiv & \langle \cos(\phi_\alpha + \phi_\beta -2{\Psi_{{\rm{RP}}}}) \rangle \\ = & \langle\cos(\phi^*_{\alpha})\cos(\phi^*_{\beta}) - \sin(\phi^*_{\alpha})\sin(\phi^*_{\beta})\rangle \\ = & (\langle v_{1,\alpha}v_{1,\beta}\rangle + B_{\rm{IN}}) -(\langle a_{1,\alpha}a_{1,\beta}\rangle + B_{\rm{OUT}}), \end{aligned} $


      where the averaging is done over all combinations of particles $ \alpha $ and $ \beta $ in an event and over all events. The trigonometric expansion exhibits the difference between in-plane and out-of-plane projections of azimuthal correlations. The third term in Eq. (3), $ \langle a_{1,\alpha}a_{1,\beta}\rangle $, represents the fluctuations of $ a_{1,\pm} $ coefficients, the main target of the CME search. Other terms are presumably irrelevant to the CME: $ \langle v_{1,\alpha}v_{1,\beta}\rangle $ is related to directed flow, and is expected to be charge-independent and unrelated to the electromagnetic field in symmetric A+A collisions; $ B_{\rm{IN}} $ and $ B_{\rm{OUT}} $ represent other possible background correlations in and out of the flow plane, respectively. When we take the difference between opposite-sign and same-sign $ \gamma_{112} $ correlators,

      $ \Delta \gamma_{112} \equiv \gamma^{\rm{OS}}_{112} - \gamma^{\rm{SS}}_{112}, $


      the $ \langle v_{1,\alpha}v_{1,\beta}\rangle $ terms cancel out, as well as a large portion of ($ B_{\rm{IN}}-B_{\rm{OUT}} $). A residual flow-plane dependent background in ($ B_{\rm{IN}}-B_{\rm{OUT}} $) still exists at a level proportional to elliptic flow, which is the major background source in the $ \Delta \gamma_{112} $ measurements. In practice, the flow plane is approximated with the "event plane" ($ \Psi_{{\rm{EP}}} $) reconstructed with detected particles, and then the measurement is corrected for the finite event plane resolution [61]. The main advantages of $ \gamma_{112} $ lie in its direct connection to $ a_1 $ and a straightforward procedure that corrects for the finite event plane resolution.

      The flow-plane-related backgrounds in $ \Delta\gamma_{112} $ can be understood with the example of resonance decays. As resonances flow with the QGP medium, their decay daughters generate random (if parents possess no global spin alignment), event-by-event charge separation across the flow plane [44, 45]. The flowing resonance picture can be generalized to transverse momentum conservation (TMC) [42, 43] and local charge conservation (LCC) [44]. Ideally, the two-particle correlator,

      $ \begin{aligned}[b] \delta \equiv& \langle \cos(\phi_\alpha -\phi_\beta) \rangle \\ =& (\langle v_{1,\alpha}v_{1,\beta}\rangle + B_{\rm{IN}}) +(\langle a_{1,\alpha}a_{1,\beta}\rangle + B_{\rm{OUT}}), \end{aligned} $


      should also reflect $ \langle a_{1,\alpha} a_{1,\beta} \rangle $, but in reality it is dominated by short-range two-particle correlations. For example, the TMC effect contributes the following pertinent correlations in $ \delta $ and $ \gamma_{112} $ [42]:

      $ \delta^{\rm{TMC}} \approx -\frac{1}{N} \frac{\langle p_T \rangle^2_{\Omega}}{\langle p_T^2 \rangle_{\rm{F}}} \frac{1+({\bar v}_{2,{\Omega}})^2-2{\bar{\bar v}}_{2,{\rm{F}}}{\bar v}_{2,{\Omega}}} {1-({\bar{\bar v}}_{2,{\rm{F}}})^2}, $


      $ \begin{aligned}[b] \gamma^{\rm{TMC}}_{112} \approx& -\frac{1}{N} \frac{\langle p_T \rangle^2_{\Omega}}{\langle p_T^2 \rangle_{\rm{F}}} \frac{2{\bar v}_{2,{\Omega}}-{\bar{\bar v}}_{2,{\rm{F}}}-{\bar{\bar v}}_{2,{\rm{F}}}({\bar v}_{2,{\Omega}})^2} {1-({\bar{\bar v}}_{2,{\rm{F}}})^2} \\ \approx& \kappa^{\rm{TMC}}_{112} \cdot v_{2,{\Omega}} \cdot \delta^{\rm{TMC}}, \end{aligned} $


      where $ N $ is the total number of particles, $ \kappa^{\rm{TMC}}_{112} = $$ (2{\bar v}_{2,{\Omega}}-{\bar{\bar v}}_{2,{\rm{F}}})/v_{2,{\Omega}} $, and $ {\bar v}_{2} $ and $ {\bar{\bar v}}_{2} $ represent the $ p_T $- and $ p_T^2 $-weighted moments of $ v_2 $, respectively. The subscript "F" denotes an average over all produced particles in the full phase space, whereas the actual measurements only sample a fraction of the full space, denoted by "$ {\Omega} $". The background contribution due to the LCC effect has a similar characteristic structure as in Eqs. (6) and (7) [43, 44]. This motivates a normalization of $ \Delta \gamma $ by $ v_2 $ and $ \Delta \delta $:

      $ \kappa_{112} \equiv \frac{\Delta \gamma_{112}}{v_2 \cdot \Delta \delta}. $


      A CME signal will make $ \kappa_{112} $ larger than the background baseline, $ \kappa^{\rm{TMC/LCC}}_{112} $. While a reliable estimate of $ \kappa^{\rm{TMC/LCC}}_{112} $ is still elusive, the comparison of $ \Delta\gamma_{112} $ (and $ \kappa_{112} $) between isobar collisions may give a more definite conclusion on the CME signal.

      It is intuitive to introduce some derivative $ \gamma $ correlators and the corresponding $ \kappa $ observables [35], for example,

      $ \gamma_{123} \equiv \langle \langle \cos(\phi_\alpha + 2\phi_\beta -3{\Psi_{3}}) \rangle\rangle,\; \ \kappa_{123} \equiv \frac{\Delta \gamma_{123}}{v_3 \cdot \Delta \delta}, $


      $ \gamma_{132} \equiv \langle \langle \cos(\phi_\alpha - 3\phi_\beta + 2{\Psi_{2}}) \rangle\rangle,\; \ \kappa_{132} \equiv \frac{\Delta \gamma_{132}}{v_2 \cdot \Delta \delta}, $


      where $ \Psi_{2} $ and $ \Psi_{3} $ represent the $ 2^{\rm{nd}} $- and $ 3^{\rm rd} $-order flow planes, respectively. However, these observables may not serve as good background estimates for $ \gamma_{112} $ or $ \kappa_{112} $. Background-only AMPT calculations indeed show that $ \kappa_{123} $ and $ \kappa_{132} $ are not equal to $ \kappa_{112} $ in Au+Au collisions at 200 GeV [65]. Therefore, in the following sections we will not extend our method study to these derivative correlators.

    • B.   $ R $ correlator

    • The $ R(\Delta S_m) $ correlator [32, 55] takes the double ratio of four event-by-event distributions,

      $ R(\Delta S_m) \equiv \frac{N(\Delta S_{m,{\rm{real}}})}{N(\Delta S_{m,{\rm{shuffled}}})} \Bigg/ \frac{N(\Delta S^{\perp}_{m,{\rm{real}}})}{N(\Delta S^{\perp}_{m,{\rm{shuffled}}})}, \; \ m = 2,3,..., $


      where in a real event the charge separation perpendicular to the $ m^{\rm{th}} $-order flow plane ($ \Psi_m $) is expressed as

      $ \Delta S_{m,{\rm{real}}} = \langle\sin\left(\frac{m}{2}\phi^*_m\right)\rangle_{N_+} - \langle\sin\left(\frac{m}{2}\phi^*_m\right)\rangle_{N_-}. $


      Here $ \phi^*_m = \phi - {\Psi}_m $, and $ N_+ $($ N_- $) is the number of positively (negatively) charged particles in a specific event. Weighted averages could be used to take into account the azimuthal acceptance of the detector. $ \Delta S^{\perp}_m $ denotes the charge separation parallel to $ \Psi_m $, and is defined similar to $ \Delta S_m $, except that $ {\Psi}_m $ is replaced with $ ({\Psi}_m + \pi/m) $ to provide a baseline unrelated to the magnetic field. The $ \Delta S^{(\perp)}_{m,{\rm{shuffled}}} $ distributions are obtained via random charge reassignment (shuffling) of the reconstructed tracks in each real event, while respecting the multiplicities of positive and negative charges. Ideally the CME will cause a concave shape in the double ratio of the $ R(\Delta S_2) $ distribution, which presumably differs from the $ R(\Delta S_3) $ shape [55]. The latter is intended as a background estimate similar to the role of $ \gamma_{123} $.

      The $ R(\Delta S_2) $ distribution is fit with a Gaussian(inverse-Gaussian) function, if it bears a convex(concave) shape, and the Gaussian width ($ \sigma_{R2} $) is used to reflect the CME signal. Since the four core components, $ \Delta S^{(\perp)}_{2,{\rm{real}}({\rm{shuffled}})} $, all roughly follow Gaussian distributions, we can establish the relation between $ \sigma_{R2} $ and the RMS values of the core-component distributions:

      $ \begin{aligned}[b] \frac{S_{{\rm{concavity}}}}{\sigma_{R2}^2} =& \frac{1}{\langle (\Delta S_{2,{\rm{real}}})^2 \rangle} - \frac{1}{\langle (\Delta S_{2,{\rm{shuffled}}})^2\rangle} \\&- \frac{1}{\langle (\Delta S_{2,{\rm{real}}}^{\perp})^2 \rangle} + \frac{1}{\langle (\Delta S_{2,{\rm{shuffled}}}^{\perp})^2 \rangle}.\end{aligned} $


      The sign of concavity, $ S_{{\rm{concavity}}} $, is positive if the double ratio is convex, and negative for the concave shape. We shall first understand the meaning of each term on the right. For simplicity, we use unity weights, and the first term can be expanded as

      $ \begin{aligned}[b] (\Delta S_{2,{\rm{real}}})^2 \equiv& \left(\frac{\displaystyle\sum\nolimits_{i = 1}^{N_+}\sin(\phi^*_i)}{N_+}-\frac{\displaystyle\sum\nolimits_{i = 1}^{N_-}\sin(\phi^*_i)}{N_-}\right)^2 \\ = & \frac{\displaystyle\sum\nolimits_{i = 1}^{N_+}\sin^2(\phi^*_i)+\displaystyle\sum\nolimits_{i\neq j}^{N_+}\sin(\phi^*_i)\sin(\phi^*_j)}{N_+^2}\\&+\frac{\displaystyle\sum\nolimits_{i = 1}^{N_-}\sin^2(\phi^*_i)+\displaystyle\sum\nolimits_{i\neq j}^{N_-}\sin(\phi^*_i)\sin(\phi^*_j)}{N_-^2}\\&-\frac{2\displaystyle\sum\nolimits_{i = 1,j = 1}^{N_+,N_-}\sin(\phi^*_i)\sin(\phi^*_j)}{N_+N_-} \\ =& \frac{\langle \sin^2(\phi^*_i)\rangle_{N_+} +(N_+-1) \langle \sin(\phi^*_i)\sin(\phi^*_j)\rangle_{N_+N_+}}{N_+}\\&+\frac{\langle \sin^2(\phi^*_i)\rangle_{N_-} +(N_–1) \langle \sin(\phi^*_i)\sin(\phi^*_j)\rangle_{N_-N_-}}{N_-} \\ &-2\langle \sin(\phi^*_i)\sin(\phi^*_j) \rangle_{N_+N_-} \end{aligned} $


      Using the trigonometric identities, $ \sin^2(x) = $$ [1-\cos(2x)]/2 $ and $ 2\sin(x)\sin(y) = \cos(x-y)-\cos(x+y) $, we have the average of Eq. (14) over all events:

      $ \begin{aligned}[b] \langle(\Delta S_{2,{\rm{real}}})^2\rangle = & \frac{1-v_2^+}{2N_+} + \frac{N_+-1}{2N_+}(\delta^{++}-\gamma^{++}_{112})\\&+\frac{1-v_2^-}{2N_-} + \frac{N_–1}{2N_-}(\delta^{–}-\gamma^{–}_{112}) - (\delta^{+-}-\gamma^{+-}_{112}) \end{aligned}$


      $\quad\quad\quad \approx \frac{2(1-v_2-\delta^{\rm{SS}}+\gamma_{112}^{\rm{SS}})}{M} - \Delta\delta + \Delta\gamma_{112}. $


      The last line assumes $ v_2^+ \approx v_2^- $ and $ N_+ \approx N_- = M/2 $. Even before the approximations are taken, it is clear that $ \langle(\Delta S_{2,{\rm{real}}})^2\rangle $ can be expressed by $ N_{+(-)} $, $ v_2 $, $ \delta $ and $ \gamma_{112} $. Similarly,

      $ \begin{aligned}[b] \langle(\Delta S_{2,{\rm{real}}}^\perp)^2\rangle =& \frac{1+v_2^+}{2N_+} + \frac{N_+-1}{2N_+}(\delta^{++}+\gamma^{++}_{112})+\frac{1+v_2^-}{2N_-} \\&+ \frac{N_–1}{2N_-}(\delta^{–}+\gamma^{–}_{112}) - (\delta^{+-}+\gamma^{+-}_{112}) \end{aligned} $


      $\quad\quad\quad \approx \frac{2(1+v_2-\delta^{\rm{SS}}-\gamma_{112}^{\rm{SS}})}{M} - \Delta\delta - \Delta\gamma_{112}. $


      For the shuffled terms, $ v_2^+ $ and $ v_2^- $ will be roughly replaced with $ (v_2^+ + v_2^-)/2 $, all the $ \delta $ correlators will be replaced with $ (\delta^{\rm{OS}}+\delta^{\rm{SS}})/2 $, and all the $ \gamma $ correlators, with $ (\gamma^{\rm{OS}}+\gamma^{\rm{SS}})/2 $. Therefore, we have

      $ \langle(\Delta S_{2,{\rm{shuffled}}})^2\rangle \approx \frac{2(1-v_2)-\delta^{\rm{SS}}-\delta^{\rm{OS}}+\gamma_{112}^{\rm{SS}}+\gamma_{112}^{\rm{OS}}}{M} $


      $ \langle(\Delta S_{2,{\rm{shuffled}}}^\perp)^2\rangle \approx \frac{2(1+v_2)-\delta^{\rm{SS}}-\delta^{\rm{OS}}-\gamma_{112}^{\rm{SS}}-\gamma_{112}^{\rm{OS}}}{M}. $


      It is appealing to construct a double-subtraction observable based on the four core components of $ R(\Delta S_2) $, which shows an apparent link to $ \Delta\gamma_{112} $:

      $ \begin{aligned}[b] \Delta_{R2} \equiv& \langle (\Delta S_{2,{\rm{real}}})^2 \rangle - \langle (\Delta S_{2,{\rm{shuffled}}})^2 \rangle - \langle (\Delta S^{\perp}_{2,{\rm{real}}})^2 \rangle \\&+ \langle (\Delta S^{\perp}_{2,{\rm{shuffled}}})^2 \rangle \approx 2\left(1-\frac{1}{M}\right)\Delta \gamma_{112}. \end{aligned} $


      This relation will be tested in Sec. IV. To make the connection between the final observable, $ \sigma^2_{R2} $ and $ \Delta\gamma $, we put Eqs. (16), (18), (19) and (20) into Eq. (13), and reach

      $ \frac{S_{{\rm{concavity}}}}{\sigma_{R2}^2} \approx -\frac{M}{2}(M-1)\Delta\gamma_{112}. $


      Here we further assume that in each of the four $ \langle (\Delta S_2)^2\rangle $ terms, $ \dfrac{2}{M} $ is much larger than other contributions, among which $ \Delta\delta $ typically has the largest magnitude. This assumption may not always hold true, depending on the collision details such as beam energy and centrality. The relation in Eq. (22) implies that if $ \Delta\gamma $ is positively (negatively) finite, the $ R(\Delta S_2) $ distribution will exhibit a concave(convex) shape, and vice versa.

      The above derivation is with respect to the known flow plane $ \Psi_m $, which is not precisely known experimentally, and is assessed by the reconstructed event plane with a finite resolution. The empirical correction for the event plane resolution used in experiment for $ \Delta_{R2} $ and $ \sigma^2_{R2} $ [66] is nontrivial. An analytical resolution correction has been derived in Ref. [67]. Eqs. (21) and (22) can provide approximate approaches to examine the resolution correction factor.

      Multiplicity fluctuations widen the $ R(\Delta S_2) $ distributions. To reduce this effect, a scaling was proposed in Ref. [55], used in experimental data analyses, and included in the frozen code of the STAR isobar blind analysis. The $ R(\Delta S_2) $ distribution is converted into the $ R(\Delta S'_2) $ distribution by dividing the horizontal axis by the RMS of the $ N(\Delta S_{2, {\rm{shuffled}}}) $ distribution, i.e., $ \Delta S_2' = \Delta S/\sqrt{\langle(\Delta S_{2,{\rm{shuffled}}})^2\rangle} $ (see Sec.5 for details). Then, the width of $ R(\Delta S'_2) $ becomes

      $ \begin{aligned}[b] \frac{S_{{\rm{concavity}}}}{\sigma_{R2'}^2} =& \frac{S_{{\rm{concavity}}}}{\sigma_{R2}^2} \langle(\Delta S_{2,{\rm{shuffled}}})^2\rangle \\ \approx& -\frac{M}{2}(M-1)\Delta\gamma_{112} \times \frac{2}{M} \approx -M \Delta\gamma_{112} . \end{aligned} $

    • C.   Signed Balance Functions

    • Another method was recently proposed and invokes signed balance functions [56],

      $ \begin{aligned}[b] \Delta B_y \equiv& \Bigg[\frac{N_{y(+-)}-N_{y(++)}}{N_+} - \frac{N_{y(-+)}-N_{y(–)}}{N_-}\Bigg] \\&- \Bigg[\frac{N_{y(-+)}-N_{y(++)}}{N_+} - \frac{N_{y(+-)}-N_{y(–)}}{N_-}\Bigg] \\ =& \frac{N_+ + N_-}{N_+N_-}\Big[N_{y(+-)} - N_{y(-+)}\Big]. \end{aligned} $


      where $ N_{y(\alpha\beta)} $ is an event-by-event quantity, and denotes the number of pairs within which particle $ \alpha $ is ahead of particle $ \beta $ in the direction perpendicular to the reaction plane ($ p_y^\alpha > p_y^\beta $). Similarly, we can construct a $ \Delta B_x $ to count the number of pairs along the in-plane direction. Then the final observable is based on the widths of the $ \Delta B_y $ and the $ \Delta B_x $ distributions:

      $ r \equiv \sigma(\Delta B_y) / \sigma(\Delta B_x). $


      Intuitively, the CME will lead to $ r>1 $, since the CME-induced charge separation will cause more fluctuations of pair ordering across the reaction plane. The ratio $ r $ can be calculated in both the laboratory frame ($ r_{{\rm{lab}}} $) and the pair's rest frame ($ r_{{\rm{rest}}} $). It is argued that the rest frame is the most appropriate frame for $ r $ to study charge separations, and the further ratio,

      $ R_{{\rm{B}}} = r_{{\rm{rest}}} / r_{{\rm{lab}}}, $


      can help differentiate the background from the real CME signal [56]. Similar to the $ \gamma_{112} $ correlator, a lower event plane resolution will lead to a lower magnitude of the observed $ r $. Extra care is also needed to correct the $ r $ observables for the event plane resolution.

      In each event, we can rewrite the core component of $ \Delta B_y $ as follows,

      $ N_{y(\alpha\beta)}-N_{y(\beta\alpha)} = \sum\limits_{\alpha,\beta} {\rm{Sign}}[p_{T,\alpha}\sin(\phi^*_\alpha)-p_{T,\beta}\sin(\phi^*_\beta)]. $


      Compared with other methods, the signed balance functions could be more sensitive to the local CME domains that move with the expanding medium, because this method takes into account the transverse-momentum ordering instead of only the azimuthal angle. For example, a pair of particles going in the same direction can still be regarded as a case of charge separation by the signed balance functions, but not by other methods that only consider the azimuthal angle. This advantage, however, probably will not make a prominent difference, if the local domains merge into a global charge separation for the whole event after a full hydrodynamic evolution. Thus, to make connection to other observables, we take the first approximation by replacing $ p_T $ with mean $ p_T $, so that $ p_T $ can be dropped out for the time being, and only the azimuthal angle is exploited as done in other methods. Next, we want to unpack the Sign() function, and directly use $ [\sin(\phi^*_\alpha) - \sin(\phi^*_\beta)] $, which requires a normalization factor, $ C_y $. In view of the event average, we have

      $\begin{aligned}[b] \langle N_{y(\alpha\beta)}-N_{y(\beta\alpha)} \rangle \approx& C_y \Big\langle \sum\limits_{\alpha,\beta} [\sin(\phi^*_\alpha) - \sin(\phi^*_\beta)] \Big\rangle \\ =& C_y \Big\langle [N_\beta \sum\limits_{\alpha}\sin(\phi^*_\alpha) - N_\alpha \sum\limits_{\beta}\sin(\phi^*_\beta)] \Big\rangle \\ =& C_y N_\alpha N_\beta \langle \langle \sin( \phi^*)\rangle _{N_\alpha}-\langle \sin(\phi^*)\rangle_{N_\beta} \rangle. \end{aligned} $


      The constant can be determined by explicitly counting the pairs, with $ \dfrac{dN}{d\phi^*} $ from Eq. (2).

      $\begin{aligned}[b] \langle N_{y(\alpha\beta)}-N_{y(\beta\alpha)} \rangle = & 2\int_{-\pi/2}^{\pi/2} \Bigg[\int_{-\pi/2}^{\phi^*_\alpha} \frac{dN}{d\phi^*_\beta}d\phi^*_\beta+\int_{\pi-\phi^*_\alpha}^{3\pi/2} \frac{dN}{d\phi^*_\beta}d\phi^*_\beta\\&-\int^{\pi-\phi^*_\alpha}_{\phi^*_\alpha} \frac{dN}{d\phi^*_\beta}d\phi^*_\beta\Bigg]\frac{dN}{d\phi^*_\alpha}d\phi^*_\alpha \\ \approx& \frac{8}{\pi^2}\left(1+\frac{2}{3}v_2\right)N_\alpha N_\beta (a_{1,\alpha}-a_{1,\beta}). \end{aligned} $


      By comparing Eqs. (28) and (29), we learn $ C_y = 8(1+2v_2/3)/\pi^2 $. Therefore, if we ignore the $ p_T $ weight, $ \langle\Delta B_y\rangle $ becomes $ \dfrac{8M}{\pi^2}\left(1+\dfrac{2}{3}v_2\right) \langle\langle \sin(\phi^*) \rangle_{N_+} - $$ \langle \sin(\phi^*)\rangle_{N_-} \rangle $, which displays a function form akin to $ \langle S_{2,{\rm{real}}}\rangle $. In a similar way, we assume

      $ \langle N_{x(\alpha\beta)}-N_{x(\beta\alpha)} \rangle \approx C_x N_\alpha N_\beta \langle\langle \cos(\phi^*) \rangle _{N_\alpha}-\langle \cos(\phi^*)\rangle_{N_\beta} \rangle, $


      and the explicit counting gives

      $ \begin{aligned}[b] \langle N_{x(\alpha\beta)}-N_{x(\beta\alpha)} \rangle = & 2\int_{0}^{\pi} \Bigg[\int_{\phi^*_\alpha}^{2\pi-\phi^*_\alpha} \frac{dN}{d\phi^*_\beta}d\phi^*_\beta\\&-\int^{\phi^*_\alpha}_{-\phi^*_\alpha} \frac{dN}{d\phi^*_\beta}d\phi^*_\beta\Bigg]\frac{dN}{d\phi^*_\alpha}d\phi^*_\alpha \\ \approx& \frac{8}{\pi^2}\left(1-\frac{2}{3}v_2\right)N_\alpha N_\beta (v_{1,\alpha}-v_{1,\beta}). \end{aligned} $


      Thus, $ C_x = 8(1-2v_2/3)/\pi^2 $, and $ \langle\Delta B_x\rangle $ becomes $ \dfrac{8M}{\pi^2}\left(1-\dfrac{2}{3}v_2\right) \langle\langle \cos(\phi^*)\rangle_{N_+} -\langle \cos(\phi^*)\rangle_{N_-} \rangle $, resembling $ \langle \Delta S_{2,{\rm{real}}}^{\perp}\rangle $.

      In reality, both $ \langle\Delta B_{y(x)}\rangle $ and $ \langle \Delta S_{2,{\rm{real}}}^{(\perp)}\rangle $ are zero, but the derivation of $ C_y $ and $ C_x $ provides insight into the meaning of the signed balance functions. Our goal is to relate the RMS values of the $ \Delta B_y $ and $ \Delta B_x $ distributions to the other observables. Here we directly give the following relations, with the details of the analytical derivations explained in Appendix A.

      $ \sigma^2(\Delta B_y) \approx \frac{4M}{3}+ \frac{64M^2}{\pi^4}\left(1+ \frac{4}{3} v_2\right)(a_{1,+}-a_{1,-})^2 $


      $ \sigma^2(\Delta B_x) \approx \frac{4M}{3}+ \frac{64M^2}{\pi^4}\left(1-\frac{4}{3}v_2\right)(v_{1,+}-v_{1,-})^2 . $


      Then we define an observable that further connects the signed balance function to the $ \gamma $ correlator:

      $ \Delta_{\rm{SBF}} \equiv \sigma^2(\Delta B_y) - \sigma^2(\Delta B_x) \approx \frac{128M^2}{\pi^4}\left(\Delta\gamma_{112}-\frac{4}{3}v_2\Delta\delta\right). $


      Note that the signed balance functions take into account not only the azimuthal angles of particles, but also their momenta. If the ratio definition in Eq. (25) is transformed to $ \sigma^2(\Delta B_y) - \sigma^2({\Delta B_x}) $, then this method is roughly equivalent to $ (\Delta \gamma_{112}-\dfrac{4}{3}v_2\Delta\delta) $ with momentum weighting.

    • We now describe a toy model that undertakes simple Monte Carlo calculations, as well as the more realistic event generator, EBE-AVFD [58-60].

    • A.   Toy Model

    • We invoke a toy model to verify the mathematical relation between different methods with simplified Monte Carlo calculations. In this setup, particle spectra, collective flow and charge-separation signals are well constrained and conveniently adjusted, such that the response of each method to the signal and the background can be understood in a controlled manner. In the simulations, each event consists of 195 $ \pi^+ $ and 195 $ \pi^- $ mesons to match the total multiplicity at midrapidities within 2 units of rapidity in 30-40% central Au+Au collisions at $ \sqrt{s_{NN}} = $ 200 GeV [68]. In the background-free case, all pions are treated as primordial ones (none from resonance decays), and their azimuthal angle distribution is configured according to Eq. (2) with $ a_1 $, $ v_2 $ and $ v_3 $. The elliptic flow ($ v_2 $) is introduced by an NCQ-inspired function [69]

      $ v_2/{\cal{N}} = a/(1+e^{-[(m_T-m_0)/{\cal{N}} -b]/c}) - d, $


      where $ {\cal{N}} = 2 $ is the number of constituent quarks in a pion, and $ m_T $ and $ m_0 $ are its transverse mass and rest mass, respectively. The parameters ($ a $, $ b $, $ c $ and $ d $) are obtained by fitting Eq. (35) to the experimental data [46]. To add the resonance background, a fraction of primordial pions (33 $ \pi^+ $ and 33 $ \pi^- $) are replaced with 33 $ \pi^+ $-$ \pi^- $ pairs from $ \rho $-meson decays, which is implemented with PYTHIA6 [70]. The $ \rho $-meson $ v_2 $ is described by Eq. (35) with the corresponding $ \rho $ masses. $ v_3 $ is fixed at $ 1/5 $ of $ v_2 $ at any given $ p_T $ for both primordial pions and $ \rho $ resonances [71]. The primordial-pion spectra follow the Bose-Einstein distribution,

      $ \frac{dN_{\pi^\pm}}{dm_T^2} \propto (e^{m_T/T_{\rm{BE}}}-1)^{-1}, $


      where $ T_{\rm{BE}} = 212 $ MeV is set to match the experimentally observed $ \langle p_T \rangle $ of 400 MeV [68]. The $ \rho $-resonance spectra obey

      $ \frac{dN_\rho}{dm_T^2} \propto \frac{e^{-(m_T-m_\rho)/T}}{T(m_\rho+T)}, $


      where $ T = 317 $ MeV is used to match its $ \langle p_T \rangle $ of 830 MeV as observed in data [71]. Pseudorapidity (Rapidity) is uniformly distributed in the range of [-1, 1] for primordial pions ($ \rho $ resonances).

    • B.   EBE-AVFD

    • The EBE-AVFD model [58-60] is a comprehensive simulation framework that dynamically describes the CME in heavy-ion collisions. This state-of-the-art tool has been developed over the past few years as an important part of the efforts within the Beam Energy Scan Theory (BEST) Collaboration, aiming to address the needs of the ongoing experimental program at RHIC collision energies. Critical to the success of the CME search is a quantitative and realistic characterization of the CME signals as well as the relevant backgrounds. Accordingly, EBE-AVFD implements the dynamical CME transport for quark currents on top of the relativistically expanding viscous QGP fluid and properly models major sources of background correlations such as LCC and resonance decays.

      More specifically, the EBE-AVFD framework starts with event-wise fluctuating initial conditions, and solves the evolution of chiral quark currents as linear perturbations in addition to the viscous bulk flow background provided by data-validated hydrodynamic simulation packages. The LCC effect is incorporated in the freeze-out process, followed by the hadron cascade simulations. This is illustrated in the following Fig. 1.

      Figure 1.  (color online) Flow chart of the EBE-AVFD framework

      The fluctuating initial conditions for entropy density profiles are generated by the Monte-Carlo Glauber model, with switching time $ \tau_0 = 0.6\; {\rm{fm}}/c $ and mixing parameter $ \alpha_{\rm{glb}} = 0.118 $. The latter factorizes the contributions of the participant density ($ n_{\rm{part}} $) and binary collision density ($ n_{\rm{coll}} $) in the local entropy density, $ s = \alpha_{\rm{glb}}\, n_{\rm{coll}} + $$ (1-\alpha_{\rm{glb}})n_{\rm{part}}/2 $. The initial axial charge density ($ n_5 $) is approximated in such a way that it is proportional to the corresponding local entropy density with a constant ratio. This ratio parameter can be varied to sensitively control the strength of the CME transport. For example, one can set $ n_5/s $ to $ 0 $, $ 0.1 $ and $ 0.2 $ in the simulations to represent scenarios of zero, modest and strong CME signals respectively. The initial electromagnetic field, $ \overrightarrow{B}_{ini}(\overrightarrow{x}) $, is computed according to the event-wise proton configuration in the Monte-Carlo Glauber initial conditions. Then, the space-time evolution is modeled as the product of the time-independent initial condition and a space-independent function characterizing its decay in time, $ \overrightarrow{B}(\overrightarrow{x},\tau) = \overrightarrow{B}_{ini}(\overrightarrow{x}) \times F(\tau) $. With the medium effects taken into account, the time-dependence is chosen to be $ F(\tau) = 1/(1+\tau^2/\tau_B^2) $ with $ \tau_B = 0.6 \;{\rm{fm}}/ c $ being the lifetime of the B-field.

      The hydrodynamic evolution is solved through two components. The bulk-matter collective flow is described by the VISH2+1 simulation package [62], with the lattice equation of state s95p-v1.2, shear-viscosity $ \eta/s = 0.08 $, and freeze-out temperature $ T_{\rm{fo}} = 160\; $MeV. Such hydrodynamic simulations of bulk flow have been extensively tested and validated with relevant experimental data. The dynamical CME transport is described by anomalous hydrodynamic equations for the quark chiral currents on top of the bulk flow background, where the magnetic-field-induced CME currents lead to a charge separation in the fireball. Additionally the conventional transport processes like diffusion and relaxation for the quark currents are consistently included, with the diffusion constant chosen to be $ \sigma = 0.1\,T $ and relaxation time $ \tau_r = 0.5/T $. More discussions of the hydrodynamics equations and relevant details can be found in Refs. [58-60].

      After the hydrodynamic stage, hadrons are locally produced in all fluid cells on the freeze-out hypersurface, using the Cooper-Frye freeze-out formula

      $ E \frac{dN}{d^3p} (x^\mu, p^\mu) = \frac{g}{(2\pi)^3} \int_{\Sigma_{\rm{fo}}} p^\mu d^3\sigma_\mu f(x,p) \,\, . $


      Here, the local distribution function automatically includes the charge separation effect due to the CME as well as non-equilibrium corrections. In the freeze-out process, the LCC effect is implemented by extending an earlier method from Ref. [63]. The approach in Ref. [63] chooses to produce all charged hadron-antihadron pairs at the same fluid cell, while their momenta are sampled independently in the local rest frame of the fluid cell. This treatment implicitly assumes the charge-correlation length to be smaller than the cell size, and hence provides an upper limit for the correlations between opposite-sign pairs. In the EBE-AVFD package, the aforementioned procedure is generalized and improved to mimic more realistically the impact of a finite charge-correlation length: a new parameter $ P_{\rm{LCC}} $ is introduced to characterize the fraction of charged hadrons that are sampled in positive-negative pairs in the same way as in Ref. [63], while the rest of the hadrons are sampled independently. Varying the parameter $ P_{\rm{LCC}} $ between 0 and 1 would tune the LCC contributions from none to its maximum. Finally, all the hadrons produced from the freeze-out hypersurface are further subject to hadron cascades through the UrQMD simulations [64], which account for various hadron resonance decay processes and automatically include their contributions to the charge-dependent correlations. The tuning of the EBE-AVFD calculations to the experimental measurements of $ \Delta\delta $ and $ \Delta\gamma_{112} $ in Au+Au collisions suggests that an optimal value of $ P_{\rm{LCC}} $ is around $ 1/3 $, and that roughly half of the background correlations come from LCC and the other half from resonance decays.

    • We will exploit the toy model and the EBE-AVFD model to simulate the core components of the experimental observables introduced in Sec. 2: $ \Delta \gamma_{112} $ for the $ \gamma $ correlator, $ \Delta_{R2} $ for the $ R $ correlator, and $ \Delta_{\rm{SBF}} $ for the signed balance functions. Our objective is to examine the responses of the core components to the CME signal and the background, and to verify the relations between these methods (Eqs. (21) and (34)). For a fair comparison with other observables, the momentum weighting is not applied in the $ \Delta_{\rm{SBF}} $ results. For simplicity, the true reaction plane is used in all the simulations in this section. The particles of interest are selected with $ |\eta|<1 $ and $ 0.2 < p_T < 2 $ GeV/c.

    • A.   Toy-model Results

    • We put the three core components on the same footing, by plotting $ 2\Delta\gamma_{112} $, $ \Delta_{R2} $ and $ \Delta'_{\rm{SBF}}\equiv\left(\dfrac{\pi^4}{64M^2}\Delta_{\rm{SBF}}+ $$ \dfrac{8}{3}v_2\Delta\delta\right) $ as function of the input $ a_1^2 $ in Fig. 2 for two scenarios: with and without resonance decays (LCC is not implemented in the toy model). In the background-free case without decays, the three observables render very similar results (open markers), all falling on the linear function of $ 4a_1^2 $. Therefore, all the three methods are sensitive to the same amount of the CME signal. When resonance decays are turned on with finite elliptic flow, sizeable background effects appear besides the pure-signal contributions for all the three approaches (solid markers), more prominent at smaller input $ a_1^2 $ values. Note that each final observable can be roughly regarded as a weighted average of the correlations due to the CME, the resonance background and the cross terms. At unrealistically large $ a_1^2 $ values, the CME contribution in an observable could be diluted by the resonance contribution, since the latter becomes smaller than the former. The three core components exhibit similar responses to the backgrounds due to flowing resonances in this toy model. In this scenario, there are some subtle differences between the results from these three approaches, probably because of some higher-order effects omitted in the derivation of Eqs. (21) and (34). Although the background contributions depend on spectra and particularly elliptic flow of resonances [43, 44, 47], to the first order, we expect the three observables to have similar responses to resonance decays for a wide range of spectra or elliptic flow. A recent study [48] has also found that the $ R $ correlator and the $ \Delta\gamma_{112} $ correlator have similar sensitivities to the CME signal and the background.

      Figure 2.  (color online) The toy-model simulations of $ 2\Delta\gamma_{112} $, $ \Delta_{R2} $ and $ \Delta'_{\rm{SBF}}\equiv(\frac{\pi^4}{64M^2}\Delta_{\rm{SBF}}+\frac{8}{3}v_2\Delta\delta) $ as function of the input $ a_1^2 $. The open markers represent the pure-signal scenario without resonances, and the solid markers denote the scenario with resonance decays. In comparison, the linear function of $ 4a^2_1 $ is also added

    • B.   EBE-AVFD Results

    • The EBE-AVFD model implements the CME and the backgrounds in a more realistic way. In the following simulations, we generate the EBE-AVFD events of 30-40% Au+Au collisions at $ \sqrt{s_{{\rm{NN}}}} = 200 $ GeV, with $ n_{5}/s $ = 0, 0.1 and 0.2. The background effects almost remain the same, whereas the CME signal is varied according to the input $ n_{5}/s $. Figure 3(a) presents the corresponding calculations of $ 2\Delta\gamma_{112} $, $ \Delta_{R2} $ and $ \Delta'_{\rm{SBF}} $ as function of $ n_{5}/s $. The three methods yield very similar results at each input $ n_{5}/s $ value, supporting the relations expressed in Eqs. (21) and (34).

      Figure 3.  (color online) (a) The EBE-AVFD simulations of $ 2\Delta\gamma_{112} $, $ \Delta_{R2} $ and $ \Delta'_{\rm{SBF}}\equiv(\frac{\pi^4}{64M^2}\Delta_{\rm{SBF}}+\frac{8}{3}v_2\Delta\delta) $ as function of $ n_5/s $ in 30-40% Au+Au collisions at 200 GeV. (b) The same results with the subtraction of the pure-background case vs $ (a_{1,+}^2 + a_{1,-}^2 - 2a_{1,+}a_{1,-}) $. In comparison, a linear function of $ y = x $ is drawn to verify the relation in Eq. (39).

      With the known reaction plane angle in each EBE-AVFD event, we can readily calculate $ a_{1,\pm} $, and check if this CME contribution explains the difference between the cases with different $ n_{5}/s $ values. $ a_{1,\pm} $ is consistent with zero for $ n_{5}/s = 0 $, and $ a_{1,+} $ and $ a_{1,-} $ are finite with opposite signs for finite $ n_{5}/s $ values. Note that $ a_{1,+} $ and $ a_{1,-} $ do not necessarily have the same magnitude, because the collision system always bears extra positive charges. Based on the expansion of the $ \gamma_{112} $ correlator in Eq. (3) and the equivalence between the three observables, we expect the following equation for any of these observables, $ O(n_{5}/s) $,

      $ O(n_{5}/s) - O(0) = a_{1,+}^2 + a_{1,-}^2 - 2a_{1,+}a_{1,-}. $


      Figure 3(b) shows that the results for each observable, after the subtraction of the pure-background case, fall on the straight line representing the relation in Eq. (39). Thus, the EBE-AVFD calculations reveal the linear superposition of the CME signal and the background contribution in the experimental observables. This property is implicitly assumed by most of the analysis techniques that attempt to separate the CME signal and the backgrounds, and it is now corroborated by the EBE-AVFD model.

      The core-component comparison using both the toy model and the EBE-AVFD model support the idea that to the first order, the three observables are equivalent to each other, with their very similar responses to the CME signal as well as the backgrounds.

    • The validity of Eq. (39) for the three experimental observables reassures the feasibility of disentangling the CME signal from the backgrounds with the isobar collisions. To eliminate unintentional biases in the analyses of the isobar-collision data, the STAR Collaboration has followed a blind-analysis procedure [52]. One important step in this procedure requires the analysis codes from all analyzers to be frozen before the mass production of the isobar-collision data. In this section, we shall apply these frozen codes to the EBE-AVFD events of isobar collisions, to investigate the realistic sensitivity of each method to the CME signal. The particles of interest come from the acceptance of $ |\eta|<1 $, with $ 0.2 < p_T < 2 $ GeV/c for the $ \gamma $ correlator and the signed balance functions, or with $ 0.35 < p_T < 2 $ GeV/c for the $ R $ correlator, as implemented in the frozen codes. Owing to different kinematic cuts used in different methods, the sensitivity study in this manner cannot guarantee apples-to-apples comparisons as was assured in the core-component comparisons in Sec 4. Nevertheless, we will obtain a reliable benchmark for interpreting experimental data.

      For each of the two isobaric collision systems, Ru+Ru and Zr+Zr at $ \sqrt{s_{{\rm{NN}}}} = 200 $ GeV, four cases of the EBE-AVFD events have been generated, with $ n_5/s = 0,\, 0.05,\, 0.1, \, {\rm{and}} \, 0.2 $, respectively. The centrality selection for all the cases focuses on 30-40% central collisions, where the potential CME signal is relatively easy to detect owing to good event plane resolutions. 200 million events are produced for each case of $ n_5/s = 0 $ and $ n_5/s = 0.2 $, and 400 million events for each of the other cases. To mimic the detection performance of the STAR Time Projection Chamber, the simulated particles in the EBE-AVFD events are randomly rejected according to a transverse-momentum dependent tracking efficiency. The event plane resolution for these EBE-AVFD events is matched to that for the real isobar data (Isobar Mixed Analysis in Ref. [53]).

      Table I lists $ a_{1,\pm} $ calculated by EBE-AVFD at different $ n_5/s $ values. Basically $ a_{1,\pm} $ displays a linear function of $ n_5/s $. $ a_{1,+} $ is about 4% larger than $ a_{1,-} $ for all the cases, reflecting the charge asymmetry in the collision systems. Meanwhile, $ a_{1,\pm} $ is roughly 7% larger in Ru+Ru than in Zr+Zr collisions, which meets the expectation from the difference in the corresponding magnetic fields [50]. For the experimental observables under study, the background levels are almost the same for the two isobaric systems, and the ratio of the Ru+Ru measurement to the Zr+Zr measurement is expected to be larger than unity, in the presence of a positive CME signal. Hence, the sensitivity of each method can be defined with the statistical significance of the deviation of this ratio from unity.

      $n_{5}/s$ $a_{1,+}$ (%) $a_{1,-}$ (%)
      Ru+Ru Zr+Zr Ru+Ru Zr+Zr
      0 0 0 0 0
      0.05 0.37 0.35 0.35 0.33
      0.10 0.74 0.69 0.71 0.66
      0.20 1.48 1.38 1.42 1.32

      Table 1.  The $a_{1,\pm}$ values calculated for the EBE-AVFD events of 30-40% isobar collisions at $\sqrt{s_{\rm{NN}}} = 200$ GeV

      Figure 4 presents the EBE-AVFD calculations of $ \gamma_{112}^{{\rm{OS}}({\rm{SS}})} $ (a) and $ \Delta \gamma_{112} $ (b) as functions of $ n_{5}/s $ for 30-40% isobar collisions at $ \sqrt{s_{\rm{NN}}} = 200 $ GeV. The ratios of $ \Delta \gamma_{112} $ between Ru+Ru and Zr+Zr is delineated in panels (c). The $ 2^{\rm{nd}} $-order event plane is reconstructed from the same kinematic region as the particles of interest, and the observed $ \gamma $ correlators and $ v_2 $ (to be shown later) have been corrected with the corresponding event plane resolution. At each $ n_5/s $ value, $ \gamma_{112}^{\rm{OS}} $ remains positive and $ \gamma_{112}^{\rm{SS}} $ stays negative, both with larger magnitudes at higher $ n_5/s $. Although the CME expects $ \gamma_{112}^{\rm{OS}} $ and $ \gamma_{112}^{\rm{SS}} $ to be symmetric around zero, there exist some charge-independent backgrounds such as momentum conservation and elliptic flow that shift both $ \gamma_{112}^{\rm{OS}} $ and $ \gamma_{112}^{\rm{SS}} $ up or down [24]. Therefore, we shall focus on $ \Delta\gamma_{112} $, which shows a finite background contribution at $ n_5/s = 0 $ and increases with the CME signal. The difference between Ru+Ru and Zr+Zr is better viewed with the ratio of $ \Delta\gamma_{112}^{{\rm{Ru+Ru}}}/\Delta\gamma_{112}^{{\rm{Zr+Zr}}} $. This ratio is consistent with unity at $ n_5/s = 0 $, and increases quadratically with $ n_5/s $ as demonstrated by the $ 2^{\rm{nd}} $-order-polynomial fit function that passes (0, 1) (dashed line). The quadratically-increasing trend is expected, because this ratio is a linear function of the CME signal fraction in $ \Delta\gamma_{112} $ in a two-component perturbative framework [50], and the latter is proportional to $ (n_5/s)^2 $ or $ a_1^2 $ as shown in Fig. 3.

      Figure 4.  (color online) EBE-AVFD calculations of $ \gamma_{112}^{ {\rm{OS}}({\rm{SS}})} $ (a) and $ \Delta \gamma_{112} $ (b) as functions of $ n_{5}/s $ for 30-40% isobar collisions at $ \sqrt{s_{\rm{NN}}} = 200 $ GeV, together with the ratio of $ \Delta \gamma_{112} $ (c) between Ru+Ru and Zr+Zr. In panel (c), the $ 2^{\rm{nd}} $-order-polynomial fit function illustrates the rising trend starting from (0, 1)

      Potentially the $ \delta $ correlator could also contain the CME signal as Eq. (5) suggests, which should be manifested in the $ \Delta\delta $ ratio between Ru+Ru and Zr+Zr [60]. Figure 5 depicts the EBE-AVFD results of $ \delta^{ {\rm{OS}}({\rm{SS}})} $ (a) and $ \Delta \delta $ (b) as functions of $ n_{5}/s $ for 30-40% isobar collisions at $ \sqrt{s_{\rm{NN}}} = 200 $ GeV. The ratio of $ \Delta \delta $ between Ru+Ru and Zr+Zr is shown in panel (c). For all the cases, $ \delta^{\rm{OS}} $ is above $ \delta^{\rm{SS}} $, leading to a positive $ \Delta\delta $, because of larger background contributions to $ \delta^{\rm{OS}} $ than $ \delta^{\rm{SS}} $. It is argued that although overshadowed by the background contributions, the CME signal in $ \Delta\delta $ could be even larger than that in $ \Delta\gamma_{112} $ measured with respect to the participant plane [60]. This is because the magnetic field difference between the two isobar systems is maximal with respect to the reaction plane, and is reduced when measured otherwise. The EBE-AVFD simulations indeed support this idea: at each $ n_5/s $ value in Fig. 5(c), the deviation of the $ \Delta\delta $ ratio from unity is at the similar level as that of the $ \Delta\gamma_{112} $ ratio. Since $ \Delta\delta $ is typically larger than $ \Delta\gamma $ by an order of magnitude, this similar deviation in their respective ratios indicates a significantly larger CME effect in $ \Delta\delta $ than in $ \Delta\gamma $. A $ 2^{\rm{nd}} $-order polynomial fit function is added to guide the eye. In view of the smaller relative statistical uncertainties of $ \Delta\delta $ than those of $ \Delta\gamma_{112} $, the former may yield even better significance levels of the CME signal than the latter in the isobar-collision data, provided equal background contributions to $ \Delta\delta $ in the two systems. The caveat is that the two-particle correlation background to $ \Delta\delta $ is significantly larger than that to $ \Delta\gamma $, so any difference in the background between the two isobaric systems would have a stronger impact on $ \Delta\delta $. However likely or unlikely this scenario is to occur, we keep these results for completeness and leave them to the test with the real data.

      Figure 5.  (color online) EBE-AVFD calculations of $ \delta^{ {\rm{OS}}({\rm{SS}})} $ (a) and $ \Delta \delta $ (b) as functions of $ n_{5}/s $ for 30-40% isobar collisions at $ \sqrt{s_{\rm{NN}}} = 200 $ GeV, together with the ratios of $ \Delta \delta $ (c) between Ru+Ru and Zr+Zr. In panels (c), the $ 2^{\rm{nd}} $-order-polynomial fit function is added to demonstrate the rising trend starting from (0, 1)

      Between Ru+Ru and Zr+Zr, the difference in the background contributions to $ \Delta\gamma_{112} $ may be small, but could still be finite owing to the possibly different $ v_2 $ values. The normalization of $ \Delta\gamma_{112} $ by $ v_2 $ would be a more robust variable for the CME search. The EBE-AVFD simulations of $ v_2 $ (a) and $ \Delta\gamma_{112}/v_2 $ (c) are presented in Fig. 6 as functions of $ n_{5}/s $ for 30-40% isobar collisions at $ \sqrt{s_{\rm{NN}}} = 200 $ GeV, with the corresponding $ v_2 $ ratio (b) and $ \Delta\gamma_{112}/v_2 $ ratio (d) between Ru+Ru and Zr+Zr. The $ v_2 $ values are very close to each other between Ru+Ru and Zr+Zr collisions, with the relative difference at the level of 0.1%. Because of this, the Ru+Ru/Zr+Zr ratio of $ \Delta\gamma_{112}/v_2 $ in Fig. 6(d) is practically identical to that of $ \Delta\gamma $ in Fig. 4(c). Besides the possible $ v_2 $ difference, the two-particle correlation strengths (which is part of the CME background) could also differ from Ru+Ru to Zr+Zr collisions. As we mentioned above, although $ \Delta\delta $ is also sensitive to the CME, it is overwhelmed by background correlations and thus may be used as an approximate gauge for the two-particle correlation strength. The additional normalization of $ \Delta\gamma_{112}/v_2 $ by $ \Delta\delta $ in $ \kappa_{112} $ could, therefore, further suppress this difference, and enhance the sensitivity to the CME signal. There are two potential scenarios where $ \kappa_{112} $ has the advantage over $ \Delta\gamma_{112} $. First, in reality the $ \Delta\delta $ ratio may disfavor the CME interpretation, and in that case, the $ \kappa_{112} $ ratio is less prone to a faked CME signal than the $ \gamma_{112} $ ratio. Second, $ \Delta\delta $ may also contain the CME signal as shown in Fig. 5(c), while the relative statistical uncertainty of $ \Delta\delta $ is smaller than that of $ \Delta\gamma_{112} $. In that case, $ \kappa_{112}^{\rm{Ru+Ru}}/\kappa_{112}^{\rm{Zr+Zr}} $ could yield a larger ratio than $ \Delta\gamma_{112}^{\rm{Ru+Ru}}/\Delta\gamma_{112}^{\rm{Zr+Zr}} $, with a similar relative statistical uncertainty, which means a better significance of the CME signal. Indeed, panel (f) also shows a quadratically increasing trend, with larger magnitudes than panel (d) at finite $ n_5/s $ values, and the calculated significance values for $ \kappa_{112}^{\rm{Ru+Ru}}/\kappa_{112}^{\rm{Zr+Zr}} $ roughly double those for $ \Delta\gamma_{112}^{\rm{Ru+Ru}}/\Delta\gamma_{112}^{\rm{Zr+Zr}} $, as documented later in Table II. In general, the EBE-AVFD simulations show good responses to the signal for $ \Delta\gamma_{112}^{\rm{Ru+Ru}} / \Delta\gamma_{112}^{\rm{Zr+Zr}} $, $ \Delta\delta^{\rm{Ru+Ru}} / \Delta\delta^{\rm{Zr+Zr}} $ and $ \kappa_{112}^{\rm{Ru+Ru}} / \kappa_{112}^{\rm{Zr+Zr}} $, and these promising features await verification by the isobar-collision data.

      Figure 6.  (color online) EBE-AVFD calculations of $ v_2 $ (a), $ \Delta\gamma_{112}/v_2 $ (c) and $ \kappa_{112} $ (e) as functions of $ n_{5}/s $ for 30-40% isobar collisions at $ \sqrt{s_{\rm{NN}}} = 200 $ GeV, together with the ratios of $ v_2 $ (b), $ \Delta\gamma_{112}/v_2 $ (d) and $ \kappa_{112} $ (f) between Ru+Ru and Zr+Zr. In panels (f), the $ 2^{\rm{nd}} $-order-polynomial fit function illustrates the rising trend starting from (0, 1)

      $ n_{5}/s $ $ N_{\rm{event}} $ $ \Delta \gamma_{112} $ $ \Delta\delta $ $ \kappa_{112} $ $ r_{\rm{lab}} $ $ \sigma_{R2}^{-1} $ $ \Delta\gamma_{112}\{{\rm{RP}}\} $ $ \sigma_{R2}^{-1}\{{\rm{RP}}\} $
      0 $ 2\times 10^8 $ −1.50 −2.89 −1.21 −0.77 1.33 0.67 0.56
      0.05 $ 4\times 10^8 $ 0.62 −6.16 1.37 0.47 0.29 2.84 3.33
      0.10 $ 4\times 10^8 $ 1.91 −16.81 3.43 3.11 0.62 11.78 10.85
      0.20 $ 2\times10^8 $ 7.73 −42.96 14.07 5.96 1.84 37.48 27.90

      Table 2.  The statistical significance of $ (O^{\rm{Ru+Ru}}/O^{\rm{Zr+Zr}} -1) $ for different experimental observables. Contrary to other observables, the $ \Delta\delta $ ratio expects negative significance values due to the CME signal. $ N_{\rm{event}} $ means the number of events used for each isobaric system in the simulation. For completeness, the last two columns list the significance values for $ \Delta\gamma_{112} $ and $ \sigma^{-1}_{R2} $ with respect to the reaction plane (not following the frozen code)

      A similar frozen-code analysis is performed for the $ R(\Delta S_2) $ correlator, and the results are presented in Figure 7. In order to minimize the influence of the particle number fluctuations, the $ R(\Delta S_2) $ distribution is converted into the $ R(\Delta S_2') $ distribution by dividing the horizontal axis by the RMS of the $ N(\Delta S_{2,{\rm{shuffled}}}) $ distribution, i.e., $ \Delta S_2' = \Delta S_2/\sqrt{\langle(\Delta S_{2,{\rm{shuffled}}})^2\rangle} $. Then $ \Delta S_2' $ is further modified to correct for the event plane resolution, i.e., $ \Delta S'' = \Delta S'/\delta_{{\rm{Res}}} $, where $ \delta_{{\rm{Res}}} $ is the correction factor whose details can be found in Ref. [55]. Panels (a) and (b) show the $ R(\Delta S_2'') $ distributions from EBE-AVFD events of 30-40% Ru+Ru and Zr+Zr collisions, respectively, at $ \sqrt{s_{\rm{NN}}} = 200 $ GeV with different $ n_{5}/s $ inputs. As $ n_{5}/s $ increases, the $ R(\Delta S_2'') $ distribution becomes more concave, qualitatively representing more CME contributions. To quantify the distribution shape, the Gaussian width ($ \sigma_{R2} $) is obtained by fitting each $ R(\Delta S_2'') $ distribution with an inverse Gaussian function, and the resultant $ \sigma^{-1}_{R2} $ values are listed in panel (c), increasing with $ n_{5}/s $. The $ \sigma^{-1}_{R2} $ ratios between Ru+Ru and Zr+Zr are shown in panel (d). According to Eq. (22), $ \sigma^{-1}_{R2} $ is proportional to $ \sqrt{\Delta\gamma_{112}} \propto [(\sigma^{-2}_{R2})^{\rm{CME}}+(\sigma^{-2}_{R2})^{{\rm{BG}}}]^{1/2} $, and hence the $ \sigma^{-1}_{R2} $ ratio is expected to follow a quadratic trend vs $ n_5/s $, assuming the CME signal is small compared with the background. We fit the $ \sigma^{-1}_{R2} $ ratios with a $ 2^{\rm{nd}} $-order polynomial function starting from (0, 1), though the statistical uncertainties are large. The significance values of these ratios are stored in Table II for later discussions.

      Figure 7.  (color online) Distributions of $ R(\Delta S_2'') $ from EBE-AVFD events of 30-40% Ru+Ru (a) and Zr+Zr (b) at 200 GeV with different $ n_{5}/s $ inputs. Panel (c) lists $ \sigma^{-1}_{R2} $ vs $ n_{5}/s $, extracted from panels (a) and (b), and the $ \sigma^{-1}_{R2} $ ratios between Ru+Ru and Zr+Zr are shown in panel (d), where the $ 2^{\rm{nd}} $-order-polynomial fit function shows the rising trend starting from (0, 1).

      Figure 8 presents the sensitivity study for the signed balance functions. This approach is not part of the STAR blind analysis, but follows the same procedure as used in the Quark Matter 2019 Conference proceedings [57]. The observables $ r_{\rm{lab}} $ and $ R_{{\rm{B}}} $ as defined in Eqs. (25) and (26) (with $ p_T $ weighting), respectively, are exhibited in panels (a) and (c) as function of $ n_{5}/s $ from the EBE-AVFD model for 30-40% Ru+Ru and Zr+Zr collisions at $ \sqrt{s_{\rm{NN}}} = 200 $ GeV. The corresponding ratios between Ru+Ru and Zr+Zr are shown in panels (b) and (d), respectively. $ r_{\rm{lab}} $ increases with the CME signal in each isobar collision. According to Eqs. (25) and (34) and the core-component comparisons, $ r_{\rm{lab}} $ is related to $ \sqrt{\Delta\gamma_{112}} $, and therefore the $ r_{\rm{lab}} $ ratio between the two systems should roughly obey a $ 2^{\rm{nd}} $-order polynomial function that starts from (0, 1). This relation is demonstrated with the corresponding fit in Fig. 8(b). Panel (d) does not show a clear trend for the ratio of $ R_{\rm{B}}^{\rm{Ru+Ru}}/R_{\rm{B}}^{\rm{Zr+Zr}} $, which is not a complete surprise: $ R_B $ looks for a higher-order effect in the difference between $ r_{\rm{lab}} $ and $ r_{\rm{rest}} $, and thus requires much more statistics than $ r_{\rm{lab}} $.

      Figure 8.  (color online) $ r_{{\rm{lab}}} $ (a) and $ R_{\rm{B}} $ (c) as function of $ n_{5}/s $ from the EBE-AVFD model for 30-40% Ru+Ru and Zr+Zr collisions at $ \sqrt{s_{\rm{NN}}} = 200 $ GeV, with their ratios between Ru+Ru and Zr+Zr in panels (b) and (d), respectively. In panel (b), the $ 2^{\rm{nd}} $-order-polynomial fit function demonstrates the rising trend starting from (0, 1).

      The sensitivity in terms of the statistical significance of $ (O^{\rm{Ru+Ru}}/O^{\rm{Zr+Zr}} -1) $ is listed in Table II for the observables, $ \Delta\gamma_{112} $, $ \Delta\delta $, $ \kappa_{112} $, $ r_{\rm{lab}} $ and $ \sigma_{R2}^{-1} $. This table serves as a reference point to interpret the STAR data of the isobar collisions. In our convention, a statistical significance within (-3, 3) is regarded as "not significant", or "consistent with no difference between the two systems". Note that opposite to other observables, the $ \Delta\delta $ ratio is supposed to be lower than unity in presence of the CME. Therefore, the more negative the statistical significance of $ (\Delta\delta^{\rm{Ru+Ru}}/\Delta\delta^{\rm{Zr+Zr}}-1) $, the more sensitive this observable is to the CME signal. The high sensitivities of the $ \Delta\delta $ ratio reported in Table II could be a special feature of EBE-AVFD, instead of a universal truth, which awaits verification/falsification of real data. In general, $ \kappa_{112} $ roughly doubles the sensitivity of $ \Delta\gamma_{112} $, which, as explained before, should be mostly due to the contribution of $ \Delta\delta $, and needs to be tested by experimental data.

      $ \Delta\gamma_{112} $ and $ r_{\rm{lab}} $ show similar significance values because of the approximate equivalence between them. Note that neither the toy model nor the EBE-AVFD model takes into account the separate CME domains that, instead of merging into a global charge separation for the whole event, still move independently from each other in the fireball. Thus $ r_{\rm{lab}} $ is expected by these models to respond to the CME signal in a similar way as $ \Delta\gamma_{112} $ that only deals with the azimuthal angle. Should the isobar-collision data show a better sensitivity of $ r_{\rm{lab}} $ than that of $ \Delta\gamma_{112} $, it may reveal the CME domains undergoing an incomplete hydrodynamic evolution due to its short duration.

      In the analysis of the $ R $ correlator using the STAR frozen code, $ \sigma_{R2}^{-1} $ yields lower significance than other observables, and this is worth to note in anticipation of the results from the STAR blind analysis. However, this is largely due to two factors in this particular implementation. First, this analysis uses the sub event plane instead of the full event plane as in the $ \Delta\gamma_{112} $ analysis, which leads to worse event plane resolutions and hence larger statistical uncertainties. Second, the particles of interest in the $ R(\Delta S_2) $ analysis come from narrower kinematic regions than other analyses, which further enlarges its statistical errors and reduces its sensitivities. When we repeat the calculations of the $ R $ correlator and the $ \gamma_{112} $ correlator both with the true reaction plane and with the same kinematic cuts (not following the frozen code anymore), $ \sigma_{R2}^{-1}\{{\rm{RP}}\} $ and $ \Delta\gamma_{112}\{{\rm{RP}}\} $ do exhibit comparable significance values, as shown in the last two columns of Table II, consistent with findings in Ref. [48]. Therefore, this result, the toy model studies and Eqs. (21) and (34) all give a coherent picture that on general grounds, the $ \gamma_{112} $ correlator, the $ R $ correlator and the signed balance functions have similar sensitivities, when used on the same set of particles.

    • Several experimental approaches have been developed to search for the CME in heavy-ion collisions. In this paper, we focus on three of them: the $ \gamma_{112} $ correlator, the $ R(\Delta S_2) $ correlator and the signed balance functions. We have established the relation between these methods via analytical derivation, and employed both simple Monte Carlo simulations and the EBE-AVFD model to verify the equivalence between the corecomponents of these observables. Our study also supports the assumption that the CME signal and the background contributions can be linearly added up in such corecomponents. For the observables of $ \Delta\gamma_{112} $, $ \Delta\delta $, $ \kappa_{112} $, $ r_{\rm{lab}} $ and $ \sigma_{R2}^{-1} $, we have extracted their sensitivities to the difference between Ru+Ru and Zr+Zr collisions at $ \sqrt{s_{\rm{NN}}} = 200 $ GeV from 30-40% central events generated by EBE-AVFD. $ \Delta\delta $ and $ \kappa_{112} $ may render better sensitivities than other observables, which could be a model-dependent feature instead of a universal truth, and needs to be further scrutinized by data. The same significance level has been corroborated for $ \Delta\gamma_{112} $, $ r_{\rm{lab}} $ and $ \sigma_{R2}^{-1} $, if put on an equal footing, but the implementation details in the STAR frozen code can cause apparent differences in their sensitivities. Therefore, this study provides a reference point to gauge the STAR isobar-collision data.

    • This effort of investigation was initiated in, and agreed upon by the CME focus group of the STAR Collaboration motivated by the desire to benchmark the various CME observables planned for use in the STAR isobar blind analysis against a well-established model EBE-AVFD. We thank the STAR Collaboration for the permitting the use of the corresponding isobar blind analysis frozen code in this simulation study, and the many collaborators who have contributed to previous STAR CME related studies besides the authors. We are especially grateful to the following people for their substantial help: Kenneth Barish, Helen Caines, Jinhui Chen, William Christie, Frank Geurts, Huanzhong Huang, Hongwei Ke, William Llope, Xiaofeng Luo, Rongrong Ma, Yugang Ma, Bedanga Mohanty, Grigory Nigmatkulov, Lijuan Ruan, Ernst Sichtermann, Yuanfang Wu, Nu Xu, Zhangbu Xu, and Zhenyu Ye.

      We thank the STAR Collaboration, the RHIC Operations Group, and RCF at RHIC for their support.

    APPENDIX A: DERIVATION OF $ \sigma(\Delta B_y) $ AND $ \sigma(\Delta B_x) $
    • To estimate the RMS of $ \Delta B_y $, $ \sigma(\Delta B_y) $, we first go through the following expansion,

      $\tag{A1}\begin{aligned}[b] (N_{x(\alpha\beta)}-N_{x(\beta\alpha)})^2 = & \Bigg\{\sum\limits_{\alpha,\beta} {\rm{Sign}}[\sin(\phi^*_\alpha)-\sin(\phi^*_\beta)] \Bigg\}^2 = \Bigg\{\sum\limits_{\alpha,\beta} {\rm{Sign}}[\sin(\phi^*_\alpha)-\sin(\phi^*_\beta)] \Bigg\} \Bigg\{\sum\limits_{\alpha',\beta'} {\rm{Sign}}[\sin(\phi^*_{\alpha'})-\sin(\phi^*_{\beta'})] \Bigg\} \\ = & \sum\limits_{\alpha\neq\alpha'}\sum\limits_{\beta\neq\beta'} {\rm{Sign}}[\sin(\phi^*_\alpha)-\sin(\phi^*_\beta)]\times{\rm{Sign}}[\sin(\phi^*_{\alpha'})-\sin(\phi^*_{\beta'})] \\& + \sum\limits_{\alpha}\sum\limits_{\beta\neq\beta'} {\rm{Sign}}[\sin(\phi^*_\alpha)-\sin(\phi^*_\beta)]\times{\rm{Sign}}[\sin(\phi^*_{\alpha})-\sin(\phi^*_{\beta'})] \\& + \sum\limits_{\alpha\neq\alpha'}\sum\limits_{\beta} {\rm{Sign}}[\sin(\phi^*_\alpha)-\sin(\phi^*_\beta)]\times{\rm{Sign}}[\sin(\phi^*_{\alpha'})-\sin(\phi^*_{\beta})] + \sum\limits_{\alpha}\sum\limits_{\beta} 1. \end{aligned} $

      Then the average of each term in Eq. (A2) is computed separately in a way similar to Eq. (29). The first one reads

      $\tag{A2}\begin{aligned}[b]& \Big\langle\sum\limits_{\alpha\neq\alpha'}\sum\limits_{\beta\neq\beta'} {\rm{Sign}}[\sin(\phi^*_\alpha)-\sin(\phi^*_\beta)]\times{\rm{Sign}}[\sin(\phi^*_{\alpha'})-\sin(\phi^*_{\beta'})]\Big\rangle \\ = & N_\alpha(N_\alpha-1)N_\beta(N_\beta-1)\times \Bigg\{ \int_{-\pi/2}^{\pi/2} \Bigg[\int_{-\pi-\phi^*_\alpha}^{\phi^*_\alpha} \frac{dN}{d\phi^*_\beta}d\phi^*_\beta-\int^{\pi-\phi^*_\alpha}_{\phi^*_\alpha} \frac{dN}{d\phi^*_\beta}d\phi^*_\beta\Bigg]\frac{dN}{d\phi^*_\alpha}d\phi^*_\alpha \\ & +\int_{\pi/2}^{3\pi/2} \Bigg[\int_{\pi-\phi^*_\alpha}^{\phi^*_\alpha} \frac{dN}{d\phi^*_\beta}d\phi^*_\beta-\int^{3\pi-\phi^*_\alpha}_{\phi^*_\alpha} \frac{dN}{d\phi^*_\beta}d\phi^*_\beta\Bigg]\frac{dN}{d\phi^*_\alpha}d\phi^*_\alpha \Bigg\}^2 = N_\alpha(N_\alpha-1)N_\beta(N_\beta-1)\Bigg[\frac{8}{\pi^2}(1+\frac{2}{3}v_2)(a_{1,\alpha}-a_{1,\beta}) \Bigg]^2. \end{aligned}$

      The second term becomes

      $\tag{A3}\begin{aligned}[b]& \Big\langle\sum\limits_{\alpha}\sum\limits_{\beta\neq\beta'} {\rm{Sign}}[\sin(\phi^*_\alpha)-\sin(\phi^*_\beta)]\times{\rm{Sign}}[\sin(\phi^*_{\alpha})-\sin(\phi^*_{\beta'})]\Big\rangle \\ = & N_\alpha N_\beta(N_\beta-1)\times \Bigg\{ \int_{-\pi/2}^{\pi/2} \Bigg[\int_{-\pi-\phi^*_\alpha}^{\phi^*_\alpha} \frac{dN}{d\phi^*_\beta}d\phi^*_\beta-\int^{\pi-\phi^*_\alpha}_{\phi^*_\alpha} \frac{dN}{d\phi^*_\beta}d\phi^*_\beta\Bigg]^2\frac{dN}{d\phi^*_\alpha}d\phi^*_\alpha \\ & +\int_{\pi/2}^{3\pi/2} \Bigg[\int_{\pi-\phi^*_\alpha}^{\phi^*_\alpha} \frac{dN}{d\phi^*_\beta}d\phi^*_\beta-\int^{3\pi-\phi^*_\alpha}_{\phi^*_\alpha} \frac{dN}{d\phi^*_\beta}d\phi^*_\beta\Bigg]^2\frac{dN}{d\phi^*_\alpha}d\phi^*_\alpha \Bigg\} = N_\alpha N_\beta(N_\beta-1)\Bigg[\frac{1}{3}-\frac{8}{\pi^2}(1+v_2)a_{1,\beta}(a_{1,\alpha}-a_{1,\beta}) \Bigg]. \end{aligned}$

      Similarly, the third term gives

      $\tag{A4} \Big\langle\sum\limits_{\alpha\neq\alpha'}\sum\limits_{\beta} {\rm{Sign}}[\sin(\phi^*_\alpha)-\sin(\phi^*_\beta)]\times{\rm{Sign}}[\sin(\phi^*_{\alpha'})-\sin(\phi^*_{\beta})]\Big\rangle = N_\alpha(N_\alpha-1) N_\beta\Bigg[\frac{1}{3}+\frac{8}{\pi^2}(1+v_2)a_{1,\alpha}(a_{1,\alpha}-a_{1,\beta}) \Bigg]. $

      The last term is simply

      $\tag{A5} \sum\limits_{\alpha}\sum\limits_{\beta} 1 = N_\alpha N_\beta. $

      According to the definition in Eq. (24), we have

      $\tag{A6} \sigma^2(\Delta B_y) = \frac{M^2}{N_+^2 N_-^2} \langle (N_{x(+-)}-N_{x(-+)})^2 \rangle. $

      For simplicity, we assume $ N_+ = N_- = M/2 \gg 1 $ and $ v_2 \ll 1 $, and put Eqs. (A3), (A4), (A5) and (A6) into Eq. (A7), so that

      $\tag{A7} \sigma^2(\Delta B_y) \approx \frac{4(M+1)}{3}+ \frac{64M^2}{\pi^4}\Bigg[\left(1+\frac{2}{3}v_2\right)^2+\frac{\pi^2(1+v_2)}{4M}\Bigg](a_{1,+}-a_{1,-})^2 \\ \approx \frac{4M}{3}+ \frac{64M^2}{\pi^4}\left(1+ \frac{4}{3}v_2\right)(a_{1,+}-a_{1,-})^2. $

      In a similar way, we also reach

      $\tag{A8} \sigma^2(\Delta B_x) \approx \frac{4(M+1)}{3}+ \frac{64M^2}{\pi^4}\Bigg[\left(1-\frac{2}{3}v_2\right)^2+\frac{\pi^2(1-v_2)}{4M}\Bigg](v_{1,+}-v_{1,-})^2 \\ \approx \frac{4M}{3}+ \frac{64M^2}{\pi^4}\left(1- \frac{4}{3}v_2\right)(v_{1,+}-v_{1,-})^2. $

      Then the difference between $ \sigma^2(\Delta B_y) $ and $ \sigma^2(\Delta B_x) $ will render

      $\tag{A9} \sigma^2(\Delta B_y) - \sigma^2(\Delta B_x) \approx \frac{128M^2}{\pi^4}\Bigg[\left(1+\frac{\pi^2}{4M}\right)\Delta\gamma_{112}-\left(\frac{4}{3}+\frac{\pi^2}{4M}\right)v_2\Delta\delta\Bigg] \approx \frac{128M^2}{\pi^4}\left(\Delta\gamma_{112}-\frac{4}{3}v_2\Delta\delta\right). $

Reference (72)



DownLoad:  Full-Size Img  PowerPoint