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

Subikash Choudhury  Fudan University, Shanghai, 200433  Xin Dong  Lawrence
Berkeley National Laboratory, Berkeley, California 94720  Jim Drachenberg
Abilene Christian University, Abilene, Texas 79699  James Dunlop  Brookhaven
National Laboratory, Upton, New York 11973  ShinIchi Esumi  University of
Tsukuba, Tsukuba, Ibaraki 305-8571, Japan  Yicheng Feng  Purdue University,
West Lafayette, Indiana 47907  Evan Finch  Southern Connecticut State
University, New Haven, Connecticut 06515  Yu Hu  Fudan University, Shanghai,
200433  Brookhaven National Laboratory, Upton, New York 11973  Jiangyong Jia
Brookhaven National Laboratory, Upton, New York 11973  State University of New
York, Stony Brook, New York 11794  Jerome Lauret  Brookhaven National
Laboratory, Upton, New York 11973  Wei Li  Rice University, Houston, Texas
77251  Jinfeng Liao  Physics Department and Center for Exploration of Energy
and Matter, Indiana University, 2401 N Milo B. Sampson Lane, Bloomington, IN
47408, USA  Yufu Lin<EMAIL_ADDRESS>College of Physics and
Technology, Guangxi Normal University, Guilin, 541004, China  Central China
Normal University, Wuhan, Hubei 430079, People’s Republic of China  Mike Lisa
Ohio State University, Columbus, Ohio 43210  Takafumi Niida  University of
Tsukuba, Tsukuba, Ibaraki 305-8571, Japan  Robert Lanny Ray  University of
Texas, Austin, Texas 78712  Masha Sergeeva  University of California, Los
Angeles, California 90095  Diyu Shen<EMAIL_ADDRESS>Fudan University,
Shanghai, 200433  Shuzhe Shi  Department of Physics, McGill University, 3600
University Street, Montreal, QC, H3A 2T8, Canada  Paul Sorensen  Brookhaven
National Laboratory, Upton, New York 11973  Aihong Tang  Brookhaven National
Laboratory, Upton, New York 11973  Prithwish Tribedy  Brookhaven National
Laboratory, Upton, New York 11973  Gene Van Buren  Brookhaven National
Laboratory, Upton, New York 11973  Sergei Voloshin  Wayne State University,
Detroit, Michigan 48201  Fuqiang Wang  Purdue University, West Lafayette,
Indiana 47907  Gang Wang  University of California, Los Angeles, California
90095  Haojie Xu  Huzhou University, Huzhou, Zhejiang 313000  Zhiwan Xu
University of California, Los Angeles, California 90095  Nanxi Yao
<EMAIL_ADDRESS>University of California, Los Angeles, California 90095
Jie Zhao  Purdue University, West Lafayette, Indiana 47907

###### 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
isobaric collisions at RHIC.

keywords

chiral magnetic effect, anisotropic flow, heavy-ion collisions, quark-gluon
plasma

##  I Introduction

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},$ (1)

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, 4, 5, 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, 12, 13, 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.

  * •

According to quantum chromodynamics (QCD), if the topological solutions of the
SU(3) gauge group are chiral, they can transfer chirality to quarks via the
chiral anomaly [15, 16], forming local chiral domains with finite $\mu_{5}$ in
a QGP [1, 2, 17, 18, 19, 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\sim
10^{14}$ 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:

$\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^{*})+...],$
(2)

where $\phi$ is the azimuthal angle of a particle, and
$\phi^{*}=\phi-{\rm\Psi_{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,+}$)
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
simultaneously advance several frontiers of fundamental physics: the evolution
of the strong magnetic field, the topological phases of QCD, and the chiral
symmetry restoration of strong interactions. A community-wide interest has
been drawn to search for the CME at RHIC [22, 23, 24, 25, 26, 27, 28, 29, 30,
31, 32] and the LHC [33, 34, 35, 36, 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, 39, 40, 41, 42, 43, 44, 45, 46, 47]. To better gauge
the background, the STAR experiment at RHIC has collected a large data sample
of isobaric collisions, and performed the corresponding CME-related analyses.

The two isobaric systems, namely ${}^{96}_{44}$Ru + ${}^{96}_{44}$Ru and
${}^{96}_{40}$Zr + ${}^{96}_{40}$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, 50, 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
currently all the analysis codes have been frozen as part of the blinding
procedure. Since the isobaric 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 isobaric-collision data coming out of the blind
analyses.

The paper is structured as follows. In Sec. II, we focus on three observables
in search of the CME: the $\gamma$ correlator [53], the $R$ correlator [32,
54] and the signed balance functions [55, 56]. 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) [57, 58, 59]. 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 isobaric collisions under study. Section VI gives the
summary and outlook.

##  II Experimental Observables

Several observables have been proposed to search for the CME in heavy-ion
collisions, including the $\gamma$ correlator [53], the $R(\Delta S_{m})$
correlator [32, 54], and the signed balance functions [55, 56]. 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.

###  II.1 $\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 [53],

$\displaystyle\gamma_{112}$ $\displaystyle\equiv$
$\displaystyle\langle\cos(\phi_{\alpha}+\phi_{\beta}-2{\rm\Psi_{RP}})\rangle$
(3) $\displaystyle=$
$\displaystyle\langle\cos(\phi^{*}_{\alpha})\cos(\phi^{*}_{\beta})-\sin(\phi^{*}_{\alpha})\sin(\phi^{*}_{\beta})\rangle$
$\displaystyle=$ $\displaystyle(\langle v_{1,\alpha}v_{1,\beta}\rangle+B_{\rm
IN})-(\langle a_{1,\alpha}a_{1,\beta}\rangle+B_{\rm OUT}),$

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},$ (4)

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” ($\rm\Psi_{EP}$) reconstructed with detected particles,
and then the measurement is corrected for the finite event plane resolution
[60]. 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 [45, 44].
The flowing resonance picture can be generalized to transverse momentum
conservation (TMC) [43, 42] and local charge conservation (LCC) [44]. Ideally,
the two-particle correlator,

$\displaystyle\delta$ $\displaystyle\equiv$
$\displaystyle\langle\cos(\phi_{\alpha}-\phi_{\beta})\rangle$ (5)
$\displaystyle=$ $\displaystyle(\langle v_{1,\alpha}v_{1,\beta}\rangle+B_{\rm
IN})+(\langle a_{1,\alpha}a_{1,\beta}\rangle+B_{\rm OUT}),$

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\delta$ and
$\Delta\gamma_{112}$ [42]:

$\displaystyle\Delta\delta^{\rm TMC}$ $\displaystyle\approx$
$\displaystyle-\frac{1}{N}\frac{\langle p_{T}\rangle^{2}_{\rm\Omega}}{\langle
p_{T}^{2}\rangle_{\rm
F}}\frac{1+({\bar{v}}_{2,{\rm\Omega}})^{2}-2{\bar{\bar{v}}}_{2,{\rm
F}}{\bar{v}}_{2,{\rm\Omega}}}{1-({\bar{\bar{v}}}_{2,{\rm F}})^{2}},$ (6)
$\displaystyle\Delta\gamma^{\rm TMC}_{112}$ $\displaystyle\approx$
$\displaystyle-\frac{1}{N}\frac{\langle p_{T}\rangle^{2}_{\rm\Omega}}{\langle
p_{T}^{2}\rangle_{\rm
F}}\frac{2{\bar{v}}_{2,{\rm\Omega}}-{\bar{\bar{v}}}_{2,{\rm
F}}-{\bar{\bar{v}}}_{2,{\rm
F}}({\bar{v}}_{2,{\rm\Omega}})^{2}}{1-({\bar{\bar{v}}}_{2,{\rm F}})^{2}}$ (7)
$\displaystyle\approx$ $\displaystyle\kappa^{\rm TMC}_{112}\cdot
v_{2,{\rm\Omega}}\cdot\Delta\delta^{\rm TMC},$

where $\kappa^{\rm
TMC}_{112}=(2{\bar{v}}_{2,{\rm\Omega}}-{\bar{\bar{v}}}_{2,{\rm
F}})/v_{2,{\rm\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 “${\rm\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}.$ (8)

A CME signal will make $\kappa_{112}$ larger than $\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 isobaric 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,

$\displaystyle\gamma_{123}\equiv\langle\langle\cos(\phi_{\alpha}+2\phi_{\beta}-3{\rm\Psi_{3}})\rangle\rangle$
$\displaystyle,~{}\ $
$\displaystyle\kappa_{123}\equiv\frac{\Delta\gamma_{123}}{v_{3}\cdot\Delta\delta},$
(9)
$\displaystyle\gamma_{132}\equiv\langle\langle\cos(\phi_{\alpha}-3\phi_{\beta}+2{\rm\Psi_{2}})\rangle\rangle$
$\displaystyle,~{}\ $
$\displaystyle\kappa_{132}\equiv\frac{\Delta\gamma_{132}}{v_{2}\cdot\Delta\delta},$
(10)

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 [64]. Therefore, in
the following sections we will not extend our method study to these derivative
correlators.

###  II.2 $R$ correlator

The $R(\Delta S_{m})$ correlator [32, 54] 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})}/\frac{N(\Delta S^{\perp}_{m,\rm real})}{N(\Delta S^{\perp}_{m,\rm
shuffled})},~{}\ m=2,3,...,$ (11)

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(\frac{m}{2}\phi^{*}_{m})\rangle_{N_{+}}-\langle\sin(\frac{m}{2}\phi^{*}_{m})\rangle_{N_{-}}.$
(12)

Here $\phi^{*}_{m}=\phi-{\rm\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
${\rm\Psi}_{m}$ is replaced with $({\rm\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 [54]. 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(shuffled)}$, all roughly follow
Gaussian distributions, we can establish the relation between $\sigma_{R2}$
and the RMS values of the core-component distributions:

$\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}.$ (13)

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

$\displaystyle(\Delta S_{2,\rm real})^{2}$ (14) $\displaystyle\equiv$
$\displaystyle(\frac{\sum_{i=1}^{N_{+}}\sin(\phi^{*}_{i})}{N_{+}}-\frac{\sum_{i=1}^{N_{-}}\sin(\phi^{*}_{i})}{N_{-}})^{2}$
$\displaystyle=$
$\displaystyle\frac{\sum_{i=1}^{N_{+}}\sin^{2}(\phi^{*}_{i})+\sum_{i\neq
j}^{N_{+}}\sin(\phi^{*}_{i})\sin(\phi^{*}_{j})}{N_{+}^{2}}+\frac{\sum_{i=1}^{N_{-}}\sin^{2}(\phi^{*}_{i})+\sum_{i\neq
j}^{N_{-}}\sin(\phi^{*}_{i})\sin(\phi^{*}_{j})}{N_{-}^{2}}-\frac{2\sum_{i=1,j=1}^{N_{+},N_{-}}\sin(\phi^{*}_{i})\sin(\phi^{*}_{j})}{N_{+}N_{-}}$
$\displaystyle=$
$\displaystyle\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_{-}}$
$\displaystyle-2\langle\sin(\phi^{*}_{i})\sin(\phi^{*}_{j})\rangle_{N_{+}N_{-}}$

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:

$\displaystyle\langle(\Delta S_{2,\rm real})^{2}\rangle$ $\displaystyle=$
$\displaystyle\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})$
(15) $\displaystyle\approx$ $\displaystyle\frac{2(1-v_{2}-\delta^{\rm
SS}+\gamma_{112}^{\rm SS})}{M}-\Delta\delta+\Delta\gamma_{112}.$ (16)

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,

$\displaystyle\langle(\Delta S_{2,\rm real}^{\perp})^{2}\rangle$
$\displaystyle=$
$\displaystyle\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})$
(17) $\displaystyle\approx$ $\displaystyle\frac{2(1+v_{2}-\delta^{\rm
SS}-\gamma_{112}^{\rm SS})}{M}-\Delta\delta-\Delta\gamma_{112}.$ (18)

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

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

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}$:

$\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(1-\frac{1}{M})\Delta\gamma_{112}.$ (21)

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}.$ (22)

Here we further assume that in each of the four $\langle(\Delta
S_{2})^{2}\rangle$ terms, $\frac{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}$
[65] is nontrivial. An analytical resolution correction has been derived in
Ref. [66]. Eqs. 21 and 22 can provide approximate approaches to examine the
resolution correction factor.

Multiplicity fluctuations widen the $R(\Delta S)$ distributions. To reduce
this effect to $R(\Delta S_{2})$, a scaling was proposed in Ref. [54], used in
experimental data analysis, 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^{{}^{\prime}}_{2})$ distribution by dividing the horizontal axis
by the RMS of the $N(\Delta S_{2,\rm shuffled})$ distribution, i.e., $\Delta
S_{2}^{{}^{\prime}}=\Delta S/\sqrt{\langle(\Delta S_{2,\rm
shuffled})^{2}\rangle}$ (see Sec.V for details). Then, the width of $R(\Delta
S^{{}^{\prime}}_{2})$ becomes

$\frac{S_{\rm concavity}}{\sigma_{R2^{\prime}}^{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}.$
(23)

###  II.3 Signed Balance Functions

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

$\displaystyle\Delta B_{y}$ $\displaystyle\equiv$
$\displaystyle\Big{[}\frac{N_{y(+-)}-N_{y(++)}}{N_{+}}-\frac{N_{y(-+)}-N_{y(--)}}{N_{-}}\Big{]}-\Big{[}\frac{N_{y(-+)}-N_{y(++)}}{N_{+}}-\frac{N_{y(+-)}-N_{y(--)}}{N_{-}}\Big{]}$
(24) $\displaystyle=$
$\displaystyle\frac{N_{+}+N_{-}}{N_{+}N_{-}}[N_{y(+-)}-N_{y(-+)}].$

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}).$ (25)

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},$ (26)

can help differentiate the background from the real CME signal [55]. 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_{\alpha,\beta}{\rm
Sign}[p_{T,\alpha}\sin(\phi^{*}_{\alpha})-p_{T,\beta}\sin(\phi^{*}_{\beta})].$
(27)

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

$\displaystyle\langle N_{y(\alpha\beta)}-N_{y(\beta\alpha)}\rangle$
$\displaystyle\approx$ $\displaystyle
C_{y}\Big{\langle}\sum_{\alpha,\beta}[\sin(\phi^{*}_{\alpha})-\sin(\phi^{*}_{\beta})]\Big{\rangle}$
(28) $\displaystyle=$ $\displaystyle
C_{y}\Big{\langle}[N_{\beta}\sum_{\alpha}\sin(\phi^{*}_{\alpha})-N_{\alpha}\sum_{\beta}\sin(\phi^{*}_{\beta})]\Big{\rangle}$
$\displaystyle=$ $\displaystyle
C_{y}N_{\alpha}N_{\beta}\langle\langle\sin(\phi^{*})\rangle_{N_{\alpha}}-\langle\sin(\phi^{*})\rangle_{N_{\beta}}\rangle.$

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

$\displaystyle\langle N_{y(\alpha\beta)}-N_{y(\beta\alpha)}\rangle$
$\displaystyle=$ $\displaystyle
2\int_{-\pi/2}^{\pi/2}\Big{[}\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}\Big{]}\frac{dN}{d\phi^{*}_{\alpha}}d\phi^{*}_{\alpha}$
(29) $\displaystyle\approx$
$\displaystyle\frac{8}{\pi^{2}}(1+\frac{2}{3}v_{2})N_{\alpha}N_{\beta}(a_{1,\alpha}-a_{1,\beta}).$

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
$\frac{8M}{\pi^{2}}(1+\frac{2}{3}v_{2})\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,$
(30)

and the explicit counting gives

$\displaystyle\langle N_{x(\alpha\beta)}-N_{x(\beta\alpha)}\rangle$
$\displaystyle=$ $\displaystyle
2\int_{0}^{\pi}\Big{[}\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}\Big{]}\frac{dN}{d\phi^{*}_{\alpha}}d\phi^{*}_{\alpha}$
(31) $\displaystyle\approx$
$\displaystyle\frac{8}{\pi^{2}}(1-\frac{2}{3}v_{2})N_{\alpha}N_{\beta}(v_{1,\alpha}-v_{1,\beta}).$

Thus, $C_{x}=8(1-2v_{2}/3)/\pi^{2}$, and $\langle\Delta B_{x}\rangle$ becomes
$\frac{8M}{\pi^{2}}(1-\frac{2}{3}v_{2})\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.

$\displaystyle\sigma^{2}(\Delta B_{y})$ $\displaystyle\approx$
$\displaystyle\frac{4M}{3}+\frac{64M^{2}}{\pi^{4}}(1+\frac{4}{3}v_{2})(a_{1,+}-a_{1,-})^{2}$
(32) $\displaystyle\sigma^{2}(\Delta B_{x})$ $\displaystyle\approx$
$\displaystyle\frac{4M}{3}+\frac{64M^{2}}{\pi^{4}}(1-\frac{4}{3}v_{2})(v_{1,+}-v_{1,-})^{2}.$
(33)

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}}(\Delta\gamma_{112}-\frac{4}{3}v_{2}\Delta\delta).$
(34)

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}-\frac{4}{3}v_{2}\Delta\delta)$ with momentum weighting.

##  III Model descriptions

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

###  III.1 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 [67]. 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 [68]

$v_{2}/\mathcal{N}=a/(1+e^{-[(m_{T}-m_{0})/\mathcal{N}-b]/c})-d,$ (35)

where $\mathcal{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 [69]. 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 [70]. 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},$ (36)

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

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

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

###  III.2 EBE-AVFD

The EBE-AVFD model [57, 58, 59] 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: 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~{}\text{fm}/c$ and mixing parameter $\alpha_{\text{glb}}=0.118$.
The latter factorizes the contributions of the participant density
($n_{\text{part}}$) and binary collision density ($n_{\text{coll}}$) in the
local entropy density,
$s=\alpha_{\text{glb}}\,n_{\text{coll}}+(1-\alpha_{\text{glb}})n_{\text{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 is computed according
to the event-wise proton configuration in the Monte-Carlo Glauber initial
conditions.

The hydrodynamic evolution is solved through two components. The bulk-matter
collective flow is described by the VISH2+1 simulation package [61], with the
lattice equation of state s95p-v1.2, shear-viscosity $\eta/s=0.08$, and
freeze-out temperature $T_{\text{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. [57, 58,
59].

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

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

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. [62]. The approach in Ref. [62] 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_{\text{LCC}}$ is introduced to
characterize the fraction of charged hadrons that are sampled in positive-
negative pairs in the same way as in Ref. [62], while the rest of the hadrons
are sampled independently. Varying the parameter $P_{\text{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 [63], 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_{\text{LCC}}$ is around
$1/3$, and that roughly half of the background correlations come from LCC and
the other half from resonance decays.

##  IV Core-Component Comparisons

We will exploit the toy model and the EBE-AVFD model to simulate the core
components of the experimental observables introduced in Sec. II:
$\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$.

###  IV.1 Toy-model Results

Figure 2: The toy-model simulations of $2\Delta\gamma_{112}$, $\Delta_{R2}$
and $\Delta^{\prime}_{\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.

We put the three core components on the same footing, by plotting
$2\Delta\gamma_{112}$, $\Delta_{R2}$ and $\Delta^{\prime}_{\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}$ 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 [47, 44, 43],
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.

###  IV.2 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^{\prime}_{\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: (a) The EBE-AVFD simulations of $2\Delta\gamma_{112}$, $\Delta_{R2}$
and $\Delta^{\prime}_{\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,-}.$ (39)

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.

##  V Sensitivity study for isobaric collisions

The validity of Eq. 39 for the three experimental observables reassures the
feasibility of disentangling the CME signal from the backgrounds with the
isobaric collisions. To eliminate unintentional biases in the analyses of the
isobaric-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
isobaric-collision data. In this section, we shall apply these frozen codes to
the EBE-AVFD events of isobaric 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 IV. 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,\,\mathrm{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.

Table 1 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.

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

$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

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.

Figure 4: EBE-AVFD calculations of $\gamma_{112}^{\rm OS(SS)}$ (a) and
$\Delta\gamma_{112}$ (b) as functions of $n_{5}/s$ for 30-40% isobaric
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).

Figure 4 presents the EBE-AVFD calculations of $\gamma_{112}^{\rm OS(SS)}$ (a)
and $\Delta\gamma_{112}$ (b) as functions of $n_{5}/s$ for 30-40% isobaric
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 5: EBE-AVFD calculations of $\delta^{\rm OS(SS)}$ (a) and
$\Delta\delta$ (b) as functions of $n_{5}/s$ for 30-40% isobaric 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).

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 [59]. Figure 5 depicts the EBE-AVFD results of $\delta^{\rm OS(SS)}$
(a) and $\Delta\delta$ (b) as functions of $n_{5}/s$ for 30-40% isobaric
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 [59]. This
is because the magnetic field difference between the two isobaric 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 isobaric-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.
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.

Figure 6: 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% isobaric
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).

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% isobaric 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 2. In
general, the EBE-AVFD simulations show good responses to the signal for
$\Delta\gamma_{112}^{\mathrm{Ru+Ru}}/\Delta\gamma_{112}^{\mathrm{Zr+Zr}}$,
$\Delta\delta^{\mathrm{Ru+Ru}}/\Delta\delta^{\mathrm{Zr+Zr}}$ and
$\kappa_{112}^{\mathrm{Ru+Ru}}/\kappa_{112}^{\mathrm{Zr+Zr}}$, and these
promising features await verification by the isobaric-collision data.

Figure 7: Distributions of $R(\Delta S_{2}^{{}^{\prime\prime}})$ 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).

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}^{{}^{\prime}})$
distribution by dividing the horizontal axis by the RMS of the $N(\Delta
S_{2,\rm shuffled})$ distribution, i.e., $\Delta S_{2}^{{}^{\prime}}=\Delta
S_{2}/\sqrt{\langle(\Delta S_{2,\rm shuffled})^{2}\rangle}$. Then $\Delta
S_{2}^{{}^{\prime}}$ is further modified to correct for the event plane
resolution, i.e., $\Delta S^{{}^{\prime\prime}}=\Delta
S^{{}^{\prime}}/\mathrm{\delta_{Res}}$, where $\mathrm{\delta_{Res}}$ is the
correction factor whose details can be found in Ref. [54]. Panels (a) and (b)
show the $R(\Delta S_{2}^{{}^{\prime\prime}})$ 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}^{{}^{\prime\prime}})$ 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}^{{}^{\prime\prime}})$ 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 2 for later discussions.

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 [56]. The
observables $r_{\mathrm{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_{\mathrm{lab}}$ increases with the CME signal in each isobaric collision.
According to Eqs. 25 and 34 and the core-component comparisons,
$r_{\mathrm{lab}}$ is related to $\sqrt{\Delta\gamma_{112}}$, and therefore
the $r_{\mathrm{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_{\mathrm{lab}}$ and $r_{\mathrm{rest}}$,
and thus requires much more statistics than $r_{\mathrm{lab}}$.

Figure 8:  $r_{\mathrm{lab}}$ (a) and $R_{\mathrm{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 2 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 isobaric collisions. 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 2 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_{\mathrm{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
isobaric-collision data show a better sensitivity of $r_{\mathrm{lab}}$ than
that of $\Delta\gamma_{112}$, it may reveal the CME domains undergoing an
incomplete hydrodynamic evolution due to its short duration.

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_{\mathrm{event}}$ means the number of events used for each isobaric
system in the simulation.

$n_{5}/s$ | $N_{\mathrm{event}}$ | $\Delta\gamma_{112}$ | $\Delta\delta$ | $\kappa_{112}$ | $r_{\mathrm{lab}}$ | $\sigma_{R2}^{-1}$
---|---|---|---|---|---|---
0 | $2\times 10^{8}$ | -1.50 | -2.89 | -1.21 | -0.77 | 1.33
0.05 | $4\times 10^{8}$ | 0.62 | -6.16 | 1.37 | 0.47 | 0.29
0.10 | $4\times 10^{8}$ | 1.91 | -16.81 | 3.43 | 3.11 | 0.62
0.20 | $2\times 10^{8}$ | 7.73 | -42.96 | 14.07 | 5.96 | 1.84

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}$ and $\Delta\gamma_{112}$ do exhibit comparable significance
values (not shown in this article), 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.

##  VI Summary

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 isobaric-
collision data.

## Acknowledgement

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. This work is funded by the US Department of Energy under
Grants No. DE-AC02-98CH10886, DE-FG02-89ER40531, DE-FG02-92ER40713, DE-
FG02-88ER40424, DE-SC0012910, DE-SC0013391, DE-SC0020651, and by the National
Natural Science Foundation of China under Contract No. 12025501, 11905059,
12075085, the Strategic Priority Research Program of Chinese Academy of
Science with Grant No. XDB34030200, the Fundamental Research Funds for the
Central Universities under Grant No. CCNU19ZN019, the Ministry of Science and
Technology (MoST) under Grant No. 2016YFE0104800 and the China Scholarship
Council (CSC), Join Large-Scale Scientific Facility Funds of NSFC and CAS
under Contracts No. U2032110, the U.S. Department of Energy, Office of
Science, Office of Nuclear Physics, within the framework of the Beam Energy
Scan Theory (BEST) Topical Collaboration, the U.S. National Science Foundation
under Grant No. PHY-1913729, the Natural Sciences and Engineering Research
Council of Canada, the Fonds de recherche du Québec - Nature et technologies
(FRQNT) through the Programmede Bourses d’Excellencepour Étudiants Étrangers
(PBEEE).

##  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,

$\displaystyle(N_{x(\alpha\beta)}-N_{x(\beta\alpha)})^{2}$ $\displaystyle=$
$\displaystyle\Big{\\{}\sum_{\alpha,\beta}{\rm
Sign}[\sin(\phi^{*}_{\alpha})-\sin(\phi^{*}_{\beta})]\Big{\\}}^{2}$ (41)
$\displaystyle=$ $\displaystyle\Big{\\{}\sum_{\alpha,\beta}{\rm
Sign}[\sin(\phi^{*}_{\alpha})-\sin(\phi^{*}_{\beta})]\Big{\\}}\Big{\\{}\sum_{\alpha^{\prime},\beta^{\prime}}{\rm
Sign}[\sin(\phi^{*}_{\alpha^{\prime}})-\sin(\phi^{*}_{\beta^{\prime}})]\Big{\\}}$
$\displaystyle=$
$\displaystyle\sum_{\alpha\neq\alpha^{\prime}}\sum_{\beta\neq\beta^{\prime}}{\rm
Sign}[\sin(\phi^{*}_{\alpha})-\sin(\phi^{*}_{\beta})]\times{\rm
Sign}[\sin(\phi^{*}_{\alpha^{\prime}})-\sin(\phi^{*}_{\beta^{\prime}})]$
$\displaystyle+\sum_{\alpha}\sum_{\beta\neq\beta^{\prime}}{\rm
Sign}[\sin(\phi^{*}_{\alpha})-\sin(\phi^{*}_{\beta})]\times{\rm
Sign}[\sin(\phi^{*}_{\alpha})-\sin(\phi^{*}_{\beta^{\prime}})]$
$\displaystyle+\sum_{\alpha\neq\alpha^{\prime}}\sum_{\beta}{\rm
Sign}[\sin(\phi^{*}_{\alpha})-\sin(\phi^{*}_{\beta})]\times{\rm
Sign}[\sin(\phi^{*}_{\alpha^{\prime}})-\sin(\phi^{*}_{\beta})]+\sum_{\alpha}\sum_{\beta}1.$

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

$\displaystyle\Big{\langle}\sum_{\alpha\neq\alpha^{\prime}}\sum_{\beta\neq\beta^{\prime}}{\rm
Sign}[\sin(\phi^{*}_{\alpha})-\sin(\phi^{*}_{\beta})]\times{\rm
Sign}[\sin(\phi^{*}_{\alpha^{\prime}})-\sin(\phi^{*}_{\beta^{\prime}})]\Big{\rangle}$
(42) $\displaystyle=$ $\displaystyle
N_{\alpha}(N_{\alpha}-1)N_{\beta}(N_{\beta}-1)\times\Big{\\{}\int_{-\pi/2}^{\pi/2}\Big{[}\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}\Big{]}\frac{dN}{d\phi^{*}_{\alpha}}d\phi^{*}_{\alpha}$
$\displaystyle\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>+\int_{\pi/2}^{3\pi/2}\Big{[}\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}\Big{]}\frac{dN}{d\phi^{*}_{\alpha}}d\phi^{*}_{\alpha}\Big{\\}}^{2}$
$\displaystyle=$ $\displaystyle
N_{\alpha}(N_{\alpha}-1)N_{\beta}(N_{\beta}-1)\Big{[}\frac{8}{\pi^{2}}(1+\frac{2}{3}v_{2})(a_{1,\alpha}-a_{1,\beta})\Big{]}^{2}.$

The second term becomes

$\displaystyle\Big{\langle}\sum_{\alpha}\sum_{\beta\neq\beta^{\prime}}{\rm
Sign}[\sin(\phi^{*}_{\alpha})-\sin(\phi^{*}_{\beta})]\times{\rm
Sign}[\sin(\phi^{*}_{\alpha})-\sin(\phi^{*}_{\beta^{\prime}})]\Big{\rangle}$
(43) $\displaystyle=$ $\displaystyle
N_{\alpha}N_{\beta}(N_{\beta}-1)\times\Big{\\{}\int_{-\pi/2}^{\pi/2}\Big{[}\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}\Big{]}^{2}\frac{dN}{d\phi^{*}_{\alpha}}d\phi^{*}_{\alpha}$
$\displaystyle\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>+\int_{\pi/2}^{3\pi/2}\Big{[}\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}\Big{]}^{2}\frac{dN}{d\phi^{*}_{\alpha}}d\phi^{*}_{\alpha}\Big{\\}}$
$\displaystyle=$ $\displaystyle
N_{\alpha}N_{\beta}(N_{\beta}-1)\Big{[}\frac{1}{3}-\frac{8}{\pi^{2}}(1+v_{2})a_{1,\beta}(a_{1,\alpha}-a_{1,\beta})\Big{]}.$

Similarly, the third term gives

$\displaystyle\Big{\langle}\sum_{\alpha\neq\alpha^{\prime}}\sum_{\beta}{\rm
Sign}[\sin(\phi^{*}_{\alpha})-\sin(\phi^{*}_{\beta})]\times{\rm
Sign}[\sin(\phi^{*}_{\alpha^{\prime}})-\sin(\phi^{*}_{\beta})]\Big{\rangle}$
(44) $\displaystyle=$ $\displaystyle
N_{\alpha}(N_{\alpha}-1)N_{\beta}\Big{[}\frac{1}{3}+\frac{8}{\pi^{2}}(1+v_{2})a_{1,\alpha}(a_{1,\alpha}-a_{1,\beta})\Big{]}.$

The last term is simply

$\sum_{\alpha}\sum_{\beta}1=N_{\alpha}N_{\beta}.$ (45)

According to the definition in Eq. 24, we have

$\sigma^{2}(\Delta
B_{y})=\frac{M^{2}}{N_{+}^{2}N_{-}^{2}}\langle(N_{x(+-)}-N_{x(-+)})^{2}\rangle.$
(46)

For simplicity, we assume $N_{+}=N_{-}=M/2\gg 1$ and $v_{2}\ll 1$, and put
Eqs. 42, 43, 44 and 45 into Eq. 46, so that

$\displaystyle\sigma^{2}(\Delta B_{y})$ $\displaystyle\approx$
$\displaystyle\frac{4(M+1)}{3}+\frac{64M^{2}}{\pi^{4}}\Big{[}(1+\frac{2}{3}v_{2})^{2}+\frac{\pi^{2}(1+v_{2})}{4M}\Big{]}(a_{1,+}-a_{1,-})^{2}$
(47) $\displaystyle\approx$
$\displaystyle\frac{4M}{3}+\frac{64M^{2}}{\pi^{4}}(1+\frac{4}{3}v_{2})(a_{1,+}-a_{1,-})^{2}.$

In a similar way, we also reach

$\displaystyle\sigma^{2}(\Delta B_{x})$ $\displaystyle\approx$
$\displaystyle\frac{4(M+1)}{3}+\frac{64M^{2}}{\pi^{4}}\Big{[}(1-\frac{2}{3}v_{2})^{2}+\frac{\pi^{2}(1-v_{2})}{4M}\Big{]}(v_{1,+}-v_{1,-})^{2}$
(48) $\displaystyle\approx$
$\displaystyle\frac{4M}{3}+\frac{64M^{2}}{\pi^{4}}(1-\frac{4}{3}v_{2})(v_{1,+}-v_{1,-})^{2}.$

Then the difference between $\sigma^{2}(\Delta B_{y})$ and $\sigma^{2}(\Delta
B_{x})$ will render

$\displaystyle\sigma^{2}(\Delta B_{y})-\sigma^{2}(\Delta B_{x})$
$\displaystyle\approx$
$\displaystyle\frac{128M^{2}}{\pi^{4}}\Big{[}(1+\frac{\pi^{2}}{4M})\Delta\gamma_{112}-(\frac{4}{3}+\frac{\pi^{2}}{4M})v_{2}\Delta\delta\Big{]}$
(49) $\displaystyle\approx$
$\displaystyle\frac{128M^{2}}{\pi^{4}}(\Delta\gamma_{112}-\frac{4}{3}v_{2}\Delta\delta).$
(50)

## References

  * [1] D. Kharzeev, Phys. Lett. B 633, 260 (2006).
  * [2] D. E. Kharzeev, L. D. McLerran and H. J. Warringa, Nucl. Phys. A 803, 227 (2008).
  * [3] D.E. Kharzeev, J. Liao, S.A. Voloshin and G. Wang, Prog. Part. Nucl. Phys. 88, 1 (2016).
  * [4] Jie Zhao and Fuqiang Wang, Prog. Part. Nucl. Phys. 107, 200 (2019).
  * [5] Wei Li and Gang Wang, Annual Review of Nuclear and Particle Science 70, no.1, 293 (2020).
  * [6] D. E. Kharzeev and J. Liao, Nature Rev. Phys. 3, no.1, 55-63 (2021).
  * [7] Q. Li, D. E. Kharzeev, C. Zhang, Y. Huang, I. Pletikosic, A. V. Fedorov, R. D. Zhong, J. A. Schneeloch, G. D. Gu and T. Valla, Nature Physics 12, 550 (2016).
  * [8] J. Xiong, S. K. Kushwaha, T. Liang, J. W. Krizan, M. Hirschberger, W. Wang, R. J. Cava, and N. P. Ong, Science 350, 413 (2015).
  * [9] X. Huang, L. Zhao, Y. Long, P. Wang, D. Chen, Z. Yang, H. Liang, M. Xue, H. Weng, Z. Fang, X. Dai, and G. Chen, Phys. Rev. X 5, 031023 (2015).
  * [10] C. Shekhar, F. Arnold, S. Wu, Y. Sun, M. Schmidt, N. Kumar, A. G. Grushin, J. H. Bardarson, R. Donizeth dos Reis, M. Naumann, M. Baenitz, H. Borrmann, M. Nicklas, E. Hassinger, C. Felser, B. Yan, Nature Communications 7, 11615 (2016).
  * [11] J. Adams et al. [STAR Collaboration], Nucl. Phys. A 757, 102 (2005).
  * [12] K. Adcox et al. [PHENIX Collaboration], Nucl. Phys. A 757, 184 (2005).
  * [13] B. B. Back et al. [PHOBOS Collaboration], Nucl. Phys. A 757, 28 (2005).
  * [14] I. Arsene et al. [BRAHMS Collaboration], Nucl. Phys. A 757, 1 (2005).
  * [15] S. L. Adler, Phys. Rev. 177, 2426 (1969).
  * [16] J. S. Bell and R. Jackiw, Nuovo Cim. A 60, 47 (1969).
  * [17] D. Kharzeev and A. Zhitnitsky, Nucl. Phys. A 797, 67 (2007).
  * [18] D. Kharzeev, A. Krasnitz and R. Venugopalan, Phys. Lett. B 545, 298 (2002).
  * [19] I. Iatrakis, S. Lin and Y. Yin, Phys. Rev. Lett. 114, 252301 (2015).
  * [20] K. Fukushima, D. E. Kharzeev and H. J. Warringa, Phys. Rev. Lett. 104, 212001 (2010).
  * [21] U. Heinz and R. Snellings, Ann. Rev. Nucl. Part. Sci. 63, 123 (2013)
  * [22] B. I. Abelev et al. [STAR Collaboration], Phys. Rev. Lett. 103, 251601 (2009).
  * [23] 7\. I. Abelev et al. [STAR Collaboration], Phys. Rev. C 81, 54908 (2010).
  * [24] L. Adamczyk et al. [STAR Collaboration], Phys. Rev. C 88, 064911 (2013).
  * [25] L. Adamczyk et al. [STAR Collaboration], Phys. Rev. C 89, 044908 (2014).
  * [26] L. Adamczyk et al. [STAR Collaboration], Phys. Rev. Lett. 113, (2014) 052302.
  * [27] J. Adam et al. [STAR Collaboration], Phys. Lett. B 798, 134975 (2019).
  * [28] Gang Wang et al. [STAR Collaboration], Nucl. Phys A 904-905, 248c (2013).
  * [29] P. Tribedy [STAR Collaboration], Nucl. Phys. A 967, 740 (2017).
  * [30] Jie Zhao et al. [STAR Collaboration], Nucl. Phys A 535-538, 982 (2019).
  * [31] N. N. Ajitanand, S. Esumi, R. A. Lacey [PHENIX Collaboration], in: Proc. of the RBRC Workshops, vol.96, 230 (2010): “P- and CP-odd ects in hot and dense matter”.
  * [32] N. N. Ajitanand, R. A. Lacey, A. Taranenko and J. M. Alexander, Phys. Rev. C 83, 011901(R) (2011).
  * [33] B. I. Abelev et al. [ALICE Collaboration], Phys. Rev. Lett. 110, 012301 (2013).
  * [34] V. Khachatryan et al. [CMS Collaboration], Phys. Rev. Lett. 118, 122301 (2017).
  * [35] A. M. Sirunyan et al. [CMS Collaboration], Phys. Rev. C 97, 044912 (2018).
  * [36] S. Acharya et al. [ALICE Collaboration], Phys. Lett. B 777, 151 (2018).
  * [37] S. Acharya et al. [ALICE Collaboration], J. High Energ. Phys. 2020, 160 (2020).
  * [38] M. Asakawa, A. Majumder and B. Müller, Phys. Rev. C 81, 064912 (2010).
  * [39] J. Liao, V. Koch and A. Bzdak, Phys. Rev. C 82, 054902 (2010).
  * [40] A. Bzdak, V. Koch and J. Liao, Phys. Rev. C 83, 014905 (2011).
  * [41] H. Petersen, T. Renk, and S. A. Bass, Phys. Rev. C 83, 014916 (2011).
  * [42] A. Bzdak, V. Koch and J. Liao, Lect. Notes Phys. 871, 503 (2013).
  * [43] S. Pratt, S. Schlichting and S. Gavin, Phys. Rev. C 84, 024909 (2011).
  * [44] S. Schlichting and S. Pratt, Phys. Rev. C 83, 014913 (2011).
  * [45] F. Wang, Phys. Rev. C 81, 064902 (2010).
  * [46] F. Wang and J. Zhao, Phys. Rev. C 95, no. 5, 051901(R) (2017).
  * [47] Y. Feng, J. Zhao and F. Wang, Phys. Rev. C 98, no. 3, 034904 (2018).
  * [48] Y. Feng, J. Zhao, H.-J. Xu and F. Wang, Phys. Rev. C 103, 034912 (2021).
  * [49] S. A. Voloshin, Phys. Rev. Lett. 105, 172301 (2010).
  * [50] W.-T. Deng, X.-G. Huang, G.-L. Ma, G. Wang, Phys. Rev. C 94, 041901 (2016).
  * [51] W.-T. Deng, X.-G. Huang, G.-L. Ma, G. Wang, Phys. Rev. C 97, 044901 (2018).
  * [52] J. Adam, et al. the STAR Collaboration, NUCL SCI TECH 32, 48 (2021).
  * [53] S.A. Voloshin, Phys.Rev.C, 70, 057901 (2004).
  * [54] N. Magdy, S.Z. Shi, J.F Liao, N.N Ajitanand and R.A. Lacey, Phys.Rev.C, 97, 061901 (2018).
  * [55] A.H. Tang, Chin.Phys.C, 44, No.5 054101 (2020).
  * [56] Y.F. Lin, [for STAR Collaboration], Nucl. Phys. A 1005, 121828 (2021).
  * [57] S. Shi, Y. Jiang, E. Lilleskov and J. Liao, Annals Phys. 394, 50 (2018) [arXiv:1711.02496 [nucl-th]].
  * [58] Y. Jiang, S. Shi, Y. Yin and J. Liao, Chin. Phys. C 42, no. 1, 011001 (2018) [arXiv:1611.04586 [nucl-th]].
  * [59] S. Shi, H. Zhang, D. Hou and J. Liao, Phys. Rev. Lett. 125, 242301 (2020) [arXiv:1910.14010 [nucl-th]].
  * [60] A. M. Poskanzer and S. A. Voloshin, Phys. Rev. C 58, 1671 (1998).
  * [61] C. Shen, Z. Qiu, H. Song, J. Bernhard, S. Bass and U. Heinz, Comput. Phys. Commun. 199, 61-85 (2016) [arXiv:1409.8164 [nucl-th]].
  * [62] B. Schenke, C. Shen and P. Tribedy, Phys. Rev. C 99, no.4, 044908 (2019) [arXiv:1901.04378 [nucl-th]].
  * [63] M. Bleicher, E. Zabrodin, C. Spieles, S. A. Bass, C. Ernst, S. Soff, L. Bravina, M. Belkacem, H. Weber and H. Stoecker, et al. J. Phys. G 25, 1859-1896 (1999) [arXiv:hep-ph/9909407 [hep-ph]].
  * [64] S. Choudhury, G. Wang, W. He, Y. Hu, H. Z. Huang, Eur. Phys. J. C 80, 383 (2020).
  * [65] N. Magdy, S. Shi, J. Liao, N. Ajitanand, and R.A. Lacey, Phys. Rev. C 97, 061901 (2018).
  * [66] Y. Feng, J. Zhao, H.-j. Xu, and F. Wang, Phys. Rev. C 103, 034912 (2021).
  * [67] B.I. Abelev et al., Phys.ReV.C, 79, 034909 (2009).
  * [68] X. Dong, S. Esumi, P.Sorensen et al., Phys.Lett.B, 597, 328 (2004).
  * [69] T. Sjostrand, S. Mrenna, and P. Z. Skands, JHEP, 0605: 026 (2006).
  * [70] J. Adams et al. Phys.Rev.Lett., 92, 092301 (2004).
  * [71] H. Xu, X. Wang, H. Li, J. Zhao, Z. Lin, C. Shen and F. Wang, Phys. Rev. Lett. 121, 022301 (2018).