[Denglu]dlredorange

#  Unbiased Estimation for Total Treatment Effect Under Interference Using
Aggregated Dyadic Data

Lu Deng  Tencent, Inc.Shenzhen, GuangdongChina<EMAIL_ADDRESS>,  Yilin
Li  Peking UniversityBeijingChina<EMAIL_ADDRESS>,  JingJing Zhang
Tencent, Inc.Shenzhen, GuangdongChina<EMAIL_ADDRESS>,  Yong Wang
Tencent, Inc.Shenzhen, GuangdongChina<EMAIL_ADDRESS>and  Chuan Chen
Tencent, Inc.Shenzhen, GuangdongChina<EMAIL_ADDRESS>
(2024)

###### Abstract.

In social media platforms, user behavior is often influenced by interactions
with other users, complicating the accurate estimation of causal effects in
traditional A/B experiments. This study investigates situations where an
individual’s outcome can be broken down into the sum of multiple pairwise
outcomes, a reflection of user interactions. These outcomes, referred to as
dyadic data, are prevalent in many social network contexts. Utilizing a
Bernoulli randomized design, we introduce a novel unbiased estimator for the
total treatment effect (TTE), which quantifies the difference in population
mean when all individuals are assigned to treatment versus control groups. We
further explore the bias of our estimator in scenarios where it is impractical
to include all individuals in the experiment, a common constraint in online
control experiments. Our numerical results reveal that our proposed estimator
consistently outperforms some commonly used estimators, underscoring its
potential for more precise causal effect estimation in social media
environments.

causal inference, interference, network effects, total treatment effect,
online control experiment, dyadic data

††copyright: acmcopyright††journalyear: 2024††doi:
XXXXXXX.XXXXXXX††conference: the ACM SIGKDD Conference on Knowledge Discovery
and Data Mining 2024; August 25–29, 2024; Barcelona, Spain††price:
15.00††isbn: 978-1-4503-XXXX-X/18/06††submissionid: 836††ccs: Probability and
Statistics Experimental design††ccs: Mathematics of computing Exploratory data
analysis

##  1\. Introduction

Online control experiments provide robust methodologies for tech companies to
assess the impact of innovative product strategies (Kohavi et al., 2013,
2014). In such experiments, the population is randomly divided into a
treatment group, which receives the treatment under investigation, and a
control group, which maintains the status quo. The treatment effect is then
estimated by comparing the difference in outcomes between these two groups
post-treatment. Classical causal inference techniques for estimating these
effects rely on the Stable Unit Treatment Value Assumption (SUTVA). SUTVA
posits that an individual’s outcome is independent of the treatment
assignments of others, ensuring that there is no interference between units.
In scenarios where SUTVA is valid, the difference-in-means estimator or
Horvitz-Thompson estimator provides an unbiased estimate of the treatment
effect (Imbens and Rubin, 2015).

However, social network settings commonly violate the SUTVA, as an
individual’s behavior can be influenced by their social connections (Eckles et
al., 2016; Pouget-Abadie et al., 2019; Karrer et al., 2021; Gui et al., 2015;
Ugander et al., 2013). This violation may introduce bias into estimates from
traditional randomized experiments. For instance, consider a social media
platform testing a feature to promote user messaging. This feature could
inadvertently stimulate friends in the control group to increase their
messaging activity. Consequently, both the treatment and control groups
exhibit increased activity, potentially leading to an underestimation of the
treatment’s actual effect. This phenomenon, known as network effects, is also
referred to as interference, spillover, or peer influence effects. Our primary
focus is to estimate the difference in population mean outcomes when all
individuals are assigned to either the treatment or control group.

Causal inference amidst interference has attracted increasing interest from
academia and industry in recent decades. There are two primary methods to
address this: the design-based approach and the model-assisted approach. The
design-based approach aims to reduce bias in effect estimation by
strategically assigning treatments. One such method is graph cluster
randomization (Gui et al., 2015; Ugander et al., 2013; Ugander and Yin, 2023;
Hayes and Moulton, 2017), which groups a network into clusters and treats each
cluster as a whole. While this can reduce bias, it may also decrease the
precision of results and is not suitable for highly interconnected networks.
Another method is LinkedIn’s ego-cluster experiment (Saint-Jacques et al.,
2019b), which focuses on individuals and their direct connections. However,
this method has low power and may still introduce bias. The model-assisted
approach involves constructing mathematical models of potential outcomes and
estimating parameters to understand effects (Särndal et al., 2003; Basse and
Airoldi, 2018; Toulis and Kao, 2013; Basse and Airoldi, 2015; Sussman and
Airoldi, 2017; Forastiere et al., 2021; Saint-Jacques et al., 2019a). The
exposure model (Manski, 2013; Aronow and Samii, 2012) is an example of this
approach, simplifying the treatment information into ”exposures” to study
causal effects. However, this approach assumes that these exposures fully
encapsulate the causal relationships, which may not always be the case. Recent
research has explored the implications of incorrect models and found that some
estimators can still be reliable under broader interference conditions, but
determining the conditions remains challenging (Chin, 2018; Sävje et al.,
2021; Sävje, 2023; Leung, 2022).

This work’s primary contribution is the development of novel estimators for
interference within the context of dyadic data. Dyadic data, capturing
interactions between pairs of users (dyads), enables the investigation of
complex dynamics and relationships among users. For example, in social network
analysis, dyads may refer to individuals within a network. Utilizing dyadic
information allows researchers to gain insights into intricate interactive
behavior. Dyadic data analysis is crucial across various scientific
disciplines, including economics (Portes and Rey, 2005), social sciences
(Apicella et al., 2012), and biomedical sciences (Luke and Harris, 2007). To
our knowledge, we are the first to propose using aggregated dyadic data to
estimate treatment effects amidst interference. We introduce and analyze
several new estimators that incorporate dyadic structures under different
experimental scenarios. The proposed framework is easy to set up and is widely
run at Weixin’s experiment platform. Our simulation results validate the
theorems and demonstrate the performance of the proposed estimators. In
summary, our contributions include:

  1. (1)

Formalizing the potential outcome model with dyadic outcomes and providing an
unbiased estimator of the total treatment effect under a simple randomized
experiment.

  2. (2)

Offering a detailed analysis of our estimator when the experiment includes
only a small fraction of the population.

  3. (3)

Proposing a novel experimental scheme that uses clustering to reduce the bias
of the estimator in subpopulation experiments.

  4. (4)

Conducting systematic simulations on two social networks and comparing our
method with traditional estimators, thereby providing a credible reference for
our method’s effectiveness.

The structure of this paper is as follows: Section 2 provides an initial
discussion on the potential outcome model and network interference. Section 3
details our proposed estimators. Section 4 presents simulation studies and a
real-world example that corroborates our theoretical findings. Section 5
concludes with a thorough discussion. Appendix A contains the proofs of some
of our theoretical results.

##  2\. Problem Setting

In this section, we’ll formalize the problem of causal inference under
interference, discuss some commonly used assumptions, and describe the
structure of dyadic data.

###  2.1. Potential outcome model

Consider a finite population of $n$ individuals. Let $W_{i}\in\\{0,1\\}$
represent the random binary treatment assignment for individuals $i$, where 1
and 0 denote treatment and control, respectively. Let
$\mathbf{W}=(W_{1},\dots,W_{n})^{\intercal}\in\Omega\overset{\Delta}{=}\\{0,1\\}^{n}$,
a vector of length $n$, be the experiment assignment vector for all
individuals. Let $Y_{i}(\mathbf{W})$ denote the potential outcome for
individual $i$ if the treatment vector $\mathbf{W}$ is implemented. Our goal
is to estimate the total treatment effect (TTE). TTE is defined as the
difference in population mean when all individuals are placed in treatment
versus control, which can be expressed as

(1)
$\tau\overset{\Delta}{=}\frac{1}{n}\sum_{i=1}^{n}\Big{[}Y_{i}(\mathbf{W}=\mathbf{1})-Y_{i}(\mathbf{W}=\mathbf{0})\Big{]}.$

where $\mathbf{1}=(1,\dots,1)^{\intercal},\mathbf{0}=(0,\dots,0)^{\intercal}$.
In what follows, expectations $\mathbb{E}$ are over the treatment assignment
only, potential outcomes are considered fixed. Under SUTVA,
$Y_{i}(\mathbf{W})=Y_{i}(W_{i})$ for all $\mathbf{W}\in\Omega$. Therefore, TTE
can be simplified as

(2)
$\tau=\frac{1}{n}\sum_{i=1}^{n}\Big{[}Y_{i}(W_{i}=1)-Y_{i}(W_{i}=0)\Big{]}.$

A commonly used estimator for $\tau$ is the difference-in-means estimator,
defined by

(3)
$\hat{\tau}_{dim}=\frac{\sum_{i=1}^{n}W_{i}Y_{i}}{\sum_{i=1}^{n}W_{i}}-\frac{\sum_{i=1}^{n}(1-W_{i})Y_{i}}{\sum_{i=1}^{n}(1-W_{i})}.$

Under Bernoulli randomization or complete randomization, $\hat{\tau}_{dim}$ is
an unbiased estimator for $\tau$ (Imbens and Rubin, 2015). However, when
interference exits, $\hat{\tau}_{dim}$ is biased and fails to provide a valid
estimate of $\tau$(Eckles et al., 2016; Gui et al., 2015).

###  2.2. Network Graph

In the following sections, we consider that the interference can be
characterized by a directed graph $G=(V,E)$, where the node set
$V=\\{1,2,\dots,n\\}$ represents the individuals, and $E=\\{(i,j)\\}$ is the
collection of edges that represent interference between individuals. We use
$|E|$ to denote the cardinality of $E$, which is the number of edges in the
graph. Each graph is represented by an adjacency matrix. Let $A$ be an
$n\times n$ matrix with entries $(A_{ij})_{1\leq i,j\leq n}$ where
$A_{ij}\in\\{0,1\\}$ with $A_{ii}=0$, and for $i\neq j$,

(4) $A_{ij}\overset{\Delta}{=}\begin{cases}0&\text{no edge between $i$ and
$j$}\\\ 1&\text{direct edge from $i$ to $j$}.\end{cases}$

We let $N_{u}(j)\overset{\Delta}{=}\\{i:A_{ij}=1\\}$ and
$N_{d}(j)\overset{\Delta}{=}\\{i:A_{ji}=1\\}$ represent the upstream neighbors
and downstream neighbors set of user $j$ with respect to the graph,
respectively. Throughout the article, we assume that the network graph $G$ is
fixed and known.

###  2.3. Neighborhood Interference

The neighborhood interference assumption (Ugander and Yin, 2023; Liu et al.,
2022; Yuan and Altenburger, 2022; Forastiere et al., 2022; Cortez et al.,
2022) posits that an individual’s outcome depends solely on the treatment
assignments of their direct neighbors within a specified network. This
assumption is frequently utilized when addressing interference, as it
simplifies the problem by focusing on local interactions.

  * (A1)

(Neighborhood Interference) The potential outcomes only depend on one’s own
and neighbors’ treatment assignments:
$Y_{i}(\mathbf{W})=Y_{i}(\mathbf{W}^{\prime})$ if $W_{i}=W_{i}^{\prime}$ and
$W_{j}=W_{j}^{\prime}$ for all $j\in N_{u}(i)\cup N_{d}(i)$.

Assuming neighborhood interference, the Horvitz-Thompson estimator combined
with exposure mapping can be employed to estimate the total treatment effect
(Aronow and Samii, 2017). However, the practical application of such
estimators is often hindered by their substantial variances or their
dependence on stringent exposure mapping assumptions.

###  2.4. Dyadic data setup

To estimate the total treatment effect, we propose an assumption that is
applicable in real-world scenarios and pertains to potential outcomes. We
denote the network graph characterizing the interference as $G=(V,E)$, and
adopt the neighborhood interference assumption (A(A1)). In a given scenario,
each direct edge $(i,j)\in E$ is associated with an outcome, represented as
$z_{i,j}$. The outcome for unit $j$, denoted as $Y_{j}$, is defined as the
summation of these associated outcomes from its upstream neighbors.

(5) $Y_{j}=\sum_{i\in N_{u}(j)}z_{i,j}.$

###### Example 2.1.

Within the framework of content sharing and forwarding, consider a scenario
where a user $i$ forwards a video to another user $j$. Upon receiving the
video, user $j$ clicks on it, leading to the sharing page, and then watches
the video. This action results in a stay duration, denoted by $z_{i,j}$. The
cumulative stay duration on the sharing page for user $j$ over a specific
period is calculated by summing up the stay durations $z_{i,j}$, each of which
is generated by clicking on videos forwarded by friends. In the absence of any
forwarding or clicking interaction between user $i$ and user $j$, we assign
$z_{i,j}$ a value of 0, signifying that no stay duration is generated under
such circumstances.

It’s important to note that $z_{i,j}$ does not necessarily need to be
symmetric for the pair $(i,j)$, which means that $z_{i,j}\neq z_{j,i}$ is
permissible. Additionally, we set $z_{i,i}=0$, indicating that we are
considering a graph without a self-loop. The value of $z_{i,i}$ will influence
the properties of our estimators, a point that will be illustrated in the
subsequent section. For each unit $j$, we denote $D_{j}$ as the summation of
paired outcomes from downstream neighbors:

(6) $D_{j}=\sum_{i\in N_{d}(j)}z_{j,i}.$

We refer to $D_{j}$ as the diffusion metric, as it characterizes the spillover
effect from individual $j$ diffusing to all of his neighbors. Moreover, we
refer to $Y_{j}$ and $D_{j}$ as aggregated dyadic data. For a more detailed
understanding of how $Y_{i}$ and $D_{i}$ are calculated, please refer to
Figure 1, which provides a visual representation of this process.

Figure 1. An illustration of the calculation of $Y_{i}$ and $D_{i}$. Suppose
node $i$ has 7 neighbor nodes, where nodes 1,2,3,4 are upstream neighbors, and
nodes 5,6,7 are downstream neighbors. Each direct edge is associated with an
outcome. $Y_{j}$ is the summation of the paired outcomes from upstream
neighbors, which is $z_{1,j}+z_{2,j}+z_{3,j}+z_{4,j}$. $D_{j}$ is the
summation of the paired outcomes from downstream neighbors, which is
$z_{j,5}+z_{j,6}+z_{j,7}$.

The forms of $Y_{i}$ and $D_{i}$ exhibit a kind of duality: one is a summation
from upstream neighbors, and the other is from downstream neighbors. In fact,
it can be shown that $\sum_{j=1}^{n}Y_{j}=\sum_{j=1}^{n}D_{j}$. Thus, the
total treatment effect defined by $Y_{i}$ and $D_{i}$ are identical. If we
denote $z_{i,j}(\mathbf{W}),Y_{j}(\mathbf{W})$ and $D_{j}(\mathbf{W})$ as the
counterfactual outcomes under treatment vector $\mathbf{W}$, then

$\tau=\frac{1}{n}\sum_{i=1}^{n}\Big{[}Y_{i}(\mathbf{W}=\mathbf{1})-Y_{i}(\mathbf{W}=\mathbf{0})\Big{]}=\frac{1}{n}\sum_{i=1}^{n}\Big{[}D_{i}(\mathbf{W}=\mathbf{1})-D_{i}(\mathbf{W}=\mathbf{0})\Big{]}.$

In past studies, researchers have only used $Y_{i}$ to estimate the total
treatment effect. Our work demonstrates that we can leverage the information
provided by $D_{j}$ to construct a more accurate estimate.

##  3\. Methodology

In this section, we will delineate the fundamental assumption that forms the
basis of our approach, introduce the novel estimator we have developed, and
conduct an in-depth analysis of the performance of this estimator under a
variety of conditions.

###  3.1. Dyadic Interference

Our overarching framework necessitates an interference assumption on the
dyadic outcome $z_{i,j}(\mathbf{W})$.

  * (A2)

For each $(i,j)$, the potential outcomes $z_{i,j}(\mathbf{W})$ only depend on
unit $i$ and unit $j$’s treatment assignments: for any
$\mathbf{W},\mathbf{W}^{\prime}\in\Omega$, if $W_{i}=W_{i}^{\prime}$ and
$W_{j}=W_{j}^{\prime}$, then
$z_{i,j}(\mathbf{W})=z_{i,j}(\mathbf{W}^{\prime})$.

Assumption (A(A2)) embodies a specific form of the neighbor interference
assumption, which is plausible considering that most dyadic data originate
from pairwise interactions between units. Under this assumption, we further
propose that $z_{i,j}(\mathbf{W})$ can be modeled as follows:

(7)
$z_{i,j}(\mathbf{W})=\alpha_{i,j}+\beta_{i,j}W_{i}+\gamma_{i,j}W_{j}+\zeta_{i,j}W_{i}W_{j}.$

The parameters $\beta_{i,j}$ and $\zeta_{i,j}$ are intended to capture network
interference, specifically addressing the potential spillover effect. On the
other hand, $\gamma_{i,j}$ is designed to capture the direct effect of
treatment.

Given (7), we can expand $Y_{j}(\mathbf{W})$ and $D_{j}(\mathbf{W})$ as

(8) $Y_{j}(\mathbf{W})=\sum_{i\in
N_{u}(j)}\Big{(}\alpha_{i,j}+\beta_{i,j}W_{i}+\gamma_{i,j}W_{j}+\zeta_{i,j}W_{i}W_{j}\Big{)},$

and

(9) $D_{j}(\mathbf{W})=\sum_{i\in
N_{d}(j)}\Big{(}\alpha_{j,i}+\beta_{j,i}W_{j}+\gamma_{j,i}W_{i}+\zeta_{j,i}W_{i}W_{j}\Big{)}.$

The total treatment effect estimand $\tau$ takes the value

$\displaystyle\ \tau=$
$\displaystyle\,\,\frac{1}{n}\sum_{i=1}^{n}\Big{[}Y_{i}(\mathbf{W}=\mathbf{1})-Y_{i}(\mathbf{W}=\mathbf{0})\Big{]}$
$\displaystyle=$ $\displaystyle\,\,\frac{1}{n}\sum_{j=1}^{n}\sum_{i\in
N_{u}(j)}\Big{(}\beta_{i,j}+\gamma_{i,j}+\zeta_{i,j}\Big{)}.$

Estimating these parameters presents a significant challenge due to the
substantial number of model parameters, which scales with the number of edges
in the network $O(|E|)$. Therefore, we do not estimate these parameters in
practice.

We can rewrite $\eqref{equa: yjmodel}$ as follows:

$Y_{j}(\mathbf{W})=\alpha_{j}+\gamma_{j}W_{j}+\sum_{i\in
N_{u}(j)}\Big{(}\beta_{i,j}W_{i}+\zeta_{i,j}W_{i}W_{j}\Big{)},$

where $\alpha_{j}=\sum_{i\in N_{u}(j)}\alpha_{i,j}$ and $\gamma_{j}=\sum_{i\in
N_{u}(j)}\gamma_{i,j}$. By setting the parameter $\zeta_{i,j}=0$, the model
simplifies to the heterogeneous additive network effects model, as discussed
by Yu et al.(Yu et al., 2022) in their 2022 paper. This demonstrates that our
model is a more general form capable of encompassing the model proposed by Yu
and colleagues when certain conditions are applied.

To simplify the formula, we assume that $z_{i,j}=0$ if $(i,j)\notin E$. Then
the TTE $\tau$ can be expressed as follows:

(10) $\tau=\frac{1}{n}\sum_{i\neq
j}\Big{(}\beta_{i,j}+\gamma_{i,j}+\zeta_{i,j}\Big{)}.$

###  3.2. Estimator

Consider a simple Bernoulli randomization experiment, where each unit in the
experiment is randomly assigned to either the treatment group or the control
group. Let $\pi=\mathbb{P}(W_{i}=1)$ represent the treatment probability. We
define two Horvitz-Thompson (HT) estimators:

(11)
$\hat{\tau}^{1}(\pi)=\frac{1}{n}\sum_{i=1}^{n}\frac{W_{i}Y_{i}}{\pi}-\frac{1}{n}\sum_{i=1}^{n}\frac{(1-W_{i})Y_{i}}{1-\pi},$

and

(12)
$\hat{\tau}^{2}(\pi)=\frac{1}{n}\sum_{i=1}^{n}\frac{W_{i}D_{i}}{\pi}-\frac{1}{n}\sum_{i=1}^{n}\frac{(1-W_{i})D_{i}}{1-\pi}.$

The parameter $\pi$ in the estimators suggests that the properties of the
estimators are dependent on $\pi$. The $\hat{\tau}^{1}(\pi)$ is a commonly
used estimator in causal inference. It is widely recognized that
$\hat{\tau}^{1}(\pi)$ is biased for the total treatment effect $\tau$ when
interference exists (Eckles et al., 2016; Gui et al., 2015).

Finally, we define a new estimator $\hat{\tau}(\pi)$ as the sum of
$\hat{\tau}^{1}(\pi)$ and $\hat{\tau}^{1}(\pi)$. That is,

(13)
$\hat{\tau}(\pi)\overset{\triangle}{=}\hat{\tau}^{1}(\pi)+\hat{\tau}^{2}(\pi).$

We use $\hat{\tau}(\pi)$ to estimate the total treatment effect $\tau$. We
will next demonstrate the properties of these three estimators under various
scenarios. All the proofs can be found in the appendix.

###  3.3. Full population experiment

We begin by considering the simplest scenario in which the entire population
participates in the experiment. The following result is provided:

###### Proposition 3.1.

$\mathbb{E}\big{[}\hat{\tau}^{1}(\pi)\big{]}=\frac{1}{n}\sum_{i\neq
j}\big{(}\gamma_{i,j}+\zeta_{i,j}\pi\big{)}.$

The estimator’s bias is given by $-\frac{1}{n}\sum_{i\neq
j}\Big{(}\beta_{i,j}+\zeta_{i,j}(1-\pi)\Big{)}$, which can scale with the
average spillover effect across the edges between the treated and control
groups. In practice, the sign of $\beta_{i,j}$ and $\zeta_{i,j}$ may not be
identical, implying that we do not even know whether our estimate is an
overestimate or an underestimate. Thankfully, with the assistance of
$\hat{\tau}^{2(\pi)}$, we can unbiasedly estimate the total treatment effect.
This additional information can then be employed to significantly simplify the
estimation of causal effects. The following result characterize the bias of
$\hat{\tau}^{2}(\pi)$.

###### Proposition 3.2.

$\mathbb{E}\big{[}\hat{\tau}^{2}(\pi)\big{]}=\frac{1}{n}\sum_{i\neq
j}\big{(}\beta_{i,j}+\zeta_{i,j}\pi\big{)}.$

$\hat{\tau}^{2}(\pi)$ is also a biased estimator for the total treatment
effect. From Proposition 3.1 and 3.2, we can infer that

$\mathbb{E}\big{[}\hat{\tau}(\pi)\big{]}=\frac{1}{n}\sum_{i\neq
j}\big{(}\beta_{i,j}+\gamma_{i,j}+2\pi\zeta_{i,j}\big{)}.$

Consequently, for $\pi=0.5$, we derive
$\mathbb{E}\big{[}\hat{\tau}(0.5)\big{]}=\frac{1}{n}\sum_{i\neq
j}\big{(}\beta_{i,j}+\gamma_{i,j}+\zeta_{i,j}\big{)}=\tau$. Therefore, when
$\pi=0.5$, $\hat{\tau}(0.5)$ serves as an unbiased estimator of the total
treatment effect. We summarize the analysis above for the estimator
$\hat{\tau}$ in the following theorem.

###### Theorem 3.1.

Assume (A(A2)) and that the experiment includes all units, the expectation of
$\hat{\tau}(\pi)$ is

$\mathbb{E}\big{[}\hat{\tau}(\pi)\big{]}=\frac{1}{n}\sum_{i\neq
j}\big{(}\beta_{i,j}+\gamma_{i,j}+2\pi\zeta_{i,j}\big{)}.$

When $\pi=0.5$, that is, when treatment and control assignments are equally
probable, $\hat{\tau}(0.5)$ is an unbiased estimator of the total treatment
effect $\tau$.

There are several important points to note in this subsection. Firstly, our
proposed estimator relies heavily on assumption (A(A2)), which states that
$z_{i,j}$ depends solely on the treatment assignments of unit $i$ and unit
$j$. This assumption is akin to the neighborhood interference assumption.
While a similar approach to the method in (Athey et al., 2018) can be used to
test assumption A(A2), it may not be as efficient. Secondly, we assume that
$z_{i,i}=0$ for all $i$. This assumption is crucial for the unbiased property
of $\hat{\tau}$. If we were to allow $z_{i,i}\neq 0$, then it can be shown
that when $\pi=0.5$, the bias of $\hat{\tau}$ is
$\frac{1}{n}\sum_{i}(\beta_{i,i}+\gamma_{i,i}+\zeta_{i,i})$. However, even
when $z_{i,i}\neq 0$, compared to $\hat{\tau}^{1}$, our proposed estimator
$\hat{\tau}$ still shows an improvement in terms of bias, which further
underscores its potential usefulness in practice. We leave the proof of this
statement to the reader.

###  3.4. Sub-population experiment

Theorem 3.1 provides an unbiased estimator for the total treatment effect.
However, this holds only when we conduct an controlled experiment with the
entire population. In reality, many industry experimental platforms run
thousands of experiments simultaneously every day using a layered structure
design (see (Tang et al., 2010; Kohavi et al., 2013; Xiong et al., 2020) for
more details about the layer and experiment structure). Typically, an
experiment encompasses only a small fraction of users, often less than 10%
according to (Kohavi et al., 2013). This leads us to a natural question: Does
the conclusion of Theorem 3.1 still holds when the experiment includes only a
subset of the population? This is a critical question, as it directly impacts
the applicability and relevance of our findings in real-world settings.

Let’s suppose we randomly select a proportion $p(0<p\leq 1)$ of the population
to participate in the experiment. For each user $i$, we use
$V_{i}\in\\{0,1\\}$ to denote whether user $i$ is included in the experiment
or not. Users who are included in the experiment are then assigned to either
the treatment or control group through a Bernoulli randomization process.
Let’s denote $\pi=\mathbb{P}(W_{i}=1|V_{i}=1)$ as the probability of users in
the experiment receiving the treatment. To simplify the problem, we’ll assume
that users who are not included in the experiment receive the control group
strategy, i.e., $\mathbb{P}(W_{i}=0|V_{i}=0)=1$. Under these assumptions, the
two HT estimators become functions of $p$ and $\pi$:

(14)
$\hat{\tau}^{1}(p,\pi)=\frac{1}{n}\sum_{i=1}^{n}\frac{V_{i}W_{i}Y_{i}}{p\pi}-\frac{1}{n}\sum_{i=1}^{n}\frac{V_{i}(1-W_{i})Y_{i}}{p(1-\pi)},$

and

(15)
$\hat{\tau}^{2}(p,\pi)=\frac{1}{n}\sum_{i=1}^{n}\frac{V_{i}W_{i}D_{i}}{p\pi}-\frac{1}{n}\sum_{i=1}^{n}\frac{V_{i}(1-W_{i})D_{i}}{p(1-\pi)}.$

This setup allows us to explore how the proportion of the population included
in the experiment and the probability of receiving treatment affect the bias
of our estimators. The bias property of these two estimators are presented
below.

###### Proposition 3.3.

$\mathbb{E}\big{[}\hat{\tau}^{1}(p,\pi)\big{]}=\frac{1}{n}\sum_{i\neq
j}\big{(}\gamma_{j,i}+p\pi\zeta_{j,i}\big{)}.$

###### Proposition 3.4.

$\mathbb{E}\big{[}\hat{\tau}^{2}(p,\pi)\big{]}=\frac{1}{n}\sum_{i\neq
j}\big{(}\beta_{i,j}+p\pi\zeta_{i,j}\big{)}.$

Figure 2. Visualization of the two-stage experiment process.

Note that $\hat{\tau}(p,\pi)=\hat{\tau}^{1}(p,\pi)+\hat{\tau}^{2}(p,\pi)$, We
have the following theorem about $\hat{\tau}(p,\pi)$.

###### Theorem 3.2.

Assume (A(A2)). The expectation of $\hat{\tau}(p,\pi)$ is

(16) $\mathbb{E}\big{[}\hat{\tau}(p,\pi)\big{]}=\frac{1}{n}\sum_{i\neq
j}\big{(}\beta_{i,j}+\gamma_{i,j}+2p\pi\zeta_{i,j}\big{)}.$

and the bias is

(17)
$\mathbb{E}\big{[}\hat{\tau}(p,\pi)-\tau\big{]}=-\frac{1-2p\pi}{n}\sum_{i\neq
j}\zeta_{i,j}.$

Theorem 3.1 can be seen as a special case of Theorem 3.2. Specifically, when
$p=1$, meaning the entire population is included in the experiment, Theorem
3.2 simplifies to Theorem 3.1. However, in general, we have $\sum_{i\neq
j}\zeta_{i,j}\neq 0$. This implies that the estimator $\hat{\tau}(p,\pi)$ in a
sub-population experiment (where $p<1$) is biased for the total treatment
effect. Furthermore, there are $O(|E|)$ terms in the summands of $\sum_{i\neq
j}\zeta_{i,j}$, which suggests that the bias is
$O\Big{(}\frac{|E|}{n}\Big{)}$. When comparing the biases of the two
estimators $\hat{\tau}(p,\pi)$ and $\hat{\tau}^{1}(p,\pi)$, we find that the
difference in their expected biases is given by:

$\mathbb{E}\big{[}\hat{\tau}^{1}(p,\pi)-\tau\big{]}-\mathbb{E}\big{[}\hat{\tau}(p,\pi)-\tau\big{]}=\frac{1}{n}\sum_{i\neq
j}(\beta_{i,j}+p\pi\zeta_{i,j}).$

This result implies that the absolute bias of $\hat{\tau}(p,\pi)$ is smaller
than the absolute bias of $\hat{\tau}^{1}(p,\pi)$. In other words, the
estimator $\hat{\tau}(p,\pi)$ is more accurate in estimating the total
treatment effect when only a subset of the population is included in the
experiment.

###  3.5. Sub-population two-stage experiment

For the remainder of this paper, we will focus on the scenario where the
experiment does not include all users in the population. Under a simple
Bernoulli randomization experiment without full population coverage, the
estimator $\hat{\tau}(p,\pi)$ is biased for the total treatment effect. This
naturally leads us to question whether it is possible to find an estimator
with a smaller bias. In this subsection, we provide some positive results by
considering a special design that is similar to two-stage randomized
experiments and the cluster network experiment randomization process (Karrer
et al., 2021). This design aims to reduce bias and improve the accuracy of our
estimates, thereby providing a more reliable understanding of the treatment
effect.

Given a network graph $G$, let’s assume we can partition the graph into
disjoint clusters using common clustering algorithm like Louvain community
detection (Blondel et al., 2008). Our two-stage experiment scheme consists of
two steps:

  1. (1)

First, we randomly select $p$ percent of the clusters to form the sub-
population in the experiment.

  2. (2)

Second, we randomly assign the units in the sub-population to either the
treatment group with probability $\pi$, or the control group with probability
$1-\pi$.

A visualization of this two-stage experiment process is provided in Figure 2.
Similar to the setting in subsection 3.4, for each user $i$, we use
$V_{i}\in\\{0,1\\}$ to indicate whether they are included in the experiment,
with $\mathbb{P}(V_{i}=1)=p$. Users who are not included in the experiment are
assigned the control strategy. To distinguish this scenario from the one
discussed in subsection 3.4, we denote the two HT estimators in (14) and (15)
as $\hat{\tau}^{1}_{c}(p,\pi)$ and $\hat{\tau}^{2}_{c}(p,\pi)$, respectively.
The HT estimator of the total treatment effect is then given by
$\hat{\tau}_{c}(p,\pi)=\hat{\tau}^{1}_{c}(p,\pi)+\hat{\tau}^{2}_{c}(p,\pi)$.
Finally, let $\bar{\sigma}$ denote the proportion of each user’s friends
within the same cluster relative to their total number of friends, i.e.,

$\bar{\sigma}=\frac{1}{n}\sum_{i}\frac{|\\{j:j\text{ and }i\text{ are in the
same cluster}\\}|}{d_{i}}$

We have the following results.

###### Proposition 3.5.

$\mathbb{E}\big{[}\hat{\tau}^{1}_{c}(p,\pi)\big{]}=\frac{1}{n}\sum_{i\neq
j}\Big{(}\gamma_{i,j}+\big{(}\bar{\sigma}\pi+(1-\bar{\sigma})p\pi\big{)}\zeta_{i,j}\Big{)}.$

###### Proposition 3.6.

$\mathbb{E}\big{[}\hat{\tau}^{2}_{c}(p,\pi)\big{]}=\frac{1}{n}\sum_{i\neq
j}\Big{(}\beta_{i,j}+\big{(}\bar{\sigma}\pi+(1-\bar{\sigma})p\pi\big{)}\zeta_{i,j}\Big{)}.$

###### Theorem 3.3.

Assume (A(A2)). Under a two-stage experiment scheme, the expectation of
$\hat{\tau}_{c}(p,\pi)$ is

(18) $\mathbb{E}\big{[}\hat{\tau}_{c}(p,\pi)\big{]}=\frac{1}{n}\sum_{i\neq
j}\Big{(}\beta_{i,j}+\gamma_{i,j}+2\pi\big{(}\bar{\sigma}+(1-\bar{\sigma})p\big{)}\zeta_{i,j}\Big{)}.$

The bias is

(19)
$\mathbb{E}\big{[}\hat{\tau}^{2}_{c}(p,\pi)-\tau\big{]}=-\frac{(1-2p\pi)(1-\bar{\sigma})}{n}\sum_{i\neq
j}\zeta_{i,j}.$

From equations (17) and (19), we can observe that the absolute bias of the
estimator $\hat{\tau}_{c}(p,\pi)$ in the two-stage experiment is reduced by
$\frac{\bar{\sigma}(1-2p\pi)}{n}\sum_{i\neq j}\zeta_{i,j}$ compared to the
absolute bias of the estimator $\hat{\tau}(p,\pi)$ in the simple randomized
experiment. This indicates that the two-stage experiment can further reduce
the bias of the estimator without sacrificing any statistical power. This is a
significant finding, as it suggests that a two-stage experimental design can
provide more accurate estimates of the treatment effect, thereby improving the
reliability of our results. In practical applications, the value of
$\bar{\sigma}$ typically falls between $30\%$ to $40\%$ (Karrer et al., 2021),
meaning that the bias can be reduced by a minimum of 30%.

##  4\. Numerical Experiments

###  4.1. Simulated data

To compare the performance of the estimators mentioned earlier under various
conditions, we conducted several simulation studies. These studies were based
on two real-world social networks: FB-Standard3 and FB-Cornell5 111Network
topology data can be found in https://networkrepository.com/socfb-
Stanford3.php and https://networkrepository.com/socfb-Cornell5.php. The FB-
Standard3 network comprises 11,586 nodes and 568,309 edges, yielding an
average degree of 98.1. On the other hand, the FB-Cornell5 network consists of
18,660 nodes and 790,777 edges, resulting in an average degree of 84.7. To
initiate our simulations, we first generated the paired outcomes based on the
model defined in equation (7):

$z_{i,j}(W_{i},W_{j})=\alpha_{i,j}+\beta_{i,j}W_{i}+\gamma_{i,j}W_{j}+\zeta_{i,j}W_{i}W_{j}.$

Next, we aggregated the paired outcomes to obtain $Y_{i}$ and $D_{i}$. The
parameters $\alpha_{i,j},\beta_{i,j},\gamma_{i,j},\zeta_{i,j}$ were generated
from either a uniform or a Bernoulli distribution. Specifically, we considered
the following two cases:

  * •

For uniform distributed coefficients, $\alpha_{i,j}=1$, $\beta_{i,j}\sim$
Uniform $(0,0.5)$, $\gamma_{i,j}\sim\text{Uniform}(0,1)$,
$\zeta_{i,j}\sim\text{Uniform}$ $(0,0.5)$ for all pair $(i,j)\in E$
independently.

  * •

For Bernoulli distributed coefficients, $\alpha_{i,j}=0$, $\beta_{i,j}\sim$
Bernoulli $(0.25)$, $\gamma_{i,j}\sim\text{Bernoulli}(0.5)$,
$\zeta_{i,j}\sim\text{Bernoulli}(0.25)$ for all pair $(i,j)\in E$
independently.

Once the parameters were generated, they remained fixed throughout the
simulation. We considered a total of four different simulation settings. For
each setting, we performed 10,000 repetitions to calculate the bias of the
estimators and assess their performance.

Firstly, we simulated a scenario where the entire population is included in
the experiment. We varied the treatment probability $\pi$ across the set
$\\{0.1,0.3,0.5,0.7,$ $0.9\\}$. The results of these simulations are presented
in Figure 3. From the figure, it is evident that both $\hat{\tau}^{1}(\pi)$
and $\hat{\tau}^{2}(\pi)$ exhibit bias under any given treatment probability.
However, $\hat{\tau}(\pi)$ remains unbiased when $\pi=0.5$. This observation
aligns perfectly with the conclusion drawn in Theorem 3.1.

Secondly, we simulated a scenario where only a subset of the population is
included in the experiment. In this case, we fixed the treatment probability
at $\pi=0.5$, varied the in-experiment probability $p$ across the set
$\\{0.05,0.1,0.25,0.5,1\\}$. The results of these simulations are presented in
Figure 4. From the figure, it is evident that when $p<1$ , all three
estimators $\hat{\tau}^{1}(p,\pi)$, $\hat{\tau}^{2}(p,\pi)$ and
$\hat{\tau}(p,\pi)$ exhibit bias. Not surprisingly, the bias of
$\hat{\tau}(p,\pi)$ decreases as $p$ approaches 1. This suggests that
increasing the traffic ratio of the experiment can help reduce the bias
introduced by interference. Furthermore, we observe that the absolute bias of
$\hat{\tau}(p,\pi)$ is consistently smaller than the absolute bias of
$\hat{\tau}^{1}(p,\pi)$ under any given in-experiment probability. This
indicates the superiority of $\hat{\tau}(p,\pi)$ over $\hat{\tau}^{1}(p,\pi)$.

Lastly, we simulated a two-stage experiment scenario. We first applied the
Louvain algorithm (Blondel et al., 2008) with a fixed random seed and a
resolution parameter set to 10 to group the units into clusters. This resulted
in 240 clusters for the FB-Cornell5 graph and 192 clusters for the FB-
Stanford3 graph. As in the previous scenario, we fixed the treatment
probability at $\pi=0.5$, and varied the in-experiment probability $p$ across
the set $\\{0.05,0.1,0.25,0.5,1\\}$. The results of these simulations are
presented in Figure 5. From the figure, it is evident that the two-stage
experiment scheme can further reduce the bias of $\hat{\tau}(p,\pi)$ compared
to the simple Bernoulli randomized experiment. This suggests that a two-stage
experimental design, where units are first grouped into clusters, can be a
more effective strategy for reducing bias in the presence of interference.

Figure 3.  Visualizations of the performance of three TTE estimators under
full population Bernoulli design on FB-Standard3 and FB-Cornell5 networks for
both Uniform and Bernoulli potential outcomes models. Each line’s height
represents the relative bias of the estimator.  Figure 4.  Visualizations of
the performance of three TTE estimators under sub-population Bernoulli design
on FB-Standard3 and FB-Cornell5 networks for both Uniform and Bernoulli
potential outcomes models. Each line’s height represents the relative bias of
the estimator.  Figure 5.  Visualizations of the performance of three TTE
estimators under sub-population two-stage experiment design on FB-Standard3
and FB-Cornell5 networks for both Uniform and Bernoulli potential outcomes
models. Each line’s height represents the relative bias of the estimator.

###  4.2. Real Data

Table 1. Estimates with a re-scaled confidence interval of the metrics in the two-stage experiment. | $\hat{\tau}^{1}$ | $\hat{\tau}^{2}$ | $\hat{\tau}$
---|---|---|---
Metric 1 | -0.605% $\pm$ 0.255% | 3.498% $\pm$ 1.198% | 0.621% $\pm$ 0.5%
Metric 2 | -0.712% $\pm$ 0.349% | 3.381% $\pm$ 1.171% | 0.678% $\pm$ 0.5%

We present an empirical experiment conducted on Weixin’s experiment platform,
utilizing aggregated dyadic data to estimate the total treatment effect.
Weixin, a popular messaging and calling application, allows users to interact
with their contacts. A built-in feature of Weixin, known as ”Channels”, is
akin to TikTok, enabling users to view and discover a plethora of personalized
short videos. Similar to Example 2.1, users typically share intriguing videos
with their friends, who can then access the ”Channels” page by clicking on the
shared video, watching it, and exploring other videos.

A two-stage experiment spanning two weeks was executed to evaluate a novel
algorithm. The experiment contains about 4% of the users. In the treatment
group, the algorithm increased the likelihood of users sharing videos. The
underlying hypothesis was that promoting video sharing would lead to users
exiting the ”Channels” page to engage in conversations with friends, thereby
reducing their consumption but increasing that of their friends. The primary
objective was to ascertain whether this new algorithm could enhance user
activity and consumption when deployed to all users.

The outcomes for two crucial metrics are detailed in Table 1. All data and
attributes are collected after user approvals and data masking to protect user
privacy. Each metric was re-scaled such that the estimator $\hat{\tau}$ has a
confidence interval (CI) width of one. The results align with our initial
expectations. The $\hat{\tau}^{1}$ indicates a significant decrease, while the
$\hat{\tau}$ exhibits a substantial increase, implying that the total
treatment effect of the new algorithm is positive. Upon deploying the new
algorithm to all users, an upward trend in the consumption of the entire user
base was observed, corroborating the results of $\hat{\tau}$.

##  5\. Discussion

The estimator $\hat{\tau}$, as developed in our study, has the potential to
estimate the total treatment effect within a dyadic data structure. This
estimator exhibits no bias when the experiment encompasses the entire
population and the probability of treatment is set at 0.5. Compared to the
conventional estimator $\hat{\tau}^{1}$, our proposed estimator $\hat{\tau}$
demonstrates reduced bias when the experiment includes only a small subset of
the population. Furthermore, we suggest the application of $\hat{\tau}$ within
a two-stage experimental scheme, a strategy that can further mitigate the
estimator’s bias.

Several potential avenues exist for extending our current work. Firstly, we
have not yet presented the variance properties of these estimators. In
practical applications, we utilize the sample variance of $Y_{i}$ and $D_{i}$
to facilitate statistical inference. Our simulation findings suggest that the
sample variance provides a reasonably accurate estimate of the true variance.
Future work could involve deriving the expression for the variance of these
estimators, demonstrating the convergence rate, and proposing alternative
estimators for variance within a dyadic data structure.

Secondly, our work relies on an important assumption (A(A2)), which posits
that the potential pair outcomes $z_{i,j}(\mathbf{W})$ depend solely on the
treatment assignments of units $i$ and $j$. This assumption could be extended
to a more general case, for instance, by assuming that the potential pair
outcomes $z_{i,j}(\mathbf{W})$ depend on the treatment assignments of units
$i,j$, and all their respective neighbors. Although an assumption of this
nature significantly increases the complexity of inferring the total treatment
effect, researchers may be able to leverage dyadic data to construct
estimators that exhibit superior properties compared to commonly used
estimators, such as the difference-in-means estimator and the Horvitz-Thompson
estimator.

Thirdly, this paper considers simple Bernoulli randomization experimentation.
There are other experimental randomization schemes, such as cluster
randomization or ego cluster randomization, which are useful for dealing with
interference. Future research could explore the integration of our estimators
with these alternative randomization schemes, potentially leading to
insightful findings for estimating the total treatment effect.

Fourthly, dyadic data represents a unique type of data structure. However, in
some scenarios, we may not have access to dyadic data, meaning that the
outcome $Y_{i}$ cannot be decomposed into the sum of edge-level outcomes.
Consider, for instance, a user activity level metric, which is binary,
represented by 1 (active) or 0 (inactive). It is not feasible to express the
user activity level metric in a summation form. For such metrics, researchers
could consider integrating our method with the approaches proposed in (Li and
Wager, 2022; Hu et al., 2022) to estimate the total treatment effect.

Exploring these potential extensions could significantly broaden the
applicability of our estimators to more general models, while also providing
the necessary theoretical justifications. These promising research directions
present exciting opportunities for future work, and we look forward to seeing
how they will further advance our understanding in this field.

Code Availability Statement: https://github.com/adamdenglu/diffusion-metric-
with-aggregated-dyadic-data

## References

  * (1)
  * Apicella et al. (2012) Coren L Apicella, Frank W Marlowe, James H Fowler, and Nicholas A Christakis. 2012.  Social networks and cooperation in hunter-gatherers.  _Nature_ 481, 7382 (2012), 497–501.
  * Aronow and Samii (2012) Peter M Aronow and Cyrus Samii. 2012.  Estimating average causal effects under general interference. In _Summer Meeting of the Society for Political Methodology, University of North Carolina, Chapel Hill, July_. 19–21.
  * Aronow and Samii (2017) Peter M Aronow and Cyrus Samii. 2017.  Estimating average causal effects under general interference, with application to a social network experiment.  (2017).
  * Athey et al. (2018) Susan Athey, Dean Eckles, and Guido W Imbens. 2018.  Exact p-values for network interference.  _J. Amer. Statist. Assoc._ 113, 521 (2018), 230–240.
  * Basse and Airoldi (2015) Guillaume W Basse and Edoardo M Airoldi. 2015.  Optimal model-assisted design of experiments for network correlated outcomes suggests new notions of network balance.  _arXiv preprint arXiv:1507.00803_ (2015).
  * Basse and Airoldi (2018) Guillaume W Basse and Edoardo M Airoldi. 2018.  Model-assisted design of experiments in the presence of network-correlated outcomes.  _Biometrika_ 105, 4 (2018), 849–858.
  * Blondel et al. (2008) Vincent D Blondel, Jean-Loup Guillaume, Renaud Lambiotte, and Etienne Lefebvre. 2008.  Fast unfolding of communities in large networks.  _Journal of statistical mechanics: theory and experiment_ 2008, 10 (2008), P10008.
  * Chin (2018) Alex Chin. 2018\.  Central limit theorems via Stein’s method for randomized experiments under interference.  _arXiv preprint arXiv:1804.03105_ (2018).
  * Cortez et al. (2022) Mayleen Cortez, Matthew Eichhorn, and Christina Yu. 2022\.  Staggered rollout designs enable causal inference under interference without network knowledge.  _Advances in Neural Information Processing Systems_ 35 (2022), 7437–7449.
  * Eckles et al. (2016) Dean Eckles, Brian Karrer, and Johan Ugander. 2016\.  Design and analysis of experiments in networks: Reducing bias from interference.  _Journal of Causal Inference_ 5, 1 (2016), 20150021\.
  * Forastiere et al. (2021) Laura Forastiere, Edoardo M Airoldi, and Fabrizia Mealli. 2021.  Identification and estimation of treatment and interference effects in observational studies on networks.  _J. Amer. Statist. Assoc._ 116, 534 (2021), 901–918.
  * Forastiere et al. (2022) Laura Forastiere, Fabrizia Mealli, Albert Wu, and Edoardo M Airoldi. 2022.  Estimating causal effects under network interference with Bayesian generalized propensity scores.  _The Journal of Machine Learning Research_ 23, 1 (2022), 13101–13161.
  * Gui et al. (2015) Huan Gui, Ya Xu, Anmol Bhasin, and Jiawei Han. 2015\.  Network a/b testing: From sampling to estimation. In _Proceedings of the 24th International Conference on World Wide Web_. 399–409.
  * Hayes and Moulton (2017) Richard J Hayes and Lawrence H Moulton. 2017.  _Cluster randomised trials_.  CRC press.
  * Hu et al. (2022) Yuchen Hu, Shuangning Li, and Stefan Wager. 2022.  Average direct and indirect causal effects under interference.  _Biometrika_ 109, 4 (2022), 1165–1172.
  * Imbens and Rubin (2015) Guido W Imbens and Donald B Rubin. 2015.  _Causal inference in statistics, social, and biomedical sciences_.  Cambridge University Press.
  * Karrer et al. (2021) Brian Karrer, Liang Shi, Monica Bhole, Matt Goldman, Tyrone Palmer, Charlie Gelman, Mikael Konutgan, and Feng Sun. 2021\.  Network experimentation at scale. In _Proceedings of the 27th acm sigkdd conference on knowledge discovery & data mining_. 3106–3116.
  * Kohavi et al. (2013) Ron Kohavi, Alex Deng, Brian Frasca, Toby Walker, Ya Xu, and Nils Pohlmann. 2013\.  Online controlled experiments at large scale. In _Proceedings of the 19th ACM SIGKDD international conference on Knowledge discovery and data mining_. 1168–1176.
  * Kohavi et al. (2014) Ron Kohavi, Alex Deng, Roger Longbotham, and Ya Xu. 2014\.  Seven rules of thumb for web site experimenters. In _Proceedings of the 20th ACM SIGKDD international conference on Knowledge discovery and data mining_. 1857–1866.
  * Leung (2022) Michael P Leung. 2022\.  Causal inference under approximate neighborhood interference.  _Econometrica_ 90, 1 (2022), 267–293.
  * Li and Wager (2022) Shuangning Li and Stefan Wager. 2022.  Random graph asymptotics for treatment effect estimation under network interference.  _The Annals of Statistics_ 50, 4 (2022), 2334–2358.
  * Liu et al. (2022) Yang Liu, Yifan Zhou, Ping Li, and Feifang Hu. 2022\.  Adaptive A/B Test on Networks with Cluster Structures. In _International Conference on Artificial Intelligence and Statistics_. PMLR, 10836–10851.
  * Luke and Harris (2007) Douglas A Luke and Jenine K Harris. 2007.  Network analysis in public health: history, methods, and applications.  _Annu. Rev. Public Health_ 28 (2007), 69–93.
  * Manski (2013) Charles F Manski. 2013\.  Identification of treatment response with social interactions.  _The Econometrics Journal_ 16, 1 (2013), S1–S23.
  * Portes and Rey (2005) Richard Portes and Helene Rey. 2005.  The determinants of cross-border equity flows.  _Journal of international Economics_ 65, 2 (2005), 269–296.
  * Pouget-Abadie et al. (2019) Jean Pouget-Abadie, Guillaume Saint-Jacques, Martin Saveski, Weitao Duan, Souvik Ghosh, Ya Xu, and Edoardo M Airoldi. 2019\.  Testing for arbitrary interference on experimentation platforms.  _Biometrika_ 106, 4 (2019), 929–940.
  * Saint-Jacques et al. (2019a) Guillaume Saint-Jacques, James Eric Sorenson, Nanyu Chen, and Ya Xu. 2019a.  A method for measuring network effects of one-to-one communication features in online a/b tests.  _arXiv preprint arXiv:1903.08766_ (2019).
  * Saint-Jacques et al. (2019b) Guillaume Saint-Jacques, Maneesh Varshney, Jeremy Simpson, and Ya Xu. 2019b.  Using ego-clusters to measure network effects at LinkedIn.  _arXiv preprint arXiv:1903.08755_ (2019).
  * Särndal et al. (2003) Carl-Erik Särndal, Bengt Swensson, and Jan Wretman. 2003\.  _Model assisted survey sampling_.  Springer Science & Business Media.
  * Sävje (2023) Fredrik Sävje. 2023\.  Causal inference with misspecified exposure mappings: separating definitions and assumptions.  _Biometrika_ (2023), asad019.
  * Sävje et al. (2021) Fredrik Sävje, Peter Aronow, and Michael Hudgens. 2021\.  Average treatment effects in the presence of unknown interference.  _Annals of statistics_ 49, 2 (2021), 673.
  * Sussman and Airoldi (2017) Daniel L Sussman and Edoardo M Airoldi. 2017.  Elements of estimation theory for causal effects in the presence of network interference.  _arXiv preprint arXiv:1702.03578_ (2017).
  * Tang et al. (2010) Diane Tang, Ashish Agarwal, Deirdre O’Brien, and Mike Meyer. 2010.  Overlapping experiment infrastructure: More, better, faster experimentation. In _Proceedings of the 16th ACM SIGKDD international conference on Knowledge discovery and data mining_. 17–26.
  * Toulis and Kao (2013) Panos Toulis and Edward Kao. 2013.  Estimation of causal peer influence effects. In _International conference on machine learning_. PMLR, 1489–1497.
  * Ugander et al. (2013) Johan Ugander, Brian Karrer, Lars Backstrom, and Jon Kleinberg. 2013.  Graph cluster randomization: Network exposure to multiple universes. In _Proceedings of the 19th ACM SIGKDD international conference on Knowledge discovery and data mining_. 329–337.
  * Ugander and Yin (2023) Johan Ugander and Hao Yin. 2023.  Randomized graph cluster randomization.  _Journal of Causal Inference_ 11, 1 (2023), 20220014\.
  * Xiong et al. (2020) Tao Xiong, Yong Wang, and Senlie Zheng. 2020.  _Orthogonal Traffic Assignment in Online Overlapping A/B Tests_.  Technical Report. EasyChair.
  * Yu et al. (2022) Christina Lee Yu, Edoardo M Airoldi, Christian Borgs, and Jennifer T Chayes. 2022.  Estimating the total treatment effect in randomized experiments with unknown network structure.  _Proceedings of the National Academy of Sciences_ 119, 44 (2022), e2208975119.
  * Yuan and Altenburger (2022) Yuan Yuan and Kristen M Altenburger. 2022.  A Two-Part Machine Learning Approach to Characterizing Network Interference in A/B Testing.

##  Appendix A Proofs

###  A.1. Proof of Proposition 3.1

###### Proof.

Let’s rewrite the estimator $\widehat{\tau}^{1}$ as

$\displaystyle\widehat{\tau}^{1}=$
$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\frac{W_{i}Y_{i}}{\pi}-\frac{1}{n}\sum_{i=1}^{n}\frac{(1-W_{i})Y_{i}}{1-\pi}$
$\displaystyle=$
$\displaystyle\frac{1}{n\pi}\sum_{i=1}^{n}W_{i}\left[\sum_{j=1}^{n}\Big{(}\alpha_{j,i}+\beta_{j,i}W_{j}+\gamma_{j,i}W_{i}+\zeta_{j,i}W_{i}W_{j}\Big{)}\right]-$
$\displaystyle\,\,\frac{1}{n(1-\pi)}\sum_{i=1}^{n}(1-W_{i})\left[\sum_{j=1}^{n}\Big{(}\alpha_{j,i}+\beta_{j,i}W_{j}+\gamma_{j,i}W_{i}+\zeta_{j,i}W_{i}W_{j}\Big{)}\right]$
$\displaystyle=$ $\displaystyle\frac{1}{n\pi}\sum_{i\neq
j}W_{i}\Big{(}\alpha_{j,i}+\beta_{j,i}W_{j}+\gamma_{j,i}+\zeta_{j,i}W_{j}\Big{)}-$
$\displaystyle\,\,\frac{1}{n(1-\pi)}\sum_{i\neq
j}(1-W_{i})\Big{(}\alpha_{j,i}+\beta_{j,i}W_{j}\Big{)}$

Note that

  1. (1)

$\mathbb{E}(W_{i})=\pi$.

  2. (2)

$\\{W_{i}\\}_{i=1}^{n}$ are independent, $\mathbb{E}(W_{i}W_{j})=\pi^{2}$ for
$i\neq j$.

So we have

$\displaystyle\mathbb{E}(\widehat{\tau}^{1})=$
$\displaystyle\frac{1}{n}\sum_{i\neq
j}\Big{(}\alpha_{j,i}+\beta_{j,i}\pi+\gamma_{j,i}+\zeta_{j,i}\pi\Big{)}-\frac{1}{n}\sum_{i\neq
j}\Big{(}\alpha_{j,i}+\beta_{j,i}\pi\Big{)}$ $\displaystyle=$
$\displaystyle\frac{1}{n}\sum_{i\neq
j}\Big{(}\gamma_{i,j}+\zeta_{i,j}\pi\Big{)}.$

∎

###  A.2. Proof of Proposition 3.2

###### Proof.

Follow the same steps,

$\displaystyle\widehat{\tau}^{2}=$
$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\frac{W_{i}D_{i}}{\pi}-\frac{1}{n}\sum_{i=1}^{n}\frac{(1-W_{i})D_{i}}{1-\pi}$
$\displaystyle=$
$\displaystyle\frac{1}{n\pi}\sum_{i=1}^{n}W_{i}\left[\sum_{j=1}^{n}\Big{(}\alpha_{i,j}+\beta_{i,j}W_{i}+\gamma_{i,j}W_{j}+\zeta_{i,j}W_{i}W_{j}\Big{)}\right]-$
$\displaystyle\,\,\frac{1}{n(1-\pi)}\sum_{i=1}^{n}(1-W_{i})\left[\sum_{j=1}^{n}\Big{(}\alpha_{i,j}+\beta_{i,j}W_{i}+\gamma_{i,j}W_{j}+\zeta_{i,j}W_{i}W_{j}\Big{)}\right]$
$\displaystyle=$ $\displaystyle\frac{1}{n\pi}\sum_{i\neq
j}^{n}W_{i}\Big{(}\alpha_{i,j}+\beta_{i,j}+\gamma_{i,j}W_{j}+\zeta_{i,j}W_{j}\Big{)}-$
$\displaystyle\,\,\frac{1}{n(1-\pi)}\sum_{i\neq
j}(1-W_{i})\Big{(}\alpha_{i,j}+\gamma_{i,j}W_{j}\Big{)}.$

Again, $\\{W_{i}\\}_{i=1}^{n}$ are independent. Hence

$\displaystyle\mathbb{E}(\widehat{\tau}^{2})=$
$\displaystyle\frac{1}{n}\sum_{i\neq
j}\Big{(}\alpha_{i,j}+\beta_{i,j}+\gamma_{i,j}\pi+\zeta_{i,j}\pi\Big{)}-\frac{1}{n}\sum_{i\neq
j}\Big{(}\alpha_{i,j}+\gamma_{i,j}\pi\Big{)}$ $\displaystyle=$
$\displaystyle\frac{1}{n}\sum_{i\neq
j}\Big{(}\beta_{i,j}+\zeta_{i,j}\pi\Big{)}.$

∎

###  A.3. Proof of Proposition 3.3

###### Proof.

Follow the same steps,

$\displaystyle\widehat{\tau}^{1}=$
$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\frac{V_{i}W_{i}Y_{i}}{p\pi}-\frac{1}{n}\sum_{i=1}^{n}\frac{V_{i}(1-W_{i})Y_{i}}{p(1-\pi)}$
$\displaystyle=$
$\displaystyle\frac{1}{np\pi}\sum_{i=1}^{n}V_{i}W_{i}\left[\sum_{j=1}^{n}\Big{(}\alpha_{j,i}+\beta_{j,i}W_{j}+\gamma_{j,i}W_{i}+\zeta_{j,i}W_{i}W_{j}\Big{)}\right]-$
$\displaystyle\,\,\frac{1}{np(1-\pi)}\sum_{i=1}^{n}V_{i}(1-W_{i})\left[\sum_{j=1}^{n}\Big{(}\alpha_{j,i}+\beta_{j,i}W_{j}+\gamma_{j,i}W_{i}+\zeta_{j,i}W_{i}W_{j}\Big{)}\right]$
$\displaystyle=$ $\displaystyle\frac{1}{np\pi}\sum_{i\neq
j}^{n}V_{i}W_{i}\Big{(}\alpha_{j,i}+\beta_{j,i}W_{j}+\gamma_{j,i}+\zeta_{j,i}W_{j}\Big{)}-$
$\displaystyle\,\,\frac{1}{np(1-\pi)}\sum_{i\neq
j}V_{i}(1-W_{i})\Big{(}\alpha_{j,i}+\beta_{j,i}W_{j}\Big{)}$

Note that

  1. (1)

$\mathbb{E}(V_{i}W_{i})=p\pi$.

  2. (2)

$\mathbb{E}(V_{i}W_{i}W_{j})=p^{2}\pi$.

  3. (3)

$\mathbb{E}\big{(}V_{i}(1-W_{i})\big{)}=p(1-\pi)$.

  4. (4)

$\mathbb{E}\big{(}V_{i}(1-W_{i})W_{j}\big{)}=p^{2}\pi(1-\pi)$.

Hence

$\displaystyle\mathbb{E}(\widehat{\tau}^{1})=$
$\displaystyle\frac{1}{n}\sum_{i\neq
j}\Big{(}\alpha_{j,i}+\beta_{j,i}p\pi+\gamma_{j,i}+\zeta_{j,i}p\pi\Big{)}-\frac{1}{n}\sum_{i\neq
j}\Big{(}\alpha_{j,i}+\beta_{j,i}p\pi\Big{)}$ $\displaystyle=$
$\displaystyle\frac{1}{n}\sum_{i\neq
j}\Big{(}\gamma_{j,i}+p\pi\zeta_{j,i}\Big{)}.$

∎

###  A.4. Proof of Proposition 3.4

###### Proof.

Follow the same steps,

$\displaystyle\widehat{\tau}^{2}=$
$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\frac{V_{i}W_{i}D_{i}}{p\pi}-\frac{1}{n}\sum_{i=1}^{n}\frac{V_{i}(1-W_{i})D_{i}}{p(1-\pi)}$
$\displaystyle=$
$\displaystyle\frac{1}{np\pi}\sum_{i=1}^{n}V_{i}W_{i}\left[\sum_{j=1}^{n}\Big{(}\alpha_{i,j}+\beta_{i,j}W_{i}+\gamma_{i,j}W_{j}+\zeta_{i,j}W_{i}W_{j}\Big{)}\right]-$
$\displaystyle\,\,\frac{1}{np(1-\pi)}\sum_{i=1}^{n}V_{i}(1-W_{i})\left[\sum_{j=1}^{n}\Big{(}\alpha_{i,j}+\beta_{i,j}W_{i}+\gamma_{i,j}W_{j}+\zeta_{i,j}W_{i}W_{j}\Big{)}\right]$
$\displaystyle=$ $\displaystyle\frac{1}{np\pi}\sum_{i\neq
j}V_{i}W_{i}\Big{(}\alpha_{i,j}+\beta_{i,j}+\gamma_{i,j}W_{j}+\zeta_{i,j}W_{j}\Big{)}-$
$\displaystyle\,\ \frac{1}{np(1-\pi)}\sum_{i\neq
j}V_{i}(1-W_{j})\Big{(}\alpha_{i,j}+\gamma_{i,j}W_{j}\Big{)}$

Hence

$\displaystyle\mathbb{E}(\widehat{\tau}^{2})=$
$\displaystyle\frac{1}{n}\sum_{i\neq
j}\Big{(}\alpha_{i,j}+\beta_{i,j}+\gamma_{i,j}p\pi+\zeta_{i,j}p\pi\Big{)}-\frac{1}{n}\sum_{i\neq
j}\Big{(}\alpha_{i,j}+\gamma_{i,j}p\pi\Big{)}$ $\displaystyle=$
$\displaystyle\frac{1}{n}\sum_{i\neq
j}\Big{(}\beta_{i,j}+p\pi\zeta_{i,j}\Big{)}$

∎

###  A.5. Proof of Proposition 3.5

###### Proof.

Follow the same steps,

$\displaystyle\widehat{\tau}^{1}=$
$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\frac{V_{i}W_{i}Y_{i}}{p\pi}-\frac{1}{n}\sum_{i=1}^{n}\frac{V_{i}(1-W_{i})Y_{i}}{p(1-\pi)}$
$\displaystyle=$
$\displaystyle\frac{1}{np\pi}\sum_{i=1}^{n}V_{i}W_{i}\left[\sum_{j=1}^{n}\Big{(}\alpha_{j,i}+\beta_{j,i}W_{j}+\gamma_{j,i}W_{i}+\zeta_{j,i}W_{i}W_{j}\Big{)}\right]-$
$\displaystyle\,\,\frac{1}{np(1-\pi)}\sum_{i=1}^{n}V_{i}(1-W_{i})\left[\sum_{j=1}^{n}\Big{(}\alpha_{j,i}+\beta_{j,i}W_{j}+\gamma_{j,i}W_{i}+\zeta_{j,i}W_{i}W_{j}\Big{)}\right]$
$\displaystyle=$ $\displaystyle\frac{1}{np\pi}\sum_{i\neq
j}V_{i}W_{i}\Big{(}\alpha_{j,i}+\beta_{j,i}W_{j}+\gamma_{j,i}+\zeta_{j,i}W_{j}\Big{)}-$
$\displaystyle\,\,\frac{1}{np(1-\pi)}\sum_{i\neq
j}V_{i}(1-W_{i})\Big{(}\alpha_{j,i}+\beta_{j,i}W_{j}\Big{)}$

Note that

  1. (1)

$\mathbb{E}(V_{i}W_{i})=p\pi,\mathbb{E}\big{(}V_{i}(1-W_{i})\big{)}=p(1-\pi).$

  2. (2)

$\mathbb{E}\big{(}V_{i}W_{i}W_{j}\big{)}=\mathbb{E}[\mathbb{E}\big{(}V_{i}W_{i}W_{j}\big{|}j\text{
and }i\text{ are in the same
cluster})]=\bar{\sigma}*\mathbb{E}\big{(}V_{i}W_{i}\big{)}+(1-\bar{\sigma})*\mathbb{E}\big{(}V_{i}W_{i}\big{)}\mathbb{E}\big{(}W_{j}\big{)}=\bar{\sigma}p\pi+(1-\bar{\sigma})p^{2}\pi^{2}.$

  3. (3)

$\mathbb{E}\big{(}V_{i}(1-W_{i})W_{j}\big{)}=\mathbb{E}[\mathbb{E}\big{(}V_{i}(1-W_{i})W_{j}\big{|}j\text{
and }i$ are not in the same cluster $)]=(1-\bar{\sigma})p^{2}\pi(1-\pi)$.

Hence

$\displaystyle\mathbb{E}(\widehat{\tau}^{1})=$
$\displaystyle\frac{1}{n}\sum_{i\neq
j}\Big{(}\alpha_{j,i}+\gamma_{j,i}+(\beta_{j,i}+\zeta_{j,i})(\bar{\sigma}+(1-\bar{\sigma})p\pi)\Big{)}-$
$\displaystyle\,\,\frac{1}{n}\sum_{i\neq
j}\Big{(}\alpha_{j,i}+\beta_{j,i}(1-\bar{\sigma})p\pi\Big{)}$ $\displaystyle=$
$\displaystyle\frac{1}{n}\sum_{i\neq
j}\Big{(}\beta_{i,j}\bar{\sigma}+\gamma_{i,j}+\big{(}\bar{\sigma}+(1-\bar{\sigma})p\pi\big{)}\zeta_{i,j}\Big{)}$

∎

###  A.6. Proof of Proposition 3.6

###### Proof.

Follow the same steps in (A.5),

$\displaystyle\widehat{\tau}^{2}=$
$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\frac{V_{i}W_{i}D_{i}}{p\pi}-\frac{1}{n}\sum_{i=1}^{n}\frac{V_{i}(1-W_{i})D_{i}}{p(1-\pi)}$
$\displaystyle=$
$\displaystyle\frac{1}{np\pi}\sum_{i=1}^{n}V_{i}W_{i}\left[\sum_{j=1}^{n}\Big{(}\alpha_{i,j}+\beta_{i,j}W_{i}+\gamma_{i,j}W_{j}+\zeta_{i,j}W_{i}W_{j}\Big{)}\right]-$
$\displaystyle\,\,\frac{1}{np(1-\pi)}\sum_{i=1}^{n}V_{i}(1-W_{i})\left[\sum_{j=1}^{n}\Big{(}\alpha_{i,j}+\beta_{i,j}W_{i}+\gamma_{i,j}W_{j}+\zeta_{i,j}W_{i}W_{j}\Big{)}\right]$
$\displaystyle=$ $\displaystyle\frac{1}{np\pi}\sum_{i\neq
j}V_{i}W_{i}\Big{(}\alpha_{i,j}+\beta_{i,j}+\gamma_{i,j}W_{j}+\zeta_{i,j}W_{j}\Big{)}-$
$\displaystyle\,\,\frac{1}{np(1-\pi)}\sum_{i\neq
j}V_{i}(1-W_{i})\Big{(}\alpha_{i,j}+\gamma_{i,j}W_{j}\Big{)}$

Hence

$\displaystyle\mathbb{E}(\widehat{\tau}^{2})=$
$\displaystyle\frac{1}{n}\sum_{i\neq
j}\Big{(}\alpha_{i,j}+\beta_{j,i}+(\gamma_{j,i}+\zeta_{i,j})(\bar{\sigma}+(1-\bar{\sigma})p\pi)\Big{)}-$
$\displaystyle\,\,\frac{1}{n}\sum_{i\neq
j}\Big{(}\alpha_{i,j}+\gamma_{i,j}(1-\bar{\sigma})p\pi\Big{)}$
$\displaystyle=$ $\displaystyle\frac{1}{n}\sum_{i\neq
j}\Big{(}\beta_{i,j}+\gamma_{i,j}\bar{\sigma}+\big{(}\bar{\sigma}+(1-\bar{\sigma})p\pi\big{)}\zeta_{i,j}\Big{)}$

∎


# Testing the Inverted Neutrino Mass Ordering with Neutrinoless Double-Beta
Decay

Matteo Agostini<EMAIL_ADDRESS>Department of Physics and
Astronomy, University College London, Gower Street, London WC1E 6BT, UK
Giovanni Benato<EMAIL_ADDRESS>INFN, Laboratori Nazionali del
Gran Sasso, 67100 Assergi, L’Aquila, Italy  Jason A. Detwiler<EMAIL_ADDRESS>Center for Experimental Nuclear Physics and Astrophysics, and Department of
Physics, University of Washington, Seattle, WA 98115 - USA  Javier Menéndez
<EMAIL_ADDRESS>Department of Quantum Physics and Astrophysics and
Institute of Cosmos Sciences, University of Barcelona, 08028 Barcelona, Spain
Francesco Vissani<EMAIL_ADDRESS>INFN, Laboratori Nazionali del Gran
Sasso, 67100 Assergi, L’Aquila, Italy  Gran Sasso Science Institute, 67100
L’Aquila, Italy

###### Abstract

We quantify the extent to which future experiments will test the existence of
neutrinoless double-beta decay mediated by light neutrinos with inverted-
ordered masses. While it remains difficult to compare measurements performed
with different isotopes, we find that future searches will fully test the
inverted ordering scenario, as a global, multi-isotope endeavor. They will
also test other possible mechanisms driving the decay, including a large
uncharted region of the allowed parameter space assuming that neutrino masses
follow the normal ordering.

Neutrino oscillation [1, 2, 3, 4] proves that the neutrino has mass and the
Standard Model of particle physics is incomplete. The unknown origin of the
neutrino mass has drawn enormous attention to neutrinoless double-beta
($0\nu\beta\beta$) decay, a matter-creating nuclear transition in which two
neutrons decay simultaneously into two protons, emitting only two new
electrons and no antineutrinos [5]. The discovery of $0\nu\beta\beta$ decay
would establish that the neutrino is its own antiparticle and has a Majorana
mass [6]. It would also mark the first observation of a lepton-creating
process, proving that neither lepton number ($L$) nor baryon minus lepton
number ($B-L$) are symmetries of the Standard Model, as predicted by leading
theories explaining the matter-antimatter asymmetry of our universe [7].
Indeed, searching for $0\nu\beta\beta$ decay is the most sensitive
experimental approach to test Majorana neutrino masses and their associated
$L$ violation. It is also a unique probe of new physics at ultrahigh energy
scales not accessible by current accelerators [8, 9].

Different physics mechanisms can lead to $0\nu\beta\beta$ decay [10]. However,
the exchange of light Majorana neutrinos interacting via standard, weak left-
handed currents plays a special role. It is the only mechanism allowed by all
theories predicting $0\nu\beta\beta$ decay, and typically dominates the rate
of the process [11, 12]. Assuming that the decay is mediated by light
neutrinos, its half-life is [13]

$T_{1/2}^{-1}=G~{}g_{A}^{4}~{}M^{2}~{}\frac{m_{\beta\beta}^{2}}{m_{e}^{2}},$
(1)

where $G$ is the phase-space integral, $g_{A}\simeq 1.276$ [14] is the axial-
vector coupling, $M$ is the nuclear matrix element, and $m_{\beta\beta}$ is
the effective Majorana mass, normalized for convenience by the electron mass
$m_{e}$. The Majorana mass captures the physics of the exchanged neutrinos and
is a function of the neutrino oscillation parameters, the neutrino mass
eigenvalues $m_{i}$, and the Majorana phases:
$m_{\beta\beta}=\bigl{|}\sum_{i}U^{2}_{ei}~{}m_{i}\bigr{|}$, where $U_{ei}$
are the elements of the full 6-parameter PMNS matrix [15].

Neutrino oscillation measurements constrain the range of allowed
$m_{\beta\beta}$ values [16] and prove that
$|U_{e3}^{2}|\ll|U_{e2}^{2}|<|U_{e1}^{2}|$ and
$m_{2}^{2}-m_{1}^{2}\ll|m_{3}^{2}-m_{2}^{2}|$ [17, 18, 19, 20, 21, 22, 23, 24,
25, 26]. This implies that the effective Majorana mass is strictly larger than
zero if neutrino masses follow the inverted ordering, i.e.
$m_{3}<m_{1}\lesssim m_{2}$. In this case, the lowest $m_{\beta\beta}$ value,
minimized with respect to the unknown Majorana phases and $m_{3}$, is given by

$(m_{\beta\beta}^{min})_{\text{IO}}=\bigl{|}U_{e1}^{2}\bigr{|}m_{1}-\bigl{|}U_{e2}^{2}\bigr{|}m_{2}-\bigl{|}U_{e3}^{2}\bigr{|}m_{3},$
(2)

with
$m_{3}=\bigl{|}U_{e3}^{2}\bigr{|}/\left(\bigl{|}U_{e1}^{2}\bigr{|}/m_{1}-\bigl{|}U_{e2}^{2}\bigr{|}/m_{2}\right)\approx
3$ meV. Using the latest values and uncertainties from the Particle Data Group
[4], we obtain

$(m_{\beta\beta}^{min})_{\text{IO}}=18.4\pm 1.3\,\text{meV},$ (3)

whose uncertainty is dominated by the uncertainty on the solar mixing angle
$\theta_{12}$. Using the latest NuFIT results [27], we obtain
$(m_{\beta\beta}^{min})_{\text{IO}}=18.6\pm 1.2$ meV. The lower bound on
$m_{\beta\beta}$ corresponds to an upper bound on $T_{1/2}$ at the scale of
$10^{27}$-$10^{28}$ years, depending on the value of the parameters in
equation (1).

In the last decades, a vast experimental program has been mounted to develop
experiments with sensitivity reaching $(m_{\beta\beta}^{min})_{\text{IO}}$,
able to exhaustively test whether $0\nu\beta\beta$ decay is mediated by the
exchange of light neutrinos with inverted-ordered masses [28, 29]. Mature
designs are now available for multiple detection techniques, and the physics
community is discussing how to proceed. As part of this process, the
Astroparticle Physics European Consortium (APPEC) is updating its
$0\nu\beta\beta$ decay roadmap [29] and the United States’ Department of
Energy has started a ton-scale-experiment portfolio review. Conceptual designs
are available for three experiments [30, 31, 32], whose construction can start
as soon as funding is available. These experiments use different
$0\nu\beta\beta$-decaying isotopes and detection technologies, and can perform
independent and complementary measurements.

As mentioned above, observing $0\nu\beta\beta$ decay would unambiguously
demonstrate matter creation and prove the Majorana nature of neutrinos.
However, the conversion between $T_{1/2}$ and $m_{\beta\beta}$ in Eq. (1) is
not trivial and requires inputs from nuclear theory. While the phase-space
integral $G$ has been calculated with negligible uncertainty [33], obtaining
reliable nuclear matrix elements $M$ is challenging, as it requires
computationally intensive many-body calculations and the evaluation of several
operators [13, 34]. Four primary many-body methods have been historically used
in the field: the nuclear shell model (NSM) [35, 36, 37], the quasiparticle
random-phase approximation (QRPA) method [38, 39, 40, 41, 42], energy-density
functional (EDF) theory [43, 44, 45], and the interacting boson model (IBM)
[46, 47]. For each of these methods, several calculations have been performed
under different assumptions and approximations. The most recent results are
listed in Table 1. They can differ by up to a factor of three for a given
isotope, and significant differences are present even within each method. The
spread of values gives a rough idea of the many-body uncertainties on $M$
(additional uncertainty contributions are discussed below). For some methods,
calculations are not available for all isotopes.

Table 1: Nuclear matrix elements $M$ for $0\nu\beta\beta$ decay mediated by light neutrinos, calculated with the NSM, QRPA, EDF, and IBM methods. The ranges correspond to the minimum and maximum values obtained with the same many-body method.  |  Ref. |  76Ge |  100Mo |  136Xe
---|---|---|---|---
NSM |  [35] |  2.89, 3.07 |  – |  2.28, 2.45
[36] |  3.37, 3.57 |  – |  1.63, 1.76
[37] |  2.66 |  – |  2.39
All |  2.66 - 3.57 |  – |  1.63 - 2.45
QRPA |  [38] |  5.09 |  – |  1.55
[39] |  5.26 |  3.90 |  2.91
[40] |  4.85 |  5.87 |  2.72
[41] |  3.12, 3.40 |  – |  1.11, 1.18
[42] |  – |  – |  3.38
All |  3.12 - 5.26 |  3.90 - 5.87 |  1.11 - 3.38
EDF |  [43] |  4.60 |  5.08 |  4.20
[44] |  5.55 |  6.59 |  4.77
[45] |  6.04 |  6.48 |  4.24
All |  4.60 - 6.04 |  5.08 - 6.59 |  4.20 - 4.77
IBM |  [46]111With the sign change in the tensor part indicated in Ref. [47]. |  5.14 |  3.84 |  3.25
[47] |  6.34 |  5.08 |  3.40
All |  5.14 - 6.34 |  3.84 - 5.08 |  3.25 - 3.40

The reach of $0\nu\beta\beta$ decay experiments is conventionally expressed in
terms of discovery and exclusion sensitivities on $m_{\beta\beta}$. The
discovery sensitivity corresponds to the smallest $m_{\beta\beta}$ value for
which an experiment has 50% probability of observing a signal at 99.7%
confidence level (CL). The exclusion sensitivity corresponds to the median
90%-CL upper limit that an experiment will set on $m_{\beta\beta}$ assuming
$0\nu\beta\beta$ decay is not observable. As stated earlier, fully testing the
inverted ordering scenario requires sensitivity to
$(m_{\beta\beta}^{min})_{\text{IO}}$, accounting for its uncertainty. For
discovery mode, this condition is met when the discovery sensitivity reaches
the central value of $(m_{\beta\beta}^{min})_{\text{IO}}$: the
$(m_{\beta\beta}^{min})_{\text{IO}}$ uncertainty is symmetric, so the
probability of lower or upper fluctuations is the same, and the 50%
probability for an observation is preserved. However, for exclusion mode, the
$(m_{\beta\beta}^{min})_{\text{IO}}$ uncertainty reduces the CL by a variable
amount that depends on the experimental parameters, mainly the background
statistical uncertainty. Therefore the exclusion sensitivity on
$m_{\beta\beta}$ cannot be used to set an experiment-independent condition
corresponding to fully covering the inverted ordering scenario. The discovery
sensitivity is anyway the most appropriate metric for searches aiming to
discover a process. Thus, reaching a discovery sensitivity of 18.4 meV is the
right concrete goal for experiments aiming to explore the full inverted
ordering parameter space.

Figure 1: Comparison of $m_{\beta\beta}$ 99.7%-CL discovery and 90%-CL median
exclusion sensitivities for different isotopes at stated half-life
sensitivities [31, 30, 48], grouped by nuclear many-body frameworks with
matrix element ranges from Table 1. The horizontal bands show the variation on
$(m_{\beta\beta}^{min})_{\text{IO}}$ under variation of the neutrino
oscillation parameters.

Figure 1 shows the discovery and exclusion sensitivities of proposed future
experiments. We converted the $T_{1/2}$ sensitivity values quoted by the
LEGEND [31], CUPID [30] and nEXO [48] collaborations to $m_{\beta\beta}$
values using the nuclear matrix elements of Table 1. We group the calculations
by many-body method to aid in the comparison of more consistent quantities.
Remarkably, these experiments all show sensitivity to measure a
$0\nu\beta\beta$ decay signal at the bottom of the inverted ordering parameter
space. Some many-body methods, such as EDF theory and IBM, give systematically
larger $M$ values, pushing the $m_{\beta\beta}$ sensitivity even lower. QRPA
calculations on the other hand give a broad range of results, partially due to
the role of nuclear deformation in this framework. The NSM provides $M$ values
which are typically smaller than for the other methods, but are not available
for 100Mo.

Comparing the performance of experiments using different isotopes is
challenging because of the large uncertainties affecting the nuclear matrix
element predictions. In particular, each many-body method uses different
approximations, which are likely to result in a common over- or
underestimation of the calculated $M$ values. Even comparisons considering one
many-body method at a time can raise concerns. The number of calculations
available for each isotope and method can be significantly different,
suggesting that not all many-body approaches are equally suitable for all
isotopes. Consequently the range of $M$ values cannot be quantitatively
interpreted as the uncertainty, which currently remains unknown. One might be
tempted to make comparisons based on the central value of the $m_{\beta\beta}$
sensitivity within a specific nuclear method, but this value can be
disproportionately affected by a single outlier matrix element. Weighted
averages are also problematic, as the weight given to each calculation would
be to some extent arbitrary. Given the lack of objective criteria to compare
experimental sensitivities in different isotopes, and the lack of a clear
estimate of the uncertainties, we advocate to refrain from ranking
experiments’ reach quantitatively, and focus instead on the fact that we have
a global, multi-isotope endeavor that will fully test the inverted ordering
scenario.

A broad effort to reduce uncertainties is ongoing within the nuclear theory
community, with significant advances made in the last few years. Ab initio
calculations that incorporate wider nuclear correlations and two-body currents
have recently succeeded in predicting single $\beta$ decay rates [49] with no
need for “quenching” – an ad hoc reduction of the value of calculated matrix
elements involving the nuclear spin required by less sophisticated
calculations [13]. The first available ab initio $0\nu\beta\beta$ matrix
element calculations in medium-sized nuclei, supported by studies in lighter
systems [50, 51], indicate a relatively mild suppression by tens of percent
with respect to the lower limit of the range given in Table 1 [52]. Efforts
are underway to improve the quality of these results, extend them to heavier
nuclei, and include two-body currents at finite momentum transfers [53]. On
the other hand, the contact term introduced in Refs. [34, 54], which until
recently went unrecognized, is a leading-order contribution to $M$. Effective
field theory and ab initio nuclear structure provide a scheme for estimating
this contribution [55, 56]. A first study in 48Ca suggests that this term can
enhance $M$ by about 40% percent [57], leading to a faster decay rate. In
heavier systems, preliminary results suggest a roughly similar impact for all
$0\nu\beta\beta$ isotopes, only slightly dependent on the nuclear many-body
method [58]. Complementary studies using e.g. lattice QCD [59, 60] will test
whether this claimed enhancement is robust. If so, it may compensate the
reduction in decay rate due to the inclusion of the “quenching” physics,
leading to a picture similar to the one represented by Fig. 1. Thus, should
the current theoretical results be confirmed, the proposed global
$0\nu\beta\beta$ decay experimental effort would still fully probe the
inverted ordering parameter space.

In this letter, we have focused on the inverted ordering scenario as a
prominent goalpost for the proposed experimental $0\nu\beta\beta$ decay
program. However, the discovery power of these experiments is high even
assuming other, equally reasonable scenarios. By reaching a sensitivity of the
order of tens of meV, these searches will probe a significant fraction of the
remaining parameter space for left-handed neutrino exchange even if neutrino
masses follow the normal ordering. Bayesian analyses suggest up to 50%
discovery probabilities for the normal ordering scenario [61, 62], and a non-
vanishing discovery probability even assuming the most unfavorable value of
the Majorana phases [63]. Significant advancement would also be made in
probing the exchange of heavy mediators. For many such models,
$0\nu\beta\beta$ decay searches probe energy scales beyond the reach of
current accelerator technology [9]. Additional physics mechanisms could
completely change the parameter space of interest, potentially even increasing
the discovery power of future experiments [64, 65, 66, 67]. In general,
pushing $0\nu\beta\beta$ decay sensitivity to increasingly large half-life
values explores uncharted territory, and new physics could manifest at any
time.

###### Acknowledgements.

We would like to thank Christoph Wiesinger and Steven R. Elliott for valuable
discussions. This work has been supported by the Science and Technology
Facilities Council (Grant No. ST/T004169/1), by the EU Horizon2020 research
and innovation program under the Marie Skłodowska-Curie Grant Agreement No.
754496, by the Spanish MICINN through the “Ramón y Cajal” program with grant
RYC-2017-22781, the Italian Research Grant Number 2017W4HA7S “NAT-NET:
Neutrino and Astroparticle Theory Network” under the program PRIN 2017 funded
by MIUR, the AEI “Unit of Excellence María de Maeztu 2020-2023” award
CEX2019-000918-M, the AEI grant FIS2017-87534-P and by the U.S. DOE Office of
Nuclear Physics under Grant Number DE-FG02-97ER41020.

## References

  * Kajita [2016] T. Kajita, Nobel Lecture: Discovery of atmospheric neutrino oscillations, Rev. Mod. Phys. 88, 030501 (2016).
  * McDonald [2016] A. B. McDonald, Nobel Lecture: The Sudbury Neutrino Observatory: Observation of flavor change for solar neutrinos, Rev. Mod. Phys. 88, 030502 (2016).
  * Eguchi  _et al._ [2003] K. Eguchi _et al._ (KamLAND), First results from KamLAND: Evidence for reactor anti-neutrino disappearance, Phys. Rev. Lett. 90, 021802 (2003), arXiv:hep-ex/0212021 .
  * Zyla  _et al._ [2020] P. A. Zyla _et al._ (Particle Data Group), Review of Particle Physics, PTEP 2020, 083C01 (2020), and 2021 update.
  * Furry [1939] W. H. Furry, On transition probabilities in double beta-disintegration, Phys. Rev. 56, 1184 (1939).
  * Schechter and Valle [1982] J. Schechter and J. Valle, Neutrinoless Double beta Decay in SU(2) x U(1) Theories, Phys.Rev. D25, 2951 (1982).
  * Fukugita and Yanagida [1986] M. Fukugita and T. Yanagida, Baryogenesis Without Grand Unification, Phys. Lett. B 174, 45 (1986).
  * Deppisch  _et al._ [2015] F. F. Deppisch, P. S. Bhupal Dev, and A. Pilaftsis, Neutrinos and Collider Physics, New J. Phys. 17, 075019 (2015), arXiv:1502.06541 [hep-ph] .
  * Peng  _et al._ [2016] T. Peng, M. J. Ramsey-Musolf, and P. Winslow, TeV lepton number violation: From neutrinoless double-$\beta$ decay to the LHC, Phys. Rev. D 93, 093002 (2016), arXiv:1508.04444 [hep-ph] .
  * Vergados  _et al._ [2012] J. D. Vergados, H. Ejiri, and F. Simkovic, Theory of Neutrinoless Double Beta Decay, Rept. Prog. Phys. 75, 106301 (2012), arXiv:1205.0649 [hep-ph] .
  * de Gouvea and Jenkins [2008] A. de Gouvea and J. Jenkins, A Survey of Lepton Number Violation Via Effective Operators, Phys. Rev. D 77, 013008 (2008), arXiv:0708.1344 [hep-ph] .
  * Mitra  _et al._ [2012] M. Mitra, G. Senjanovic, and F. Vissani, Neutrinoless Double Beta Decay and Heavy Sterile Neutrinos, Nucl. Phys. B 856, 26 (2012), arXiv:1108.0004 [hep-ph] .
  * Engel and Menéndez [2017] J. Engel and J. Menéndez, Status and Future of Nuclear Matrix Elements for Neutrinoless Double-Beta Decay: A Review, Rept. Prog. Phys. 80, 046301 (2017), arXiv:1610.06548 [nucl-th] .
  * Märkisch  _et al._ [2019] B. Märkisch _et al._ , Measurement of the Weak Axial-Vector Coupling Constant in the Decay of Free Neutrons Using a Pulsed Cold Neutron Beam, Phys. Rev. Lett. 122, 242501 (2019), arXiv:1812.04666 [nucl-ex] .
  * Bilenky and Giunti [2015] S. M. Bilenky and C. Giunti, Neutrinoless Double-Beta Decay: a Probe of Physics Beyond the Standard Model, Int. J. Mod. Phys. A 30, 1530001 (2015), arXiv:1411.4791 [hep-ph] .
  * Vissani [1999] F. Vissani, Signal of neutrinoless double beta decay, neutrino spectrum and oscillation scenarios, J. High Energy Phys. 06, 022 (1999).
  * Abe  _et al._ [2016] K. Abe _et al._ (Super-Kamiokande), Solar Neutrino Measurements in Super-Kamiokande-IV, Phys. Rev. D 94, 052010 (2016), arXiv:1606.07538 [hep-ex] .
  * Abe  _et al._ [2018] K. Abe _et al._ (Super-Kamiokande), Atmospheric neutrino oscillation analysis with external constraints in Super-Kamiokande I-IV, Phys. Rev. D 97, 072001 (2018), arXiv:1710.09126 [hep-ex] .
  * Abe  _et al._ [2020] K. Abe _et al._ (T2K), Constraint on the matter–antimatter symmetry-violating phase in neutrino oscillations, Nature 580, 339 (2020), [Erratum: Nature 583, E16 (2020)], arXiv:1910.03887 [hep-ex] .
  * Acero  _et al._ [2019] M. A. Acero _et al._ (NOvA), First Measurement of Neutrino Oscillation Parameters using Neutrinos and Antineutrinos by NOvA, Phys. Rev. Lett. 123, 151803 (2019), arXiv:1906.04907 [hep-ex] .
  * Adamson  _et al._ [2020] P. Adamson _et al._ (MINOS+), Precision Constraints for Three-Flavor Neutrino Oscillations from the Full MINOS+ and MINOS Dataset, Phys. Rev. Lett. 125, 131802 (2020), arXiv:2006.15208 [hep-ex] .
  * Aartsen  _et al._ [2018] M. G. Aartsen _et al._ (IceCube), Measurement of Atmospheric Neutrino Oscillations at 6–56 GeV with IceCube DeepCore, Phys. Rev. Lett. 120, 071801 (2018), arXiv:1707.07081 [hep-ex] .
  * Adey  _et al._ [2018] D. Adey _et al._ (Daya Bay), Measurement of the Electron Antineutrino Oscillation with 1958 Days of Operation at Daya Bay, Phys. Rev. Lett. 121, 241805 (2018), arXiv:1809.02261 [hep-ex] .
  * Bak  _et al._ [2018] G. Bak _et al._ (RENO), Measurement of Reactor Antineutrino Oscillation Amplitude and Frequency at RENO, Phys. Rev. Lett. 121, 201801 (2018), arXiv:1806.00248 [hep-ex] .
  * de Kerret  _et al._ [2020] H. de Kerret _et al._ (Double Chooz), Double Chooz $\theta$13 measurement via total neutron capture detection, Nature Phys. 16, 558 (2020), arXiv:1901.09445 [hep-ex] .
  * Gando  _et al._ [2013] A. Gando _et al._ (KamLAND), Reactor On-Off Antineutrino Measurement with KamLAND, Phys. Rev. D 88, 033001 (2013), arXiv:1303.4667 [hep-ex] .
  * Esteban  _et al._ [2020] I. Esteban, M. C. Gonzalez-Garcia, M. Maltoni, T. Schwetz, and A. Zhou, The fate of hints: updated global analysis of three-flavor neutrino oscillations, JHEP 09, 178, arXiv:2007.14792 [hep-ph] .
  * NSAC NLDBD Subcommittee [2015] NSAC NLDBD Subcommittee, Report to the nuclear science advisory committee: Neutrinoless double beta decay (2015), http://science.energy.gov/np/nsac/reports.
  * Giuliani  _et al._ [2019] A. Giuliani, J. J. Gomez Cadenas, S. Pascoli, E. Previtali, R. Saakyan, K. Schäffner, and S. Schönert (APPEC Committee), Double Beta Decay APPEC Committee Report (2019), arXiv:1910.04688 [hep-ex] .
  * Armstrong  _et al._ [2019] W. R. Armstrong _et al._ (CUPID), CUPID pre-CDR (2019), arXiv:1907.09376 [physics.ins-det] .
  * Abgrall  _et al._ [2021] N. Abgrall _et al._ (LEGEND), LEGEND-1000 Preconceptual Design Report (2021), arXiv:2107.11462 [physics.ins-det] .
  * Kharusi  _et al._ [2018] S. A. Kharusi _et al._ (nEXO), nEXO Pre-Conceptual Design Report (2018), arXiv:1805.11142 [physics.ins-det] .
  * Kotila and Iachello [2012] J. Kotila and F. Iachello, Phase space factors for double-$\beta$ decay, Phys. Rev. C 85, 034316 (2012), arXiv:1209.5722 [nucl-th] .
  * Cirigliano  _et al._ [2018a] V. Cirigliano, W. Dekens, E. Mereghetti, and A. Walker-Loud, Neutrinoless double-$\beta$ decay in effective field theory: The light-Majorana neutrino-exchange mechanism, Phys. Rev. C 97, 065501 (2018a), [Erratum: Phys.Rev.C 100, 019903 (2019)], arXiv:1710.01729 [hep-ph] .
  * Menéndez [2018] J. Menéndez, Neutrinoless $\beta\beta$ decay mediated by the exchange of light and heavy neutrinos: The role of nuclear structure correlations, J. Phys. G 45, 014003 (2018).
  * Horoi and Neacsu [2016] M. Horoi and A. Neacsu, Shell model predictions for 124Sn double-$\beta$ decay, Phys. Rev. C 93, 024308 (2016).
  * Coraggio  _et al._ [2020] L. Coraggio, A. Gargano, N. Itaco, R. Mancino, and F. Nowacki, The calculation of the neutrinoless double-beta decay matrix element within the realistic shell model, Phys. Rev. C 101, 044315 (2020), arXiv:2001.00890 [nucl-th] .
  * Mustonen and Engel [2013] M. T. Mustonen and J. Engel, Large-scale calculations of the double-$\beta$ decay of ${}^{76}Ge,^{130}Te,^{136}Xe$, and ${}^{150}Nd$ in the deformed self-consistent Skyrme quasiparticle random-phase approximation, Phys. Rev. C 87, 064302 (2013).
  * Hyvarinen and Suhonen [2015] J. Hyvarinen and J. Suhonen, Nuclear matrix elements for $0\nu\beta\beta$ decays with light or heavy Majorana-neutrino exchange, Phys. Rev. C 91, 024613 (2015).
  * Šimkovic  _et al._ [2018] F. Šimkovic, A. Smetana, and P. Vogel, $0\nu\beta\beta$ nuclear matrix elements, neutrino potentials and $\mathrm{SU}(4)$ symmetry, Phys. Rev. C 98, 064325 (2018).
  * Fang  _et al._ [2018] D.-L. Fang, A. Faessler, and F. Šimkovic, 0$\nu\beta\beta$-decay nuclear matrix element for light and heavy neutrino mass mechanisms from deformed quasiparticle random-phase approximation calculations for 76Ge, 82Se, 130Te, 136Xe, and 150Nd with isospin restoration, Phys. Rev. C 97, 045503 (2018).
  * Terasaki [2020] J. Terasaki, Strength of the isoscalar pairing interaction determined by a relation between double-charge change and double-pair transfer for double- $\beta$ decay, Phys. Rev. C 102, 044303 (2020), arXiv:2003.03542 [nucl-th] .
  * Rodriguez and Martinez-Pinedo [2010] T. R. Rodriguez and G. Martinez-Pinedo, Energy density functional study of nuclear matrix elements for neutrinoless $\beta\beta$ decay, Phys. Rev. Lett. 105, 252503 (2010).
  * López Vaquero  _et al._ [2013] N. López Vaquero, T. R. Rodríguez, and J. L. Egido, Shape and pairing fluctuations effects on neutrinoless double beta decay nuclear matrix elements, Phys. Rev. Lett. 111, 142501 (2013).
  * Song  _et al._ [2017] L. S. Song, J. M. Yao, P. Ring, and J. Meng, Nuclear matrix element of neutrinoless double-$\beta$ decay: Relativity and short-range correlations, Phys. Rev. C 95, 024305 (2017).
  * Barea  _et al._ [2015] J. Barea, J. Kotila, and F. Iachello, $0\nu\beta\beta$ and $2\nu\beta\beta$ nuclear matrix elements in the interacting boson model with isospin restoration, Phys. Rev. C 91, 034304 (2015).
  * Deppisch  _et al._ [2020] F. F. Deppisch, L. Graf, F. Iachello, and J. Kotila, Analysis of light neutrino exchange and short-range mechanisms in $0\nu\beta\beta$ decay, Phys. Rev. D 102, 095016 (2020), arXiv:2009.10119 [hep-ph] .
  * Adhikari  _et al._ [2021] G. Adhikari _et al._ (nEXO), nEXO: Neutrinoless double beta decay search beyond $10^{28}$ year half-life sensitivity (2021), arXiv:2106.16243 [nucl-ex] .
  * Gysbers  _et al._ [2019] P. Gysbers _et al._ , Discrepancy between experimental and theoretical $\beta$-decay rates resolved from first principles, Nature Phys. 15, 428 (2019), arXiv:1903.00047 [nucl-th] .
  * Yao  _et al._ [2020] J. M. Yao, B. Bally, J. Engel, R. Wirth, T. R. Rodríguez, and H. Hergert, $AbInitio$ Treatment of Collective Correlations and the Neutrinoless Double Beta Decay of 48Ca, Phys. Rev. Lett. 124, 232501 (2020), arXiv:1908.05424 [nucl-th] .
  * Novario  _et al._ [2021] S. Novario, P. Gysbers, J. Engel, G. Hagen, G. R. Jansen, T. D. Morris, P. Navrátil, T. Papenbrock, and S. Quaglioni, Coupled-Cluster Calculations of Neutrinoless Double-$\beta$ Decay in 48Ca, Phys. Rev. Lett. 126, 182502 (2021), arXiv:2008.09696 [nucl-th] .
  * Belley  _et al._ [2021] A. Belley, C. G. Payne, S. R. Stroberg, T. Miyagi, and J. D. Holt, $AbInitio$ Neutrinoless Double-Beta Decay Matrix Elements for 48Ca , 76Ge , and 82Se, Phys. Rev. Lett. 126, 042502 (2021), arXiv:2008.06588 [nucl-th] .
  * Menéndez  _et al._ [2011] J. Menéndez, D. Gazit, and A. Schwenk, Chiral two-body currents in nuclei: Gamow-Teller transitions and neutrinoless double-beta decay, Phys. Rev. Lett. 107, 062501 (2011), arXiv:1103.3622 [nucl-th] .
  * Cirigliano  _et al._ [2018b] V. Cirigliano, W. Dekens, J. De Vries, M. L. Graesser, E. Mereghetti, S. Pastore, and U. Van Kolck, New Leading Contribution to Neutrinoless Double-$\beta$ Decay, Phys. Rev. Lett. 120, 202001 (2018b), arXiv:1802.10097 [hep-ph] .
  * Cirigliano  _et al._ [2021a] V. Cirigliano, W. Dekens, J. de Vries, M. Hoferichter, and E. Mereghetti, Toward Complete Leading-Order Predictions for Neutrinoless Double $\beta$ Decay, Phys. Rev. Lett. 126, 172002 (2021a), arXiv:2012.11602 [nucl-th] .
  * Cirigliano  _et al._ [2021b] V. Cirigliano, W. Dekens, J. de Vries, M. Hoferichter, and E. Mereghetti, Determining the leading-order contact term in neutrinoless double $\beta$ decay, JHEP 05, 289, arXiv:2102.03371 [nucl-th] .
  * Wirth  _et al._ [2021] R. Wirth, J. M. Yao, and H. Hergert, Ab initio calculation of the contact operator contribution in the standard mechanism for neutrinoless double beta decay (2021), arXiv:2105.05415 [nucl-th] .
  * Jokiniemi  _et al._ [2021] L. Jokiniemi, P. Soriano, and J. Menéndez, Impact of the leading-order short-range nuclear matrix element on the neutrinoless double-beta decay of heavy nuclei (2021), arXiv:2107.13354 [nucl-th] .
  * Cirigliano  _et al._ [2020] V. Cirigliano, W. Detmold, A. Nicholson, and P. Shanahan, Lattice QCD Inputs for Nuclear Double Beta Decay (2020), arXiv:2003.08493 [nucl-th] .
  * Davoudi and Kadam [2021] Z. Davoudi and S. V. Kadam, Path from Lattice QCD to the Short-Distance Contribution to 0$\nu\beta\beta$ Decay with a Light Majorana Neutrino, Phys. Rev. Lett. 126, 152003 (2021), arXiv:2012.02083 [hep-lat] .
  * Agostini  _et al._ [2017] M. Agostini, G. Benato, and J. Detwiler, Discovery probability of next-generation neutrinoless double- $\beta$ decay experiments, Phys. Rev. D 96, 053001 (2017), arXiv:1705.02996 [hep-ex] .
  * Caldwell  _et al._ [2017] A. Caldwell, A. Merle, O. Schulz, and M. Totzauer, Global Bayesian analysis of neutrino mass data, Phys. Rev. D 96, 073001 (2017), arXiv:1705.01945 [hep-ph] .
  * Agostini  _et al._ [2021] M. Agostini, G. Benato, S. Dell’Oro, S. Pirro, and F. Vissani, Discovery probabilities of Majorana neutrinos based on cosmological data, Phys. Rev. D 103, 033008 (2021), arXiv:2012.13938 [hep-ph] .
  * Cirigliano  _et al._ [2018c] V. Cirigliano, W. Dekens, J. de Vries, M. L. Graesser, and E. Mereghetti, A neutrinoless double beta decay master formula from effective field theory, JHEP 12, 097, arXiv:1806.02780 [hep-ph] .
  * Tello  _et al._ [2011] V. Tello, M. Nemevsek, F. Nesti, G. Senjanovic, and F. Vissani, Left-Right Symmetry: from LHC to Neutrinoless Double Beta Decay, Phys. Rev. Lett. 106, 151801 (2011), arXiv:1011.3522 [hep-ph] .
  * King  _et al._ [2013] S. F. King, A. Merle, and A. J. Stuart, The Power of Neutrino Mass Sum Rules for Neutrinoless Double Beta Decay Experiments, JHEP 12, 005, arXiv:1307.2901 [hep-ph] .
  * Rodejohann [2011] W. Rodejohann, Neutrino-less Double Beta Decay and Particle Physics, Int. J. Mod. Phys. E 20, 1833 (2011), arXiv:1106.1334 [hep-ph] .