# Sampling an Edge in Sublinear Time Exactly and Optimally

Talya Eden Bar Ilan University<EMAIL_ADDRESS>supported by the NSF
TRIPODS program, award CCF-1740751. This work was partially done while
affiliated with MIT and Boston University. Shyam Narayanan MIT,
<EMAIL_ADDRESS>supported by the NSF GRFP Fellowship and the NSF TRIPODS
Program (award DMS-2022448). Jakub Tětek BARC, Univ. of Copenhagen,
<EMAIL_ADDRESS>supported by the VILLUM Foundation grant 16582. This work
was partially done while visiting MIT.

###### Abstract

Sampling edges from a graph in sublinear time is a fundamental problem and a
powerful subroutine for designing sublinear-time algorithms. Suppose we have
access to the vertices of the graph and know a constant-factor approximation
to the number of edges. An algorithm for pointwise $\varepsilon$-approximate
edge sampling with complexity $O(n/\sqrt{\varepsilon m})$ has been given by
Eden and Rosenbaum [SOSA 2018]. This has been later improved by Tětek and
Thorup [STOC 2022] to $O(n\log(\varepsilon^{-1})/\sqrt{m})$. At the same time,
$\Omega(n/\sqrt{m})$ time is necessary. We close the problem, by giving an
algorithm with complexity $O(n/\sqrt{m})$ for the task of sampling an edge
exactly uniformly.

##  1 Introduction

Suppose we have a graph too big to even read the whole input. We then need an
algorithm running in time sublinear in the input size. Such algorithms have
recently received a lot of attention. In the sublinear-time settings, one
usually has direct access to the vertices of the input graph, but not to the
edges. Because of this, one tool commonly used for designing sublinear-time
graph algorithms is an algorithm for sampling edges. This allows us to design
an algorithm that uses random edge queries, as we can simulate these queries
by the edge sampling algorithm.

##### The task and query access.

Our goal is to sample an edge uniformly, i.e., to return an edge so that each
edge is returned with exactly equal probability. We assume that the algorithm
may $(i)$ ask for the $i$-th vertex of the input graph, $(ii)$ ask for the
degree of a given vertex, and $(iii)$ ask for the $j$-th neighbor of a given
vertex. We assume the algorithm has (approximate) knowledge of the number of
edges $m$. This assumption of knowing $m$ was not made in the previous work,
and we think getting rid of this assumption is a very interesting open
problem. This assumption is, however, not a barrier to using our algorithm as
a subroutine for implementing random edge queries, as we discuss below.

##### Ours and previous results.

In past work only algorithms for sampling $\varepsilon$-_approximately_
uniformly in the pointwise distance (or equivalently approximately in the
$\ell_{\infty}$ metric) were given, where the state-of-the-art complexity for
that problem was $O(n\log(\varepsilon^{-1})/\sqrt{m})$ given by Tětek and
Thorup [13]. Our algorithm is not only exact, but also more efficient.
Specifically, the expected complexity of our algorithm is $O(n/\sqrt{m})$.
This is known to be the best possible. If we have a graph with
$n-\Theta(\sqrt{m})$ isolated vertices and a clique over $\Theta(\sqrt{m})$
vertices with $m$ edges, we need to sample $\Omega(n/\sqrt{m})$ vertices
before we expect to see a single edge, giving us a simple matching lower
bound.

##### Using our algorithm as a subroutine and the necessity of knowing $m$.

As we said above, our algorithm needs to have a constant-factor approximation
of the number of edges. This was not necessary in previous works. The reason
is that the previous state-of-the-art has complexity in which one can also
afford to independently estimate the number of edges; the possibility of the
estimate being incorrect is then added to the bias of the edge sampling
algorithm. However, as our algorithm is more efficient, we are not able to
estimate $m$ in that complexity, and since we want to sample exactly, we
cannot accept that the estimate could be wrong.

Suppose we have an algorithm $A$ that performs random edge queries. We may
then use our algorithm in a black box manner to implement these queries
(unlike, for example, the algorithm for sampling multiple edges from [8] which
has polynomial dependency on $\varepsilon$). Specifically, if $A$ uses $q$
random edge queries, then we may set $\varepsilon=1/(10q)$ and it will only
decrease the success probability of $A$ by at most $1/10$. 111This holds
because the total variation distance from uniform of each query is at most
$1/(10q)$, so the total variation distance from uniform of the sequence of $q$
queries is at most $1/10$, meaning that the output from the algorithm has
total variation distance at most $1/10$ from the distribution the output would
have if the queries were answered exactly.

If the goal is to get an algorithm with a constant success probability (which
can then be amplified) that uses our edge sampling algorithm as a subroutine,
then we may remove the need for having an a priori constant-factor
approximation $\tilde{m}$ of $m$ by computing it using the algorithm from [9],
only adding a constant to the failure probability. The algorithm from [9] has
expected complexity $O(n/\sqrt{m})$. The complexity of our algorithm will also
still be as desired: it follows from our analysis that the complexity is
$O(\frac{n\sqrt{\tilde{m}}}{m})$. It holds by the Jensen inequality that
$\mathbb{E}[\frac{n\sqrt{\tilde{m}}}{m}]\leq\frac{n\sqrt{\mathbb{E}[\tilde{m}}]}{m}+\log\varepsilon^{-1}=O(n/\sqrt{m})$
since it holds $\mathbb{E}[\tilde{m}]=O(m)$.

To summarize, we may remove the assumption of a priori knowledge of $m$ when
sampling multiple edges at the cost of adding a constant failure probability.
This means that we may use our algorithm as a subroutine in an algorithm with
constant probability of error, even without knowing $m$ a priori.

###  1.1 Technical overview

The starting point of our algorithm is the algorithm by Eden and Rosenbaum
[6], which we now shortly recall.

#### The algorithm by Eden and Rosenbaum [6].

Consider each undirected edge as two directed edges, and let $\theta$ be a
degree threshold. We refer to vertices with degree at most $\theta$ as _light_
vertices, and to all other vertices as _heavy_. We refer to edges originating
in light vertices as _light edges_ , and to all other edges as _heavy edges_.
Using rejection sampling, light edges can be sampled with probability
_exactly_ $\frac{1}{n\theta}$: by sampling a uniform vertex $v$, then sampling
one of its incident edges u.a.r., and then returning that edge with
probability $\frac{d(v)}{\theta}$. Sampling heavy vertices is done by first
sampling a light edge $uv$ as described above, and if the second endpoint $v$
of the sampled light edge is heavy, sampling one of its incident edges. This
procedure results in every heavy edge $vw$ being sampled with probability
$\frac{d_{\ell}(v)}{n\theta}\cdot\frac{1}{d(v)}$, where $d_{\ell}(v)$ is the
number of light neighbors of $v$. In Eden and Rosenbaum [6], $\theta$ is set
to $\sqrt{2m/\varepsilon}$ which implies that for every heavy vertex $v$,
$d_{\ell}(v)\in[(1-\varepsilon)d(v),d(v)]$.222Since, denoting by $\cal{H}$ the
set of heavy vertices, and by $d_{h}(v)$ the number of heavy neighbors of
vertex $v$, we have the following. For every $v\in\cal{H}$,
$d_{h}(v)\leq|H|\leq\frac{2m}{\theta}=\sqrt{2\varepsilon
m}=\varepsilon\theta\leq\varepsilon d(v)$. Therefore,
$d_{\ell}(v)>(1-\varepsilon)d(v)$. Hence, each (heavy) edge is sampled with
probability in $[\frac{(1-\varepsilon)}{n\theta},\frac{1}{n\theta}]$. The
total probability of sampling some edge (with the algorithm failing otherwise)
is thus at least
$\frac{(1-\varepsilon)m}{n\theta}\approx\frac{\sqrt{\varepsilon m}}{n}$.
Therefore, the number of attempts needed before we expect to sample an edge is
$O(\frac{n}{\sqrt{\varepsilon m}})$, implying a multiplicative dependence on
$\varepsilon$.

#### Improving the dependency on $\varepsilon$.

In order to avoid the multiplicative dependency in $\varepsilon$, we instead
set the threshold $\theta$ to $\sqrt{cm}$ for some constant $c$. Considering
the same sampling procedures as before, light edges can still be sampled with
probability exactly $\frac{1}{n\theta}$. For heavy edges, however, the values
$d_{\ell}(v)/d(v)$ can vary up to a constant factor between the different
heavy vertices, leading to a large bias towards heavy edges originating in
vertices $v$ with higher values of $\frac{d_{\ell}(v)}{d(v)}$. If for each
vertex $v$, we knew the value of $d_{\ell}(v)$, we could use rejection
sampling with probability $q$ that is inversely proportional to
$p=\frac{d_{\ell}(v)}{d(v)}$, e.g., $q=\frac{d(v)}{2d_{\ell(v)}}=\frac{1}{2p}$
(we may assume that, say, $p\geq 2/3$ and thus $q<1$, by making $c$ large
enough). This would result in each heavy edge being sampled with exactly equal
probability
$\frac{d_{\ell(v)}}{n\theta}\cdot\frac{d(v)}{2d_{\ell}(v)}=\frac{1}{2n\theta}$.

While we do not know the exact value of $d_{\ell}(v)$, we can approximate it
up to a $(1\pm\Theta(\varepsilon))$-multiplicative factor using
$O(1/\varepsilon^{2})$ neighbor queries. This results in
$(1\pm\Theta(\varepsilon))$-approximation of $q$ and thus leads to a
distribution that is $\varepsilon$-close to uniform. Note that we only need to
approximate $q$ when the algorithm samples a heavy edge. Moreover, when we do
that, we return the edge with constant probability. Thus, in expectation, we
only need to approximate $q$ a constant number of times. This means that the
total expected time complexity is $O(n/\sqrt{m}+1/\varepsilon^{2})$.

To remove the dependence on $1/\varepsilon^{2}$, our main observation is that
we do not actually need to (approximately) learn the value of $p$, in order to
reject with probability proportional to $q=1/(2p)$. Rather, we can “simulate”
a Bernoulli trial with probability exactly $\frac{1}{2p}$ by using the results
of only $O(1)$ many $Bern(p)$ trials in expectation (though possibly more in
the worst case), using the Bernoulli Factory technique of Nacu and Peres [12].

When we sample a uniform neighbor of a heavy vertex $v$, we see a light
neighbor of $v$ with probability exactly $p=\frac{d_{\ell}(v)}{d(v)}$, where,
as discussed above, we can set $\theta$ so that $p>2/3$. Therefore, we have
access to a Bernoulli trial that succeeds with probability $Bern(p)$ for
$p>2/3$. As previously explained, in order to achieve uniformity, we need to
perform rejection sampling (corresponding to a Bernoulli trial) that succeeds
with probability $=\frac{1}{2p}$. We can simulate $Bern(1/(2p))$ by relying on
the results of an expected $O(1)$ independent copies of $Bern(p)$. Namely, we
perform an expected $O(1)$ neighbor queries where each results in a light
neighbor with probability $p$, giving us an independent copy of a random
variable distributed as $Bern(p)$. If we knew that $p$ is bounded away from
$1$, we could directly use the result of Nacu and Peres [12]. We can ensure
that this is the case by first rejecting with probability, say, $1/2$. The
probability of getting a light neighbor and not rejecting is then
$p/2\in[1/3,2/3]$. We then use the result of [12] with the function $1/(4x)$,
thus simulating $Bern(1/(4p/2))=Bern(1/(2p))$.

###  1.2 Related work

Using uniform edge samples as a basic query in the sublinear time setting was
first suggested by Aliakbarpour et al. [2] in the context of estimating the
number of $s$-stars in a graph, where they showed that this access allows to
circumvent lower bounds that hold in the standard adjacency list access. It
was later used for the more general tasks of estimating and uniformly sampling
arbitrary subgraphs in sublinear time [3, 10, 5, 13].

As mentioned in the introduction, sampling edges from a distribution that is
pointwise close to uniform in sublinear time was first suggested by Eden and
Rosenbaum [6] who gave an algorithm with complexity
$O\left(n/\sqrt{\varepsilon m}\right)$. This was later improved by Tětek and
Thorup [13] to an algorithm with complexity
$O(n\log(\varepsilon^{-1})/\sqrt{m})$. They also considered two additional
access models (full neighborhood access and hash-ordered access) and gave new
lower and upper bounds for these settings. In [8], Eden, Mossel, and Rubinfeld
gave an upper bound for the problem of sampling $k$ edges from a pointwise
close to uniform distribution. The complexity of their algorithm is
$O\left(\sqrt{k}\cdot\frac{n}{\sqrt{m}}\cdot\frac{\log^{2}n}{\varepsilon^{2.5}}+k\right)$.
This was later shown to be essentially optimal (i.e., up to the dependencies
on $\varepsilon$ and $\log n$) by [13]. In [7], Eden, Ron and Rosenbaum gave
an $O\left(\frac{n\alpha}{m}\cdot\frac{\log^{3}n}{\varepsilon}\right)$
algorithm for sampling edges in graphs with arboricity at most $\alpha$. They
also showed their algorithm is optimal up to the $\text{poly}(\log
n,\varepsilon^{-1})$ dependencies.

The task of sampling close to uniform edges was also recently considered in
the setting where the access to the graph is given via Bipartite Independent
Set (BIS) queries [11, 4, 1].

##  2 Preliminaries

##### The query model:

Since we do not have time to read the whole input (and thus to change its
representation), it matters how exactly we are able to query it. Throughout
this paper, we assume that the graphs’ vertices are labeled arbitrarily in
$[n]$, the edges of every vertex $v$ are labeled arbitrarily by $[d(v)]$, and
that the algorithm knows $n$. We then assume the following standard set of
queries:

  * •

Uniform vertex queries: given $i\in[n]$, return the $i$-th vertex

  * •

Degree queries: given a vertex $v$, return its degree $d(v)$

  * •

Neighbor queries: given a vertex $v$ and $j\in[d(v)]$, return the $j$-th
neighbor of $v$

This setting has been previously called the adjacency list or indexed neighbor
access model and is among the most-studied settings for sublinear-time
algorithms.

A model with an additional query

  * •

Uniform edge queries: given $i\in[m]$, return vertices $u,v$ such that $uv$ is
the $i$-th edge of the graph

has also been considered. Our algorithm can be thought of as a reduction
between the two models. One can show (by randomly permuting the edges) that up
to a logarithmic factor this setting is equivalent to just assuming random
edge queries (with replacement). Our algorithm then allows us to simulate
random edge queries in the indexed neighborhood access model.

##  3 Sampling an edge

In this section, we give our algorithm for sampling an edge. We start by
stating a result of Nacu and Peres [12] about “Bernoulli factories”. We recall
that a function $f$ from a closed interval $I$ to $\mathbb{R}$ is said to be
_real analytic_ if for any point $x_{0}$ in the interior of $I$,
$f(x)=\sum_{n=0}^{\infty}\frac{f^{(n)}(x_{0})}{n!}\cdot(x-x_{0})^{n}$ where
$f^{(n)}(x_{0})$ represents the $n$th derivative of $f$ at $x_{0}$.
Equivalently, $f$ matches its Taylor series about $x_{0}$ for all of $I$.

###### Theorem 3.1

Let $I\subset(0,1)$ be a closed interval and $f:I\rightarrow(0,1)$ be a real
analytic function. Let $p$ be a number in $I$. Then, there exists an algorithm
independent of $p$ that performs in expectation $O(1)$ independent trials from
$Bern(p)$ and returns one Bernoulli trial with distribution $Bern(f(p))$. In
addition, the probability of using more than $k$ independent trials is at most
$C\rho^{k}$, for some constants $C\geq 1,0<\rho<1$ that only depend on $I$.

Theorem 3.1 gives us the following corollary.

###### Corollary 3.1

There is an algorithm that for any $p\in[2/3,1]$ performs in expectation
$O(1)$ trials from $Bern(p)$ and returns one Bernoulli trial distributed as
$Bern(1/(2p))$. (Note that the algorithm may also use its own randomness,
independent of the $Bern(p)$’s given.)

  * Proof.

First, note that for any $p$ we can simulate $Bern(p/2)$ from a single
$Bern(p)$, by simulating an independent $Bern(1/2)$ and considering the event
where both $Bern(p)$ and $Bern(1/2)$ equal $1$. Then, it is well-known that
$\frac{1}{4x}$ is real analytic on the interval $[1/2,2/3]$, and
$\frac{1}{4x}$ is contained in $[3/8,1/2]$ for $x\in[1/2,2/3]$. Since we can
generate $Bern(p/2)$, we can therefore apply Theorem 3.1 to get a trial from
$Bern\left(\frac{1}{4(p/2)}\right)=Bern\left(\frac{1}{2p}\right)$.

We now give an algorithm for sampling an edge. The algorithm closely follows
the approach from [6] but uses the Bernoulli factory of [12] to reduce the
sampling error in a significantly more efficient way than in [6]. We first
give an algorithm for $\Theta(1)$-approximate edge sampling.

Throughout this section, we assume for sake of simplicity that we know the
number of edges exactly. The analysis of correctness only uses that we have an
upper bound, while the analysis of the complexity needs that we have a lower
bound up to a constant factor. Putting this together, it is in fact sufficient
to have a constant-factor approximation.

1 $u\leftarrow$ uniformly random vertex

2 $j\leftarrow Unif([\theta]),$ where $\theta=\lceil\sqrt{6m}\rceil$.

3 Fail if $d(v)>\theta$ or $d(v)\leq j$

4 $v\leftarrow j$-th neighbor of $u$

5 $B\sim Bern(1/3)$

6 if  _$B=1$_ then

7 return _$uv$_

8

9else if  _$B=0$ and $v$ is heavy_ then

10 $w\leftarrow$ random neighbor of $v$

11 return _$vw$_

12 end if

return _Fail_

Algorithm 1 Sample an edge pointwise $\Theta(1)$-close to uniform

###### Lemma 3.1

Let $e$ be the edge returned by Algorithm 1 if successful. Then for any light
edge $e^{\prime}$, it holds that $\mathbb{P}(e=e^{\prime})=1/(3n\theta)$, and
for any heavy edge, it holds that
$\mathbb{P}(e=e^{\prime})=\frac{2}{3}\cdot\frac{d_{\ell}(v)}{d(v)}\cdot\frac{1}{n\theta}$.

  * Proof.

Fix a light edge $e^{\prime}=uv$. Recall that by definition, $uv$ is light iff
$d(u)\leq\theta$ for $\theta=\lceil\sqrt{6m}\rceil$. The edge $uv$ is returned
only in the case that (1) $u$ is sampled in Step 1, (2) the chosen index $j$
in Step 1 is the label of $v$, and (3) $B=1$ in Step 1. Therefore,
$\Pr[e=e^{\prime}]=\frac{1}{n}\cdot\frac{1}{\lceil\sqrt{6m}\rceil}\cdot\frac{1}{3}=\frac{1}{3n\theta}$.

Now fix a heavy edge $e^{\prime}=vw$. The edge $vw$ is returned in the event
that (1) the sampled vertex $u$ in Step 1 is a light neighbor of $v$, (2) the
chosen index $j$ in Step 1 is the label of $v$, (3) $B=0$ in Step 1, and (4)
$w$ is the sampled neighbor in Step 1. Therefore, if we define $\Gamma_{L}(v)$
to be the set of light neighbors of $v$, then
$\Pr[e=e^{\prime}]=\sum_{u\in\Gamma_{L}(v)}\frac{1}{n}\cdot\frac{1}{\lceil\sqrt{6m}\rceil}\cdot\frac{2}{3}\cdot\frac{1}{d(v)}=\frac{2}{3}\cdot\frac{d_{\ell}(v)}{d(v)}\cdot\frac{1}{n\theta}$.

We are now able to give an algorithm for sampling an edge perfectly uniformly.
Simply re-running the above algorithm until it succeeds would result in
$\Theta(1)$-pointwise close to uniform sampling. The algorithm below differs
in that if Algorithm 1 returns a heavy edge (which has some bias), we use
rejection sampling based on Bernoulli factories to reduce the bias.

1 repeat  __

2 $vw\leftarrow$ Algorithm 1

3 if  _$v$ is light_ then

4 return _vw_

5 end if

6 if  _$v$ is heavy_ then

7 Let $w_{1},\dots,$ be random neighbors of $v$

8 $Y\leftarrow$ use Corollary 3.1 on Bernoulli trials defined as
$B_{i}=[d(w_{i})\leq\sqrt{6m}]$

9 if  _$Y=0$_ then

10 return _vw_

11 end if

12

13 end if

14

15 end

Algorithm 2 Sample an edge $1\pm\varepsilon$-pointwise-close to uniform

###### Theorem 3.2

Algorithm 2 returns a perfectly uniform edge. Its expected complexity is
$O(n/\sqrt{m})$.

  * Proof.

We start with proving the correctness of the algorithm. By Lemma 3.1, each
invocation of Algorithm 1 returns each light edge with probability
$\frac{1}{3n\theta}$, and each heavy edge with probability
$\frac{2}{3}\cdot\frac{d_{\ell}(v)}{d(v)}\cdot\frac{1}{n\theta}$. If Algorithm
1 returns a heavy edge $vw$, then for every $w_{i}$ sampled in Step 2 in
Algorithm 2, it holds that the indicator of the event
$[d(w_{i})\leq\lceil\sqrt{6m}\rceil]$ is the result of a Bernoulli trial
$Bern(p)$ with $p=\frac{d_{\ell}(v)}{d(v)}$. Let $H$ denote the set of
vertices with degree greater than $\lceil\sqrt{6m}\rceil$. Then
$deg_{H}(v)\leq|H|\leq\frac{2m}{\lceil\sqrt{6m}\rceil}\leq\sqrt{\frac{2}{3}m}\leq\frac{1}{3}d(v)$,
where the last is since $v$ is heavy (so $d(v)>\lceil\sqrt{6m}\rceil)$.
Therefore, $p=\frac{d_{\ell}(v)}{d(v)}\in[\frac{2}{3},1]$. Hence, by Corollary
3.1, the value $Y$ returned by Algorithm 2 has distribution
$Y\sim\text{Bern}(\frac{1}{2p})$. Therefore, in a single iteration of the
repeat loop, every fixed heavy edge $vw$ is returned with probability
$\frac{2}{3}\cdot\frac{d_{\ell}(v)}{d(v)}\cdot\frac{1}{n\theta}\cdot\left(\frac{1}{2p}\right)\in\frac{1}{3n\theta}$,
as $p=\frac{d_{\ell}(v)}{d(v)}$.

Therefore, every edge is sampled with probability exactly $\frac{1}{3n},$ so
conditioning on some edge being returned, each edge is returned with
probability in $\frac{1}{m}$, as claimed.

We turn to analyze the complexity of the algorithm. By the above analysis,
every invocation of the loop returns an edge with probability at least
$m\cdot\frac{1}{3n\theta}\geq\frac{\sqrt{m}}{10n}.$ Also, note that each
invocation is independent. Therefore, the expected number of iterations until
an edge is returned is $O(n/\sqrt{m})$. Furthermore, each invocation of the
loop invokes Algorithm 1 once, and Corollary 3.1 at most once. Algorithm 1
clearly takes a constant number of queries. If Algorithm 1 returns a heavy
edge, then sampling the $w_{i}$ neighbors in Step 2 takes $O(1)$ queries in
expectation. Therefore, the expected number of queries in each loop is
constant. Hence, the expected query complexity is
$\Theta\left(n/\sqrt{m}\right)$.

## Acknowledgments

We thank Nima Anari and Peter Occil for informing us about the existence of
Bernoulli factories in the literature, and Nima for pointing us to the
reference [12].

## References

  * Addanki et al. [2022] Raghavendra Addanki, Andrew McGregor, and Cameron Musco.  Non-adaptive edge counting and sampling via bipartite independent set queries.  _arXiv preprint arXiv:2207.02817_ , 2022.
  * Aliakbarpour et al. [2018] Maryam Aliakbarpour, Amartya Shankha Biswas, Themis Gouleakis, John Peebles, Ronitt Rubinfeld, and Anak Yodpinyanee.  Sublinear-time algorithms for counting star subgraphs via edge sampling.  _Algorithmica_ , 80(2):668–697, 2018.
  * Assadi et al. [2019] Sepehr Assadi, Michael Kapralov, and Sanjeev Khanna.  A simple sublinear-time algorithm for counting arbitrary subgraphs via edge sampling.  In _Innovations in Theoretical Computer Science Conference ITCS_ , volume 124 of _LIPIcs_ , pages 6:1–6:20. Schloss Dagstuhl - Leibniz-Zentrum fuer Informatik, 2019.
  * Bhattacharya et al. [2022] Anup Bhattacharya, Arijit Bishnu, Arijit Ghosh, and Gopinath Mishra.  Faster counting and sampling algorithms using colorful decision oracle.  In _39th International Symposium on Theoretical Aspects of Computer Science (STACS 2022)_. Schloss Dagstuhl-Leibniz-Zentrum für Informatik, 2022.
  * Biswas et al. [2021] Amartya Shankha Biswas, Talya Eden, and Ronitt Rubinfeld.  Towards a decomposition-optimal algorithm for counting and sampling arbitrary motifs in sublinear time.  In _Approximation, Randomization, and Combinatorial Optimization. Algorithms and Techniques, APPROX/RANDOM 2021, to appear_ , 2021\.
  * Eden and Rosenbaum [2018] Talya Eden and Will Rosenbaum.  On Sampling Edges Almost Uniformly.  In Raimund Seidel, editor, _1st Symposium on Simplicity in Algorithms (SOSA 2018)_ , volume 61 of _OpenAccess Series in Informatics (OASIcs)_ , pages 7:1–7:9, Dagstuhl, Germany, 2018. Schloss Dagstuhl–Leibniz-Zentrum fuer Informatik.  ISBN 978-3-95977-064-4.  doi: 10.4230/OASIcs.SOSA.2018.7.  URL http://drops.dagstuhl.de/opus/volltexte/2018/8300.
  * Eden et al. [2019] Talya Eden, Dana Ron, and Will Rosenbaum.  The Arboricity Captures the Complexity of Sampling Edges.  In Christel Baier, Ioannis Chatzigiannakis, Paola Flocchini, and Stefano Leonardi, editors, _46th International Colloquium on Automata, Languages, and Programming (ICALP 2019)_ , volume 132 of _Leibniz International Proceedings in Informatics (LIPIcs)_ , pages 52:1–52:14, Dagstuhl, Germany, 2019. Schloss Dagstuhl–Leibniz-Zentrum fuer Informatik.  ISBN 978-3-95977-109-2.  doi: 10.4230/LIPIcs.ICALP.2019.52.  URL http://drops.dagstuhl.de/opus/volltexte/2019/10628.
  * Eden et al. [2021] Talya Eden, Saleet Mossel, and Ronitt Rubinfeld.  Sampling multiple edges efficiently.  In _Approximation, Randomization, and Combinatorial Optimization. Algorithms and Techniques (APPROX/RANDOM 2021)_. Schloss Dagstuhl-Leibniz-Zentrum für Informatik, 2021.
  * Feige [2006] Uriel Feige.  On sums of independent random variables with unbounded variance and estimating the average degree in a graph.  _SIAM Journal on Computing_ , 35(4):964–984, 2006\.
  * Fichtenberger et al. [2020] Hendrik Fichtenberger, Mingze Gao, and Pan Peng.  Sampling arbitrary subgraphs exactly uniformly in sublinear time.  In Artur Czumaj, Anuj Dawar, and Emanuela Merelli, editors, _47th International Colloquium on Automata, Languages, and Programming, ICALP 2020, July 8-11, 2020, Saarbrücken, Germany (Virtual Conference)_ , volume 168 of _LIPIcs_ , pages 45:1–45:13. Schloss Dagstuhl - Leibniz-Zentrum für Informatik, 2020.  doi: 10.4230/LIPIcs.ICALP.2020.45.  URL https://doi.org/10.4230/LIPIcs.ICALP.2020.45.
  * Lapinskas et al. [2019] John A Lapinskas, Holger Dell, and Kitty Meeks.  Approximately counting and sampling small witnesses using a colourful decision oracle.  In _ACM-SIAM Symposium on Discrete Algorithms (SODA20)_ , 2019.
  * Nacu and Peres [2006] Şerban Nacu and Yuval Peres.  Fast simulation of new coins from old.  _SIAM Journal on Computing_ , 35(4):964–984, 2006\.
  * Tětek and Thorup [2022] Jakub Tětek and Mikkel Thorup.  Edge sampling and graph parameter estimation via vertex neighborhood accesses.  In _Proceedings of the 54th Annual ACM SIGACT Symposium on Theory of Computing_ , STOC 2022, page 1116–1129, New York, NY, USA, 2022. Association for Computing Machinery.  ISBN 9781450392648.  doi: 10.1145/3519935.3520059.  URL https://doi.org/10.1145/3519935.3520059.


# Progress Towards Two-loop QCD Corrections to $pp\to t\bar{t}j$

Simon Badger  Dipartimento di Fisica and Arnold-Regge Center, Università di
Torino, and INFN, Sezione di Torino, Via P. Giuria 1, I-10125 Torino, Italy

###### Abstract

We present the status of on-going efforts to compute the two-loop virtual
corrections to $pp\to t\bar{t}j$ in the leading colour approximation. We
review the recent study of the master integrals and their differential
equations and present some ideas on modern techniques that can be applied to
provide complete evaluation of helicity amplitudes containing information on
top quark decays in the narrow width approximation.

##  1 Introduction

Top quark pair production in association with an additional jet is a high
priority process with current and future LHC experiments. This signal forms a
significant fraction, around 50%, of all top-quark pair events and has been
studied extensively by both ATLAS and CMS. [1, 2, 3, 4] The normalized
distributions for $t\bar{t}j$ have been shown to be highly sensitive to the
top-quark mass [5] and therefore precision predictions are in high demand.

While next-to-leading order (NLO) QCD corrections have been known for many
years, [6] the next-to-next-to-leading order (NNLO) QCD corrections are
currently unavailable due to the missing two-loop virtual contributions. The
computational complexity of these amplitudes presents a formidable challenge
to current theoretical methods and new techniques are required. New methods
using finite field arithmetic, [7, 8] highly optimized integration-by-parts
reduction [9, 10, 11, 12] and improved analytic understanding of multi-scale
Feynman integrals [13, 14, 15, 16, 17] have made considerable impact on $2\to
3$ scattering processes. Recent highlights have been the completion of full
colour corrections to $pp\to 3j$ and $pp\to\gamma+2j$ at NNLO. [18, 19, 20,
21] The first steps have been taken to develop these ideas for applications to
$pp\to t{\bar{t}}j$ amplitudes [22, 23] where promising results have been
obtained. Here we report on the status of the leading colour two-loop
amplitudes where differential equations for the master integrals have recently
been derived. [24]

##  2 Amplitudes in the leading colour limit

Since the amplitudes are of an extremely high complexity, we start by
considering the limit of a large number of colour charges $N_{c}\to\infty$.
There is evidence to suggest that such an approximation to the two-loop
virtual correction would be sufficient for next-to-next-to-leading order since
the size with respect to double-real and real-virtual contributions is often
small using standard infra-red subtraction schemes.

Helicity amplitude techniques have been extremely successful for predictions
with massless particles. For massive particles helicity is not well defined
since one can always boost to a frame in which helicity is flipped. We can
therefore only define helicity with respect to a reference direction. Helicity
amplitude techniques can be extremely useful in the case of top-quark pair
production since it is simple to include decays in the narrow width
approximation. [25]

The leading colour decomposition of the loop amplitudes, $\mathcal{A}^{(L)}$
can be written as

$\displaystyle\mathcal{A}^{(L)}$
$\displaystyle(1_{\bar{t}},2_{t},3_{g},4_{g},5_{g})=g_{s}^{3+2L}N_{c}^{L}$
$\displaystyle\sum_{\sigma\in
S_{3}}(t^{a_{\sigma(3)}}t^{a_{\sigma(4)}}t^{a_{\sigma(5)}})_{i_{2}}^{\bar{i}_{1}}A^{(L)}_{1}(1_{\bar{t}},2_{t},\sigma(3)_{g},\sigma(4)_{g},\sigma(5)_{g})$
(1)

in the gluon channel and

$\displaystyle\mathcal{A}^{(L)}$
$\displaystyle(1_{\bar{t}},2_{t},3_{q},4_{\bar{q}},5_{g})=g_{s}^{3+2L}N_{c}^{L}\bigg{\\{}$
$\displaystyle\delta^{\bar{i}_{4}}_{i_{1}}(t^{a_{5}})^{\bar{i}_{2}}_{i_{3}}A^{(L)}_{1}(1_{\bar{t}},2_{t},3_{\bar{q}},4_{q},5_{g})+\delta^{\bar{i}_{3}}_{i_{2}}(t^{a_{5}})^{\bar{i}_{4}}_{i_{1}}A^{(L)}_{2}(1_{\bar{t}},2_{t},3_{\bar{q}},4_{q},5_{g})\bigg{\\}}$
(2)

in the (light) quark channel. Here we have used $g_{s}$ as the strong
coupling, $N_{c}$ as the number of colours and $t^{a}$ are the $SU(N_{c})$
generators. Due to the high gluon luminosity at the LHC, the gluon channel
will dominate the contributions. At leading colour the external particles
always appear with a cyclic ordering and so the top-quark pairs remain
adjacent. This has the effect of minimising the number of massive internal
propagators and therefore leads to considerably simpler Feynman integrals.
Feynman integrals with internal masses can give rise to elliptic integrals
from two-loops onwards. These structures are currently of great interest to
mathematicians and there is an active field of research aiming to obtain an
understanding suitable for collider physics applications. Unfortunately, for
complicated kinematics such as $pp\to t\bar{t}j$, the technology for numerical
evaluation of such functions is not available at the same level as it is for
more standard situations that can be written in terms of multiple
polylogarithms.

From analysis of the Feynman diagrams of the two-loop leading colour
amplitudes, one can find nine different integral families. Three topologies
are formed from bubble insertions into one-loop pentagon graphs and can be
trivially mapped into the three ‘hexagon-triangle’ topologies. After solving
integration-by-parts identities to determine the master integrals, one can
determine that there are no master integrals in the top sector of the hexagon-
triangles and so these master integrals can also be mapped into the ‘pentagon-
box’ families. The families are shown in Figure 1.

(a) Topology ${\rm PB}_{A}$.

(b) Topology ${\rm PB}_{B}$.

(c) Topology ${\rm PB}_{C}$.

(d) Topology ${\rm HT}_{A}$.

(e) Topology ${\rm HT}_{B}$.

(f) Topology ${\rm HT}_{C}$.

Figure 1: The six integral families appearing in the tensor integral
representation of $pp\rightarrow t\bar{t}j$ at two-loops in the leading colour
limit. Black lines denote massless particles and red double-lines denote
massive particles. The master integrals of all three hexagon-triangle
topologies can be mapped to the pentagon-box topologies.

##  3 Evaluation of the pentagon-box master integrals

In the recent preprint article [24], differential equations (DEs) for all
pentagon-box master integrals have been presented. Bases of master integrals
were found for topologies ${\rm PB}_{A}$ and ${\rm PB}_{C}$ where the DEs can
be written in terms of logarithmic one-forms. Contrary to all other pentagon-
box configurations considered so far, topology ${\rm PB}_{B}$ was shown to
have two problematic sectors. One could only be rotated in to
$\epsilon$-factorized form by introducing a nested square-root while another
was shown to contain elliptic structures. As a result the strategy of
expanding the integrals into a basis of stable pentagon functions which can be
evaluated efficiently is not straightforward. As an alternative, we
demonstrated that numerical values could be obtained in a reasonable time from
generalized series expansions of the DEs. [26, 27] Topology ${\rm PB}_{B}$ was
not obtained in $\epsilon$-factorized form but as a quadratic polynomial in
$\epsilon$ which could be represented compactly in terms of linearly
independent one-forms.

We find that topologies ${\rm PB}_{A}$ and ${\rm PB}_{C}$ have 71 and 79
independent one-forms respectively, all of which can be written as
$\mathrm{d}\log$ forms. Topology ${\rm PB}_{B}$ has 135 independent one-forms
of which 72 can be written as $\mathrm{d}\log$ forms. Evaluation times per
segment of the path between points in the generalized series expansion are
around 5 times longer for ${\rm PB}_{B}$ than for ${\rm PB}_{A}$ and ${\rm
PB}_{C}$ but, at around 1 minute per point, feasible for phenomenological
applications. Considerations such as numerical stability and efficient
algorithms for spanning the phasespace while minimizing the number of segments
remain for further study. For the two topologies admitting $\mathrm{d}\log$
representations, it is straightfoward to derive analytic representations in
terms of multiple polylogarithms or pentagon functions that would provide fast
and stable numerics.

##  4 Outlook

Having established a set of master integrals suitable for the complete
amplitude in the planar limit that can be reliably evaluated numerically, the
next step is to obtain their rational coefficients for each amplitude. This
could be attempted though analytic reconstruction from finite field
evaluations (modulo a prime number), for example using the FiniteFlow package,
[28] as has been successful for other $2\to 3$ scattering processes. Analytic
forms for rational coefficients may well be extremely complicated and it may
be useful to also consider numerical approaches. While this step will require
the reduction of high rank tensor integrals, the finite field strategy
combined with optimized IBP systems does appear to be viable. The unknown
aspect of the evaluations is the average time per phase-space point which will
certainly be considerably more expensive than massless $2\to 3$ processes.
There are many unexplored optimizations for efficient evaluation of the
integrals and numerical interpolation techniques that lead us to believe the
two-loop virtual corrections can be obtained via this strategy in the not too
distant future.

## Acknowledgments

I would like to thank my collaborators Matteo Becchetti, Colomba Branccacio,
Nicolò Giraudo, Heribertus Bayu Hartanto and Simone Zoia. This project has
received funding from the European Union’s Horizon Europe research and
innovation programme under the Marie Skłodowska-Curie grant agreement No.
101105486, and ERC Starting Grant No. 101040760 _FFHiggsTop_. This work has
received funding from the Italian Ministry of Universities and Research
through FARE grant R207777C4R. This research was supported in part by the
Swiss National Science Foundation (SNF) under contract 200021_212729. SB has
been partially supported by the Italian Ministry of Universities and Research
(MUR) through grant PRIN 2022BCXSW9.

## References

## References

  * [1] V. Khachatryan et al. [CMS], Phys. Rev. D 95 (2017) no.9, 092001 doi:10.1103/PhysRevD.95.092001
  * [2] M. Aaboud et al. [ATLAS], JHEP 10 (2018), 159 doi:10.1007/JHEP10(2018)159
  * [3] A. M. Sirunyan et al. [CMS], JHEP 07 (2020), 125 doi:10.1007/JHEP07(2020)125
  * [4] [CMS], [arXiv:2402.08486 [hep-ex]].
  * [5] S. Alioli, P. Fernandez, J. Fuster, A. Irles, S. O. Moch, P. Uwer and M. Vos, Eur. Phys. J. C 73 (2013), 2438 doi:10.1140/epjc/s10052-013-2438-2
  * [6] S. Dittmaier, P. Uwer and S. Weinzierl, Phys. Rev. Lett. 98 (2007), 262002 doi:10.1103/PhysRevLett.98.262002
  * [7] A. von Manteuffel and R. M. Schabinger, Phys. Lett. B 744 (2015), 101-104 doi:10.1016/j.physletb.2015.03.029
  * [8] T. Peraro, JHEP 12 (2016), 030 doi:10.1007/JHEP12(2016)030
  * [9] J. Gluza, K. Kajda and D. A. Kosower, Phys. Rev. D 83 (2011), 045012 doi:10.1103/PhysRevD.83.045012
  * [10] H. Ita, Phys. Rev. D 94 (2016) no.11, 116015 doi:10.1103/PhysRevD.94.116015
  * [11] K. J. Larsen and Y. Zhang, Phys. Rev. D 93 (2016) no.4, 041701 doi:10.1103/PhysRevD.93.041701
  * [12] Z. Wu, J. Boehm, R. Ma, H. Xu and Y. Zhang, Comput. Phys. Commun. 295 (2024), 108999 doi:10.1016/j.cpc.2023.108999
  * [13] T. Gehrmann, J. M. Henn and N. A. Lo Presti, Phys. Rev. Lett. 116 (2016) no.6, 062001 [erratum: Phys. Rev. Lett. 116 (2016) no.18, 189903] doi:10.1103/PhysRevLett.116.062001
  * [14] T. Gehrmann, J. M. Henn and N. A. Lo Presti, JHEP 10 (2018), 103 doi:10.1007/JHEP10(2018)103
  * [15] D. Chicherin and V. Sotnikov, JHEP 20 (2020), 167 doi:10.1007/JHEP12(2020)167
  * [16] D. Chicherin, V. Sotnikov and S. Zoia, JHEP 01 (2022), 096 doi:10.1007/JHEP01(2022)096
  * [17] S. Abreu, D. Chicherin, H. Ita, B. Page, V. Sotnikov, W. Tschernow and S. Zoia, Phys. Rev. Lett. 132 (2024) no.14, 141601 doi:10.1103/PhysRevLett.132.141601
  * [18] S. Badger, M. Czakon, H. B. Hartanto, R. Moodie, T. Peraro, R. Poncelet and S. Zoia, JHEP 10 (2023), 071 doi:10.1007/JHEP10(2023)071
  * [19] B. Agarwal, F. Buccioni, F. Devoto, G. Gambuti, A. von Manteuffel and L. Tancredi, Phys. Rev. D 109 (2024) no.9, 094025 doi:10.1103/PhysRevD.109.094025
  * [20] G. De Laurentis, H. Ita, M. Klinkert and V. Sotnikov, Phys. Rev. D 109 (2024) no.9, 094023 doi:10.1103/PhysRevD.109.094023
  * [21] G. De Laurentis, H. Ita and V. Sotnikov, Phys. Rev. D 109 (2024) no.9, 094024 doi:10.1103/PhysRevD.109.094024
  * [22] S. Badger, M. Becchetti, E. Chaubey, R. Marzucca and F. Sarandrea, JHEP 06 (2022), 066 doi:10.1007/JHEP06(2022)066
  * [23] S. Badger, M. Becchetti, E. Chaubey and R. Marzucca, JHEP 01 (2023), 156 doi:10.1007/JHEP01(2023)156
  * [24] S. Badger, M. Becchetti, N. Giraudo and S. Zoia, [arXiv:2404.12325 [hep-ph]].
  * [25] K. Melnikov and M. Schulze, Nucl. Phys. B 840 (2010), 129-159 doi:10.1016/j.nuclphysb.2010.07.003
  * [26] F. Moriello, JHEP 01 (2020), 150 doi:10.1007/JHEP01(2020)150
  * [27] M. Hidding, Comput. Phys. Commun. 269 (2021), 108125 doi:10.1016/j.cpc.2021.108125
  * [28] T. Peraro, JHEP 07 (2019), 031 doi:10.1007/JHEP07(2019)031