# On the Probability of Immunity

Jose M. Peña1 1Linköping University, Sweden<EMAIL_ADDRESS>
###### Abstract.

This work is devoted to the study of the probability of immunity, i.e. the
effect occurs whether exposed or not. We derive necessary and sufficient
conditions for non-immunity and $\epsilon$-bounded immunity, i.e. the
probability of immunity is zero and $\epsilon$-bounded, respectively. The
former allows us to estimate the probability of benefit (i.e., the effect
occurs if and only if exposed) from a randomized controlled trial, and the
latter allows us to produce bounds of the probability of benefit that are
tighter than the existing ones. We also introduce the concept of indirect
immunity (i.e., through a mediator) and repeat our previous analysis for it.
Finally, we propose a method for sensitivity analysis of the probability of
immunity under unmeasured confounding.

##  1\. Introduction

Let $X$ and $Y$ denote an exposure and its outcome, respectively. Let $X$ and
$Y$ be binary taking values in $\\{x,x^{\prime}\\}$ and $\\{y,y^{\prime}\\}$.
Let $Y_{x}$ and $Y_{x^{\prime}}$ denote the counterfactual outcome when the
exposure is set to level $X=x$ and $X=x^{\prime}$. Let $y_{x}$,
$y^{\prime}_{x}$, $y_{x^{\prime}}$ and $y^{\prime}_{x^{\prime}}$ denote the
events $Y_{x}=y$, $Y_{x}=y^{\prime}$, $Y_{x^{\prime}}=y$ and
$Y_{x^{\prime}}=y^{\prime}$. For instance, let $X$ represent whether a patient
gets treated or not for a deadly disease, and $Y$ represent whether she
survives it or not. Individual patients can be classified into immune (they
survive whether they are treated or not, i.e. $y_{x}\land y_{x^{\prime}}$),
doomed (they die whether they are treated or not, i.e. $y^{\prime}_{x}\land
y^{\prime}_{x^{\prime}}$), benefited (they survive if and only if treated,
i.e. $y_{x}\land y^{\prime}_{x^{\prime}}$), and harmed (they die if and only
if treated, i.e. $y^{\prime}_{x}\land y_{x^{\prime}}$).

In general, the average treatment effect (ATE) estimated from a randomized
controlled trial (RCT) does not inform about the probability of benefit (or of
any of the other response types, i.e. harm, immunity, and doom). However, it
may do it under certain conditions. For instance,

$\displaystyle ATE=p(y_{x})-p(y_{x^{\prime}})$
$\displaystyle=p(y_{x},y_{x^{\prime}})+p(y_{x},y^{\prime}_{x^{\prime}})-[p(y_{x},y_{x^{\prime}})+p(y^{\prime}_{x},y_{x^{\prime}})]$
$\displaystyle=p(\text{benefit})-p(\text{harm})$ (1)

and thus $p(\text{benefit})=ATE$ if $p(\text{harm})=0$ (a.k.a. monotonicity).
Necessary and sufficient conditions are derived by Mueller and Pearl [1] to
determine from observational and experimental data if monotonicity holds. In
this work, we derive similar conditions for non-immunity, i.e.
$p(\text{immunity})=p(y_{x},y_{x^{\prime}})=0$. These are interesting because
under non-monotonicity, they turn an RCT informative about the probabilities
of benefit and harm. To see it, consider

$ATE=p(y_{x})-p(y_{x^{\prime}})$

where the terms on the right-hand side of the equation are estimated from an
RCT. Moreover,

$\displaystyle p(y_{x})$
$\displaystyle=p(y_{x},y_{x^{\prime}})+p(y_{x},y^{\prime}_{x^{\prime}})=p(\text{immunity})+p(\text{benefit})$
(2) $\displaystyle p(y_{x^{\prime}})$
$\displaystyle=p(y_{x},y_{x^{\prime}})+p(y^{\prime}_{x},y_{x^{\prime}})=p(\text{immunity})+p(\text{harm})$
(3)

and thus $p(\text{benefit})=p(y_{x})$ and $p(\text{harm})=p(y_{x^{\prime}})$
if $p(\text{immunity})=0$.

In some cases, non-immunity is assured. For instance, when evaluating the
effect of advertising on the purchase of a new product. The control group not
being exposed to the ad has no way of purchasing the product, i.e.
$p(y_{x^{\prime}})=0$ and thus $p(y_{x},y_{x^{\prime}})=0$. In other cases,
non-immunity cannot be assured. For instance, when evaluating the effect of a
drug. An individual may carry a gene variant that makes her recover from the
disease regardless of whether she takes the drug or not, i.e.
$p(y_{x},y_{x^{\prime}})\geq 0$. However, it may still be bounded as
$p(y_{x},y_{x^{\prime}})\leq\epsilon$ from expert knowledge. We show that our
necessary and sufficient conditions for non-immunity can trivially be adapted
to $\epsilon$-bounded immunity. Moreover, we show that the knowledge of
$\epsilon$-bounded immunity may tighten the bounds of the probabilities of
benefit and harm by Tian and Pearl [2]. We also introduce the concepts of
indirect benefit and harm (i.e., through a mediator) and repeat our previous
analysis for them. Finally, we propose a method for sensitivity analysis of
immunity under unmeasured confounding. We illustrate our results with concrete
examples.

##  2\. Conditions for Non-Immunity

Consider the bounds of $p(\text{benefit})$ derived by Tian and Pearl [2]:

$\max\left\\{\begin{array}[]{cc}0,\\\ p(y_{x})-p(y_{x^{\prime}}),\\\
p(y)-p(y_{x^{\prime}}),\\\ p(y_{x})-p(y)\end{array}\right\\}\leq
p(\text{benefit})\leq\min\left\\{\begin{array}[]{cc}p(y_{x}),\\\
p(y^{\prime}_{x^{\prime}}),\\\ p(x,y)+p(x^{\prime},y^{\prime}),\\\
p(y_{x})-p(y_{x^{\prime}})+\\\
p(x,y^{\prime})+p(x^{\prime},y)\end{array}\right\\}.$ (4)

Then, combining Equations 2 or 3 with 4 gives

$\max\left\\{\begin{array}[]{cc}0,\\\ p(y_{x})-p(y^{\prime}_{x^{\prime}}),\\\
p(y_{x})-p(x,y)-\\\ p(x^{\prime},y^{\prime}),\\\
p(y_{x^{\prime}})-p(x,y^{\prime})-\\\ p(x^{\prime},y)\end{array}\right\\}\leq
p(\text{immunity})\leq\min\left\\{\begin{array}[]{cc}p(y_{x}),\\\
p(y_{x^{\prime}}),\\\ p(y_{x})-p(y)+\\\ p(y_{x^{\prime}}),\\\
p(y)\end{array}\right\\}.$ (5)

A sufficient condition for $p(\text{immunity})=0$ to hold is that some
argument to the min function in Equation 5 is equal to 0, that is

$p(y_{x})=0\text{ or }p(y_{x^{\prime}})=0\text{ or
}p(y_{x})+p(y_{x^{\prime}})=p(y)\text{ or }p(y)=0.$ (6)

Likewise, a necessary condition for $p(\text{immunity})=0$ to hold is that all
the arguments to the max function are non-positive, that is

$\displaystyle p(y_{x})+p(y_{x^{\prime}})\leq 1\text{ and }$ $\displaystyle
p(y_{x})\leq p(x,y)+p(x^{\prime},y^{\prime})\text{ and }$ $\displaystyle
p(y_{x^{\prime}})\leq p(x,y^{\prime})+p(x^{\prime},y).$ (7)

###  2.1. Conditions for $\epsilon$-Bounded Immunity

The conditions in the previous section can be relaxed to allow certain degree
of immunity (e.g., based on expert knowledge), making them more applicable in
practice. Specifically, a sufficient condition for
$p(\text{immunity})\leq\epsilon$ to hold is

$p(y_{x})\leq\epsilon\text{ or }p(y_{x^{\prime}})\leq\epsilon\text{ or
}p(y_{x})+p(y_{x^{\prime}})\leq p(y)+\epsilon\text{ or }p(y)\leq\epsilon.$

Likewise, a necessary condition for $p(\text{immunity})\leq\epsilon$ to hold
is

$\displaystyle p(y_{x})+p(y_{x^{\prime}})\leq 1+\epsilon\text{ and }$
$\displaystyle p(y_{x})\leq p(x,y)+p(x^{\prime},y^{\prime})+\epsilon\text{ and
}$ $\displaystyle p(y_{x^{\prime}})\leq
p(x,y^{\prime})+p(x^{\prime},y)+\epsilon.$ (8)

###  2.2. $\epsilon$-Bounds on Benefit and Harm

Assuming $\epsilon$-bounded immunity (e.g., based on expert knowledge) can
help narrowing the bounds on $p(\text{benefit})$ and $p(\text{harm})$.
Specifically, if $p(\text{immunity})\leq\epsilon$ then Equation 2 gives

$p(y_{x})-\epsilon\leq p(\text{benefit})\leq p(y_{x}).$

Incorporating this into Equation 4 gives

$\max\left\\{\begin{array}[]{cc}0,\\\ p(y_{x})-p(y_{x^{\prime}}),\\\
p(y)-p(y_{x^{\prime}}),\\\ p(y_{x})-p(y),\\\
p(y_{x})-\epsilon\end{array}\right\\}\leq
p(\text{benefit})\leq\min\left\\{\begin{array}[]{cc}p(y_{x}),\\\
p(y^{\prime}_{x^{\prime}}),\\\ p(x,y)+p(x^{\prime},y^{\prime}),\\\
p(y_{x})-p(y_{x^{\prime}})+\\\
p(x,y^{\prime})+p(x^{\prime},y)\end{array}\right\\}$ (9)

which can potentially return a tighter lower bound than Equation 4, i.e. if
$\epsilon<\min(p(y_{x^{\prime}}),p(y))$. Although the value of $\epsilon$ is
typically determined from expert knowledge and not from data, the experimental
and observational data available do restrict the values that are valid, as
indicated by Equation 2.1. In short, $\epsilon$ can take any value as long as
the lower bound is not greater than the upper bound in Equation 9. Moreover,
$p(\text{harm})$ can likewise be bounded by simply swapping $x$ and
$x^{\prime}$ in Equation 9.

###  2.3. Examples

This section illustrates the results above with two concrete examples.111R
code for the examples can be found at https://tinyurl.com/2s3bxmyu.

####  2.3.1. Example 1

A pharmaceutical company wants to market their drug to cure a disease by
claiming that no one is immune. The RCT they conducted for the drug approval
yielded the following:

$\displaystyle p(y_{x})$ $\displaystyle=0.76$ $\displaystyle
p(y_{x^{\prime}})$ $\displaystyle=0.31$

which correspond to the following unknown data generation model:

$\displaystyle p(u)=0.3$ $\displaystyle p(x|u)=0.2$ $\displaystyle p(y|x,u)$
$\displaystyle=0.9$ $\displaystyle p(y|x,u^{\prime})$ $\displaystyle=0.7$
$\displaystyle p(x|u^{\prime})=0.9$ $\displaystyle p(y|x^{\prime},u)$
$\displaystyle=0.8$ $\displaystyle p(y|x^{\prime},u^{\prime})$
$\displaystyle=0.1.$

Therefore, the necessary condition for non-immunity in Equation 2 does not
hold, and thus the company is not entitled to make the claim they intended to
make. The company changes strategy and now wishes to market their drug as
having a minimum of 50 % efficacy, i.e. benefit. To do so, they first conduct
an observational study that yields the following:

$\displaystyle p(x,y)=0.5$ $\displaystyle p(x,y^{\prime})$ $\displaystyle=0.2$
$\displaystyle p(x^{\prime},y)=0.2$ $\displaystyle p(x^{\prime},y^{\prime})$
$\displaystyle=0.1.$

Then, they apply Equation 4 to the RCT and observational results to conclude
that $0.45\leq p(\text{benefit})\leq 0.61$. Again, the company cannot proceed
with their marketing strategy. A few months later, a research publication
reports that no more than 25 % of the population is immune. The company
realizes that this value is compatible with their RCT and observational
results, by checking the necessary condition for $\epsilon$-bounded immunity
in Equation 2.1. More importantly, the company realizes that Equation 9 with
$\epsilon=0.25$ allows to conclude that $0.51\leq p(\text{benefit})\leq 0.61$,
and thus they can resume their marketing strategy.

####  2.3.2. Example 2

The previous example has shown that expert knowledge on immunity may
complement experimental and observational data. While data alone rarely
provide precise information on immunity (or on any other response type, for
that matter), there are cases where data alone provide enough actionable
information. The following example illustrates this.

A pharmaceutical company is concerned by the poor sales of a drug to cure a
disease. The RCT conducted for the drug approval and a subsequent
observational study yielded the following:

$\displaystyle p(y_{x})$ $\displaystyle=0.48$ $\displaystyle
p(y_{x^{\prime}})$ $\displaystyle=0.36$

and

$\displaystyle p(x,y)=0.08$ $\displaystyle p(x,y^{\prime})$
$\displaystyle=0.2$ $\displaystyle p(x^{\prime},y)=0.25$ $\displaystyle
p(x^{\prime},y^{\prime})$ $\displaystyle=0.47.$

which correspond to the following unknown data generation model:

$\displaystyle p(u)=0.4$ $\displaystyle p(x|u)=0.1$ $\displaystyle p(y|x,u)$
$\displaystyle=0.9$ $\displaystyle p(y|x,u^{\prime})$ $\displaystyle=0.2$
$\displaystyle p(x|u^{\prime})=0.4$ $\displaystyle p(y|x^{\prime},u)$
$\displaystyle=0.3$ $\displaystyle p(y|x^{\prime},u^{\prime})$
$\displaystyle=0.4.$

The fact that 36 % of the untreated recover from the disease makes the company
suspect that the low sales are due to a large part of the population being
immune. Equation 5 allows to conclude that $0\leq p(\text{immunity})\leq
0.34$, which suggests that the explanation offered by the company is rather
unlikely. A more plausible explanation for the low sales may be that the
efficacy or benefit of the drug is not very high, as $0.14\leq
p(\text{benefit})\leq 0.48$ by Equation 4.

##  3\. Indirect Benefit and Harm

In the previous sections, the causal graph of the domain under study was
unknown. In this section, we assume that the graph is available (e.g., from
expert knowledge) and discuss two advantages that follow with it.
Specifically, suppose that the domain under study corresponds to the following
causal graph:

$X$$Y$$Z$

and thus $p(y_{x})=p(y|x)$ and $p(y_{x^{\prime}})=p(y|x^{\prime})$. Then,
$p(y_{x})$ and $p(y_{x^{\prime}})$ can be estimated from observational data
and thus, unlike in the previous sections, no RCT is required. A further
advantage is that we can now compute the probabilities of benefit and harm
mediated by $Z$. We elaborate on this below.

The effect of $X$ on $Y$ mediated by $Z$ (a.k.a. indirect effect) corresponds
to the effect due to the indirect path $X\rightarrow Z\rightarrow Y$, i.e.
after deactivating the direct path $X\rightarrow Y$. Different ways of
deactivating the direct path have resulted in different indirect effect
measures in the literature. Pearl [3] proposes deactivating the direct path by
setting $X$ to non-exposure and comparing the expected outcome when $Z$ takes
the value it would under exposure and non-exposure:

$NIE=E[Y_{x^{\prime},Z_{x}}]-E[Y_{x^{\prime}}]$

which is known as the average natural (or pure) indirect effect. Geneletti [4]
also proposes deactivating the direct path by setting $X$ to non-exposure but
instead, she proposes comparing the expected outcome when $Z$ is drawn from
the distributions $\mathcal{Z}_{x}$ and $\mathcal{Z}_{x^{\prime}}$ of $Z_{x}$
and $Z_{x^{\prime}}$:

$IIE=E[Y_{x^{\prime},\mathcal{Z}_{x}}]-E[Y_{x^{\prime},\mathcal{Z}_{x^{\prime}}}]$

which is known as the interventional indirect effect. Although $NIE$ and $IIE$
do not coincide in general, they coincide for the causal graph above [5].
Finally, Fulcher et al. [6] proposes deactivating the direct path by setting
$X$ to its natural (observed) value and comparing the expected outcome when
$Z$ takes its natural value and the value it would under no exposure:

$PIIE=E[Y_{X,Z_{X}}]-E[Y_{X,Z_{x^{\prime}}}]$

which is also known as the population intervention indirect effect. This
measure is suitable when the exposure is harmful (e.g., smoking), and thus one
may be more interested in elucidating the effect (e.g., disease prevalence) of
eliminating the exposure rather than in contrasting the effects of exposure
and non-exposure.

We propose an alternative way of deactivating the direct path $X\rightarrow Y$
and measuring the indirect effect of $X$ on $Y$ through $Z$. Specifically, we
assume that the direct path $X\rightarrow Y$ is actually mediated by an
unmeasured random variable $U$ that is left unmodelled. This arguably holds in
most domains. The identity of $U$ is irrelevant. Let $G$ denote the causal
graph below, i.e. the original causal graph refined with the addition of $U$.

$X$$Y$$U$$Z$

Now, deactivating the direct path $X\rightarrow Y$ in the original causal
graph can be achieved by adjusting for $U$ in $G$, i.e.
$\sum_{u}E[Y|x,u]p(u)$. Unfortunately, $U$ is unmeasured. Instead, we propose
the following way of deactivating $X\rightarrow Y$. Let $H$ denote the causal
graph below, i.e. the result of reversing the edge $X\rightarrow U$ in $G$.

$X$$Y$$U$$Z$

The average total effect of $X$ on $Y$ in $H$ can be computed by the front-
door criterion [7]:

$\displaystyle TE$ $\displaystyle=E[Y_{x}]-E[Y_{x^{\prime}}]$
$\displaystyle=\sum_{z}p(z|x)\sum_{\dot{x}}E[Y|\dot{x},z]p(\dot{x})-\sum_{z}p(z|x^{\prime})\sum_{\dot{x}}E[Y|\dot{x},z]p(\dot{x}).$

Note that $G$ and $H$ are distribution equivalent, i.e. every probability
distribution that is representable by $G$ is representable by $H$ and vice
versa [7]. Then, evaluating the second line of the equation above in $G$ or
$H$ gives the same result. If we evaluate it in $H$, then it corresponds to
the part of association between $X$ and $Y$ that is attributable to the path
$X\rightarrow Z\rightarrow Y$. If we evaluate it in $G$, then it corresponds
to the part of $TE$ in $G$ that is attributable to the path $X\rightarrow
Z\rightarrow Y$, because $TE$ in $G$ equals the association between $X$ and
$Y$, since $G$ has only directed paths from $X$ to $Y$. Therefore, the second
line in the equation above corresponds to the part of $TE$ in the original
causal graph that is attributable to the path $X\rightarrow Z\rightarrow Y$,
thereby deactivating the direct path $X\rightarrow Y$. We propose to use the
second line in the equation above as a measure of the indirect effect of $X$
on $Y$ in the original causal graph.

The reasoning above can be extended to the probabilities of benefit and harm,
and thereby measure the benefit and harm mediated by $Z$. As mentioned above,
the causal graphs $G$ and $H$ represent different data generation mechanisms
but the same probability distribution over $X$, $Y$ and $Z$. Therefore, the
mechanisms agree on observational probabilities but may disagree on
counterfactual probabilities. We use $p()$ to denote observational
probabilities obtained from either mechanism, and $q()$ to denote
counterfactual probabilities obtained from the mechanism corresponding to $H$.
The probabilities of benefit and harm of $X$ on $Y$ mediated by $Z$ in $G$ and
thus in the original causal graph (henceforth indirect benefit and harm, or
$IB$ and $IH$) can be computed by applying Equation 2 to $H$. That is,

$IB=q(\text{benefit})=q(y_{x})=\sum_{z}p(z|x)\sum_{\dot{x}}p(y|\dot{x},z)p(\dot{x})$

where the second equality holds if $q(\text{immunity})=0$, and the third is
due to the front-door criterion on $H$. Likewise for $IH$ simply replacing $x$
by $x^{\prime}$. Applying Equation 5 to $H$ yields necessary and sufficient
conditions for $q(\text{immunity})=0$. That is,

$\displaystyle\sum_{z}p(z|x)\sum_{\dot{x}}p(y|\dot{x},z)p(\dot{x})=0\text{ or
}$
$\displaystyle\sum_{z}p(z|x^{\prime})\sum_{\dot{x}}p(y|\dot{x},z)p(\dot{x})=0\text{
or }$
$\displaystyle\sum_{z}[p(z|x)+p(z|x^{\prime})]\sum_{\dot{x}}p(y|\dot{x},z)p(\dot{x})=p(y)\text{
or }$ $\displaystyle p(y)=0$ (10)

is a sufficient condition, whereas

$\displaystyle\sum_{z}[p(z|x)+p(z|x^{\prime})]\sum_{\dot{x}}p(y|\dot{x},z)p(\dot{x})\leq
1\text{ and }$
$\displaystyle\sum_{z}p(z|x)\sum_{\dot{x}}p(y|\dot{x},z)p(\dot{x})\leq
p(x,y)+p(x^{\prime},y^{\prime})\text{ and }$
$\displaystyle\sum_{z}p(z|x^{\prime})\sum_{\dot{x}}p(y|\dot{x},z)p(\dot{x})\leq
p(x,y^{\prime})+p(x^{\prime},y)$ (11)

is a necessary condition. Necessary and sufficient conditions for
$\epsilon$-bounded immunity on $H$ (i.e., $q(\text{immunity})\leq\epsilon$)
can be obtained much like in Section 2.1. That is, it suffices to add
$\epsilon$ to the right-hand sides of the conditions above and replace $=$
with $\leq$. Finally, we can adapt accordingly the equations in Section 2.2 to
obtain $\epsilon$-bounds on $IB$ and $IH$. Note that the analysis of indirect
benefit and harm presented here does not require an RCT, i.e. all the
expressions involved can be estimated from just observational data.

###  3.1. Example

This section illustrates the results above with a concrete example borrowed
from Pearl [8]. It concerns the following causal graph:

$X$$Y$$Z$

where $X$ represents a drug treatment, $Z$ the presence of a certain enzyme in
a patient’s blood, and $Y$ recovery. Moreover, we have that

$\displaystyle p(z|x)=0.75$ $\displaystyle p(y|x,z)=0.8$ $\displaystyle
p(y|x,z^{\prime})=0.4$ $\displaystyle p(z|x^{\prime})=0.4$ $\displaystyle
p(y|x^{\prime},z)=0.3$ $\displaystyle p(y|x^{\prime},z^{\prime})=0.2.$

Since $p(x)$ is not given in the original example, we take $p(x)=0.6$.

Pearl imagines a scenario where the pharmaceutical company plans to develop a
cheaper drug that is equal to the existing one except for the lack of direct
effect on recovery, i.e. it just stimulates enzyme production as much as the
existing drug. Therefore, the probability of benefit of the planned drug is
the probability of benefit of the existing drug that is mediated by the
enzyme. The company wants to market their drugs by claiming that no one is
immune. The sufficient conditions for non-immunity in Equations 6 and 3 do not
hold for the drugs. However, while the existing drug satisfies the necessary
condition for non-immunity in Equation 2, the planned drug does not satisfy
the corresponding condition in Equation 3. Therefore, the company should
either abandon their marketing strategy or abandon the plan to develop the new
drug and instead focus on trying to confirm non-immunity for the existing
drug.

##  4\. Sensitivity Analysis of Immunity

In this section, like in the previous section, we assume that the causal graph
of the domain under study is available, e.g. from expert knowledge. We also
assume that we only have access to observational data, i.e. no RCT is
available. Specifically, consider the following causal graph:

$X$$Y$$Z$

which includes potential unmeasured exposure-outcome confounding. Since
$p(y_{x})=\sum_{z}p(z|x)\sum_{\dot{x}}E[Y|\dot{x},z]p(\dot{x})$ by the front-
door criterion, we can proceed as in the previous section to derive necessary
and sufficient conditions for non-immunity. Suppose now that $Z$ is unmeasured
or that the effect of $X$ on $Y$ is direct rather than mediated by $Z$. Then,
$p(y_{x})$ is unidentifiable from observational data [7], and thus we cannot
proceed as in the previous section. We therefore take an alternative approach
to inform the analyst about the probability of immunity and thereby help her
in decision making. In particular, we propose a sensitivity analysis method to
bound the probability of immunity as a function of the observed data
distribution and some intuitive sensitivity parameters. Our method is an
straightforward adaption of the method by Peña [9], originally developed to
bound the probabilities of benefit and harm.

Let $U$ denote the unmeasured exposure-outcome confounders. For simplicity, we
assume that all these confounders are categorical, but our results also hold
for ordinal and continuous confounders.222If $U$ is continuous then
sums/maxima/minimima over $u$ should be replaced by integrals/suprema/infima.
For simplicity, we treat $U$ as a categorical random variable whose levels are
the Cartesian product of the levels of the elements in the original $U$.

Note that

$p(y_{x})=p(y_{x}|x)p(x)+p(y_{x}|x^{\prime})p(x^{\prime})=p(y|x)p(x)+p(y_{x}|x^{\prime})p(x^{\prime})$

where the second equality follows from counterfactual consistency, i.e.
$X=x\Rightarrow Y_{x}=Y$. Moreover,

$p(y_{x}|x^{\prime})=\sum_{u}p(y_{x}|x^{\prime},u)p(u|x^{\prime})=\sum_{u}p(y|x,u)p(u|x^{\prime})\leq\max_{u}p(y|x,u)$

where the second equality follows from $Y_{x}\\!\perp\\!X|U$ for all $x$, and
counterfactual consistency. Likewise,

$p(y_{x}|x^{\prime})\geq\min_{u}p(y|x,u).$

Now, let us define

$M_{x}=\max_{u}p(y|x,u)$

and

$m_{x}=\min_{u}p(y|x,u)$

and likewise $M_{x^{\prime}}$ and $m_{x^{\prime}}$. Then,

$p(x,y)+p(x^{\prime})m_{x}\leq p(y_{x})\leq p(x,y)+p(x^{\prime})M_{x}$

and likewise

$p(x^{\prime},y)+p(x)m_{x^{\prime}}\leq p(y_{x^{\prime}})\leq
p(x^{\prime},y)+p(x)M_{x^{\prime}}.$

These equations together with Equation 5 give

$\max\left\\{\begin{array}[]{cc}0,\\\
p(x^{\prime})m_{x}+p(x)m_{x^{\prime}}-p(y^{\prime}),\\\
p(x^{\prime})m_{x}-p(x^{\prime},y^{\prime}),\\\
p(x)m_{x^{\prime}}-p(x,y^{\prime})\end{array}\right\\}\leq
p(\text{immunity})\leq\min\left\\{\begin{array}[]{cc}p(x,y)+p(x^{\prime})M_{x},\\\
p(x^{\prime},y)+p(x)M_{x^{\prime}},\\\
p(x^{\prime})M_{x}+p(x)M_{x^{\prime}},\\\ p(y)\end{array}\right\\}$ (12)

where $m_{x}$, $M_{x}$, $m_{x^{\prime}}$ and $M_{x^{\prime}}$ are sensitivity
parameters. The possible regions for $m_{x}$ and $M_{x}$ are

$0\leq m_{x}\leq p(y|x)\leq M_{x}\leq 1$ (13)

and likewise for $m_{x^{\prime}}$ and $M_{x^{\prime}}$.

Our lower bound in Equation 12 is informative if and only if333Note that the
second row in the maximum equals the third plus the fourth rows.

$0<p(x^{\prime})m_{x}-p(x^{\prime},y^{\prime})$

or

$0<p(x)m_{x^{\prime}}-p(x,y^{\prime}).$

Then, the informative regions for $m_{x}$ and $m_{x^{\prime}}$ are

$p(y^{\prime}|x^{\prime})<m_{x}\leq p(y|x)$

and

$p(y^{\prime}|x)\leq m_{x^{\prime}}<p(y|x^{\prime}).$

On the other hand, our upper bound in Equation 12 is informative444Note that
we already know that $p(\text{immunity})\leq p(y)$ by Equation 5. if and only
if555Note that the third row in the minimum equals the first plus the second
minus the fourth rows.

$p(x,y)+p(x^{\prime})M_{x}<p(y)$

or

$p(x^{\prime},y)+p(x)M_{x^{\prime}}<p(y)$

which occurs if and only if $p(y|x)<p(y|x^{\prime})$ or
$p(y|x^{\prime})<p(y|x)$.666To see it, rewrite $p(y)=p(x,y)+p(x^{\prime},y)$
and recall Equation 13. Therefore, our upper bound is always informative, and
thus the informative regions for $M_{x}$ and $M_{x^{\prime}}$ coincide with
their possible regions.

Figure 1. Lower and upper bounds of $p(\text{immunity})$ in the example in
Section 4.1 as functions of the sensitivity parameters $m_{x}$,
$m_{x^{\prime}}$, $M_{x}$ and $M_{x^{\prime}}$.

###  4.1. Example

We illustrate our method for sensitivity analysis of $p(\text{immunity})$ with
the following fictitious epidemiological example. Consider a population
consisting of a majority and a minority group. Let the binary random variable
$U$ represent the group an individual belongs to. Let $X$ represent whether
the individual gets treated or not for a certain disease. Let $Y$ represent
whether the individual survives the disease. Assume that the scientific
community agrees that $U$ is a confounder for $X$ and $Y$. Assume also that it
is illegal to store the values of $U$, to avoid discrimination complaints. In
other words, the identity of the confounder is known but its values are not.
More specifically, consider the following unknown data generation model:

$\displaystyle p(u)=0.2$ $\displaystyle p(x|u)=0.4$ $\displaystyle
p(y|x,u)=0.9$ $\displaystyle p(y|x,u^{\prime})=0.8$ $\displaystyle
p(x|u^{\prime})=0.2$ $\displaystyle p(y|x^{\prime},u)=0.2$ $\displaystyle
p(y|x^{\prime},u^{\prime})=0.7.$

Since this model does not specify the functional forms of the causal
mechanisms, we cannot compute the true $p(\text{immunity})$ [7]. However, we
can bound it by Equation 5 and the fact that $p(y_{x})=\sum_{u}p(y|x,u)p(u)$
[7], which yields $p(\text{immunity})\in[0.42,0.6]$. Note that these bounds
cannot be computed in practice because $U$ is unmeasured.

Figure 1 (top) shows the lower bound of $p(\text{immunity})$ in Equation 12 as
a function of the sensitivity parameters $m_{x}$ and $m_{x^{\prime}}$. The
axes span the possible regions of the parameters. The dashed lines indicate
the informative regions of the parameters. Specifically, the bottom left
quadrant corresponds to the non-informative region, i.e. the lower bound is
zero. In the data generation model considered, $m_{x}=0.8$ and
$m_{x^{\prime}}=0.2$. These values are unknown to the epidemiologist, because
$U$ is unobserved. However, the figure reveals that the epidemiologist only
needs to have some rough idea of these values to confidently conclude that
$p(\text{immunity})$ is lower bounded by 0.2. Figure 1 (bottom) shows our
upper bound of $p(\text{immunity})$ in Equation 12 as a function of the
sensitivity parameters $M_{x}$ and $M_{x^{\prime}}$. Likewise, having some
rough idea of the unknown values $M_{x}=0.9$ and $M_{x^{\prime}}=0.7$ enables
the epidemiologist to confidently conclude that the $p(\text{immunity})$ is
upper bounded by 0.65. Applying Equation 5 with just observational data
produces looser bounds, namely 0 and 0.67. Recall that
$p(\text{immunity})\in[0.42,0.6]$ in truth.

##  5\. Discussion

The analysis in this work can be repeated for $p(\text{doom})$ instead of
$p(\text{immunity})$ by simply swapping $y$ and $y^{\prime}$, and
$p(\text{benefit})$ and $p(\text{harm})$. Additionally, the analysis of
indirect benefit can be repeated for $p(\text{harm})$ instead of
$p(\text{immunity})$ due to Equation 1, and thereby extend the analysis by
Mueller and Pearl [1].

## References

  * Mueller and Pearl [2023] S. Mueller and J. Pearl.  Monotonicity: Detection, Refutation, and Ramification.  _UCLA Cognitive Systems Laboratory, Technical Report (R-529)_ , 2023\.
  * Tian and Pearl [2000] J. Tian and J. Pearl.  Probabilities of Causation: Bounds and Identification.  _Annals of Mathematics and Artificial Intelligence_ , 28:287–313, 2000.
  * Pearl [2001] J. Pearl.  Direct and Indirect Effects.  In _Proceedings of the 17th Conference on Uncertainty in Artificial Intelligence_ , pages 411–420, 2001.
  * Geneletti [2007] S. Geneletti.  Identifying Direct and Indirect Effects in a Non-Counterfactual Framework.  _Journal of the Royal Statistical Society Series B_ , 69:199–215, 2007.
  * VanderWeele et al. [2014] T. J. VanderWeele, S. Vansteelandt, and J. M. Robins.  Effect Decomposition in the Presence of an Exposure-Induced Mediator-Outcome Confounder.  _Epidemiology_ , 25:300–306, 2014.
  * Fulcher et al. [2020] I. R. Fulcher, I. Shpitser, S. Marealle, and E. J. Tchetgen Tchetgen.  Robust Inference on Population Indirect Causal Effects: The Generalized Front Door Criterion.  _Journal of the Royal Statistical Society Series B_ , 82:199–214, 2020.
  * Pearl [2009] J. Pearl.  _Causality: Models, Reasoning, and Inference_.  Cambridge University Press, 2009.
  * Pearl [2012] J. Pearl.  The Causal Mediation Formula - A Guide to the Assessment of Pathways and Mechanisms.  _Prevention Science_ , 13:426–436, 2012.
  * Peña [2023] J. M. Peña.  Bounding the Probabilities of Benefit and Harm Through Sensitivity Parameters and Proxies.  _Journal of Causal Inference_ , 11:20230012, 2023.


# Evidence for the Cusp Effect in $\eta^{\prime}$ Decays into
$\eta\pi^{0}\pi^{0}$

M. Ablikim1, M. N. Achasov11,b, P. Adlarson70, M. Albrecht4, R. Aliberti31, A.
Amoroso69A,69C, M. R. An35, Q. An66,53, X. H. Bai61, Y. Bai52, O. Bakina32, R.
Baldini Ferroli26A, I. Balossino1,27A, Y. Ban42,g, V. Batozskaya1,40, D.
Becker31, K. Begzsuren29, N. Berger31, M. Bertani26A, D. Bettoni27A, F.
Bianchi69A,69C, J. Bloms63, A. Bortone69A,69C, I. Boyko32, R. A. Briere5, A.
Brueggemann63, H. Cai71, X. Cai1,53, A. Calcaterra26A, G. F. Cao1,58, N.
Cao1,58, S. A. Cetin57A, J. F. Chang1,53, W. L. Chang1,58, G. Chelkov32,a, C.
Chen39, Chao Chen50, G. Chen1, H. S. Chen1,58, M. L. Chen1,53, S. J. Chen38,
S. M. Chen56, T. Chen1, X. R. Chen28,58, X. T. Chen1, Y. B. Chen1,53, Z. J.
Chen23,h, W. S. Cheng69C, X. Chu39, G. Cibinetto27A, F. Cossio69C, J. J.
Cui45, H. L. Dai1,53, J. P. Dai73, A. Dbeyssi17, R. E. de Boer4, D.
Dedovich32, Z. Y. Deng1, A. Denig31, I. Denysenko32, M. Destefanis69A,69C, F.
De Mori69A,69C, Y. Ding36, J. Dong1,53, L. Y. Dong1,58, M. Y. Dong1,53,58, X.
Dong71, S. X. Du75, P. Egorov32,a, Y. L. Fan71, J. Fang1,53, S. S. Fang1,58,
W. X. Fang1, Y. Fang1, R. Farinelli27A, L. Fava69B,69C, F. Feldbauer4, G.
Felici26A, C. Q. Feng66,53, J. H. Feng54, K Fischer64, M. Fritsch4, C.
Fritzsch63, C. D. Fu1, H. Gao58, Y. N. Gao42,g, Yang Gao66,53, S.
Garbolino69C, I. Garzia27A,27B, P. T. Ge71, Z. W. Ge38, C. Geng54, E. M.
Gersabeck62, A Gilman64, K. Goetzen12, L. Gong36, W. X. Gong1,53, W. Gradl31,
M. Greco69A,69C, L. M. Gu38, M. H. Gu1,53, Y. T. Gu14, C. Y Guan1,58, A. Q.
Guo28,58, L. B. Guo37, R. P. Guo44, Y. P. Guo10,f, A. Guskov32,a, T. T. Han45,
W. Y. Han35, X. Q. Hao18, F. A. Harris60, K. K. He50, K. L. He1,58, F. H.
Heinsius4, C. H. Heinz31, Y. K. Heng1,53,58, C. Herold55, M. Himmelreich12,d,
G. Y. Hou1,58, Y. R. Hou58, Z. L. Hou1, H. M. Hu1,58, J. F. Hu51,i, T.
Hu1,53,58, Y. Hu1, G. S. Huang66,53, K. X. Huang54, L. Q. Huang67, L. Q.
Huang28,58, X. T. Huang45, Y. P. Huang1, Z. Huang42,g, T. Hussain68, N
Hüsken25,31, W. Imoehl25, M. Irshad66,53, J. Jackson25, S. Jaeger4, S.
Janchiv29, Q. Ji1, Q. P. Ji18, X. B. Ji1,58, X. L. Ji1,53, Y. Y. Ji45, Z. K.
Jia66,53, H. B. Jiang45, S. S. Jiang35, X. S. Jiang1,53,58, Y. Jiang58, J. B.
Jiao45, Z. Jiao21, S. Jin38, Y. Jin61, M. Q. Jing1,58, T. Johansson70, N.
Kalantar-Nayestanaki59, X. L. Kang76, X. S. Kang36, R. Kappert59, M.
Kavatsyuk59, B. C. Ke75, I. K. Keshk4, A. Khoukaz63, P. Kiese31, R. Kiuchi1,
R. Kliemt12, L. Koch33, O. B. Kolcu57A, B. Kopf4, M. Kuemmel4, M. Kuessner4,
A. Kupsc40,70, W. Kühn33, J. J. Lane62, J. S. Lange33, P. Larin17, A.
Lavania24, L. Lavezzi69A,69C, Z. H. Lei66,53, H. Leithoff31, M. Lellmann31, T.
Lenz31, C. Li43, C. Li39, C. H. Li35, Cheng Li66,53, D. M. Li75, F. Li1,53, G.
Li1, H. Li47, H. Li66,53, H. B. Li1,58, H. J. Li18, H. N. Li51,i, J. Q. Li4,
J. S. Li54, J. W. Li45, Ke Li1, L. J Li1, L. K. Li1, Lei Li3, M. H. Li39, P.
R. Li34,j,k, S. X. Li10, S. Y. Li56, T. Li45, W. D. Li1,58, W. G. Li1, X. H.
Li66,53, X. L. Li45, Xiaoyu Li1,58, H. Liang66,53, H. Liang1,58, H. Liang30,
Y. F. Liang49, Y. T. Liang28,58, G. R. Liao13, L. Z. Liao45, J. Libby24, A.
Limphirat55, C. X. Lin54, D. X. Lin28,58, T. Lin1, B. J. Liu1, C. X. Liu1, D.
Liu17,66, F. H. Liu48, Fang Liu1, Feng Liu6, G. M. Liu51,i, H. Liu34,j,k, H.
B. Liu14, H. M. Liu1,58, Huanhuan Liu1, Huihui Liu19, J. B. Liu66,53, J. L.
Liu67, J. Y. Liu1,58, K. Liu1, K. Y. Liu36, Ke Liu20, L. Liu66,53, M. H.
Liu10,f, P. L. Liu1, Q. Liu58, S. B. Liu66,53, T. Liu10,f, W. K. Liu39, W. M.
Liu66,53, X. Liu34,j,k, Y. Liu34,j,k, Y. B. Liu39, Z. A. Liu1,53,58, Z. Q.
Liu45, X. C. Lou1,53,58, F. X. Lu54, H. J. Lu21, J. G. Lu1,53, X. L. Lu1, Y.
Lu7, Y. P. Lu1,53, Z. H. Lu1, C. L. Luo37, M. X. Luo74, T. Luo10,f, X. L.
Luo1,53, X. R. Lyu58, Y. F. Lyu39, F. C. Ma36, H. L. Ma1, L. L. Ma45, M. M.
Ma1,58, Q. M. Ma1, R. Q. Ma1,58, R. T. Ma58, X. Y. Ma1,53, Y. Ma42,g, F. E.
Maas17, M. Maggiora69A,69C, S. Maldaner4, S. Malde64, Q. A. Malik68, A.
Mangoni26B, Y. J. Mao42,g,g, Z. P. Mao1, S. Marcello69A,69C, Z. X. Meng61, J.
G. Messchendorp59,12, G. Mezzadri1,27A, H. Miao1, T. J. Min38, R. E.
Mitchell25, X. H. Mo1,53,58, N. Yu. Muchnoi11,b, Y. Nefedov32, F. Nerling17,
I. B. Nikolaev11, Z. Ning1,53, S. Nisar9,l, Y. Niu45, S. L. Olsen58, Q.
Ouyang1,53,58, S. Pacetti26B,26C, X. Pan10,f, Y. Pan52, A. Pathak1, A.
Pathak30, M. Pelizaeus4, H. P. Peng66,53, K. Peters12,d, J. Pettersson70, J.
L. Ping37, R. G. Ping1,58, S. Plura31, S. Pogodin32, V. Prasad66,53, F. Z.
Qi1, H. Qi66,53, H. R. Qi56, M. Qi38, T. Y. Qi10,f, S. Qian1,53, W. B. Qian58,
Z. Qian54, C. F. Qiao58, J. J. Qin67, L. Q. Qin13, X. P. Qin10,f, X. S. Qin45,
Z. H. Qin1,53, J. F. Qiu1, S. Q. Qu39, S. Q. Qu56, K. H. Rashid68, C. F.
Redmer31, K. J. Ren35, A. Rivetti69C, V. Rodin59, M. Rolo69C, G. Rong1,58, Ch.
Rosner17, S. N. Ruan39, H. S. Sang66, A. Sarantsev32,c, Y. Schelhaas31, C.
Schnier4, K. Schönning70, M. Scodeggio27A,27B, K. Y. Shan10,f, W. Shan22, X.
Y. Shan66,53, J. F. Shangguan50, L. G. Shao1,58, M. Shao66,53, C. P. Shen10,f,
H. F. Shen1,58, X. Y. Shen1,58, B.-A. Shi58, H. C. Shi66,53, J. Y. Shi1, q. q.
Shi50, R. S. Shi1,58, X. Shi1,53, X. D Shi66,53, J. J. Song18, W. M. Song1,30,
Y. X. Song42,g, S. Sosio69A,69C, S. Spataro69A,69C, F. Stieler31, K. X. Su71,
P. P. Su50, Y.-J. Su58, G. X. Sun1, H. Sun58, H. K. Sun1, J. F. Sun18, L.
Sun71, S. S. Sun1,58, T. Sun1,58, W. Y. Sun30, X Sun23,h, Y. J. Sun66,53, Y.
Z. Sun1, Z. T. Sun45, Y. H. Tan71, Y. X. Tan66,53, C. J. Tang49, G. Y. Tang1,
J. Tang54, L. Y Tao67, Q. T. Tao23,h, M. Tat64, J. X. Teng66,53, V. Thoren70,
W. H. Tian47, Y. Tian28,58, I. Uman57B, B. Wang1, B. L. Wang58, C. W. Wang38,
D. Y. Wang42,g, F. Wang67, H. J. Wang34,j,k, H. P. Wang1,58, K. Wang1,53, L.
L. Wang1, M. Wang45, M. Z. Wang42,g, Meng Wang1,58, S. Wang10,f, T. Wang10,f,
T. J. Wang39, W. Wang54, W. H. Wang71, W. P. Wang66,53, X. Wang42,g, X. F.
Wang34,j,k, X. L. Wang10,f, Y. D. Wang41, Y. F. Wang1,53,58, Y. H. Wang43, Y.
Q. Wang1, Y. Q. Wang16, Y. Wang56, Z. Wang1,53, Z. Y. Wang1,58, Ziyi Wang58,
D. H. Wei13, F. Weidner63, S. P. Wen1, D. J. White62, U. Wiedner4, G.
Wilkinson64, M. Wolke70, L. Wollenberg4, J. F. Wu1,58, L. H. Wu1, L. J.
Wu1,58, X. Wu10,f, X. H. Wu30, Y. Wu66, Z. Wu1,53, L. Xia66,53, T. Xiang42,g,
D. Xiao34,j,k, G. Y. Xiao38, H. Xiao10,f, S. Y. Xiao1, Y. L. Xiao10,f, Z. J.
Xiao37, C. Xie38, X. H. Xie42,g, Y. Xie45, Y. G. Xie1,53, Y. H. Xie6, Z. P.
Xie66,53, T. Y. Xing1,58, C. F. Xu1, C. J. Xu54, G. F. Xu1, H. Y. Xu61, Q. J.
Xu15, S. Y. Xu65, X. P. Xu50, Y. C. Xu58, Z. P. Xu38, F. Yan10,f, L. Yan10,f,
W. B. Yan66,53, W. C. Yan75, H. J. Yang46,e, H. L. Yang30, H. X. Yang1, L.
Yang47, S. L. Yang58, Tao Yang1, Y. X. Yang1,58, Yifan Yang1,58, M. Ye1,53, M.
H. Ye8, J. H. Yin1, Z. Y. You54, B. X. Yu1,53,58, C. X. Yu39, G. Yu1,58, T.
Yu67, C. Z. Yuan1,58, L. Yuan2, S. C. Yuan1, X. Q. Yuan1, Y. Yuan1,58, Z. Y.
Yuan54, C. X. Yue35, A. A. Zafar68, F. R. Zeng45, X. Zeng6, Y. Zeng23,h, Y. H.
Zhan54, A. Q. Zhang1, B. L. Zhang1, B. X. Zhang1, D. H. Zhang39, G. Y.
Zhang18, H. Zhang66, H. H. Zhang54, H. H. Zhang30, H. Y. Zhang1,53, J. L.
Zhang72, J. Q. Zhang37, J. W. Zhang1,53,58, J. X. Zhang34,j,k, J. Y. Zhang1,
J. Z. Zhang1,58, Jianyu Zhang1,58, Jiawei Zhang1,58, L. M. Zhang56, L. Q.
Zhang54, Lei Zhang38, P. Zhang1, Q. Y. Zhang35,75, Shulei Zhang23,h, X. D.
Zhang41, X. M. Zhang1, X. Y. Zhang45, X. Y. Zhang50, Y. Zhang64, Y. T.
Zhang75, Y. H. Zhang1,53, Yan Zhang66,53, Yao Zhang1, Z. H. Zhang1, Z. Y.
Zhang71, Z. Y. Zhang39, G. Zhao1, J. Zhao35, J. Y. Zhao1,58, J. Z. Zhao1,53,
Lei Zhao66,53, Ling Zhao1, M. G. Zhao39, Q. Zhao1, S. J. Zhao75, Y. B.
Zhao1,53, Y. X. Zhao28,58, Z. G. Zhao66,53, A. Zhemchugov32,a, B. Zheng67, J.
P. Zheng1,53, Y. H. Zheng58, B. Zhong37, C. Zhong67, X. Zhong54, H. Zhou45, L.
P. Zhou1,58, X. Zhou71, X. K. Zhou58, X. R. Zhou66,53, X. Y. Zhou35, Y. Z.
Zhou10,f, J. Zhu39, K. Zhu1, K. J. Zhu1,53,58, L. X. Zhu58, S. H. Zhu65, S. Q.
Zhu38, T. J. Zhu72, W. J. Zhu10,f, Y. C. Zhu66,53, Z. A. Zhu1,58, B. S. Zou1,
J. H. Zou1 (BESIII Collaboration) 1 Institute of High Energy Physics, Beijing
100049, People’s Republic of China
2 Beihang University, Beijing 100191, People’s Republic of China
3 Beijing Institute of Petrochemical Technology, Beijing 102617, People’s
Republic of China
4 Bochum Ruhr-University, D-44780 Bochum, Germany
5 Carnegie Mellon University, Pittsburgh, Pennsylvania 15213, USA
6 Central China Normal University, Wuhan 430079, People’s Republic of China
7 Central South University, Changsha 410083, People’s Republic of China
8 China Center of Advanced Science and Technology, Beijing 100190, People’s
Republic of China
9 COMSATS University Islamabad, Lahore Campus, Defence Road, Off Raiwind Road,
54000 Lahore, Pakistan
10 Fudan University, Shanghai 200433, People’s Republic of China
11 G.I. Budker Institute of Nuclear Physics SB RAS (BINP), Novosibirsk 630090,
Russia
12 GSI Helmholtzcentre for Heavy Ion Research GmbH, D-64291 Darmstadt, Germany
13 Guangxi Normal University, Guilin 541004, People’s Republic of China
14 Guangxi University, Nanning 530004, People’s Republic of China
15 Hangzhou Normal University, Hangzhou 310036, People’s Republic of China
16 Hebei University, Baoding 071002, People’s Republic of China
17 Helmholtz Institute Mainz, Staudinger Weg 18, D-55099 Mainz, Germany
18 Henan Normal University, Xinxiang 453007, People’s Republic of China
19 Henan University of Science and Technology, Luoyang 471003, People’s
Republic of China
20 Henan University of Technology, Zhengzhou 450001, People’s Republic of
China
21 Huangshan College, Huangshan 245000, People’s Republic of China
22 Hunan Normal University, Changsha 410081, People’s Republic of China
23 Hunan University, Changsha 410082, People’s Republic of China
24 Indian Institute of Technology Madras, Chennai 600036, India
25 Indiana University, Bloomington, Indiana 47405, USA
26 INFN Laboratori Nazionali di Frascati, (A)INFN Laboratori Nazionali di
Frascati, I-00044, Frascati, Italy; (B)INFN Sezione di Perugia, I-06100,
Perugia, Italy; (C)University of Perugia, I-06100, Perugia, Italy
27 INFN Sezione di Ferrara, (A)INFN Sezione di Ferrara, I-44122, Ferrara,
Italy; (B)University of Ferrara, I-44122, Ferrara, Italy
28 Institute of Modern Physics, Lanzhou 730000, People’s Republic of China
29 Institute of Physics and Technology, Peace Avenue 54B, Ulaanbaatar 13330,
Mongolia
30 Jilin University, Changchun 130012, People’s Republic of China
31 Johannes Gutenberg University of Mainz, Johann-Joachim-Becher-Weg 45,
D-55099 Mainz, Germany
32 Joint Institute for Nuclear Research, 141980 Dubna, Moscow Region, Russia
33 Justus-Liebig-Universitaet Giessen, II. Physikalisches Institut, Heinrich-
Buff-Ring 16, D-35392 Giessen, Germany
34 Lanzhou University, Lanzhou 730000, People’s Republic of China
35 Liaoning Normal University, Dalian 116029, People’s Republic of China
36 Liaoning University, Shenyang 110036, People’s Republic of China
37 Nanjing Normal University, Nanjing 210023, People’s Republic of China
38 Nanjing University, Nanjing 210093, People’s Republic of China
39 Nankai University, Tianjin 300071, People’s Republic of China
40 National Centre for Nuclear Research, Warsaw 02-093, Poland
41 North China Electric Power University, Beijing 102206, People’s Republic of
China
42 Peking University, Beijing 100871, People’s Republic of China
43 Qufu Normal University, Qufu 273165, People’s Republic of China
44 Shandong Normal University, Jinan 250014, People’s Republic of China
45 Shandong University, Jinan 250100, People’s Republic of China
46 Shanghai Jiao Tong University, Shanghai 200240, People’s Republic of China
47 Shanxi Normal University, Linfen 041004, People’s Republic of China
48 Shanxi University, Taiyuan 030006, People’s Republic of China
49 Sichuan University, Chengdu 610064, People’s Republic of China
50 Soochow University, Suzhou 215006, People’s Republic of China
51 South China Normal University, Guangzhou 510006, People’s Republic of China
52 Southeast University, Nanjing 211100, People’s Republic of China
53 State Key Laboratory of Particle Detection and Electronics, Beijing 100049,
Hefei 230026, People’s Republic of China
54 Sun Yat-Sen University, Guangzhou 510275, People’s Republic of China
55 Suranaree University of Technology, University Avenue 111, Nakhon
Ratchasima 30000, Thailand
56 Tsinghua University, Beijing 100084, People’s Republic of China
57 Turkish Accelerator Center Particle Factory Group, (A)Istinye University,
34010, Istanbul, Turkey; (B)Near East University, Nicosia, North Cyprus,
Mersin 10, Turkey
58 University of Chinese Academy of Sciences, Beijing 100049, People’s
Republic of China
59 University of Groningen, NL-9747 AA Groningen, Netherlands
60 University of Hawaii, Honolulu, Hawaii 96822, USA
61 University of Jinan, Jinan 250022, People’s Republic of China
62 University of Manchester, Oxford Road, Manchester, M13 9PL, United Kingdom
63 University of Muenster, Wilhelm-Klemm-Strasse 9, 48149 Muenster, Germany
64 University of Oxford, Keble Road, Oxford, OX13RH, United Kingdom
65 University of Science and Technology Liaoning, Anshan 114051, People’s
Republic of China
66 University of Science and Technology of China, Hefei 230026, People’s
Republic of China
67 University of South China, Hengyang 421001, People’s Republic of China
68 University of the Punjab, Lahore-54590, Pakistan
69 University of Turin and INFN, (A)University of Turin, I-10125, Turin,
Italy; (B)University of Eastern Piedmont, I-15121, Alessandria, Italy;
(C)INFN, I-10125, Turin, Italy
70 Uppsala University, Box 516, SE-75120 Uppsala, Sweden
71 Wuhan University, Wuhan 430072, People’s Republic of China
72 Xinyang Normal University, Xinyang 464000, People’s Republic of China
73 Yunnan University, Kunming 650500, People’s Republic of China
74 Zhejiang University, Hangzhou 310027, People’s Republic of China
75 Zhengzhou University, Zhengzhou 450001, People’s Republic of China
76 China University of Geosciences, Wuhan 430074, People’s Republic of China
a Also at the Moscow Institute of Physics and Technology, Moscow 141700,
Russia
b Also at the Novosibirsk State University, Novosibirsk, 630090, Russia
c Also at the NRC ”Kurchatov Institute”, PNPI, 188300, Gatchina, Russia
d Also at Goethe University Frankfurt, 60323 Frankfurt am Main, Germany
e Also at Key Laboratory for Particle Physics, Astrophysics and Cosmology,
Ministry of Education; Shanghai Key Laboratory for Particle Physics and
Cosmology; Institute of Nuclear and Particle Physics, Shanghai 200240,
People’s Republic of China
f Also at Key Laboratory of Nuclear Physics and Ion-beam Application (MOE) and
Institute of Modern Physics, Fudan University, Shanghai 200443, People’s
Republic of China
g Also at State Key Laboratory of Nuclear Physics and Technology, Peking
University, Beijing 100871, People’s Republic of China
h Also at School of Physics and Electronics, Hunan University, Changsha
410082, China
i Also at Guangdong Provincial Key Laboratory of Nuclear Science, Institute of
Quantum Matter, South China Normal University, Guangzhou 510006, China
j Also at Frontiers Science Center for Rare Isotopes, Lanzhou University,
Lanzhou 730000, People’s Republic of China
k Also at Lanzhou Center for Theoretical Physics, Lanzhou University, Lanzhou
730000, People’s Republic of China
l Also at the Department of Mathematical Sciences, IBA, Karachi , Pakistan

(February 26, 2023)

###### Abstract

Using a sample of $4.3\times 10^{5}$ $\eta^{\prime}\to\eta\pi^{0}\pi^{0}$
events selected from the ten billion $J/\psi$ event dataset collected with the
BESIII detector, we study the decay $\eta^{\prime}\to\eta\pi^{0}\pi^{0}$
within the framework of nonrelativistic effective field theory. Evidence for a
structure at $\pi^{+}\pi^{-}$ mass threshold is observed in the invariant mass
spectrum of $\pi^{0}\pi^{0}$ with a statistical significance of around
$3.5\sigma$, which is consistent with the cusp effect as predicted by the
nonrelativistic effective field theory. After introducing the amplitude for
describing the cusp effect, the $\pi\pi$ scattering length combination
$a_{0}-a_{2}$ is determined to be $\rm 0.226\pm 0.060_{stat}\pm 0.013_{syst}$,
which is in good agreement with theoretical calculation of $0.2644\pm 0.0051$.

Experimental studies of light meson decays are important guides to our
understanding of how QCD works in the nonperturbative regime. In this contest,
the $\pi\pi$ and $\pi K$ interactions at low energies have been the subject of
investigations for a few decades. In $\pi\pi$ interaction, one of the
prominent features is the loop contribution to the $\pi\pi$ scattering: the
$S$-wave charge-exchange rescattering $\pi^{+}\pi^{-}\to\pi^{0}\pi^{0}$ (as
shown in Fig. 1) causes a prominent cusp at the center of mass energy
corresponding to the summed mass of two charged pions. The cusp effect can
shed light on the fundamental properties of QCD at low energies, by
determining the strength of the $S$-wave $\pi\pi$ interaction [1, 2, 3, 4, 5,
6]. Six decades ago this effect was predicted to be seen in
$K^{+}\to\pi^{0}\pi^{0}\pi^{+}$ [7], and it was finally observed in 2006 [8]
by the NA48/2 experiment and studied further [9]. These results inspired
theoretical predictions for the cusp in other decays, such as $K_{L}\to
3\pi^{0}$ [3, 5, 10] and $\eta\to 3\pi^{0}$ [11, 12], which were
experimentally investigated: it was observed in the decay of $K_{L}\to
3\pi^{0}$ by KTeV [13], while no clear evidence was seen in $\eta\to 3\pi^{0}$
decay [14, 15, 16, 17].

Figure 1:  One-loop contribution in $P\to P^{\prime}\pi^{0}\pi^{0}$ decay,
where $P$ and $P^{\prime}$ denote pseudoscalar particles in initial and final
states, respectively. Different behaviors below and above the charged pion
mass threshold cause the cusp effect.

Another process where the cusp effect is expected to have a sizable
contribution is the hadronic decay $\eta^{\prime}\to\eta\pi^{0}\pi^{0}$ [18];
this has been experimentally investigated by BESIII [19], with 5.6$\times$104
$\eta^{\prime}\to\eta\pi^{0}\pi^{0}$ events, and no evidence was seen, while
the A2 experiment [20] accumulated about 1.24$\times$105
$\eta^{\prime}\to\eta\pi^{0}\pi^{0}$ decays and reported a deviation with a
significance of about $2.5\sigma$ and it is also studied in chiral
perturbation theory [21]. Therefore, it is essential to further investigate
this decay with higher precision.

The recently available data of ten billion $J/\psi$ events [22] at BESIII
imply an increased data sample of $\eta^{\prime}$ decays by nearly an order of
magnitude, offering a unique opportunity for further investigations of the
cusp effect. In this Letter, we present the first evidence of the cusp effect
in $\eta^{\prime}\to\eta\pi^{0}\pi^{0}$ and the corresponding measurement of
the $\pi\pi$ scattering length based on the nonrelativistic effective field
theory (NREFT) [18].

The BESIII detector [23, 24] records symmetric $e^{+}e^{-}$ collisions
provided by the BEPCII storage ring [25]. The cylindrical core of the BESIII
detector covers 93% of the full solid angle and consists of a helium-based
multilayer drift chamber (MDC), a plastic scintillator time-of-flight (TOF)
system, and a CsI(Tl) electromagnetic calorimeter (EMC), which are all
enclosed in a superconducting solenoidal magnet providing a 1.0 T (0.9 T in
2012) magnetic field. The end cap TOF system was upgraded in 2015 using
multigap resistive plate chamber technology [26, 27, 28].

To reconstruct events of $J/\psi\to\gamma\eta^{\prime}$ with
$\eta^{\prime}\to\eta\pi^{0}\pi^{0}$, the $\pi^{0}$ and $\eta$ are selected by
$\pi^{0}/\eta\to 2\gamma$ process. The charged tracks are reconstructed from
hits in the MDC. The polar angle with respect to the MDC symmetry axis should
be in the range $|\cos\theta|<0.93$. The distance away from the interaction
point should be less than 10.0 cm in the beam direction and 1.0 cm in the
radial direction. The photon candidates are reconstructed using clusters of
energy deposited in the EMC. The energy deposited in the nearby TOF system is
included in EMC measurements to improve the reconstruction efficiency and the
energy resolution. Photon candidates are required to have a deposited energy
larger than 25 MeV in the barrel region ($|\cos\theta|<0.80$) and 50 MeV in
the end cap regions ($0.86<|\cos\theta|<0.92$). A requirement on the EMC
cluster timing with respect to the most energetic photon, $-500<T<500$ ns, is
used to suppress electronic noise and energy deposits unrelated to the event.
The events with at least seven photon candidates and no charged tracks are
kept for further analysis.

For each candidate event, the photon with the maximum energy is assumed to be
the radiative photon originating from the decay of $J/\psi$, while the
remaining photons are used to reconstruct $\pi^{0}/\eta$ candidates. A one-
constraint (1C) kinematic fit is performed by constraining the invariant mass
of photon pairs to the $\pi^{0}$ or $\eta$ mass, and the $\chi^{2}$ for this
fit is required to be less than 25. Since the $\pi^{0}$ decays into two
photons isotropically in its rest frame, the angle of one photon in the
$\pi^{0}$ rest frame with respect to the $\pi^{0}$ momentum direction is
required to satisfy $|\cos\theta_{\pi^{0}}|<0.95$. Afterward, an eight-
constraint (8C) kinematic fit is performed for the $\gamma\eta\pi^{0}\pi^{0}$
combinations, requiring energy-momentum conservation and constraining the
invariant masses of the three photon pairs to the nominal $\pi^{0}/\eta$
masses and of the $\eta\pi^{0}\pi^{0}$ combinations to the $\eta^{\prime}$
mass. If more than one $\gamma\eta\pi^{0}\pi^{0}$ combination is found, only
the one with the least $\rm\chi^{2}_{8C}$ is retained. After the requirement
of $\rm\chi^{2}_{8C}<100$, 432295 candidate events are accepted for further
analysis; the corresponding Dalitz plot is shown in Fig. 2.

Figure 2:  Dalitz plot of $\eta^{\prime}\to\eta\pi^{0}\pi^{0}$.

To investigate the background contamination, a 6C kinematic fit, instead of
the 8C fit, is performed on candidate events, in which the constraints on the
masses of $\eta$ and $\eta^{\prime}$ are removed. Figure 3 shows the
$\eta\pi^{0}\pi^{0}$ invariant mass distribution of the data sample, after
requiring $\rm\chi^{2}_{6C}<100$ and $\eta$ mass window cut
$|M(\gamma\gamma)-M_{\eta}|<30$ MeV/$c^{2}$ on the unconstrained photon pair,
a clear $\eta^{\prime}$ peak is observed. In addition, a ten billion $J/\psi$
inclusive decay Monte Carlo (MC) sample generated with LUNDCHARM [29, 30] is
used to check possible background sources; the surviving events mainly consist
of the peaking background $\eta^{\prime}\to 3\pi^{0}$ decay channel and the
flat contribution from $J/\psi\to\omega\eta$, with $\omega\to\gamma\pi^{0}$
and $\eta\to 3\pi^{0}$. The background contamination rate is estimated about
0.82%, and its shape on $\pi^{0}\pi^{0}$ and $\eta\pi^{0}$ mass spectrum is
smooth; therefore, it is neglected in the further analysis.

Figure 3:  The $\eta\pi^{0}\pi^{0}$ invariant mass distribution after 6C
kinematic fit without the $\eta$ and $\eta^{\prime}$ mass constraints. The
dots with error bars are experimental data, and the blue histogram is the
$\eta^{\prime}\to 3\pi^{0}$ peaking background from MC sample.

Using an unbinned maximum likelihood method, we fit the Dalitz plot of
$M^{2}(\pi^{0}\pi^{0})$ versus $M^{2}(\eta\pi^{0})$ within the framework of
NREFT. (More details are given in the Supplemental Material [31].) The
resolution effect and detection efficiency are studied by MC simulation and
taken into account in the fit.

\begin{overpic}[width=214.63977pt]{General_X.eps} \put(80.0,85.0){(a)}
\end{overpic}\begin{overpic}[width=214.63977pt]{General_Y.eps}
\put(80.0,85.0){(b)}
\end{overpic}\begin{overpic}[width=214.63977pt]{cusp_tree_s1.eps}
\put(80.0,85.0){(c)}
\end{overpic}\begin{overpic}[width=214.63977pt]{cusp_tree_s3.eps}
\put(80.0,85.0){(d)} \end{overpic}

Figure 4:  The fit result of fit I. The projections to $X$ and $Y$ are shown in (a) and (b), and the mass spectrums of $M^{2}(\eta\pi^{0})$ and $M^{2}(\pi^{0}\pi^{0})$ divided by phase space are shown in (c) and (d). The black dots with error bars are from data and the blue line is the fit result of the tree level amplitude. The red dashed line indicates the charged pion mass threshold.  Table 1: Experimental values of the matrix element parameters for $\eta^{\prime}\to\eta\pi^{0}\pi^{0}$. Parameters | Fit I | Fit II | Fit III | Fit IV
---|---|---|---|---
$a$ | $-0.075\pm 0.003\pm 0.001$ | $-0.207\pm 0.013$ | $-0.143\pm 0.010$ | $-0.077\pm 0.003\pm 0.001$
$b$ | $-0.073\pm 0.005\pm 0.001$ | $-0.051\pm 0.014$ | $-0.038\pm 0.006$ | $-0.066\pm 0.006\pm 0.001$
$d$ | $-0.066\pm 0.003\pm 0.001$ | $-0.068\pm 0.004$ | $-0.067\pm 0.003$ | $-0.068\pm 0.004\pm 0.001$
$a_{0}-a_{2}$ | - | $0.174\pm 0.066$ | $0.225\pm 0.062$ | $0.226\pm 0.060\pm 0.013$
$a_{0}$ | - | $0.497\pm 0.094$ | - | -
$a_{2}$ | - | $0.322\pm 0.129$ | - | -
Statistical significance | - | $3.4\sigma$ | $3.7\sigma$ | $3.6\sigma$

In the simplest case (fit I), only the tree level contribution is included in
which the final state interaction effect is ignored. In this case, the
amplitude is the same as the general parametrization used in Ref. [19]. The
projections of the fit result to the Dalitz coordinates $X$ and $Y$ are shown
in Figs. 4(a) and 4(b), and they indicate a good description of data. The
fitted parameter values, shown in Table 1, are consistent with the previous
BESIII measurement and the statistical uncertainties are about one third of
the previous results [19]. In Figs. 4(c) and 4(d) the comparisons between data
and the fit projections of the $\eta\pi^{0}$ and $\pi^{0}\pi^{0}$ invariant
mass distributions divided by the phase space are presented. The discrepancy
between data and fit result below the charged pion mass threshold corresponds
to the cusp effect. Therefore, we perform alternative fits by including the
loop contributions within the framework of NREFT to evaluate this effect, and
fit I is taken as the baseline for the further loop level fits.

At the loop level amplitude, only $\pi\pi$ scattering is considered while
$\eta\pi$ scattering is ignored; the $S$-wave $\pi\pi$ scattering lengths
$a_{0}$ and $a_{2}$ are included in the loop level amplitude by matching
between NREFT amplitude and partial wave decomposition,

$\begin{split}C_{00}&=\frac{16\pi}{3}(a_{0}+2a_{2})(1-\xi),\\\
C_{x}&=\frac{16\pi}{3}(a_{2}-a_{0})(1+\frac{\xi}{3}),\\\
C_{+-}&=\frac{8\pi}{3}(2a_{0}+a_{2})(1+\xi),\\\
\xi&=\frac{M^{2}_{\pi^{\pm}}-M^{2}_{\pi^{0}}}{M^{2}_{\pi^{\pm}}}.\end{split}$
(1)

where $C_{x}$ denotes the coupling coefficient of the cusp term
$\pi^{+}\pi^{-}\to\pi^{0}\pi^{0}$, and $C_{00}$ and $C_{+-}$ are the coupling
coefficients of noncusp terms $\pi^{0}\pi^{0}\to\pi^{0}\pi^{0}$ and
$\pi^{+}\pi^{-}\to\pi^{+}\pi^{-}$, which are defined in Ref. [18].

The distribution of $M^{2}(\pi^{0}\pi^{0})$ is determined by the whole
amplitude and all five parameters $a$, $b$, $d$, $a_{0}$, and $a_{2}$, while
the distribution of $M^{2}(\eta\pi^{0})$ is mainly determined by parameter
$d$, where $a$, $b$ and $d$ are coefficients in tree level amplitude.

To verify the prediction of NREFT and evaluate the scattering length
combination $a_{0}-a_{2}$, we perform many unbinned maximum likelihood fits in
different cases after including the contributions from the amplitudes at one-
and two-loop levels.

In the case when all the parameters are free (fit II), the fit quality is
improved, and we obtain a statistical significance of $3.4\sigma$ compared to
fit I. In Fig. 5, the comparison between the fit and data for the projections
in different variables, as well as the pull distributions, shows that the fit
provides a good description, in particular, for the region below the charged
pion mass threshold. However, the correlation between the four parameters $a$,
$b$, $a_{0}$, and $a_{2}$ is very large, as shown in Eq. (2). This strong
correlation between $a$, $b$, $a_{0}$, and $a_{2}$ may be caused by the loop
level amplitude contribution to noncusp terms. The scattering length
combination is calculated to be $a_{0}-a_{2}=0.174\pm 0.066$,

$\left(\begin{array}[]{c|cccc}&b&d&a_{0}&a_{2}\\\ \hline\cr a&\hskip
7.0pt0.831&\hskip 7.0pt0.189&-0.966&-0.789\\\ b&&\hskip
7.0pt0.348&-0.918&-0.839\\\ d&&&-0.257&-0.210\\\ a_{0}&&&&\hskip 7.0pt0.872\\\
\end{array}\right)$ (2)

To reduce the correlations between parameters, we also made an attempt (fit
III) by fixing $a_{0}+2a_{2}=0.1312$ according to the theoretical values
$a_{0}=0.220\pm 0.005$ and $a_{2}=-0.0444\pm 0.0010$ [18] and setting
$a_{0}-a_{2}$ as a free parameter, since only $C_{x}$ contributes to the cusp
effect. The fit result presented in Fig. 5 shows a good agreement with data,
also in the region below the charged pion mass threshold. The fitted parameter
values are summarized in Table 1, and the corresponding correlations are shown
in Eq. (3); the obtained value of $a_{0}-a_{2}=0.225\pm 0.062$ is in agreement
with the theoretical value $0.2644\pm 0.0051$ [18]. We also test by changing
the value of $a_{0}+2a_{2}$ or fixing both $a_{0}+2a_{2}$ and $2a_{0}+a_{2}$
to theoretical values, the fit results are consistent with the result of fit
III, and $a_{0}-a_{2}$ is not sensitive to fixed value,

$\left(\begin{array}[]{c|ccc}&b&d&a_{0}-a_{2}\\\ \hline\cr
a&-0.560&-0.046&-0.955\\\ b&&\hskip 7.0pt0.249&\hskip 7.0pt0.457\\\
d&&&-0.032\\\ \end{array}\right)$ (3)

Comparing to the tree level amplitude, the loop contributions with $C_{00}$
and $C_{+-}$ are expected to be small. Additionally, we performed an
alternative fit (fit IV) by ignoring noncusp terms with $C_{00}$ and $C_{+-}$
and only introducing the decay amplitude with $C_{x}$ for the description of
the cusp effect. In this case, the fitted values of different parameters,
summarized in Table 1, are in agreement with those of fit I, the correlations
shown in Eq. (4) are reduced and the statistical significance of the cusp
effect is $3.6\sigma$, while the scattering length combination
$a_{0}-a_{2}=0.226\pm 0.060$ is consistent with that in Ref. [18]. In
addition, we found that the log-likelihood value of fit IV is very close to
those of fit II fit III, which implies that the introduction of the loop
contributions with $C_{00}$ and $C_{+-}$ has little impact on the improvement
of the fit quality and the cusp effect, but significantly increases the
correlations between the different parameters. Therefore, in this analysis, it
is reasonable to ignore these loop contributions in fitting data,

$\left(\begin{array}[]{c|ccc}&b&d&a_{0}-a_{2}\\\ \hline\cr
a&-0.363&-0.253&\hskip 7.0pt0.126\\\ b&&\hskip 7.0pt0.257&\hskip 7.0pt0.237\\\
d&&&-0.107\\\ \end{array}\right)$ (4)

\begin{overpic}[width=195.12767pt]{cusp_s1.eps} \put(80.0,65.0){(a)}
\end{overpic}\begin{overpic}[width=195.12767pt]{cusp_s3.eps}
\put(80.0,65.0){(b)} \end{overpic}

Figure 5:  The fit result projections divided by phase space of different
models to variable (a)$M^{2}(\eta\pi^{0})$ and (b)$M^{2}(\pi^{0}\pi^{0})$. The
black dots with error bars are from data. The solid lines are fit results from
the corresponding models. The red dashed line indicates the charged pion mass
threshold. The cusp region is also shown in the inset.

The systematic uncertainties for the Dalitz plots analysis are listed in Table
2. We calculate the total systematic uncertainty by assuming that all the
contributions are independent and adding them in quadrature.

The photon detection efficiency is studied with the control sample of
$J/\psi\to\rho^{0}\pi^{0}$ events. To evaluate the impact from the slight
discrepancy between data and MC simulation, we perform a correction on the
photon detection and the change of the fit results is considered as the
systematic uncertainty.

To estimate the uncertainties from the 1C kinematic fit for $\pi^{0}$ and
$\eta$, we selected as control samples $J/\psi\to\pi^{+}\pi^{-}\pi^{0}$ and
$J/\psi\to\gamma\eta^{\prime}$ with $\eta^{\prime}\to\eta\pi^{+}\pi^{-}$,
without kinematic fit. After taking into account the discrepancy between data
and MC simulation, repeating the fit with the weighted events leads to changes
of the parameter values, which are assigned as the systematic uncertainties.

To check if the photon miscombinations can effect the fitted parameters, we
generate a MC sample based on NREFT amplitude and tag miscombination events by
matching the truth and the reconstructed value of photon momentum. Two fits
are performed to the sample with and without miscombination events, and the
change of the results is taken as the systematic uncertainty.

To evaluate the uncertainty associated with the efficiency parametrization, we
change the Dalitz plot variables to $M^{2}(\eta\pi)$ and $\cos\theta$, where
$\theta$ is the angle between the directions of the two $\pi$s in the rest
frame of $\eta\pi$. We repeat the fit based on the newly defined Dalitz plot
variables, and the change of the resulting parameters with respect to the
nominal results is assigned as the systematic uncertainty.

The uncertainty of the 8C kinematic fit mainly comes from the inconsistency of
the photon resolution between data and MC simulation. We adjust the energy
resolution in the reconstructed photon error matrix to ensure that the MC
simulation provides a good description of data. Afterward, an alternative fit
is performed and the change of the fitted parameters with respect to the
nominal result is taken as the systematic uncertainty.

To estimate the uncertainty from resolution effect, we vary the resolution by
$\pm 10\%$ and perform alternative fits. The maximum change with respect to
the nominal result is taken as the systematic uncertainty.

Table 2: Summary of the systematic uncertainty sources and their corresponding contributions (%). Parametrization | Fit I | Fit IV
---|---|---
Source | $a$ | $b$ | $d$ | $a$ | $b$ | $d$ | $a_{0}-a_{2}$
Photon detection | 0.7 | 0.4 | 1.0 | 0.6 | 0.4 | 0.9 | 1.8
$\eta$ 1C kinematic fit | 0.1 | 0.6 | 0.0 | 0.1 | 0.7 | 0.0 | 0.2
$\pi^{0}$ 1C kinematic fit | 0.1 | 0.2 | 1.0 | 0.1 | 0.2 | 0.9 | 0.3
Photon miscombination | 0.0 | 0.2 | 1.1 | 0.0 | 0.2 | 1.1 | 1.6
Efficiency presentation | 0.7 | 1.0 | 0.4 | 0.7 | 0.9 | 0.4 | 1.9
Kinematic fit | 0.5 | 1.3 | 0.7 | 0.4 | 0.9 | 0.8 | 4.2
Resolution | 0.0 | 0.0 | 0.0 | 0.1 | 0.3 | 0.0 | 2.0
Total | 1.1 | 1.8 | 2.0 | 1.0 | 1.6 | 1.9 | 5.6

In summary, using ten billion $J/\psi$ events collected with the BESIII
detector, we select a $\eta^{\prime}\to\eta\pi^{0}\pi^{0}$ sample 8 times
larger than that previously analyzed by BESIII, and perform a Dalitz plot
analysis within the framework of nonrelativistic effective field theory. The
fit with tree level amplitude shows a discrepancy below the charged pion mass
threshold, which implies the existence of the cusp effect. To describe the
data in this region, the contributions at one- and two-loop level are
introduced in the decay amplitude. We perform alternative analyses by taking
into account the cusp effect and the results are summarized in Table 1. For
each case, the amplitude provides a good description of the structure around
the charged pion mass threshold and the statistical significance is found to
be around $3.5\sigma$. The scattering length combination $a_{0}-a_{2}$ is
measured to be $0.226\pm 0.060\pm 0.013$, which is in good agreement with the
theoretical value of $0.2644\pm 0.0051$ [18] within the uncertainties. The
observation of the evidence of the cusp effect in
$\eta^{\prime}\to\eta\pi^{0}\pi^{0}$ decay demonstrates the excellent
potential to investigate the underlying dynamics of light mesons at the BESIII
experiment. The prospects [32] for the precise measurements are very promising
at the planned Super Tau-Charm Factories [33, 34].

The BESIII Collaboration thanks the staff of BEPCII and the IHEP computing
center for their strong support. This work is supported in part by National
Key R&D Program of China under Contracts No. 2020YFA0406300, No.
2020YFA0406400; National Natural Science Foundation of China (NSFC) under
Contracts No. 11635010, No. 11735014, No. 11835012, No. 11935015, No.
11935016, No. 11935018, No. 11961141012, No. 12005195, No. 12022510, No.
12025502, No. 12035009, No. 12035013, No. 12192260, No. 12192261, No.
12192262, No. 12192263, No. 12192264, No. 12192265, No. 12225509; the Chinese
Academy of Sciences (CAS) Large-Scale Scientific Facility Program; Joint
Large-Scale Scientific Facility Funds of the NSFC and CAS under Contract No.
U1832207; CAS Key Research Program of Frontier Sciences under Contract No.
QYZDJ-SSW-SLH040; 100 Talents Program of CAS; INPAC and Shanghai Key
Laboratory for Particle Physics and Cosmology; ERC under Contract No. 758462;
European Union’s Horizon 2020 research and innovation program under Marie
Sklodowska-Curie Grant Agreement under Contract No. 894790; German Research
Foundation DFG under Contracts No. 443159800, Collaborative Research Center
CRC 1044, GRK 2149; Istituto Nazionale di Fisica Nucleare, Italy; Ministry of
Development of Turkey under Contract No. DPT2006K-120470; National Science and
Technology fund; STFC (United Kingdom); The Royal Society, UK under Contracts
No. DH140054, No. DH160214; The Swedish Research Council; U. S. Department of
Energy under Award No. DE-FG02-05ER41374.

Supplemental Material:A Brief Description of NREFT amplitude of
$\eta^{\prime}\to\eta\pi^{0}\pi^{0}$ Decay

This supplemental material is based on Ref. [18]. In the
$\eta^{\prime}\to\eta\pi^{0}\pi^{0}$ decay

$\eta^{\prime}(P_{\eta^{\prime}})\to\pi^{0}(p_{1})\pi^{0}(p_{2})\eta(p_{3}),$
(5)

the kinematical variables $s_{i}$ are defined as
$s_{i}=(P_{\eta^{\prime}}-p_{i})^{2},i=1,2,3$, and
$s_{1}+s_{2}+s_{3}=M^{2}_{\eta^{\prime}}+M^{2}_{\eta}+2M^{2}_{\pi^{0}}$. The
Dalitz plot distribution of this decay can also described by kinematical
variables $X$ and $Y$

$\begin{split}X&=\frac{\sqrt{3}|s_{1}-s_{2}|}{2M_{\eta^{\prime}}Q_{\eta^{\prime}}}=\frac{\sqrt{3}|T_{\pi^{0}_{1}}-T_{\pi^{0}_{2}}|}{Q_{\eta^{\prime}}},\\\
Y&=\frac{(M_{\eta}+2M_{\pi^{0}})[(M_{\eta^{\prime}}-M_{\eta})^{2}-s_{3}]}{2M_{\eta^{\prime}}M_{\pi^{0}}Q_{\eta^{\prime}}}-1\\\
&=\frac{(M_{\eta}+2M_{\pi^{0}})T_{\eta}}{M_{\pi^{0}}Q_{\eta^{\prime}}}-1,\end{split}$
(6)

where $T_{i}$ denote kinetic energy of mesons in the rest frame of
$\eta^{\prime}$, and
$Q_{\eta^{\prime}}=M_{\eta^{\prime}}-M_{\eta}-2M_{\pi^{0}}$. The Dalitz plot
distribution can be expanded by $X$ and $Y$ around the center of the Dalitz
plot

$|\mathcal{M}(X,Y)|^{2}=|\mathcal{N}|^{2}(1+aY+bY^{2}+cX+dX^{2}+\cdots),$ (7)

which is known as general parameterization. Here $\mathcal{N}$ is a
normalization factor and parameter $c$ is fixed at 0 since two $\pi^{0}$s are
identical bosons. The general parameterization can be also expressed as

$\mathcal{M}(X,Y)=\mathcal{N}\\{1+\frac{a}{2}Y+\frac{1}{2}(b-\frac{a^{2}}{4})Y^{2}+\frac{d}{2}X^{2}+\cdots\\}.$
(8)

The NREFT amplitude of $\eta^{\prime}\to\eta\pi\pi$ can be decomposed to

$\begin{split}\mathcal{M}_{\eta^{\prime}\to\eta\pi^{0}\pi^{0}}&=\mathcal{M}^{tree}_{N}+\mathcal{M}^{one-
loop}_{N}+\mathcal{M}^{two-loop}_{N}+\cdots,\\\
\mathcal{M}_{\eta^{\prime}\to\eta\pi^{+}\pi^{-}}&=\mathcal{M}^{tree}_{C}+\mathcal{M}^{one-
loop}_{C}+\mathcal{M}^{two-loop}_{C}+\cdots.\end{split}$ (9)

The tree level amplitudes are

$\begin{split}\mathcal{M}^{tree}_{N}(s_{1},s_{2},s_{3})&=\sum_{i=0}^{2}G_{i}X^{i}_{3}+G_{3}(X_{1}-X_{2})^{2},\\\
\mathcal{M}^{tree}_{C}(s_{1},s_{2},s_{3})&=\sum_{i=0}^{2}H_{i}X^{i}_{3}+H_{3}(X_{1}-X_{2})^{2},\end{split}$
(10)

where $X_{k}=p^{0}_{k}-M_{\eta},k=1,2,3$, $p^{0}_{i}$ is the energy of
particle $i$ in the $\eta^{\prime}$ rest frame, and parameters $G_{i}$ are the
low-energy coupling coefficients of $\eta^{\prime}\to\eta\pi^{0}\pi^{0}$ decay
and $H_{i}$ for $\eta^{\prime}\to\eta\pi^{+}\pi^{-}$ decay. The charged decay
mode is introduced for the further description of loop level amplitude and we
assume $H_{i}=-G_{i}$ according to the isospin limit [18, 20, 19]. $G_{i}$ can
be evaulated by matching to the general parameterization

$\begin{split}G_{0}&=\mathcal{N}\\{1-\frac{a}{2}+\frac{1}{2}(b-\frac{a^{2}}{4})\\},\\\
G_{1}&=\mathcal{N}\\{\frac{a}{2}-(b-\frac{a^{2}}{4})\\}\frac{M_{\eta}+2M_{\pi^{0}}}{M_{\pi^{0}}Q_{\eta^{\prime}}},\\\
G_{2}&=\mathcal{N}(b-\frac{a^{2}}{4})\frac{(M_{\eta}+2M_{\pi^{0}})^{2}}{2M^{2}_{\pi^{0}}Q^{2}_{\eta^{\prime}}},\\\
G_{3}&=\mathcal{N}\frac{3d}{2Q^{2}_{\eta^{\prime}}}.\end{split}$ (11)

The loop level amplitude of $\eta^{\prime}\to\eta\pi^{0}\pi^{0}$ decay are

$\begin{split}\mathcal{M}^{one-
loop}_{N}(s_{1},s_{2},s_{3})&=\mathcal{B}_{N1}(s_{3})J_{+-}(s_{3})\\\
&+\mathcal{B}_{N2}(s_{3})J_{00}(s_{3}),\\\ \mathcal{M}^{two-
loop}_{N}(s_{1},s_{2},s_{3})&=C_{00}(s_{3})\mathcal{B}_{N2}(s_{3})J^{2}_{00}(s_{3})\\\
&+C_{00}(s_{3})\mathcal{B}_{N1}(s_{3})J_{00}(s_{3})J_{+-}(s_{3})\\\
&+2C_{x}(s_{3})\mathcal{B}_{C2}(s_{3})J_{00}(s_{3})J_{+-}(s_{3})\\\
&+2C_{x}(s_{3})\mathcal{B}_{C1}(s_{3})J^{2}_{+-}(s_{3}),\end{split}$ (12)

with one-loop function

$\begin{split}J_{ab}&=\frac{iq_{ab}(s_{k})}{8\pi\sqrt{s_{k}}},\\\
q^{2}_{ab}(s)&=\frac{\lambda(s,M^{2}_{a},M^{2}_{b})}{4s},\\\
\lambda(a,b,c)&=a^{2}+b^{2}+c^{2}-2(ab+ac+bc),\end{split}$ (13)

and neutral channel polynomials

$\begin{split}\mathcal{B}_{N1}(s_{3})&=2C_{x}(s_{3})[\sum_{i=0}^{2}H_{i}X^{i}_{3}+H_{3}\frac{4Q^{2}_{3}}{3s_{3}}q^{2}_{+-}(s_{3})],\\\
\mathcal{B}_{N2}(s_{3})&=C_{00}(s_{3})[\sum_{i=0}^{2}G_{i}X^{i}_{3}+G_{3}\frac{4Q^{2}_{3}}{3s_{3}}q^{2}_{00}(s_{3})],\end{split}$
(14)

and charged channel polynominals

$\begin{split}\mathcal{B}_{C1}(s_{3})&=2C_{+-}(s_{3})[\sum_{i=0}^{2}H_{i}X^{i}_{3}+H_{3}\frac{4Q^{2}_{3}}{3s_{3}}q^{2}_{+-}(s_{3})],\\\
\mathcal{B}_{C2}(s_{3})&=C_{x}(s_{3})[\sum_{i=0}^{2}G_{i}X^{i}_{3}+G_{3}\frac{4Q^{2}_{3}}{3s_{3}}q^{2}_{00}(s_{3})],\\\
\end{split}$ (15)

where

$\begin{split}C_{bc}(s_{a})&=C_{bc}+4D_{bc}q^{2}_{bc}(s_{a})+16F_{bc}q^{4}_{bc}(s_{a}),\\\
Q^{2}_{a}&=\frac{\lambda(M^{2}_{\eta^{\prime}},M^{2}_{a},s_{a})}{4M^{2}_{\eta^{\prime}}}.\end{split}$
(16)

The parameters $C_{i}$, $D_{i}$ and $F_{i}$ are coupling coefficients of
$\pi\pi$ interaction and are evaulated by matching to the effective range
expansion of $\pi\pi$ scattering, where $i$ represent differnet $\pi\pi$
rescattering channels:
($00$)$\pi^{0}\pi^{0}\to\pi^{0}\pi^{0}$;($x$)$\pi^{+}\pi^{-}\to\pi^{0}\pi^{0}$;($+-$)$\pi^{+}\pi^{-}\to\pi^{+}\pi^{-}$,

$\begin{split}C_{00}&=\frac{16\pi}{3}(a_{0}+2a_{2})(1-\xi),\\\
C_{x}&=\frac{16\pi}{3}(a_{2}-a_{0})(1+\frac{\xi}{3}),\\\
C_{+-}&=\frac{8\pi}{3}(2a_{0}+a_{2})(1+\xi),\\\
\xi&=\frac{M^{2}_{\pi^{\pm}}-M^{2}_{\pi^{0}}}{M^{2}_{\pi^{\pm}}},\\\
D_{00}&=\frac{4\pi}{3}(b_{0}+2b_{2}),\\\
D_{x}&=\frac{4\pi}{3}(b_{2}-b_{0}),\\\
D_{+-}&=\frac{2\pi}{3}(2b_{0}+b_{2}),\\\
F_{00}&=\frac{\pi}{3}(f_{0}+2f_{2}),\\\ F_{x}&=\frac{\pi}{3}(f_{2}-f_{0}),\\\
F_{+-}&=\frac{\pi}{6}(2f_{0}+f_{2}),\end{split}$ (17)

where $a_{i}$, $b_{i}$ and $f_{i}$ are S-wave scattering length, effective
ranges and shape parameters of isospin 0 and 2, respectively. $a_{0}$ and
$a_{2}$ are taken as free or fixed parameters in differnet cases in our study,
and $b_{i}$ are fixed to theoretical value $b_{0}=(0.276\pm 0.006)\times
M^{-2}_{\pi}$, $b_{2}=(-0.0803\pm 0.0012)\times M^{-2}_{\pi}$, and $f_{i}$ are
fixed to $0$. The $\pi\eta$ scattering terms are ignored because the $\pi\eta$
scattering is much weaker than the $\pi\pi$ scattering.

## References

  * Meissner  _et al._ [1997] U.-G. Meissner, G. Muller, and S. Steininger, Phys. Lett. B 406, 154 (1997), [Erratum: Phys.Lett.B 407, 454 (1997)].
  * Cabibbo [2004] N. Cabibbo, Phys. Rev. Lett. 93, 121801 (2004).
  * Cabibbo and Isidori [2005] N. Cabibbo and G. Isidori, JHEP 03 (2005), 021.
  * Colangelo  _et al._ [2006] G. Colangelo, J. Gasser, B. Kubis, and A. Rusetsky, Phys. Lett. B 638, 187 (2006).
  * Gamiz  _et al._ [2007] E. Gamiz, J. Prades, and I. Scimemi, Eur. Phys. J. C 50, 405 (2007).
  * Bissegger  _et al._ [2009] M. Bissegger, A. Fuhrer, J. Gasser, B. Kubis, and A. Rusetsky, Nucl. Phys. B 806, 178 (2009).
  * Budini and Fonda [1961] P. Budini and L. Fonda, Phys. Rev. Lett. 6, 419 (1961).
  * Batley  _et al._ [2006] J. R. Batley _et al._ (NA48/2 Collaboration), Phys. Lett. B 633, 173 (2006).
  * Batley  _et al._ [2009] J. R. Batley _et al._ (NA48/2 Collaboration), Eur. Phys. J. C 64, 589 (2009).
  * Bissegger  _et al._ [2008] M. Bissegger, A. Fuhrer, J. Gasser, B. Kubis, and A. Rusetsky, Phys. Lett. B 659, 576 (2008).
  * Ditsche  _et al._ [2009] C. Ditsche, B. Kubis, and U.-G. Meissner, Eur. Phys. J. C 60, 83 (2009).
  * Gullstrom  _et al._ [2009] C. O. Gullstrom, A. Kupsc, and A. Rusetsky, Phys. Rev. C 79, 028201 (2009).
  * Abouzaid  _et al._ [2008] E. Abouzaid _et al._ (KTeV Collaboration), Phys. Rev. D 78, 032009 (2008).
  * Adolph  _et al._ [2009] C. Adolph _et al._ (WASA-at-COSY Collaboration), Phys. Lett. B 677, 24 (2009).
  * Unverzagt  _et al._ [2009] M. Unverzagt _et al._ (Crystal Ball at MAMI, TAPS, A2 Collaborations), Eur. Phys. J. A 39, 169 (2009).
  * Prakhov  _et al._ [2009] S. Prakhov _et al._ (Crystal Ball at MAMI, A2 Collaborations), Phys. Rev. C 79, 035204 (2009).
  * Prakhov  _et al._ [2018] S. Prakhov _et al._ (A2 Collaboration), Phys. Rev. C 97, 065203 (2018).
  * Kubis and Schneider [2009] B. Kubis and S. P. Schneider, Eur. Phys. J. C 62, 511 (2009).
  * Ablikim  _et al._ [2018] M. Ablikim _et al._ (BESIII Collaboration), Phys. Rev. D 97, 012003 (2018).
  * Adlarson  _et al._ [2018] P. Adlarson _et al._ (A2 Collaboration), Phys. Rev. D 98, 012001 (2018).
  * Gonzàlez-Solís and Passemar [2018] S. Gonzàlez-Solís and E. Passemar, Eur. Phys. J. C 78, 758 (2018).
  * Ablikim  _et al._ [2022] M. Ablikim _et al._ (BESIII), Chin. Phys. C 46, 074001 (2022).
  * Ablikim  _et al._ [2010] M. Ablikim _et al._ (BESIII Collaboration), Nucl. Instrum. Meth. A 614, 345 (2010).
  * Ablikim  _et al._ [2020] M. Ablikim _et al._ (BESIII Collaboration), Chin. Phys. C 44, 040001 (2020).
  * Yu  _et al._ [2016] C. Yu _et al._ , in  _7th International Particle Accelerator Conference_ (2016) p. TUYA01.
  * Li  _et al._ [2017] X. Li _et al._ , Radiat. Detect. Technol. Methods 1, 13 (2017).
  * Guo  _et al._ [2017] Y. X. Guo _et al._ , Radiat. Detect. Technol. Methods 1, 15 (2017).
  * Cao  _et al._ [2020] P. Cao _et al._ , Nucl. Instrum. Meth. A 953, 163053 (2020).
  * Chen  _et al._ [2000] J. C. Chen, G. S. Huang, X. R. Qi, D. H. Zhang, and Y. S. Zhu, Phys. Rev. D 62, 034003 (2000).
  * Yang  _et al._ [2014] R. L. Yang, R. G. Ping, and H. Chen, Chin. Phys. Lett. 31, 061301 (2014).
  * [31] See supplemental material at [URL will be inserted by publisher] for additional information.
  * Fang  _et al._ [2021] S. S. Fang, B. Kubis, and A. Kupsc, Prog. Part. Nucl. Phys. 120, 103884 (2021).
  * Bondar  _et al._ [2013] A. E. Bondar _et al._ (Charm-Tau Factory Collaboration), Phys. Atom. Nucl. 76, 1072 (2013).
  * Peng  _et al._ [2020] H. P. Peng, Y. H. Zheng, and X. R. Zhou, Physics 49, 513 (2020).