Estimate of Background Baseline and Upper Limit on the Chiral Magnetic Effect in Isobar Collisions at $\sqrt{s_{{}_{\textsc{NN}}}}=200$ GeV at the Relativistic Heavy-Ion Collider

Abstract

For the search of the chiral magnetic effect (CME), STAR previously presented the results from isobar collisions (${{}^{96}_{44}\text{Ru}}+{{}^{96}_{44}\text{Ru}}$, ${{}^{96}_{40}\text{Zr}}+{{}^{96}_{40}\text{Zr}}$) obtained through a blind analysis. The ratio of results in Ru+Ru to Zr+Zr collisions for the CME-sensitive charge-dependent azimuthal correlator ($\Delta\gamma$), normalized by elliptic anisotropy ($v_{2}$), was observed to be close to but systematically larger than the inverse multiplicity ratio. The background baseline for the isobar ratio, $Y=\frac{(\Delta\gamma/v_{2})^{\text{Ru}}}{(\Delta\gamma/v_{2})^{\text{Zr}}}$, is naively expected to be $\frac{(1/N)^{{\rm Ru}}}{(1/N)^{{\rm Zr}}}$; however, genuine two- and three-particle correlations are expected to alter it. We estimate the contributions to $Y$ from those correlations, utilizing both the isobar data and hijing simulations. After including those contributions, we arrive at a final background baseline for $Y$, which is consistent with the isobar data. We extract an upper limit for the CME fraction in the $\Delta\gamma$ measurement of approximately $10\%$ at a $95\%$ confidence level on in isobar collisions at $\sqrt{s_{\textsc{nn}}}=$ 200 GeV, with an expected 15% difference in their squared magnetic fields.

1 Introduction

The quantum chromodynamics (QCD) predicts vacuum fluctuations, rendering nonzero topological charges in a local domain with odd $\mathcal{P}$ and $\mathcal{CP}$ symmetries [1, 2, 3, 4], which may be pertinent to the matter-antimatter asymmetry of our universe [5]. As a result, there would be more particles with one certain chirality than the other, which is called chirality imbalance. If there is a strong magnetic field, the spins of particles would be locked either parallel or anti-parallel to the magnetic field direction, depending on their charges. Then, with the same chirality imbalance and opposite spin directions, the positive and negative charged particles would have opposite momentum directions. This charge separation phenomenon is called the chiral magnetic effect (CME) [6].

In heavy-ion collisions, we refer to the collided nucleons as participants and the others as spectators. From the collision geometry, the magnetic field created by the spectator protons is generally perpendicular to the reaction plane (rp, spanned by the impact parameter and beam direction), so the CME searches usually take rp as the reference direction. The particle distribution can be expanded w.r.t. rp in azimuth ($\psi_{{\textsc{rp}}}$) into a Fourier series

$$\frac{2\pi}{N}\frac{\text{d}N}{\text{d}\phi}=1+2v_{1}\cos(\phi-\psi_{{\textsc{rp}}})+2v_{2}\cos 2(\phi-\psi_{{\textsc{rp}}})+\ldots\,,$$ (1)

where $N$ is the number of particles and $\phi$ is the azimuthal angle in the plane perpendicular to the beam axis.

The charge-dependent azimuthal correlator $\Delta\gamma$ is used to measure the charge separation of CME [7]. Its definition is as follows:

$$\gamma_{\alpha\beta}\equiv\langle\cos(\phi_{\alpha}+\phi_{\beta}-2\psi_{{\textsc{rp}}})\rangle\,.$$ (2)

Here, the subscripts $\alpha$, $\beta$ represent two different particles of interest (poi) in the same event, and their electric charge signs determine whether the pair ($\alpha\beta$) is opposite-sign (os) or same-sign (ss). For a given pair sign (os or ss), the angle brackets $\langle\cdots\rangle$ represent the average over those pairs within one event and then further averaged over multiple events in each centrality. The signature of CME is that $\gamma_{\textsc{os}}>0$ and $\gamma_{\textsc{ss}}<0$ with the same magnitude [7]. However, there are correlation backgrounds that can also cause the $\gamma$ correlators to deviate from zero [7]. Many of these backgrounds are charge independent, such as those arising from momentum conservation. In order to remove those charge-independent backgrounds, the difference between os and ss is computed as follows:

$$\Delta\gamma\equiv\gamma_{{\textsc{os}}}-\gamma_{{\textsc{ss}}}\,.$$ (3)

Experiments have focused on measuring the $\Delta\gamma$ observable. Large signals of $\Delta\gamma$ have been measured at the Relativistic Heavy-Ion Collider (RHIC) [8, 9, 10, 11, 12, 13, 14, 15, 16] and the Large Hadron Collider (LHC) [17, 18, 19, 20].

Although charge-independent backgrounds are canceled in $\Delta\gamma$, backgrounds remain from charge-dependent two-particle (2p) correlations coupled with elliptic flow of those correlation sources such as resonance decays and jets [7, 21, 22, 23, 24]. These backgrounds can be expressed as

$$\gamma_{\textsc{os}}^{\rm bkgd}=\frac{N_{\rm 2p}}{N_{\textsc{os}}}\left\langle{\cos(\phi_{\alpha}+\phi_{\beta}-2\phi_{{\rm 2p}})}\right\rangle_{{\rm 2p}}v_{2,{\rm 2p}}\,.$$ (4)

Here, $N_{\rm 2p}$, $\phi_{{\rm 2p}}$, and $v_{2,{\rm 2p}}$ represent the number, azimuthal angle, and elliptic flow parameter of those os ${\rm 2p}$ correlation sources (e.g., $\rho$ resonances), respectively. Elliptic flow $v_{2}$ is defined as

$$v_{2}=\left\langle{\cos(2\phi-2\psi_{{\textsc{rp}}})}\right\rangle\,,$$ (5)

which is also a coefficient in Eq. (1). It has been found that those backgrounds dominate charge-separation measurements [16, 24].

To mitigate these backgrounds, STAR conducted experiments of isobar collisions ${{}^{96}_{44}\text{Ru}}+{{}^{96}_{44}\text{Ru}}$ and ${{}^{96}_{40}\text{Zr}}+{{}^{96}_{40}\text{Zr}}$ at $\sqrt{s_{{}_{\textsc{NN}}}}=200$ GeV in 2018 [25]. The choice of these two isobaric species was based on the strategy of keeping the background constant while changing the signal. It was anticipated that the CME-related signal would be larger in Ru+Ru due to the larger number of protons creating a stronger magnetic field. It was also initially expected that with the same number of nucleons in Zr and Ru, flow-related backgrounds would be the same in the two species’ collisions, ideally with charge-independent quantities like multiplicity and elliptic flow ($v_{2}$) the same between the two isobars. However, the STAR data shows a few percent difference between the two isobaric species for both quantities. This discrepancy arises because Ru and Zr have different nuclear sizes and structures, which were predicted by energy density functional theory (DFT) calculations [26, 27, 28]. For the charge-dependent and CME-sensitive quantity, the STAR data show the isobar ratio (Ru+Ru/Zr+Zr) of $\Delta\gamma/v_{2}$ below unity [25]. This seems opposite to the initial expectation of larger CME in Ru+Ru but, in fact, results from the different multiplicities. If the number of background correlation sources is proportional to multiplicity ($N$), then background $\Delta\gamma$ would be diluted by $N$ because $\Delta\gamma$ is a pair-wide average and the pair multiplicity $N_{\textsc{os}},N_{\textsc{ss}}\propto N^{2}$. Then, the naive background baseline of $\frac{(\Delta\gamma/v_{2})^{{\rm Ru}}}{(\Delta\gamma/v_{2})^{{\rm Zr}}}$ would be $\frac{(1/N)^{{\rm Ru}}}{(1/N)^{{\rm Zr}}}$, and Ref. [25] finds that the former is larger than the latter, though with large uncertainties, suggesting a small CME signal [29].

However, this $1/N$ scaling for $\Delta\gamma/v_{2}$ is approximate. It assumes that the number of background sources is proportional to multiplicity, $N_{\rm 2p}\propto N$. This assumption may not hold for the two isobar collision systems and may be violated differently in the two due to their slightly different energy densities. A more precise way of scaling may be to divide $\Delta\gamma/v_{2}$ by $r=N_{\rm 2p}/N_{\textsc{os}}$ (Sec. 3.1). Then, the STAR data show the isobar ratio of $\frac{1}{r}\Delta\gamma/v_{2}$ below unity [25], indicating more complications to be considered [30]. The fact that scaling by $r$ does not give the same result as scaling by multiplicity indicates that the assumption of proportionality is not necessarily valid–there are residual backgrounds in addition to the naive baseline of unity.

Consequently, it is critical to pin down the exact background baseline in order to access information on the CME from the $\Delta\gamma$ measurements in isobar collisions. It is the goal of this paper and Ref. [31] to arrive at a rigorous estimate of the background baseline and examine what the isobar data entail with respect to the possible CME signal. Reference [31] summarizes the essential findings, and this paper provides necessary details of the analysis work.

2 Background baseline

In heavy-ion experiments, rp is unknown. One often reconstructs an event plane (ep) from the particle momentum distribution, taking this reconstructed ep as a proxy for rp. More elegantly, one exploits particle cumulants; instead of reconstructing an ep, one can calculate $v_{2}$ from the two-particle correlator

$$v_{2}^{*}=\sqrt{\left\langle{\cos 2(\phi_{\alpha}-\phi_{\beta})}\right\rangle}\,.$$ (6)

This $v_{2}$ measurement includes nonflow backgrounds, such as two-particle correlations from jets or resonance decays, whose contribution we quantify by

$$\epsilon_{\rm nf}=\frac{v_{2,{\text{nf}}}^{2}}{v_{2}^{2}}=\frac{{v_{2}^{*}}^{2}}{v_{2}^{2}}-1\,.$$ (7)

Here, $v_{2}^{*}$ stands for the measurement that includes nonflow correlations unrelated to the global collision geometry, while $v_{2}$ refers to the true flow. Note that the cumulant measurement of $v_{2}^{*}$ is similar to that obtained by reconstructing an ep, but with subtle differences [32, 33]. It should also be noted that the elliptic flow measured by the ep method is also contaminated by nonflow, as the ep is reconstructed using final-state particles [34]. However, the decomposition of nonflow from flow in the ep method is less straightforward than that in the cumulant method [35].

Similarly, instead of calculating the $\gamma$ correlation using Eq. 2, one can measure a three-particle (3p) correlator

$$C_{3,\alpha\beta}=\langle\cos(\phi_{\alpha}+\phi_{\beta}-2\phi_{c})\rangle\,.$$ (8)

Here, $c$ represents a third particle that is different from $\alpha$, $\beta$. The $\gamma$ correlator can be calculated using $\gamma_{\alpha\beta}=C_{3,\alpha\beta}/v_{2,c}$, where $v_{2,c}$ is the elliptic flow of particle of type $c$, given by Eq. (5). Effectively, particle $c$ serves as the event plane, and its resolution is simply equal to the particle’s elliptic flow $v_{2,c}$. In practice, $\gamma$ is determined by dividing by the measured quantity $v_{2}^{*}$,

$$\gamma_{\alpha\beta}=C_{3,\alpha\beta}/v_{2}^{*}\,.$$ (9)

The main background in $C_{3}$ ($=C_{3,{\textsc{os}}}-C_{3,{\textsc{ss}}}$) arises from the flow-induced background [16, 24]. In this scenario, some of the poi’s are correlated with one another via a 2p source and these 2p sources are all correlated with each particle $c$ through the global flow correlation. This flow-induced background is described by Eq. (4). In addition to this flow-induced background, there is contamination in $C_{3}$ from genuine 3p correlations, where the three particles $\alpha$, $\beta$ and $c$ are intrinsically correlated. Thus, the background contributions to the 3p correlators can be expressed as:

$$\begin{array}[]{llllll}C_{3,{\textsc{os}}}^{\text{bkgd}}=&\frac{N_{\textsc{ss}}}{N_{\textsc{os}}}&\gamma_{\textsc{ss}}v_{2,c}&+\frac{N_{\rm 2p}}{N_{\textsc{os}}}C_{{\rm 2p},{\textsc{os}}}v_{2,{\rm 2p}}v_{2,c}&+&\frac{N_{{\rm 3p},{\textsc{os}}}}{N_{\textsc{os}}N_{c}}C_{{\rm 3p},{\textsc{os}}}\,,\\ C_{3,{\textsc{ss}}}^{\text{bkgd}}=&&\gamma_{\textsc{ss}}v_{2,c}&&+&\frac{N_{{\rm 3p},{\textsc{ss}}}}{N_{\textsc{ss}}N_{c}}C_{{\rm 3p},{\textsc{ss}}}\,,\end{array}$$ (10)

where the first terms in both lines are the charge-independent backgrounds, such as momentum conservation, which will be largely canceled out when taking ${\textsc{os}}-{\textsc{ss}}$. The shorthand notations stand for

$$\begin{array}[]{lcll}C_{{\rm 2p},{\textsc{os}}}&=&\left\langle{\cos(\phi_{\alpha}+\phi_{\beta}-2\phi_{{\rm 2p}})}\right\rangle_{{\rm 2p},{\textsc{os}}}&,\\ C_{{\rm 3p},{\textsc{os}}}&=&\left\langle{\cos(\phi_{\alpha}+\phi_{\beta}-2\phi_{c})}\right\rangle_{{\rm 3p},{\textsc{os}}}&,\\ C_{{\rm 3p},{\textsc{ss}}}&=&\left\langle{\cos(\phi_{\alpha}+\phi_{\beta}-2\phi_{c})}\right\rangle_{{\rm 3p},{\textsc{ss}}}&\end{array}$$ (11)

where the average $\langle\cdots\rangle_{{\rm 2p},{\textsc{os}}}$ runs only over the correlated background pairs with parent cluster azimuth $\phi_{\rm 2p}$. The averages $\langle\cdots\rangle_{{\rm 3p},{\textsc{os}}}$ and $\langle\cdots\rangle_{{\rm 3p},{\textsc{ss}}}$ run over only the correlated background triplets, whose multiplicities are $N_{{\rm 3p},{\textsc{os}}}$ and $N_{{\rm 3p},{\textsc{ss}}}$, respectively. In other words, these quantities characterize the angular properties of the correlated clusters and are not diluted by combinatorial multiplicities. In Eq. (10), the set of os pairs consists of two components: the first component is the correlated 2p os pairs, while the second component comprises the remaining os pairs that are identical to the ss pairs (in terms of the $\gamma$ quantity). The number of correlated 2p pairs is denoted by $N_{\rm 2p}=N_{\textsc{os}}-N_{\textsc{ss}}$. The correlated 3p triplets also contribute nonflow background to $\Delta\gamma$ and are treated as a separate term in Eq. (10). It is worth noting that the $v_{2}$’s in Eq. (10) are the true elliptic flows, as they arise from the correlations between the common global symmetry of the correlated 2p source ($v_{2,{\rm 2p}}$) and the particle $c$ ($v_{2,c}$), giving rise to the contribution to $C_{3}$. In this study, $v_{2,c}=v_{2}$ because the $c$ particles use the same cut as poi.

The $\gamma$ correlators are calculated from $C_{3}$ by Eq. (9). The backgrounds in $\Delta\gamma/v_{2}^{*}$ can be expressed using Eq. (10) as

$$\frac{\Delta\gamma_{\rm bkgd}}{v_{2}^{*}}=\frac{C_{\rm 2p}}{N}\cdot\frac{v_{2}^{2}}{v_{2}^{*2}}+\frac{C_{\rm 3p}}{N}\cdot\frac{1}{N_{c}v_{2}^{*2}}.$$ (12)

The first term comes from the two-particle (2p) nonflow backgrounds (like resonance decay daughter pairs, and os pairs from intra-jet correlation)

$$\frac{C_{{\rm 2p}}}{N}=\frac{N_{{\rm 2p}}}{N_{\textsc{os}}}\left(C_{{\rm 2p},{\textsc{os}}}\frac{v_{2,{\rm 2p}}}{v_{2}}-\frac{\gamma_{{\textsc{ss}}}}{v_{2}}\right)\,,$$ (13)

where $v_{2,{\rm 2p}}$ is the elliptic flow parameter of the correlated pairs (like resonance decays). The second term arises from the three-particle (3p) nonflow backgrounds, such as jets

$$\frac{C_{\rm 3p}}{N}=\frac{N_{{\rm 3p},{\textsc{os}}}}{N_{\textsc{os}}}C_{{\rm 3p},{\textsc{os}}}-\frac{N_{{\rm 3p},{\textsc{ss}}}}{N_{\textsc{ss}}}C_{{\rm 3p},{\textsc{ss}}}\,,$$ (14)

where the $N$ is the multiplicity of poi, and $N_{c}$ is that of particle $c$ (in this analysis $N=N_{c}$). Both the Eqs. (13) and (14) have their corresponding ss contributions subtracted from their os components, as how $\Delta\gamma$ is defined.

With this decomposition, the isobar ratio due to these background sources can be calculated for $\Delta\gamma/v_{2}^{*}$. After approximation to the leading order, the expression becomes

$$Y_{\text{bkgd}}\equiv\frac{\left(\Delta\gamma_{\rm bkgd}/v_{2}^{*}\right)^{{\rm Ru}}}{\left(\Delta\gamma_{\rm bkgd}/v_{2}^{*}\right)^{{\rm Zr}}}\approx 1+\frac{\delta(C_{\rm 2p}/N)}{C_{\rm 2p}/N}-\frac{\delta\epsilon_{\rm nf}}{1+\epsilon_{\rm nf}}\\ +\frac{1}{1+\frac{Nv_{2}^{2}}{C_{\rm 3p}/C_{\rm 2p}}}\left(\frac{\delta C_{\rm 3p}}{C_{\rm 3p}}-\frac{\delta C_{\rm 2p}}{C_{\rm 2p}}-\frac{\delta N}{N}-\frac{\delta v_{2}^{2}}{v_{2}^{2}}\right)\,,$$ (15)

where $\delta X\equiv X^{{\rm Ru}}-X^{{\rm Zr}}$ for any $X=C_{{\rm 3p}}$, $C_{{\rm 2p}}$, etc., while all other quantities without “$\delta$” refer to those in Zr+Zr.

We note that the details of 2p and 3p correlation sources can be complicated. For example, besides 2- and 3-body decays of resonances, they can also come from multi-particle decays of clusters or jets. The os and ss background pairs, besides the extra os pairs from decays of neutral objects, may still not be strictly symmetric (e.g., decays from the $\Delta$ resonances). In addition, the decay kinematics themselves can be altered in heavy-ion collisions, such as by the possible global spin alignment of $\rho$ mesons [36]. Our analysis formalisms, however, do not rely on those details, but only on the overall os-ss difference in the $\gamma$ correlators and the overall nonflow contribution to $v_{2}^{*}$.

3 Analysis

In order to estimate the background baseline, we need to analyze the isobar data to assess the quantities required by Eq. (15). Since the purpose is to estimate the background in the measurements of the STAR isobar blind analysis, we use the same datasets, follow the same event and track selections, and apply identical analysis cuts as described in Ref. [25].

The majority of the analysis cuts are the same for all $\Delta\gamma/v_{2}$ measurements in the blind analysis. These are as follows: On the event level, the minimum-bias trigger is used–correlated hits in both Vertex Position Detectors (VPD) within a time window. The primary vertex reconstructed by the Time Projection Chamber (TPC) [37, 38] is also required to have a longitudinal $z$ position between $-35$ cm and 25 cm ($-35<V_{z}<25$ cm) and a transverse position to be within 2 cm ($V_{\perp}<2$ cm) with respect to the center of TPC. In addition, the interaction position measured online by the Vertex Position Detector (VPD) [39] is required to be within 5 cm from the reconstructed primary vertex along the beam ($|V_{z}-V_{z}^{\textsc{vpd}}|<5$ cm). On the track level, the reconstructed particle tracks must have more than 15 space points measured in the TPC. As the CME is a signal in primordial particles, the track’s distance of closest approach to the primary vertex (DCA) is required to be less than 3 cm ($\text{DCA}<3$ cm). The particle transverse momentum is restricted within $0.2<p_{T}<2$ GeV/$c$. The pseudorapidity range is limited within $|\eta|<1$ for the full-event analyses. For the subevent analyses, where the event is divided into two subevents based on their $\eta$ ranges, the measurements in the STAR blind analysis used slightly different ranges. However, for the estimates reported in this paper, all subevent measurements used the same $\eta$ ranges; specifically, $-1<\eta<-0.1$ was considered as the East subevent, and $0.1<\eta<1$ was considered as the West subevent. The centrality is defined by the number of tracks in the rapidity range $-0.5<\eta<0.5$.

In the STAR blind analysis of the isobar data  [25], a total of seven $\Delta\gamma/v_{2}$ measurements were performed by four research Groups. Among these measurements, four utilized the two-particle cumulant method for the $v_{2}$ measurement and the three-particle cumulant method for the $\Delta\gamma$ measurement. The remaining three measurements used the event-plane method. While both the cumulant and event-plane methods produced similar results, the event-plane method posed significantly more challenges in distinguishing and accounting for nonflow contributions compared to the cumulant method, as aforementioned. For this reason, we concentrate on the former four measurements, namely, the full-event measurements from Group-2 and Group-3, and the subevent measurements from Group-2 and Group-4. The analysis details of the four measurements differ slightly, and those details are tabulated in the first block of Table 2 for the corresponding measurements. In our background baseline estimates for the four measurements, we maintain consistency by using the identical analysis cuts for each corresponding measurement. The block in Table 2 also includes the results for the average over the 20-50% centrality range for the isobar ratios of $\Delta\gamma/v_{2}$, denoted as $Y\equiv\left\langle{\frac{(\Delta\gamma/v_{2})^{\rm Ru}}{(\Delta\gamma/v_{2})^{\rm Zr}}}\right\rangle$, which were obtained from the STAR blind analysis [25]. The average of each term in Eq. 15 over centrality bins is calculated separately using the inverse statistical uncertainty squared as weight, and then summed together to get the background baseline.

Equation (15) suggests categorizing the nonflow contributions to the background into three ingredients: (1) $\delta(C_{{\rm 2p}}/N)/(C_{{\rm 2p}}/N)$ which characterizes the relative difference of flowing clusters between the two isobars, (2) Differences that arise from using $v_{2}^{*}$ rather than true flow in the calculation of $\Delta\gamma$, characterized by $\epsilon_{\rm nf}$, and (3) differences in the relative amounts (or character) of three particle clusters between the isobars. In the next section we will discuss each of these three in turn.

3.1 Pair versus single multiplicity difference

Figure 1: The relative excess of os over ss pion pair multiplicity, $r=(N_{{\textsc{os}}}-N_{{\textsc{ss}}})/N_{{\textsc{os}}}$, as a function of the pair invariant mass ($m_{\rm inv}$) for the 30-40% centrality in Ru+Ru and Zr+Zr collisions (upper panels). The lower panels show the ratio of $r$ in Ru+Ru to that in Zr+Zr collisions. Full-event results are shown on the left, corresponding to the Group-3 analysis cuts. Subevent results are shown on the right, corresponding to the Group-2 analysis cuts. Only statistical uncertainties are shown. The pion pair is required to have an $\eta$ gap $\Delta\eta_{\alpha\beta}>0.05$. For subevent method, the pion pair needs to come from the same subevent.

Figure 2: The Ru+Ru to Zr+Zr ratio of $r=(N_{{\textsc{os}}}-N_{{\textsc{ss}}})/N_{{\textsc{os}}}$ and that of the efficiency-corrected inverse multiplicity $1/N$ as functions of centrality for full-event (left panel, with Group-3 cuts) and subevent (right panel, with Group-2 cuts) methods. Vertical bars are statistical uncertainties; hollow boxes on $r$ data points are systematic uncertainties.

The main background in a $\Delta\gamma$ measurement comes from the charge-independent correlated pairs, such as resonance decays and intra-jet correlations. The number of those pairs is denoted as $N_{{\rm 2p}}=N_{{\textsc{os}}}-N_{{\textsc{ss}}}$, the excess of os pairs over ss, so the charge-independent part is removed (similar to $\Delta\gamma=\gamma_{{\textsc{os}}}-\gamma_{{\textsc{ss}}}$). The relative excess is

$$r=\frac{N_{{\rm 2p}}}{N_{{\textsc{os}}}}=\frac{N_{{\textsc{os}}}-N_{{\textsc{ss}}}}{N_{{\textsc{os}}}}\,.$$ (16)

This is reflected in the background contribution to $\Delta\gamma$ of Eq. (15).

The $C_{\rm 2p}$ in Eq. (15) refers to the background contribution, and the most relevant quantity is $\delta C_{\rm 2p}$. The ZDC (Zero Degree Calorimeter) [40] measurement of $N\Delta\gamma/v_{2}$ is close to the background $C_{\rm 2p}$ because the possible CME signal is small and ZDC has a large $\eta$-gap with TPC (Sec. 3.3). However, the difference in CME signal contributions in $\delta C_{\rm 2p}$ between the isobar systems may not be negligible compared to that in the background contributions, which is what we need for the baseline estimation. Therefore, we cannot simply take the difference between the ZDC measurements in Ru+Ru and Zr+Zr collisions. Even without this complication, due to the low resolution of the ZDC event plane, the statistical uncertainties of the ZDC $C_{{\rm 2p}}$ measurements are too large for the precision needed to achieve the background baseline estimate. Therefore, we go to the $C_{\rm 2p}/N$ term itself as defined in Eq. (13). The $C_{{\rm 2p},{\textsc{os}}}$ represents average angular correlations per 2p cluster (determined by decay kinematics in the case of a resonance decay), and is therefore insensitive to the collision species. The $v_{2}$’s should well scale between various particle/cluster types. The $\gamma_{{\textsc{os}}}/v_{2}$ is relatively insignificant compared to $C_{{\rm 2p},{\textsc{os}}}$. It is therefore reasonably safe to assume that the quantity in the parentheses of Eq. (13), i.e., $\left(C_{{\rm 2p},{\textsc{os}}}\frac{v_{2,{\rm 2p}}}{v_{2}}-\frac{\gamma_{\textsc{ss}}}{v_{2}}\right)$, is the same for Ru+Ru and Zr+Zr. Thus, we have

$$\frac{\delta(C_{{\rm 2p}}/N)}{C_{{\rm 2p}}/N}\approx\frac{\delta r}{r}\,.$$ (17)

In this analysis, we only use pions to calculate $r$, because the main source of this background is the $\rho$ meson decay and the baryon stopping effect does not affect the produced pions. To identify a track as a pion, the TPC reconstructed track is required to match with a TOF (Time Of Flight detector) [39] hit, and additional selections are the TPC energy loss deviation $n\sigma_{\pi}<3$, TOF mass $m^{2}<0.1\text{ GeV}^{2}/c^{4}$, and transverse momentum $0.2<p_{T}<1.8\text{ GeV}/c$. To study the invariant pair mass ($m_{\rm inv}$) dependence, $r$ is calculated as a function of $m_{\rm inv}$ (Fig. 1) for Ru+Ru and Zr+Zr separately, and then the isobar ratio is taken between the two. A constant fit is used to extract the average of the ratio. By applying the same procedure to each centrality bin, we can get their centrality dependence as shown in Fig. 2.

The quantity in the parentheses of Eq. (13), assumed to be equal between the two isobar systems, may have a $m_{\rm inv}$ dependence. $Y_{\text{bkgd}}$ therefore can be described by the average $\delta r/r$ weighted by a $m_{\rm inv}$-dependent function. We estimate the systematic uncertainty by the variation in $r$ obtained by changing the fit range to $m_{\rm inv}<1\text{ GeV}/c^{2}$. The resultant difference from the default is expanded to be symmetric and assigned as part of the systematic uncertainty. This is listed in Table 1 together with other main sources of systematic uncertainties.

The isobar ratio of $1/N$ is also plotted in Fig. 2, which is different from the isobar ratio of $r$ as mentioned in Sec. 1 and observed in Ref. [25]. This difference arises from the fact that the background source multiplicity (such as the $\rho$ mesons) does not scale identically with multiplicity for Ru+Ru and Zr+Zr collisions. As a result, the background does not strictly scale as the inverse multiplicity, and thus $Y_{\text{bkgd}}$ deviates from unity. The amount of deviation is the difference between the two curves in each panel of Fig. 2.

Table 1: The main sources of systematic uncertainties in $Y_{\text{bkgd}}$. Those contributions from one-sided variations are expanded to be symmetric in the calculation of systematic uncertainties.

$Y_{\text{bkgd}}$ syst. source/variation	affected quantity	Group-2 FE	Group-3 FE	Group-2 SE	Group-4 SE
$m_{\text{inv}}<1$ GeV/$c^{2}$ fit range	$r$	$\pm 0.00061$	$\pm 0.00061$	$\pm 0.00055$	$\pm 0.00052$
flow decorrelation $\pm 3\%$	$v_{2}$, $\epsilon_{\rm nf}$	$\pm 0.00040$	$\pm 0.00040$	$\pm 0.00030$	$\pm 0.00035$
alternative 2D fits	$v_{2}$, $\epsilon_{\rm nf}$	$\pm 0.00016$	$\pm 0.00107$	$\pm 0.00058$	$\pm 0.00065$
possible $5\%$ CME in $C_{{\rm 2p}}$	$C_{{\rm 2p}}$	$\pm 0.00070$	$\pm 0.00070$	$\pm 0.00053$	$\pm 0.00060$
HIJING jet-quenching off	$C_{{\rm 3p}}$	$\pm 0.00064$	$\pm 0.00064$	$\pm 0.00051$	$\pm 0.00140$

3.2 Nonflow contamination in $v_{2}^{*}$ measurements

Figure 3: The left plots show the $(\Delta\eta,\Delta\phi)$ 2D distributions (colored) for ss pairs after acceptance and kink corrections with fit result (black mesh). The middle and right plots are the projections of data (black markers) and fit results (red for total fit function, blue for flow component) to $\Delta\eta$, $\Delta\phi$ directions respectively. The upper row is for Ru+Ru and lower row for Zr+Zr. Only the full-event method and ss pairs are shown here. The fitted flow is assumed to be the same between os and ss pairs, and the same between full-event and subevent methods. The centrality range $30$–$40\%$ is used here, and other centrality ranges are similar.

Figure 4: The ${v_{2}^{*}}^{2}$ components as functions of centrality for Ru+Ru (left panel) and Zr+Zr (right panel). The blue curves are the fitted $v_{2}^{2}$. The black and red solid markers are the inclusive ${v_{2}^{*}}^{2}$ for ss pairs from data and fit function respectively, and they agree with each other. The black and orange open squares are the inclusive ${v_{2}^{*}}^{2}$ for all pairs from this measurement and the cited Ref. [25] Group-3 (see the plot legends), and they also agree with each other. Vertical bars and hollow boxes show the statistical and systematic uncertainties, respectively.

The two-particle correlation is usually used to calculate the elliptic flow in the TPC. However, this method cannot avoid nonflow backgrounds that also correlate with those two particles. To separate the true flow from this inclusive measurement, a data-driven approach is used in this analysis. We fit the 2D two-particle $(\Delta\eta,\Delta\phi)$ distribution, where $\Delta\eta$ is the pseudorapidity difference between the two particles, and $\Delta\phi$ azimuth difference. Since the global anisotropy, flow, is supposed to be charge-independent, the os and ss pairs should have the same true flow. We only focus on the ss pairs to avoid the charge-dependent backgrounds that are stronger than charge-independent ones. We assume the true flow is $\eta$-independent at $|\eta|<1$, so we only consider the full event method in this section, and use the same fitted flow value to treat the subevent method. To avoid confusion, we will use “fitted flow” or “fitted $v_{2}$” to refer to our estimate, which should closely reflect the underlying true flow.

In the range, $|\eta|<1$, the single particle $\eta$ distribution is roughly uniform. Therefore, the $\Delta\eta$ projection of the 2D distribution is mainly a triangle due to the finite $\eta$ range. This acceptance effect also happens in mixed events. Each pair has one particle from the current event and the other from another similar event (in the same centrality bin, same $V_{z}$ bin of width 1 cm). However, all other correlations do not exist in those mixed events. Thus, we use the 2D distribution from mixed events, with its peak at $\Delta\eta=0$ scaled to one, to correct the acceptance effect in real events by taking the ratio of the real over mixed. We assume that the fitted flow does not have a $\Delta\eta$ dependence, so this operation does not affect the fitted flow measurement.

After this acceptance correction, the $\Delta\eta$ distribution generally appears flat, revealing the fine structures of nonflow (shown in Figs. 3a, 3d). For each centrality bin, we can get such a corrected 2D distribution, where we see some nonphysical kinks at $|\Delta\eta|=1\pm 0.5$. The centrality of this STAR dataset is defined by the charged particle multiplicity within $|\eta|<0.5$ [25, 41]. On the other hand, the poi in this analysis are within $|\eta|<1$. This distinction means that the centrality bin implicitly imposes an additional constraint on the particle number within $|\eta|<0.5$ but not for $0.5<|\eta|<1$. Consequently, these differences cause those artificial kinks. To eliminate this effect, we project the corrected 2D distribution to $\Delta\eta$ from the away side range of $0.8\pi<|\Delta\phi|<\pi$, where the fine structure of nonflow is small and can be ignored. We then subtract its integral from this projection to remove the pedestal, effectively making it a 2D distribution, independent of $\Delta\phi$. By taking the difference between the acceptance-corrected 2D distribution and this projection, we can eliminate the kinks and any $\phi$-independent detector effects. It is important to note that this operation does not affect the true flow, as the fitted flow is assumed to be independent of $\Delta\eta$.

With all those operations above, the 2D fit can be conducted, and the fit function from observation and tuning is

$$\begin{split}f(\Delta\eta,\Delta\phi)=&A_{1}G\left(\frac{|\Delta\eta|-\mu}{\sigma_{1}}\right)G\left(\frac{\Delta\phi}{\rho_{1}}\right)+A_{2}G\left(\frac{\Delta\eta}{\sigma_{2}}\right)G\left(\frac{\Delta\phi}{\rho_{2}}\right)+\left[A_{3}G\left(\frac{\Delta\eta}{\sigma_{3}}\right)+A_{4}G\left(\frac{\Delta\eta}{\sigma_{4}}\right)\right]G\left(\frac{\Delta\phi}{\rho_{3}}\right)\\ &+C\big[1+2V_{1}\cos(\Delta\phi)+2V_{2}\cos(2\Delta\phi)+2V_{3}\cos(3\Delta\phi)\big]\,,\end{split}$$ (18)

where $G(x)=e^{-x^{2}/2}$ is the Gaussian function. The second line represents the flow pedestal, and the parameter $V_{n}$ corresponds to the squared fitted $v_{n}$ ($n=1,2,3$) assuming the “true flow” does not depend on $\eta$. The 2D Gaussians are empirical models for nonflow corrections, guided by the data shape (Fig. 3a, 3d): track merging effect and Coulomb effect for SS pairs could result in the dip at $\Delta\eta=0$; HBT, resonance decays, and intra-jet correlations are short-range shown by the narrow peak; inter-jet correlations are long-range characterized by a wide $\Delta\eta$ Gaussian. In Eq. (18), all the parameters ($A$, $C$, $V$, $\mu$, $\sigma$, $\rho$) are free in the fitting, starting from reasonable initial values. Figure 3 shows the 2D fit results (left) along with their $\Delta\eta$ (middle), $\Delta\phi$ (right) projections in the centrality bin 30-40% for Ru+Ru (upper) and Zr+Zr (lower) separately.

We note that correlation studies by 2-dimensional $(\Delta\eta,\Delta\phi)$ distributions have been performed previously [42, 43, 44, 45, 46, 47, 48, 49]. The analysis procedure is well established and produces consistent nonflow contributions. Other analyses to quantify nonflow contributions have also been carried out, for example, by utilizing reflection symmetry in $\eta$ in symmetric heavy ion collisions [50], by varying two-particle or subevent $\eta$ gap [51, 52, 53, 54], and by extrapolating from proton-proton, proton-nucleus, and peripheral heavy ion collisions to more central collisions assuming inverse multiplicity scaling of nonflow [55].

Figure 4 shows the ${v_{2}^{*}}^{2}$ components as functions of centrality for Ru+Ru and Zr+Zr separately. As listed in the legend, this study closely reproduced the previous STAR measurement for ${v_{2}^{*}}^{2}$ (Ref. [25], Group-3) using all pairs (os+ss), ; the difference is negligible and the numerical difference can be attributed to nonidentical datasets dynamically accessed at run time. The fitted flow is extracted from the ss pair correlations. Thus, the plot also includes the ${v_{2}^{*}}^{2}$ measured only from ss pairs, together with the ${v_{2}^{*}}^{2}$ calculated from the total fit function, both of which are found to be consistent with each other. Since multiple corrections are applied before conducting the fits, the fit results have been folded back to be comparable to the data. As also shown in the projection plots in Fig. 3, the fitted flow is one component of the total fit function, and all the rest are regarded as nonflow components in this study.

The nonflow fraction $\epsilon_{\rm nf}$ can therefore be calculated by Eq. (7) from the fitted flow $v_{2}^{2}$ and the inclusive ${v_{2}^{*}}^{2}$ measurements from the STAR isobar blind analysis [25]. Figure 5 presents the $\epsilon_{\rm nf}$ for the two isobars and the isobar difference $\frac{-\delta\epsilon_{\rm nf}}{1+\epsilon_{\rm nf}}$ for both full event and subevent cases. The boxes represent the systematic uncertainties.

Figure 5: The relative nonflow strength, $\epsilon_{\text{nf}}=v_{2,{\text{nf}}}^{2}/v_{2}^{2}$, in isobar collisions as functions of centrality for full-event data (upper left panel, with Group-3 cuts) and subevent data (upper right panel, with Group-2 cuts). The lower panels depict the corresponding $-\delta\epsilon_{\rm nf}/(1+\epsilon_{\rm nf})$, one of the contributions to the deviation of $Y_{\text{bkgd}}$ from unity; see Eq. (15). Vertical bars and hollow boxes show the statistical and systematic uncertainties, respectively.

Since the STAR data in Ref. [25] already include the systematics from selection variations, this study does not duplicate them in the baseline estimation to avoid double counting. Instead, it focuses on considering the systematic uncertainties arising from sources specific to the background baseline estimation of this work. These sources include fit uncertainties and model dependencies. The full-event inclusive ${v_{2}^{*}}^{2}$ can be calculated from the 2D distribution (os+ss) in this study, represented by the open black curves in Fig. 4. These results are essentially a replication of the measurement conducted by Group-3 in the STAR isobar blind analysis [25], shown as open orange curves in Fig. 4. In addition, a $3\%$ flow decorrelation over one unit of pseudorapidity has been observed [56]. As a result, $\pm 3\%$ variation of $v_{2}^{2}$ is also considered as a systematic uncertainty. Additionally, an alternative 2D fit is used with a different functional form (and corresponding kink correction). The fitted flow deviation from the default is considered as part of the systematic uncertainty. Those systematic sources and their contributions are listed in Table 1, where the contributions from one-sided variations are expanded to be symmetric in the calculation of systematic uncertainties.

Figure 6: (Left panel) The isobar Ru+Ru/Zr+Zr ratio of the fitted $v_{2}$ flow parameters as a function of centrality in black solid circles, compared to those of the inclusive $v_{2}^{*}$ measurements from the cited STAR isobar blind analysis [25] in colored squares. Full-event (FE) measurements are shown in solid squares, and subevent (SE) measurements are shown in open squares where the two particles come from two different subevents. Group-2 (FE) exploited Gaussian fits to reduce nonflow contributions in their measurement, resulting in larger statistical uncertainties and isobar ratios closer to the SE measurements. (Right panel) The isobar ratio of the extracted nonflow components $v_{2,{\text{nf}}}$ scaled by the square root of multiplicity as a function of centrality. Vertical bars and hollow boxes show the statistical and systematic uncertainties, respectively. The line connecting one set of the data points is to guide the eye.

It is of interest to examine the relative strength of the fitted flows of the two isobar systems. Figure 6 displays the isobar ratio of the fitted $v_{2}$ parameters, where the systematic uncertainties are represented as boxes. The fitted $v_{2}$ values averaged over the 20-50% centrality range are $0.0561\pm 0.0010$ and $0.0548\pm 0.0010$ for Ru+Ru and Zr+Zr collisions, respectively, where the quoted uncertainties are dominated by systematic uncertainties. The average ratio within the 20-50% centrality range is $1.0215\pm 0.0004({\rm stat.})\pm 0.0009({\rm syst.})$. The difference in $v_{2}$ between the isobar systems originates from variations in the initial collision geometries. These differences can be attributed to distinct nuclear structures, as predicted by DFT calculations [26, 57, 28]. Both hydrodynamic calculations [58] and transport models [27] can produce the isobar $v_{2}$ ratios similar to data, including the subtle hump structure of the ratio in the medium centrality region, once the DFT calculated densities are implemented in those models. The significant difference in the most central collisions arises predominantly from nuclear deformations. For comparison, the isobar ratio of the inclusive ${v_{2}^{*}}$ from STAR data [25] are also plotted in Fig. 6. The measured $v_{2}^{*}$ ratios are significantly smaller than those of the fitted, presumably true $v_{2}$. This is because nonflow contamination in Ru+Ru is smaller than that in Zr+Zr due to the higher charged particle multiplicity dilution. This is shown in Fig. 6 by the triangles where most data points are below unity. Factoring out the multiplicity dilutions, Fig. 6(b) shows the ratio of $\sqrt{N}v_{\rm 2,nf}$, which reflects the genuine difference in nonflow correlations between the isobar systems. The ratio is in fact larger than unity for most centralities, which is consistent with the presumably larger energy density achieved in Ru+Ru than Zr+Zr collisions.

Azimuthal anisotropies in central heavy-ion collisions are particularly sensitive to nuclear deformations [59, 60]. It is noteworthy that the difference between the measured $v_{2}^{*}$ ratio and the fitted one is significant also in the most central collisions, as shown in Fig. 6(a). This suggests that using comparisons of the measured $v_{2}^{*}$ ratio to hydrodynamic or other model calculations, which often fail to describe nonflow contributions, to infer nuclear deformations should be taken with caution [61, 62]. It is interesting to observe that there is no significant difference between full-event and subevent $v_{2}^{*}$ results in the most central collisions. This similarity arises because the relative nonflow contribution to $v_{2}^{*}$ is similar between full-event and subevent for those central collisions, as shown in Fig. 5.

3.3 Three-particle correlation background

Figure 7: The $C_{{\rm 3p}}$ in Ru+Ru and Zr+Zr collisions as functions of centrality for full-event (left panel, with Group-3 cuts) and subevent (right panel, with Group-2 cuts) analysis, obtained from hijing simulations of 7.0 billion minimum-bias events for each system. The lower panels show the Ru+Ru over Zr+Zr ratio of $C_{{\rm 3p}}$. Vertical bars indicate statistical uncertainties. The hollow boxes show systematic uncertainties, estimated from hijing quenching-off simulations. The two poi’s have an $\eta$ gap, $\Delta\eta_{\alpha\beta}>0.05$. For the subevent method, the two poi’s come from the same subevent, while the reference particle comes from the other subevent.

The three-particle (3p) background correlation is the last piece needed to form the background baseline estimate, but it is challenging to measure due to the significant combinatorial background. We resort to the hijing (Heavy Ion Jet INteraction Generator) model [63, 64], which simulates parton-parton hard scatterings based on perturbative QCD and gives a reasonable description of partonic energy loss in the QGP medium (jet quenching). Since hijing does not have collective flow, the correlator $C_{3}\equiv C_{3,{\textsc{os}}}-C_{3,{\textsc{ss}}}$ in hijing is entirely composed of genuine 3-particle correlations, $C_{3,{\textsc{os}}}=\frac{N_{{\rm 3p},{\textsc{os}}}}{N_{\textsc{os}}N}C_{{\rm 3p},{\textsc{os}}},C_{3,{\textsc{ss}}}=\frac{N_{{\rm 3p},{\textsc{ss}}}}{N_{\textsc{ss}}N}C_{{\rm 3p},{\textsc{ss}}}$ (c.f. Eq. (10)). The pathlength-dependent jet quenching does produce some degree of anisotropy in the final-state particle azimuthal distribution, sensitive to the initial geometry and indistinguishable from collective anisotropy. This anisotropy, however, is negligible compared to the effect of the 3p correlations [65]. The 3p correlation strength, $C_{{\rm 3p}}$, can be readily obtained from the 3p correlator in hijing,

$$C_{\rm 3p}=N^{2}C_{3}\,.$$ (19)

We use hijing version v1.411 and simulate 7.0 billion events each for Ru+Ru and Zr+Zr collisions at $\sqrt{s_{{}_{\textsc{NN}}}}=200$ GeV. The nuclear structure density distributions are given by energy density functional theory calculations [26, 27, 28], and they are implemented in the initial geometry setup in hijing. Only the final-state charged pions, kaons, protons and their antiparticles are used for centrality definition and poi. The centrality is defined by the multiplicity distribution of particles with $|\eta|<0.5$. The same analysis cuts as the STAR isobar data analysis [25] (see Table 2) are used to process hijing simulation data. Figure 7 shows the $C_{\rm 3p}$ in isobar collisions as functions of centrality, as obtained from hijing simulations, and the Ru+Ru over Zr+Zr ratios of $C_{\rm 3p}$. The centrality dependence is weak, as one would expect for $C_{\rm 3p}$, which is defined for the correlated triplets with the dilution effect factored out as in Eq. (14). The $C_{\rm 3p}$ is larger in full-event than subevent because the number of triplets drops with acceptance more rapidly than single multiplicity.

Figure 8: Comparison of $C_{{\rm 3p}}$ between hijing (jet quenching on) and the cited STAR measurements [25] in Ru+Ru collisions for full-event (left panel, with Group-3 cuts) and subevent (right panel, with Group-2 cuts) analyses (Zr+Zr is similar). The data measurements on the plots use Eq. (19), which are inclusive. The peripheral data are similar to the hijing result, suggesting that the peripheral data are dominated by 3p correlation background. The non-peripheral data are dominated by flow-induced background, absent from hijing. The $C_{\rm 3p}$ results from hijing simulation with quenching turned off are also shown, which is taken as the maximum systematic uncertainty. The two poi’s have an $\eta$ gap, $\Delta\eta_{\alpha\beta}>0.05$. For the subevent method, the two poi’s come from the same subevent, while the reference particle comes from the other subevent. To give a magnitude assessment, the hijing (quenching on) results without the $\Delta\eta_{\alpha\beta}>0.05$ cut (as in Group-4 data analysis [25]) are also displayed.

An important question is how well hijing describes data in terms of 3p correlations? The question cannot be directly answered because of the difficulties to access the 3p correlations in real data as aforementioned. However, there are a few checks one can make. We first examine how well hijing describes the $C_{3}$ in peripheral collisions, where the nonflow effects dominate due to smaller multiplicity dilution. We show the measured $C_{3}$ in data and in hijing in Fig. 8 for full-event and subevent analysis. As shown in Fig. 8, the 70-80% peripheral data are reasonably well described by hijing. This suggests that the peripheral $C_{3}$ data are dominated by 3p correlations, the flow-induced background is small in peripheral collisions. hijing is a reasonable model for (mini-)jet production as well as soft physics via string fragmentation, so it is considered as suitable description for 3p correlations.

The default setup of our hijing simulations include jet quenching. We also simulate hijing with jet quenching turned off. The quenching-off $C_{{\rm 3p}}$ is about 20% higher than the quenching-on result, as shown in Fig. 8. We take this difference as the maximum systematic uncertainty on 3p correlation estimate. As another check, we also show in Fig. 8 the hijing $C_{3}$ without cutting on $\Delta\eta_{\alpha\beta}$ (as in the Group-4 data analysis [25]). The difference is small, and we conclude that hijing with quenching-on and -off provide a safe estimate of the one-side maximum systematic uncertainty. We thus assign their difference divided by $\sqrt{3}$, assuming a uniform probability for the systematic uncertainty, to be the one standard deviation systematic uncertainty on our 3p correlation background estimate, and expand it to be symmetric, as listed in Table 1.

Figure 9: The $C_{{\rm 2p}}$ in isobar collisions for the full-event (left panel, from the Group-3 isobar result) and subevent (right panel, with Group-2 cuts) analyses, obtained from the ZDC measurements of $\Delta\gamma/v_{2}$ multiplied by the efficiency-corrected multiplicity. Vertical bars indicate statistical uncertainties. The asymmetric systematic uncertainty, shown by boxes, is composed of the data systematic uncertainty from the cited Ref. [25] and a one-sided systematic uncertainty of $-5\%$ to account for possible CME contributions in the ZDC measurements. The two poi’s in $\Delta\gamma$ have an $\eta$ gap, $|\Delta\eta_{\alpha\beta}|>0.05$. For subevent method, the poi pair comes from the same subevent.

Figure 10: The ratio $C_{{\rm 3p}}/C_{{\rm 2p}}$ as functions of centrality for full-event (left panel, with Group-3 cuts) and subevent methods (right panel, with Group-2 cuts). Some of the central and peripheral data points are off the plots with large error bars. Vertical bars and hollow boxes show the statistical and systematic uncertainties, respectively.

In order to estimate the baseline contribution from 3p correlations, we also need the quantity $C_{\rm 2p}$ (Eq. (15)). As shown in Eq. (13), $C_{{\rm 2p}}$ represents the 2p nonflow background correlations in $C_{3}$ measurements. If the azimuth of rp is known, then Eq. (5) gives the true elliptic flow and Eq. (2) gives $\Delta\gamma$ without 3p nonflow by definition. If so, Eq. (12) simply becomes

$$N\frac{\Delta\gamma_{\text{bkgd}}\{{\textsc{rp}}\}}{v_{2}\{{\textsc{rp}}\}}=C_{{\rm 2p}}\,.$$ (20)

In STAR, the ZDC is separated from TPC by a large $\eta$ gap and at high collision energies measures only spectator neutrons. Thus, the ZDC measurements of the event plane are not correlated with poi, and the $N\Delta\gamma/v_{2}$ w.r.t. ZDC can be used to estimate $C_{{\rm 2p}}$. Figure 9(a) shows the full-event ZDC measurement from the STAR isobar blind analysis (Group-3) [25]. Figure 9(b) shows the subevent ZDC measurement from this analysis because the subevent ZDC measurement in the STAR blind analysis [25] used $0.05<|\eta|<1$ instead of $0.1<|\eta|<1$.

The poi multiplicities have been corrected for the centrality- and $p_{T}$-dependent track reconstruction efficiencies of the STAR TPC. The average efficiency is on the order of $85\%$. The efficiency is obtained from Monte Carlo tracks simulated by geant in the STAR detector and embedded into the isobar data on the pixel level with proper detector response simulations. In this embedding sample, the input tracks and reconstructed tracks are distributed as functions of centrality, particle species, charge, and $p_{T}$, $\eta$, which are slightly different between the two isobars. For each centrality bin, we obtain the integral (total number of $p\bar{p}$, $K^{\pm}$, $\pi^{\pm}$) inside our cuts (poi) for the reconstructed and input tracks, and their ratio gives us the efficiency. Due to the $\eta$ gap, the subevent has slightly different efficiency compared to the full event.

The systematic uncertainties shown in Fig. 9 include those on the ZDC measurements from the blind analysis [25] for both full-event and subevent analyses. The ZDC measurements could contain some CME signal, possibly on the order of a few percent [16, 66]. Because $C_{\rm 2p}$ refers to 2p background correlation, we additionally assign a one-sided systematic uncertainty of $-5\%$ on $C_{{\rm 2p}}$. This one-sided uncertainty is expanded to be symmetric in the calculation of systematic uncertainties on $Y_{\text{bkgd}}$, as listed in Table 1. The uncertainties on the ZDC measurement of $C_{{\rm 2p}}$ are not included in those on $Y_{\text{bkgd}}$ but only on $Y$ [25], to avoid double counting.

Figure 10 shows the ratios of $C_{\rm 3p}/C_{\rm 2p}$ as a function of centrality in both full-event and subevent analyses. The ratios are on the order of a few percent. The ratio in the full event is larger than in the subevent due to a simple acceptance effect–the smaller the acceptance, the smaller the high-order correlations.

With all the ingredients ready, as in Figs. 2, 9, 6, and 7, one can easily calculate the background contribution from 3p correlations. The results are shown in Fig. 11. The prefactor is depicted in the left column, the sum of the various isobar differences is depicted in the middle column, and the final 3p background difference is depicted in the right column.

Figure 11: The nonflow components related to 3p nonflow as functions of centrality for full-event (upper row, with Group-3 cuts) and subevent (lower row, with Group-2 cuts) methods. The multiplicity $N$ has been corrected for tracking efficiency. Vertical bars and hollow boxes show the statistical and systematic uncertainties, respectively.

4 Result

The previous section discusses all the ingredients needed for background estimation. Besides the unity, there are three terms in Eq. (15). These terms are presented in Fig. 2, Fig. 5, and Fig. 11, separately. The sum of all those terms by Eq. (15) is the background baseline estimate. It is shown in Fig. 12 as a function of centrality, together with the STAR data from Group-3 full-event and Group-2 subevent measurements. We also apply the same procedure for the other two measurements in the STAR isobar blind analysis [25], namely Group-2 full-event and Group-4 subevent measurements.

Figure 12: The $Y\equiv\frac{(\Delta\gamma/v_{2}^{*})^{{\rm Ru}}}{(\Delta\gamma/v_{2}^{*})^{{\rm Zr}}}$ measurements from the cited Ref. [25] (orange curves) and background baseline $Y_{\text{bkgd}}$ from this study (black curves) as functions of centrality for full-event (left panel, with Group-3 cuts) and subevent (right panel, with Group-2 cuts) methods. Vertical bars and hollow boxes show the statistical and systematic uncertainties, respectively.

We compute an average background baseline over the centrality range of 20-50%. Each of the three background terms in Eq. 15 is averaged first, weighted by the corresponding inverse squared statistical uncertainty. Then, the three terms are added to yield the average $Y_{\text{bkgd}}$ baseline. The average $Y_{\text{bkgd}}$ baselines are tabulated in Table 2 for the four measurements and plotted on the summary plot in Fig. 13.

Table 2: Isobar measurements in the 20-50% centrality range from the STAR blind analysis [25] and the corresponding background baseline estimates from this work, for the four measurements using cumulant analysis techniques (two-particle cumulant for $v_{2}^{*}$ and three-particle cumulant for $\Delta\gamma$). Kinematic cuts [25] are $|\eta|<1$ for full event (FE) and $0.1<|\eta|<1$ for subevent (SE), both with $0.2<p_{T}<2$ GeV/$c$. The first block lists the important analysis cuts [25], slightly differing among the four measurements, along with the isobar ratios of the $\left\langle{\Delta\gamma/v_{2}}\right\rangle$ measurements [25] (denoted as $Y$). The second block tabulates the ingredients used for the background baseline estimates, either from data or hijing simulations with the corresponding analysis cuts, along with final background baseline ($Y_{\rm bkgd}$) as well as the magnitudes of its three components in Eq. (15). The last block lists the background subtracted signal ($Y_{\text{signal}}$) and the CME upper limit at the 95% confidence level assuming the CME signal difference between Ru+Ru and Zr+Zr is 15%. The first quoted uncertainty is statistical and the second systematic.

	Group-2 FE	Group-3 FE	Group-2 SE	Group-4 SE
$\Delta\gamma$ cuts	$|\Delta\eta_{\alpha\beta}|>0.05$		$|\Delta\eta_{\alpha\beta}|>0.05$	–
$v_{2}^{*}$ cuts	$|\Delta\eta_{c}|>0.05$ & Gaus. fit	–	–	Same-sign only
$Y\equiv\left\langle{\frac{(\Delta\gamma/v_{2}^{*})^{\rm Ru}}{(\Delta\gamma/v_{2}^{*})^{\rm Zr}}}\right\rangle$	$0.9658\pm 0.0050\pm 0.0007$	$0.9733\pm 0.0040\pm 0.0010$	$0.9611\pm 0.0070\pm 0.0016$	$0.9629\pm 0.0050\pm 0.0003$
$\left\langle{\epsilon_{\rm nf}}\right\rangle$	$0.2528\pm 0.0027\pm 0.0489$	$0.3419\pm 0.0008\pm 0.0555$	$0.1812\pm 0.0007\pm 0.0464$	$0.1948\pm 0.0007\pm 0.0465$
$\left\langle{C_{\rm 2p}}\right\rangle$ (ZDC)	$0.8302\pm 0.0511\pm 0.0415$		$0.5236\pm 0.0375\pm 0.0262$	$0.6097\pm 0.0540\pm 0.0305$
$\left\langle{C_{\rm 3p}}\right\rangle$ (hijing)	$0.0707\pm 0.0005\pm 0.0084$		$0.0133\pm 0.0002\pm 0.0013$	$0.0135\pm 0.0002\pm 0.0013$
$\left\langle{\delta C_{\rm 3p}/C_{\rm 3p}}\right\rangle$ (hijing)	$-0.0054\pm 0.0102\pm 0.0049$		$-0.0035\pm 0.0262\pm 0.0065$	$-0.0180\pm 0.0241\pm 0.0073$
$\left\langle{\delta r/r}\right\rangle$	$-0.0329\pm 0.0003\pm 0.0007$		$-0.0308\pm 0.0004\pm 0.0006$	$-0.0323\pm 0.0003\pm 0.0006$
$\left\langle{-\delta\epsilon_{\rm nf}/(1+\epsilon_{\rm nf})}\right\rangle$	$0.0097\pm 0.0028\pm 0.0001$	$0.0162\pm 0.0008\pm 0.0013$	$0.0080\pm 0.0008\pm 0.0007$	$0.0075\pm 0.0008\pm 0.0008$
the third term	$-0.0144\pm 0.0017\pm 0.0011$		$-0.0108\pm 0.0028\pm 0.0008$	$-0.0125\pm 0.0025\pm 0.0016$
$Y_{\text{bkgd}}$	$0.9625\pm 0.0033\pm 0.0013$	$0.9689\pm 0.0019\pm 0.0016$	$0.9664\pm 0.0029\pm 0.0011$	$0.9628\pm 0.0027\pm 0.0018$
$Y_{\text{signal}}$	$0.0033\pm 0.0060\pm 0.0014$	$0.0044\pm 0.0043\pm 0.0019$	$-0.0052\pm 0.0075\pm 0.0020$	$0.0001\pm 0.0061\pm 0.0018$
$f_{\textsc{cme}}^{{\rm Ru}}$ upper limit	$11.5\%$	$10.3\%$	$8.3\%$	$9.8\%$

Figure 13: Compilation of the Ru+Ru to Zr+Zr isobar ratios of the $\Delta\gamma/v_{2}$ measurements, $Y$ (black squares, with statistical uncertainties indicated by the vertical bars and systematic uncertainties by hollow boxes), and the $r$ and $1/N$ measurements (purple diamonds), with the left coordinate labeling, from the cited STAR blind analysis [25]. The estimated background baselines from this analysis, $Y_{\text{bkgd}}$, for the four cumulant measurements of the isobar $\Delta\gamma/v_{2}$ ratios are shown by the horizontal bars (central values) and the shaded areas (total uncertainties). The total uncertainties are the quadratic sum of the statistical and systematic uncertainties on the background baseline estimates [25].

The difference between the STAR data from Ref. [25] and the baseline from this study, $Y_{\text{signal}}=Y-Y_{\text{bkgd}}$, reflects the relative isobar difference of the possible CME signals in the inclusive $\Delta\gamma$ measurements. Simple algebra indicates

$$\begin{split}Y_{\rm signal}=&Y_{\text{bkgd}}\frac{\delta f_{\textsc{cme}}}{1-f_{\textsc{cme}}^{{\rm Ru}}}\end{split}$$ (21)

where $f_{\textsc{cme}}\equiv\Delta\gamma_{{\textsc{cme}}}/\Delta\gamma=1-\Delta\gamma_{\text{bkgd}}/\Delta\gamma$ is the fraction of CME signal in the $\Delta\gamma$ measurement, and $\delta f_{\textsc{cme}}\equiv f_{\textsc{cme}}^{{\rm Ru}}-f_{\textsc{cme}}^{{\rm Zr}}$ is its difference between Ru+Ru and Zr+Zr collisions. Assuming the CME signal is proportional to the squared magnetic field, $\Delta\gamma_{{\textsc{cme}}}\propto B^{2}$, then after a few steps of algebra we obtain

$$\begin{split}f_{\textsc{cme}}^{{\rm Ru}}=&\frac{Y_{\text{signal}}}{Y}\left/\left[1-Y_{\text{bkgd}}\left/\frac{(B^{2}/v_{2}^{*})^{{\rm Ru}}}{(B^{2}/v_{2}^{*})^{{\rm Zr}}}\right.\right]\right.\\ \approx&\frac{Y_{\text{signal}}}{Y}\left/\left[1-Y_{\text{bkgd}}\left(1+\frac{\delta v_{2}^{*}}{v_{2}^{*}}-\frac{\delta B^{2}}{B^{2}}\right)\right]\right.\end{split}$$ (22)

where $\delta B^{2}/B^{2}\equiv(B^{2}_{{\rm Ru}}-B^{2}_{{\rm Zr}})/B^{2}_{{\rm Zr}}$ is the relative squared magnetic field difference between the isobar collisions. It reduces to the simple relationship $f_{\textsc{cme}}^{{\rm Ru}}=Y_{\text{signal}}\left/\frac{\delta B^{2}}{B^{2}}\right.$ if neglecting small quantities.

Our results for $Y_{\text{signal}}$ and $f_{\textsc{cme}}$ are consistent with zero. Assuming $\delta B^{2}/B^{2}=15\%$ [67, 68] and that $f_{\textsc{cme}}$ has a Gaussian probability distribution with a lower bound at the origin [69] ($f_{\textsc{cme}}\geq 0$), we extract an upper limit on $f_{\textsc{cme}}^{{\rm Ru}}$ of roughly 10% for all four measurements at 95% confidence level. (The upper limit on $f_{\textsc{cme}}^{{\rm Zr}}$ is accordingly smaller.) The upper limits are listed in Table 2 and illustrated in Fig. 14(a). Figure 14(b) shows the $f_{\textsc{cme}}^{{\rm Ru}}$ upper limits extracted for a range of values of $\delta B^{2}/B^{2}$.

We note that in Fig. 12 the difference between data and baseline in the 50-80% centrality range is $0.0305\pm 0.0100\pm 0.0041$ ($2.8\sigma$) for full events and $0.0355\pm 0.0215\pm 0.0046$ ($1.6\sigma$) for subevents. In the blind analysis [25] and in this work, we have concentrated on the mid-central 20-50% centrality range where the CME is predicted to be more probable than peripheral or central collisions [3, 4, 6]. We speculate that the peripheral collision results are likely due to fluctuations.

Figure 14: (Left panel) Upper limits at 95% confidence level on the CME fraction $f_{\textsc{cme}}^{\rm Ru}$ in the inclusive $\Delta\gamma$ measurement in 20–50% centrality Ru+Ru collisions extracted from the isobar blind analyses [25] with the background baseline estimates for four measurements from this work, assuming a squared magnetic field difference of $\delta B^{2}/B^{2}=15\%$ between the two isobar systems [68, 67]. (Right panel) Extracted $f_{\textsc{cme}}^{\rm Ru}$ upper limits as a function of $\delta B^{2}/B^{2}$. The statistical and systematic uncertainties are added in quadrature in extracting the upper limits.

5 Summary

In this study, we have estimated nonflow contributions in $v_{2}$ by fitting two-particle $(\Delta\eta,\Delta\phi)$ distributions and analyzing the deviations from simple multiplicity scaling of the three-particle correlator using the STAR isobar data. The 3-particle nonflow correlation contributions to $C_{3}$ are evaluated using hijing simulations. With these inputs, we have obtained an improved background estimate of the Ru+Ru to Zr+Zr ratio of the $\Delta\gamma/v_{2}$ variable. The estimated background baselines are found to be consistent with the STAR measurements for both the full-event and subevent methods. We have also extracted an upper limit of the CME fraction of approximately $10\%$ with a $95\%$ confidence level in isobar collisions at 200 GeV.

This paper focuses on the STAR isobar experiments. On the other hand, the Au+Au collision data from STAR indicate a possible finite CME signal [16]. This is consistent with the expectation that the signal to background ratio is approximately a factor of three larger in Au+Au collisions than in isobar collisions [66]. To outlook, an order of magnitude increase in statistics is expected from future data taking of Au+Au collisions at 200 GeV [70]. STAR will continue the CME search with those data as well as data from the beam energy scan [71].