# Triangle Counting Through Cover-Edges

David A. Bader∗, Fuhuan Li, Anya Ganeshan+, Ahmet Gundogdu#, Jason Lew, Oliver
Alvarado Rodriguez, Zhihui Du  Department of Data Science
New Jersey Institute of Technology
Newark, New Jersey, USA
<EMAIL_ADDRESS>
###### Abstract

Counting and finding triangles in graphs is often used in real-world analytics
to characterize cohesiveness and identify communities in graphs. In this
paper, we propose the novel concept of a cover-edge set that can be used to
find triangles more efficiently. We use a breadth-first search (BFS) to
quickly generate a compact cover-edge set. Novel sequential and parallel
triangle counting algorithms are presented that employ cover-edge sets. The
sequential algorithm avoids unnecessary triangle-checking operations, and the
parallel algorithm is communication-efficient. The parallel algorithm can
asymptotically reduce communication on massive graphs such as from real social
networks and synthetic graphs from the Graph500 Benchmark. In our estimate
from massive-scale Graph500 graphs, our new parallel algorithm can reduce the
communication on a scale 36 graph by 1156x and on a scale 42 graph by 2368x.

###### Index Terms:

Graph Algorithms, Triangle Counting, Parallel Algorithms, High Performance
Data Analytics

${}^{+}$${}^{+}$footnotetext: Anya Ganeshan attends Bergen County Academies,
Hackensack, NJ${}^{\\#}$${}^{\\#}$footnotetext: Ahmet Gundogdu attends Paramus
High School, Paramus, NJ${}^{*}$${}^{*}$footnotetext: This research was
partially supported by NSF grant number CCF-2109988.

##  I Introduction

_Triangle counting_ [1] is a fundamental problem in graph analytics, which
involves finding the number of unique triangles in a graph. It plays a crucial
role in various graph analysis techniques such as clustering coefficients [2],
k-truss [3], and triangle centrality [4]. The significance of triangle
counting is evident in its application in high-performance computing
benchmarks like Graph500 [5] and the MIT/Amazon/IEEE Graph Challenge [6], as
well as in the design of future architecture systems (e.g., IARPA AGILE [7]).

Both sequential and parallel triangle counting algorithms have been studied
extensively since 1977 [8]. Latapy [9] provides a comprehensive overview of
sequential triangle counting and various finding algorithms. Existing
techniques, including list intersection, matrix multiplication, and subgraph
matching [10], are techniques used to count triangles.

To enhance the performance of triangle counting, Cohen [11] introduced a novel
map-reduce parallelization technique that generates _open wedges_ between
triples of vertices in the graph. It determines whether a closing edge exists
to complete a triangle, thus avoiding the redundant counting of the same
triangle while maintaining load balancing. Many parallel approaches for
triangle counting [12, 13] partition the sparse graph data structure across
multiple compute nodes and adopt the strategy of generating open wedges, which
are sent to other compute nodes to determine the presence of a closing edge.
Consequently, the communication time for these open wedges often dominates the
running time of parallel triangle counting.

In traditional edge-based triangle counting methods, all triangles are
identified by accumulating the sizes of intersections between pairs of
endpoints for each edge. _Direction-oriented_ approaches can avoid counting
the same triangle multiple times. However, in this paper, we propose a novel
approach that efficiently identifies all triangles using a reduced set of
edges known as a cover-edge set. By leveraging the cover-edge-based triangle
counting method, unnecessary edge checks can be skipped while ensuring that no
triangles are missed. This significantly reduces the number of computational
operations compared to existing methods. Furthermore, for distributed parallel
algorithms, the cover-edge-based method can greatly reduce overall
communication requirements. As a result, our proposed method offers improved
efficiency and scalability for triangle counting.

Our contributions include:

  * •

A novel concept, _Cover-Edge Set_ , is proposed to support efficient triangle
counting. The essential idea is that we can identify all triangles from a
significantly reduced cover-edge set instead of the complete edge set. A
simple breadth-first search (BFS) is used to orient the graph’s vertices into
levels and to generate the cover-edge set.

  * •

A novel triangle counting and finding algorithm, _CETC_ , is developed based
on the concept of _Cover-Edge Set_. _CETC_ runs in $\mathcal{O}(m\cdot
d_{\text{max}})$ time and $\mathcal{O}(n+m)$ space, where $d_{\text{max}}$ is
the maximal degree of a vertex $v\in V$.

  * •

A novel communication-efficient distributed parallel algorithm for triangle
counting and finding, _Comm-CETC_ , is also developed based on the concept of
_Cover-Edge Set_. _Comm-CETC_ can asymptotically reduce the communication to
improve total performance.

The remainder of the paper is organized as follows. Section II presents our
new approach for triangle counting. In Section III, we employ our idea in
distributed parallel triangle counting to reduce communication. Section V
discusses related work. Lastly, in Section VI, we conclude the paper.

##  II Cover-Edge Set Based Triangle Counting

###  II-A Notations and Basic Idea

Let $G=(V,E)$ be an undirected graph with $n=|V|$ vertices and $m=|E|$ edges.
A _triangle_ in the graph is a set of three vertices
$\\{v_{a},v_{b},v_{c}\\}\subseteq V$ such that
$\\{(v_{a},v_{b}),(v_{a},v_{c}),(v_{b},v_{c})\\}\subseteq E$. We will use
$N(v)=\\{u|u\in V\wedge((v,u)\in E)\\}$ to denote the _neighbor set_ of vertex
$v\in V$. The degree of vertex $v\in V$ is $d(v)=|N(v)|$, and $d_{\text{max}}$
is the maximal degree of a vertex in graph $G$.

######  Definition 1 (Cover-Edge and Cover-Edge Set).

For any edge $e$ of a triangle $\Delta$ in graph $G$, $e$ is referred to as a
cover-edge of $\Delta$. For a given graph $G$, an edge set $S\subseteq E$ is
called a cover-edge set if it contains at least one cover-edge for every
triangle in $G$.

Based on the given definition, it is evident that the entire edge set $E$ can
serve as a cover-edge set $S$ for graph $G$. However, our proposed method aims
to efficiently count all triangles using a smaller subset of edges instead of
$E$. Thus, the primary challenge lies in generating a compact cover-edge set,
which forms the initial problem to be addressed in our approach. Let
$k=|S|/|E|$. Our goal is to identify cover-edge sets with the smallest $k$. In
this paper, we propose using breadth-first search (BFS) to generate a compact
cover-edge set.

######  Definition 2 (BFS-Edge).

Let $r$ be the root vertex of an undirected graph $G$. The level $L(v)$ of a
vertex $v$ is defined as the shortest distance from $r$ to $v$ obtained
through a breadth-first search (BFS). From the BFS, we classify the edges into
three types:

  * •

_Tree-Edges_ : These edges belong to the BFS tree.

  * •

_Strut-Edges_ : These are non-tree edges with endpoints on two adjacent levels
in the BFS traversal.

  * •

_Horizontal-Edges_ : These are non-tree edges with endpoints on the same level
in the BFS traversal.

Figure 1: Example BFS tree of Graph $G$. The tree-edges are black, strut-edges
are blue, and horizontal-edges are red.

Fig. 1 gives an example of these different edge types.

######  Lemma 1.

Each triangle $\\{u,v,w\\}$ in a graph contains at least one horizontal-edge.

###### Proof.

(Proof by contradiction) A triangle is a path of length 3 that starts and ends
at the same vertex. Suppose there are no horizontal-edges in the triangle. In
that case, every edge in the path (i.e., a tree-edge or strut-edge) either
increases or decreases the level by one.

Since the path must end on the same level as the starting vertex, the number
of edges in the path that decrease the level must be equal to the number of
edges that increase the level. Consequently, the length of the path must be
even to maintain level parity. However, this contradicts the fact that a
triangle has an odd path length of 3.

Therefore, we conclude that there must be at least one horizontal-edge in
every triangle. ∎

######  Theorem 3 (Cover-Edge Set Generation).

All horizontal-edges form a valid cover-edge set.

###### Proof.

According to Definition 1, for any triangle $\Delta$ in graph $G$, we can
always find at least one horizontal-edge that serves as a cover-edge for
$\Delta$. Thus, the set of all horizontal-edges constitutes a cover-edge set.
∎

Therefore, we can construct a cover-edge set, denoted as _BFS-CES_ , by
selecting all the horizontal-edges obtained during a breadth-first search
(_BFS_). It is evident that _BFS-CES_ is a subset of $E$ and is typically much
smaller than the complete edge set $E$.

###  II-B Cover-Edge based Triangle Counting

In this subsection, we provide a comprehensive description of the process
involved in identifying all triangles using a cover-edge set generated through
a breadth-first search.

######  Lemma 2.

Each triangle $\\{u,v,w\\}$ must contain either one or three horizontal-edges.

###### Proof.

By referring to the proof of Lemma 1, we know that the path corresponding to
the triangle’s three edges consists of an even number of tree-edges and strut-
edges. This implies that there can be either 0 or 2 tree- or strut-edges
within each triangle.

In the case where there are 0 tree- or strut-edges, all three edges of the
triangle must be horizontal-edges. This is because the absence of tree- or
strut-edges implies that the entire path is composed of horizontal-edges.

In the case where there are 2 tree- or strut-edges, the triangle contains
exactly one horizontal-edge. This is because having two tree- or strut-edges
in the path means that there is one horizontal-edge connecting the remaining
two vertices.

Therefore, we conclude that each triangle $\\{u,v,w\\}$ must contain either
one or three horizontal-edges. ∎

Our triangle counting approach, described in Alg. 1, efficiently counts
triangles using a cover-edge set. In line 3, we initialize the counter $T$ to
0, which will store the total number of triangles. To generate the cover-edge
set, we perform a breadth-first search (BFS) starting from any unvisited
vertex, identifying the level ($L(v)$) of each vertex $v$ in its respective
component, as shown in lines 4 to 5. In lines 6 to 10 the algorithm iterates
over each edge, selecting the cover-set of horizontal edges $(u,v)$ in a
direction-oriented fashion in line 7. For each vertex $w$ in the intersection
of $u$ and $v$’s neighborhoods (line 8), we check the following two conditions
to determine if $(u,v,w)$ is a unique triangle to be counted (line 9). If
$L(u)\neq L(w)$ then the edge $(u,v)$ is the only horizontal-edge in the
triangle $(u,v,w)$. If $L(u)\equiv L(w)$, then the edge $(u,v)$ is one of
three horizontal-edges in the triangle $(u,v,w)$. To ensure uniqueness, the
algorithm then checks the added constraint that $v<w$. If the constraints are
satisfied, we increment the triangle counter $T$ in line 10.

This approach effectively counts the triangles in the graph while avoiding
redundant counting.

Algorithm 1 CETC:Cover-Edge Triangle Counting

1:Graph $G=(V,E)$

2:Triangle Count $T$

3:$T\leftarrow 0$

4:$\forall v\in V$

5: if $v$ unvisited, then BFS($G$, $v$)

6:$\forall(u,v)\in E$

7: if $(L(u)\equiv L(v))\land(u<v)$ $\triangleright$ $(u,v)$ is horizontal

8: $\forall w\in N(u)\cap N(v)$

9: if $(L(u)\neq L(w))\lor\left((L(u)\equiv L(w))\land(v<w)\right)$ then

10: $T\leftarrow T+1$

######  Theorem 4 (Correctness).

Alg. 1 can accurately count all triangles in a graph $G$.

###### Proof.

Lemma 2 establishes that a triangle in the graph falls into one of two cases:
1) the two endpoint vertices of the horizontal-edge are on the same level
while the apex vertex is on a different level, or 2) all three vertices of the
triangle are at the same level.

Consider a triangle $\\{v_{a},v_{b},v_{c}\\}$ in $G$. Without loss of
generality, assume that $(v_{a},v_{b})$ is a horizontal-edge, implying
$L(v_{a})\equiv L(v_{b})$. Let $v_{c}$ be the apex vertex. The two cases can
be distinguished as follows:

For the first case, each triangle is uniquely defined by a horizontal-edge and
an apex vertex from the common neighbors of the horizontal-edge’s endpoint
vertices. Whenever Alg. 1 identifies such a triangle
$\\{v_{a},v_{b},v_{c}\\}$, it increments the total triangle count $T$ by 1.

In the second case, where all three vertices are at the same level
($L(v_{c})\equiv L(v_{a})\equiv L(v_{b})$), Alg. 1 ensures that $T$ is
increased by 1 only when $v_{a}<v_{b}<v_{c}$. This condition ensures that
triangle $\\{v_{a},v_{b},v_{c}\\}$ is counted only once, preventing triple-
counting and ensuring the correctness of the triangle count.

Hence, Alg. 1 is proven to accurately count all triangles in the graph $G$. ∎

The time complexity of Alg. 1 can be analyzed as follows. The computation of
breadth-first search, including determining the level of each vertex and
marking horizontal-edges, requires $\mathcal{O}(n+m)$ time.

Since there are at most $\mathcal{O}(m)$ horizontal-edges, finding the common
neighbors of each horizontal-edge individually can be done in
$\mathcal{O}(d_{\text{max}})$ time. Here, $d_{\text{max}}$ represents the
maximal degree of a vertex in the graph.

Therefore, the overall time complexity of Alg. 1 is $\mathcal{O}(m\cdot
d_{\text{max}})$.

##  III Communication Efficient Triangle Counting Algorithm

This section presents our communication-efficient parallel algorithm for
counting triangles in massive graphs on a $p$-processor distributed-memory
parallel computer. We will take advantage of the concept of _Cover-Edge Set_
to significantly improve the communication performance of our triangle
counting method. Since distributed triangle counting is communication-bound
[12], this algorithm is expected to improve the overall running time. The
input graph $G$ is stored in a compressed sparse row (CSR) format. The
vertices are partitioned non-uniformly to the $p$ processors such that each
processor stores approximately $2m/p$ edge endpoints. This graph input follows
the format used by the majority of parallel graph algorithm implementations
and benchmarks such as Graph500 and Graph Challenge.

###  III-A Parallel Algorithm Description

Our communication-efficient parallel algorithm (see Alg. 2) is based on the
same cover-edge approach proposed in section II. The binary operator $\oplus$
used in line 11 is bitwise exclusive OR (XOR).

Algorithm 2 Comm-CETC: Communication Efficient Triangle Counting

1:Graph $G=(V,E)$

2:Triangle Count $T$

3:Run parallel BFS($G$) and build partial cover-edge set $S_{i}$ on $p_{i}$

4:For all $p_{i},i\in\\{0\ldots p-1\\}$ in parallel do:

5: $t_{i}\leftarrow 0$

6: $\forall(u,v)\in S_{i}$ with $u<v$ on $p_{i}$

7: $\forall w\in V_{i}$ such that $w\in N(u),N(v)$

8: if $(L(u)\neq L(w))\lor\left((L(u)\equiv L(w))\land(v<w)\right)$ then

9: $t_{i}=t_{i}+1$

10: For $j\leftarrow 1$ to $p-1$ do:

11: Processors $i$ and $i\oplus j$ swap edge sets $S_{i}$ and $S_{j}$.

12: $\forall(u,v)\in S_{j}$ with $u<v$ on $p_{i}$

13: $\forall w\in V_{i}$ such that $w\in N(u),N(v)$

14: if $(L(u)\neq L(w))\lor\left((L(u)\equiv L(w))\land(v<w)\right)$ then

15: $t_{i}=t_{i}+1$

16:$T\leftarrow$ Reduce$(t_{i},+)$

Similar to the baseline _CETC_ algorithm, the cover-edge set
$S=\cup_{i=0}^{p-1}S_{i}$ is determined in line 3 by labeling the horizontal
edges from a parallel BFS.

Each processor runs lines 4 to 15 in parallel that consists of two main
substeps. Local triangles are counted in lines 6 to 9 and a total exchange of
cover-edges between each pair of processors to count triangles is performed in
lines 10 to 15. Note at the end of each iteration of the _for_ loop, processor
$p_{i}$ can discard the cover-edge set $S_{j}$. In lines 7 and 13, processor
$p_{i}$ determines for each cover edge $(u,v)$ all the apex vertices $w$ held
locally that are adjacent to both $u$ and $v$. The logic for counting
triangles in lines 8 and 14 is similar to Alg. 1 as to only count unique
triangles. Finally, a reduction operation in line 16 calculates the total
number of triangles by accumulating the triangle counters across the system,
i.e., $T=\sum_{i=0}^{p-1}t_{i}$.

###  III-B Cost Analysis

####  III-B1 Space

In addition to the input graph data structure, an additional bit is needed per
edge (for marking a horizontal-edge) and $\mathcal{O}\left(\lceil\log
D\rceil\right)$ bits per vertex to store its level, where $D$ is the diameter
of the graph. This is a total of at most $m+n\lceil\log D\rceil$ bits across
the $p$ processors. Preserving the graph requires additional
$\mathcal{O}\left(n+m\right)$ space for the graph.

####  III-B2 Compute

The BFS costs $\mathcal{O}\left((n+m)/p\right)$ [14], the modified neighbor
sets take $\mathcal{O}\left(m/p\right)$. The search corresponding to one
cover-edge in a vertex’s adjacency list takes at most
$\mathcal{O}\left(\log(d_{max})\right)$ time using binary search, and only
$\mathcal{O}\left(1\right)$ expected time using a hash table. Let $d_{i}$ be
the degree of vertex $v_{i}$ where $0\leq i<n$. Searching $km$ edges in all
vertices’ adjacency lists takes
$\mathcal{O}(km\sum_{i=0}^{n-1}\log(d_{i}))=\mathcal{O}(km\log(\Pi_{i=0}^{n-1}d_{i}))$
time. Since $\sum_{i=0}^{n-1}d_{i}=2m$, we know that
$\log(\Pi_{i=0}^{n-1}d_{i})$ reaches its maximum value when $d_{i}=2m/n$ for
$0\leq i<n$. Thus,
$\mathcal{O}(km\log(\Pi_{i=0}^{n-1}d_{i}))\leq\mathcal{O}(km\log((2m/n)^{n}))\leq\mathcal{O}(kmn\log(2n^{2}/n))=\mathcal{O}(mn\log(n))$.

####  III-B3 Total Communication

In our analysis of communication cost for BFS, we measure the total
communication volume independent of the number of processors. Thus, this is a
conservative overestimate of communication since a fraction (e.g., $1/p$) of
accesses will be on the same compute node versus message traffic between
nodes. At the same time, we do not consider the savings from overlapping with
the computation cost.

The cost of the breadth-first search is $m$ edge traversals with $\lceil\log
D\rceil+3\lceil\log n\rceil$ bits communicated per edge traversal for the
level information, pair of vertex ids, and vertex degree, yielding
$m\cdot(\lceil\log D\rceil+3\lceil\log n\rceil)$ bits for the BFS.
Transferring $km$ horizontal-edges requires $kmp\lceil\log n\rceil$ bits,
where $p$ is the number of processors. The final reduction to find the total
number of triangles requires $(p-1)\lceil\log n\rceil$ bits.

Hence, the total communication volume is $m\cdot(\lceil\log
D\rceil+3\lceil\log n\rceil)+kmp\lceil\log n\rceil+(p-1)\lceil\log
n\rceil=m\cdot(\lceil\log D\rceil+(kp+3)\lceil\log n\rceil)+(p-1)\lceil\log
n\rceil$ bits. Hence, since the word size is $\Theta$$\left(\log n\right)$ and
$D\leq n$, the communication is $\mathcal{O}\left(pm\right)$ words.

##  IV Communication Analysis on Real and Synthetic Graphs

TABLE I: Communication costs for real and synthetic graph. The synthetic graphs are Graph500 RMAT graphs of scale 36 and 42. The column ‘Previous’ represents the communication volume of the best prior parallel algorithms [15, 16, 17], that use wedge-checking based algorithms and ‘This paper’ represents the communication cost of our new approach. ‘Reduction’ represents the communication reduction between these two, and thus, the expected speedup of the parallel algorithm. Entries in _italics_ are estimated values. Graph | n | m | # Triangles | # Wedges | $k$ | $p$ | Previous | This paper | Reduction
---|---|---|---|---|---|---|---|---|---
ca-GrQc | 5242 | 14484 | 48260 | 165798 | 0.522 | 4 | 526KB | 122KB | 4.31
ca-HepTh | 9877 | 25973 | 28339 | 277389 | 0.423 | 4 | 948KB | 218KB | 4.35
as-caida20071105 | 26475 | 53381 | 36365 | 776895 | 0.225 | 4 | 2.78MB | 401KB | 7.10
facebook_combined | 4039 | 88234 | 1612010 | 17051688 | 0.914 | 4 | 48.8MB | 893KB | 56.0
ca-CondMat | 23133 | 93439 | 173361 | 1567373 | 0.511 | 4 | 5.61MB | 897KB | 6.40
ca-HepPh | 12008 | 118489 | 3358499 | 5081984 | 0.621 | 4 | 17.0MB | 1.13MB | 15.1
email-Enron | 36692 | 183831 | 727044 | 5933045 | 0.478 | 4 | 22.6MB | 1.79MB | 12.7
ca-AstroPh | 18772 | 198050 | 1351441 | 8451765 | 0.667 | 4 | 30.2MB | 2.08MB | 14.6
loc-brightkite_edges | 58228 | 214078 | 494728 | 6956250 | 0.441 | 4 | 26.5MB | 2.02MB | 20.4
soc-Epinions1 | 75879 | 405740 | 1624481 | 21377935 | 0.498 | 4 | 86.7MB | 4.25MB | 10.7
amazon0601 | 403394 | 2443408 | 3986507 | 96348699 | 0.529 | 8 | 436MB | 40.9MB | 10.7
com-Youtube | 1134890 | 2987624 | 3056386 | 209811585 | 0.347 | 8 | 1.03GB | 44.3MB | 23.7
RMAT-36 | 68719476736 | 1099511627776 | _1.2E+14_ | _2.73E+16_ | _0.311_ | 128 | 218PB | _192TB_ | _1156_
RMAT-42 | 4398046511104 | 70368744177664 | _1.3E+16_ | _5.79E+18_ | _0.260_ | 256 | 52.8EB | _22.8PB_ | _2368_

In this section, we analyze the performance of the parallel triangle counting
algorithm on both real and synthetic graphs. We implemented our new triangle
counting algorithm using Python to accurately compute the exact communication
volume and determine an analytic model based on the size of the graph and
number of processors, and the ratio or percentage ($k$) of cover-edges from
the BFS. The results given in Table I are exact communication volumes from our
new algorithm on all of the graphs except the two large RMAT graphs where we
compute the communication volume from the validated analytic model. For the
comparison with prior approaches [15, 16, 17], we estimate the communication
volume from the number of wedges which is exact for all graphs other than the
last two large RMAT graphs where we estimate the number of wedges using graph
theory.

For the real graphs, we find the actual value of $k$, the percentage of graph
edges that are cover-edges, for an arbitrary breadth-first search, and set the
number $p$ of processors to a reasonable number given the size of the graph.
For the synthetic graphs, we use large Graph500 RMAT graphs [18] with
parameters $a=0.57$, $b=0.19$, $c=0.19$, and $d=0.05$, for scale 36 and 42
with $n=2^{\mbox{scale}}$ and $m=16n$, similar with the IARPA AGILE benchmark
graphs, and set $p$ according to estimates of potential system sizes with
sufficient memory to hold these large instances.

For comparison, most prior parallel algorithms for triangle counting operate
on the graph as follows. A parallel loop over the vertices $v\in V$ produces
all 2-paths (_wedges_) where $(v,v_{1}),(v,v_{2})\in E$ and (w.l.o.g.)
$v_{1}<v_{2}$. The processor that produces this wedge will send an open wedge
query message containing the vertex ids of $v_{1}$ and $v_{2}$ to the
processor that owns vertex $v_{1}$. If the consumer processor that receives
this query message finds an edge $(v_{1},v_{2})\in E$, then a local triangle
counter is incremented. After producers and consumers complete all work, a
global reduction over the $p$ triangle counts computes the total number of
triangles in $G$.

###  IV-A Graph500 RMAT Graphs

Figure 2: Estimate of $k$ using an exponential model, based on observations of
$k$ for RMAT graph scale 6 to 23 graphs.

For the large Graph500 RMAT graphs, the number of triangles is estimated from
our model based on the number of triangles found in RMAT graphs up to scale 29
in the literature [19, 20, 18, 21]. The fitting equation is
$\mbox{\\#Triangles}=77.422n^{1.125}$ with $R^{2}=1.0$, where $n$ is the total
number of vertices. The number of triangles estimated for scale 36 and 42 RMAT
graphs are $1.20\times 10^{14}$ and $1.30\times 10^{16}$, respectively.

We estimate the number of wedges for the scale 36 and 42 Graph500 RMAT graphs
based on the theorem given by Seshadhri _et al._ in [22]. According to their
formula, we can estimate the expected number of vertices $N(d)$ for a given
out-degree $d$. The number of wedges that can be formed by vertices with such
a degree is calculated as $\binom{d}{2}\times N(d)$, where $\binom{d}{2}$
means choosing two from $d$.

By summing all such wedges generated from the minimum ($e\ln n$) to the
maximum degree ($\sqrt{n}$), which is the assumption of the formula, we can
approximate the total number of wedges in the given graph, where $n$ is the
total number of vertices. This is a conservative estimate because it only
considers the out-degree instead of the sum of out and in-degrees. Employing
the formula, we calculate the number of wedges to be $2.73\times 10^{16}$ for
scale 36 and $5.8\times 10^{18}$ for scale 42. With $2\log n$ bits/wedge, the
total volume of wedge checks is 218PB and 52.8EB for RMAT graphs of scales 36
and 42, respectively111Throughout this paper, a petabyte (PB) is $2^{50}$
bytes and an exabyte (EB) is $2^{60}$ bytes..

Beamer  _et al._ [23] find a typical BFS on a scale 27 Graph500 RMAT graph has
7 levels, so 4 bits is a reasonable estimate for $\log D$ in our analyses of
scale 36 and 42 graphs.

The methodology for estimating the value of $k$ for RMAT graphs is as follows.
RMAT graphs from scale 6 to 23 are generated, and the exact value of $k$ is
determined for each by counting the horizontal-edges after a breadth-first
search. The data fit to an exponential model
$k=1.1773e^{-0.036\cdot\mbox{scale}}$ with very high $R^{2}=0.9956$ (see Fig.
2). For scale 36, $k$ is estimated to be 0.311 and for scale 42, $k$ is
estimated to be 0.260.

In our new approach for scale 36, where the communication cost is
$m\cdot(\lceil\log D\rceil+(kp+3)\lceil\log n\rceil)+(p-1)\lceil\log n\rceil$
bits. With $\lceil\log D\rceil=4$, and assuming $p=128$ processors, we have a
total communication volume of 192TB, for a communication reduction of
$1156\times$. For scale 42, and assuming $p=256$ processors, we estimate the
communication of our new triangle counting algorithm as 22.8PB, for a
communication reduction of $2368\times$.

##  V Related Work

###  V-A Sequential Algorithms

The naïve approach for triangle counting uses brute-force: find all the
triplets $\\{v_{a},v_{b},v_{c}\\}$, that is, permutations of three arbitrary
vertices in the graph, and check whether each edge in the triplet exists. The
time complexity is $\Omega$$\left(n^{3}\right)$. Latapy [9] and Schank and
Wagner [24] provide surveys of faster sequential algorithms. Triangle counting
generally can be formulated as three kinds of problems: set (list)
intersection, matrix multiplication and subgraph (cycle) query.

The three main intersection-based triangle counting algorithms are: 1) the
node-iterator algorithm iterates over all vertices and tests for each pair of
neighbors whether they are connected by an edge, 2) the edge-iterator
algorithm iterates over all edges and searches for common neighbors of the two
endpoints of each edge, and 3) the forward algorithm is a refinement of the
edge-iterator algorithm that computes the intersection of a subset of
neighborhoods by using an orientation of the graph. The time complexity of
node-iterator and edge-iterator are both $\mathcal{O}\left(m\cdot
d_{\mbox{max}}\right)$ and the forward algorithm is
$\mathcal{O}\left(m^{\frac{3}{2}}\right)$, which has significantly better
performance when $d_{\mbox{max}}\gg\sqrt{m}$ [9].

When performing the intersection of two lists, the commonly used techniques
are _merge-path_ , _binary search_ and _hashing-based_ algorithms.

Merge-path algorithms (e.g., [25, 26]) use two pointers to scan through
neighbor lists of two endpoints from beginning to end in order to find the
list intersection. During the scan, the pointer that points to a smaller value
will be incremented. A triangle is enumerated if both pointers are incremented
(i.e., they both point to the same vertex). Binary-search algorithms (e.g.,
[27, 19]) organize the longer list as a binary tree and use the shorter list
as search keys. For each search key, it descends through the binary-search
tree in order to find the equal entry, which is a triangle. Hashing-based
algorithms (e.g., [10, 26]) construct a hash table for one list and use the
other list as search keys to find the common elements in the hash table. The
hash table is used here to find the intersection of two adjacency lists, so it
is not necessary to sort all the adjacency lists to find all the triangles.
The running time is proportional to the size of the two adjacency lists.

Triangle counting using matrix multiplication [28] relies on a linear algebra
formulation for triangle counting. This approach can be optimized [29] using
matrix decomposition by decomposing $A$ into lower and upper triangular
matrices $L$ and $U$, and then computing $(L\times U)\odot L$, or $(L\times
L)\odot L$ to determine the number of triangles. The binary operator $\odot$
denotes the Hadamard product.

A subgraph-based approach for triangle counting searches for all occurrences
of a query graph, which is a triangle, in the input graph. Wang and Owens [30]
use breadth-first search to update the subgraph matching approach by pruning
more invalid vertices based on neighborhood encoding information, and using
optimizations like $k$-step look-ahead to reduce unwanted intermediate
results. Alon _et al._ [31] proposed a $\mathcal{O}\left(m^{1.41}\right)$
algorithm to find length 3 cycles (triangle) in a graph, which is an
improvement over the Itai and Rodeh sequential
$\mathcal{O}\left(m^{\frac{3}{2}}\right)$ algorithm [8].

###  V-B Parallel Algorithms

Map-reduce is a standard platform for large scale distributed computation.
Cohen [11] first demonstrated the capability of map-reduce to solve triangle
counting in an approach that generates _open wedges_ between triples of
vertices in the graph and determines if a closing edge exists that completes a
triangle. Suri _et al._ [32] implemented triangle counting using map-reduce
that ranks vertices by degree and distributes them across hosts. Pearce [12]
developed an algorithm that is based on creating an augmented degree-ordered
directed graph, where the original undirected edges are directed from low-
degree to high degree, and implemented this approach in the distributed
asynchronous graph processing framework HavoqGT. DistTC [19] is a distributed
triangle counting implementation for multiple machines that uses mirror proxy
on each partition to eliminate almost all the inner-host communication.
TriCore [27] partitions the graph held in a compressed-sparse row (CSR) data
structure for multiple GPUs and uses stream buffers to load edge lists from
CPU memory to GPU memory on-the-fly and then uses binary search to find the
intersection. Hu _et al._ [33] employed a “copy-synchronize-search” pattern to
improve the parallel threads efficiency of GPU and mixed the computing and
memory intensive workloads together to improve the resource efficiency. Pandey
_et al._ [10] employed an vertex-centric hash-based design to scale triangle
counting to over 1,000 GPUs. TriC [13] exploits the vertex-based distributed
triangle counting and sends vertices rather than edges (vertex pairs), and
then the remote processor could translate the sequence of vertex IDs to
correct combination of vertices as edges to reduce communication. An
enhancement is then presented to TriC [34] that added a user-defined buffer to
improve the flexibility of controlling the memory usage for large data sets
and used a probabilistic data structure to optimize the edge lookups by
trading off the accuracy. Strausz _et al._ [35] use CLaMPI, a software caching
layer that caches data retrieved through MPI remote memory access operations,
to reduce the overall communication cost. Zeng _et al._ [36] proposed a
triangle counting algorithm that adaptively selects vertex-parallel and edge-
parallel paradigm.

Panduranga _et al._ [37] and Dolev _et al._ [15]’s work focused on the
communication cost. Compared with our work, there are two major differences.
First, they use the number of communication rounds to measure the total
communication with a bandwidth restriction. However, we use the total volume
of messages to evaluate the communication. Second, they are probabilistic
algorithms, but our algorithm is a deterministic algorithm (Dolev  _et al._
[15] also contains a deterministic version). Probabilistic methods cannot be
used under scenarios with an exact result requirement. Uhl [38, 17] also
focuses on reducing the communication cost of triangle counting. The paper’s
basic idea is only requiring communication for counting triangles consisting
of cut edges. If the partition generates many cut edges, the proposed method
cannot significantly reduce communication. In contrast, our method identifies
a subset of the total set of edges independent of the partitioning and only
transfers this smaller set of edges during the triangle counting to
significantly reduce the total communication.

##  VI Conclusions

In this paper, we present novel sequential and parallel algorithms for
counting and finding triangles in graphs based on a compact cover-edge set.
The parallel algorithm is the first communication-efficient triangle counting
algorithm by exploiting BFS horizontal-edges to significantly reduce the
communication volume on massive graphs of practical interest. Our approach
uses the breadth-first search to significantly reduce the number of edges
examined and minimize the communication required for triangle checking. The
parallel algorithm achieves an order of magnitude or more reduction of
communication volume for large graphs as communication is the main bottleneck
for triangle counting on distributed memory systems.

##  VII Reproducibility

The sequential triangle counting source code and the Python code for
determining the communication volume of the parallel algorithm are open source
and available on GitHub at https://github.com/Bader-Research/triangle-
counting. The input graphs are from the Stanford Network Analysis Project
(SNAP) available from http://snap.stanford.edu/.

## References

  * [1] M. Al Hasan and V. S. Dave, “Triangle counting in large networks: a review,” _Wiley Interdisciplinary Reviews: Data Mining and Knowledge Discovery_ , vol. 8, no. 2, p. e1226, 2018.
  * [2] D. J. Watts and S. H. Strogatz, “Collective dynamics of ‘small-world’ networks,” _Nature_ , vol. 393, no. 6684, pp. 440–442, 1998.
  * [3] J. Cohen, “Trusses: Cohesive subgraphs for social network analysis,” _National Security Agency Technical Report_ , vol. 16, no. 3.1, 2008.
  * [4] P. Burkhardt, “Triangle centrality,” _CoRR_ , vol. abs/2105.00110, 2021. [Online]. Available: https://arxiv.org/abs/2105.00110
  * [5] Graph 500 Steering Committee, “The Graph500 benchmark,” 2010. [Online]. Available: https://www.graph500.org
  * [6] S. Samsi, V. Gadepally, M. Hurley, M. Jones, E. Kao, S. Mohindra, P. Monticciolo, A. Reuther, S. Smith, W. Song, D. Staheli, and J. Kepner, “GraphChallenge.org: Raising the bar on graph analytic performance,” in _2018 IEEE High Performance Extreme Computing Conference (HPEC)_ , 2018, pp. 1–7.
  * [7] W. Harrod, “Advanced graphical intelligence logical computing environment (AGILE),” https://www.iarpa.gov/images/PropsersDayPDFs/AGILE/AGILE_Proposers_Day_sm.pdf, 2020, accessed: 2022–05-24.
  * [8] A. Itai and M. Rodeh, “Finding a minimum circuit in a graph,” in _Proceedings of the Ninth Annual ACM Symposium on Theory of Computing_ , ser. STOC ’77. New York, NY, USA: Association for Computing Machinery, 1977, p. 1–10.
  * [9] M. Latapy, “Main-memory triangle computations for very large (sparse (power-law)) graphs,” _Theoretical Computer Science_ , vol. 407, no. 1, pp. 458–473, 2008.
  * [10] S. Pandey, Z. Wang, S. Zhong, C. Tian, B. Zheng, X. Li, L. Li, A. Hoisie, C. Ding, D. Li, and H. Liu, “Trust: Triangle counting reloaded on GPUs,” _IEEE Transactions on Parallel and Distributed Systems_ , vol. 32, no. 11, pp. 2646–2660, 2021.
  * [11] J. Cohen, “Graph twiddling in a MapReduce world,” _Computing in Science & Engineering_, vol. 11, no. 4, pp. 29–41, 2009.
  * [12] R. Pearce, “Triangle counting for scale-free graphs at scale in distributed memory,” in _IEEE High Performance Extreme Computing Conference (HPEC)_ , 2017, pp. 1–4.
  * [13] S. Ghosh and M. Halappanavar, “TriC: Distributed-memory triangle counting by exploiting the graph structure,” in _IEEE High Performance Extreme Computing Conference (HPEC)_ , 2020, pp. 1–6.
  * [14] T. H. Cormen, C. E. Leiserson, R. L. Rivest, and C. Stein, _Introduction to Algorithms, Fourth Edition_. MIT Press, 2022.
  * [15] D. Dolev, C. Lenzen, and S. Peled, “’Tri, tri again’: finding triangles and small subgraphs in a distributed setting,” in _International Symposium on Distributed Computing_. Springer, 2012, pp. 195–209.
  * [16] R. Pearce and G. Sanders, “K-truss decomposition for scale-free graphs at scale in distributed memory,” in _2018 IEEE high performance extreme computing conference (HPEC)_. IEEE, 2018, pp. 1–6.
  * [17] P. Sanders and T. N. Uhl, “Engineering a distributed-memory triangle counting algorithm,” _arXiv preprint arXiv:2302.11443_ , 2023.
  * [18] D. Chakrabarti, Y. Zhan, and C. Faloutsos, “R-MAT: A recursive model for graph mining,” in _Proceedings of the 2004 SIAM International Conference on Data Mining_. SIAM, 2004, pp. 442–446.
  * [19] L. Hoang, V. Jatala, X. Chen, U. Agarwal, R. Dathathri, G. Gill, and K. Pingali, “DistTC: High performance distributed triangle counting,” in _IEEE High Performance Extreme Computing Conference (HPEC)_ , 2019, pp. 1–7.
  * [20] I. Giechaskiel, G. Panagopoulos, and E. Yoneki, “PDTL: Parallel and distributed triangle listing for massive graphs,” in _2015 44th International Conference on Parallel Processing_. IEEE, 2015, pp. 370–379.
  * [21] P. Burkhardt, “Graphing trillions of triangles,” _Information Visualization_ , vol. 16, no. 3, pp. 157–166, 2017.
  * [22] C. Seshadhri, A. Pinar, and T. G. Kolda, “A hitchhiker’s guide to choosing parameters of stochastic kronecker graphs,” _CoRR, abs/1102.5046_ , vol. 1, 2011.
  * [23] S. Beamer, K. Asanovic, and D. Patterson, “Searching for a parent instead of fighting over children: A fast breadth-first search implementation for Graph500,” _EECS Department, University of California, Berkeley, Tech. Rep. UCB/EECS-2011-117_ , 2011.
  * [24] T. Schank and D. Wagner, “Finding, counting and listing all triangles in large graphs, an experimental study,” in _International workshop on experimental and efficient algorithms_. Springer, 2005, pp. 606–609.
  * [25] O. Green, P. Yalamanchili, and L.-M. Munguía, “Fast triangle counting on the GPU,” in _Proceedings of the 4th Workshop on Irregular Applications: Architectures and Algorithms_ , 2014, pp. 1–8.
  * [26] J. Shun and K. Tangwongsan, “Multicore triangle computations without tuning,” in _2015 IEEE 31st International Conference on Data Engineering_. IEEE, 2015, pp. 149–160.
  * [27] Y. Hu, H. Liu, and H. H. Huang, “TriCore: Parallel triangle counting on GPUs,” in _SC18: International Conference for High Performance Computing, Networking, Storage and Analysis_. IEEE, 2018, pp. 171–182.
  * [28] A. Azad, A. Buluç, and J. Gilbert, “Parallel triangle counting and enumeration using matrix algebra,” in _2015 IEEE International Parallel and Distributed Processing Symposium Workshop_. IEEE, 2015, pp. 804–811.
  * [29] S. Acer, A. Yaşar, S. Rajamanickam, M. Wolf, and Ü. V. Catalyürek, “Scalable triangle counting on distributed-memory systems,” in _IEEE High Performance Extreme Computing Conference (HPEC)_ , 2019, pp. 1–5.
  * [30] L. Wang and J. D. Owens, “Fast BFS-based triangle counting on GPUs,” in _IEEE High Performance Extreme Computing Conference (HPEC)_ , 2019, pp. 1–6.
  * [31] N. Alon, R. Yuster, and U. Zwick, “Finding and counting given length cycles,” _Algorithmica_ , vol. 17, no. 3, pp. 209–223, 1997.
  * [32] S. Suri and S. Vassilvitskii, “Counting triangles and the curse of the last reducer,” in _Proceedings of the 20th International Conference on World Wide Web_ , 2011, pp. 607–614.
  * [33] L. Hu, L. Zou, and Y. Liu, “Accelerating triangle counting on GPU,” in _Proceedings of the 2021 International Conference on Management of Data_ , 2021, pp. 736–748.
  * [34] S. Ghosh, “Improved distributed-memory triangle counting by exploiting the graph structure,” in _IEEE High Performance Extreme Computing Conference (HPEC)_ , 2022, pp. 1–6.
  * [35] A. Strausz, F. Vella, S. Di Girolamo, M. Besta, and T. Hoefler, “Asynchronous distributed-memory triangle counting and LCC with RMA caching,” in _36th IEEE International Parallel and Distributed Processing Symposium, IPDPS_. IEEE, 2022, pp. 291–301.
  * [36] L. Zeng, K. Yang, H. Cai, J. Zhou, R. Zhao, and X. Chen, “HTC: Hybrid vertex-parallel and edge-parallel triangle counting,” in _IEEE High Performance Extreme Computing Conference (HPEC)_ , 2022, pp. 1–7.
  * [37] G. Pandurangan, P. Robinson, and M. Scquizzato, “On the distributed complexity of large-scale graph computations,” _ACM Transactions on Parallel Computing (TOPC)_ , vol. 8, no. 2, pp. 1–28, 2021.
  * [38] T. N. Uhl, “Communication efficient triangle counting,” Ph.D. dissertation, Karlsruhe Institute of Technology, 2021.


# Getting in shape with minimal energy. A variational principle for
protohaloes

Marcello Musso1 and Ravi K. Sheth,2,3
1Departamento de Física Fundamental and IUFFyM, Universidad de Salamanca,
E-37008 Salamanca, Spain
2Center for Particle Cosmology, University of Pennsylvania, 209 S. 33rd St.,
Philadelphia, PA 19104, USA
3The Abdus Salam International Center for Theoretical Physics, Strada
Costiera, 11, Trieste 34151, Italy E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

In analytical models of structure formation, protohalos are routinely assumed
to be peaks of the smoothed initial density field, with the smoothing filter
being spherically symmetric. This works reasonably well for identifying a
protohalo’s center of mass, but not its shape. To provide a more realistic
description of protohalo boundaries, one must go beyond the spherical picture.
We suggest that this can be done by looking for regions of fixed volume, but
arbitrary shape, that minimize the enclosed energy. Such regions are
surrounded by surfaces over which (a slightly modified version of) the
gravitational potential is constant. We show that these equipotential surfaces
provide an excellent description of protohalo shapes, orientations and
associated torques.

###### keywords:

large-scale structure of Universe

††pubyear: 2022††pagerange: Getting in shape with minimal energy. A
variational principle for protohaloes–Getting in shape with minimal energy. A
variational principle for protohaloes

##  1 Introduction

The key ingredient of most analytical models of structure formation is the
initial (over-)density field $\delta$ smoothed with a suitable spherical
filter. The spherical collapse model (Gunn & Gott, 1972) establishes a mapping
between the initial enclosed density and the collapse time, leading to a
critical initial mean density $\delta_{c}$ that a region must have for it to
collapse at the present time. Motivated by this, the earliest models tried to
characterize protohalo regions as spheres of a prescribed initial mean density
(Press & Schechter, 1974). Subsequent work noticed that if a sphere of
critical mean density is enclosed within a larger one of the same density,
then the protohalo should be identified with the larger volume. Hence, to
avoid double counting, a protohalo should be assigned the mass contained in
the _largest_ sphere of critical mean density centered on that location (the
so-called _cloud-in-cloud_ problem of Bond et al., 1991).

Other models, focusing more on the locations of protohalos rather than their
formation times, emphasized that collapse does not take place around random
locations (Sheth et al., 2001). For instance, spheres that are denser than
their neighbors collapse faster, so it is reasonable to identify protohalo
centers with the centers of spheres that are local maxima of the mean density
(the _peaks theory_ of Bardeen et al., 1986). The two approaches have been
combined together, leading to the description of protohalos as peaks of the
smoothed density field that have critical height, and are not contained within
other peaks of the same height but larger smoothing scale (the _excursion set
peaks_ approach: Bond, 1989; Appel & Jones, 1990; Paranjape & Sheth, 2012;
Paranjape et al., 2013).

Considerable attention has also been devoted to the precise value of the
critical value for collapse and the functional form of the filter. The former,
leading to more realistic thresholds that are both stochastic (i.e. dependent
on variables other than the mean overdensity) and scale dependent, and can
include the anisotropy of the environment (Bond & Myers, 1996; Sheth et al.,
2001; Chiueh & Lee, 2001). The latter, in order to avoid the mathematical
drawback of divergences arising with the usual top-hat filter (e.g. the _CUSP_
approach of Salvador-Solé & Manrique, 2021, and references therein), but also
to deal with more physically motivated quantities such as the protohalo’s
energy, rather than mass density (Musso & Sheth, 2021).

Although halos in simulations are not usually spherically symmetric (Bonamigo
et al., 2015, and references therein), and the protohalo patches from which
they formed are not either (Despali et al., 2013; Ludlow et al., 2014;
Nikakhtar et al., 2022), to the best of our knowledge, most analytical studies
to date have used a spherical volume to model protohalos (but see Betancort-
Rijo & López-Corredoira, 2002; Lam & Sheth, 2008; Borzyszkowski et al., 2017,
for initially nonspherical regions that evolve to become spheres). Because
protohalos are not spherical, this can, at best, identify the center of mass.
Although suitable averages with this filter can offer clues to the actual
protohalo shape (Bond & Myers, 1996; Monaco, 1997) and its subsequent
evolution (Rossi et al., 2011; Salvador-Solé & Manrique, 2021), their accuracy
is limited (Ludlow et al., 2014), in part because they do not provide a
physical handle to describe the protohalo boundary.

In this Letter, we assume no _a priori_ shape, but suggest that both the
location _and_ the shape of a protohalo are associated with a region of
minimal initial energy, whose boundary is asurface on which (a slightly
modified form of) the gravitational potential perturbation remains constant.
This equipotential surface is the locus of points where the initial infall
velocity (the component of the velocity directed towards the center of mass)
is constant. Section 2 describes the minimum energy principle in detail.
Section 3 tests this _ansatz_ against protohaloes in an N-body simulation. A
final section summarizes.

##  2 Conceptual foundation

This section makes two conceptual points: (i) the shape that best describes a
protohalo of a given mass is the one that, while preserving the volume,
minimises the enclosed energy, and (ii) the boundary of this special region is
an ‘equipotential’ surface.

###  2.1 The miniumum energy principle

In an earlier paper on this subject, Musso & Sheth (2021) showed that the
center of mass of a protohalo typically coincides with the center of a sphere
which contains the same mass $M$ and maximises the _potential energy
overdensity_

$\displaystyle\epsilon$
$\displaystyle\equiv\frac{5}{MR_{I}^{2}}\int_{V}\mathrm{d}\bm{r}\,\rho(\bm{r})\,(\bm{r}-\bm{r}_{\mathrm{cm}})\cdot[\nabla\phi(\bm{r})-\nabla\phi_{\mathrm{cm}}]$
(1)

where $\phi$ is the potential perturbation, normalized so that
$\nabla^{2}\phi=\delta$, $\bm{r}_{\mathrm{cm}}$ and
$-\nabla\phi_{\mathrm{cm}}$ are the center of mass position and acceleration
of the volume $V$, and $R_{I}$ is its _inertial radius_ , defined as

$R_{I}^{2}\equiv\frac{5}{3M}\int_{V}\mathrm{d}\bm{r}\,\rho(\bm{r})\,|\bm{r}-\bm{r}_{\mathrm{cm}}|^{2}\,.$
(2)

Although the definitions above are general, they simplify in the initial
conditions, when $\rho(\bm{r})\simeq\bar{\rho}$, $M\simeq\bar{\rho}V$ and
$\nabla\phi\simeq-\bm{v}/fDH$, where $\bm{v}$ is the initial velocity and
$fHD\equiv(d\ln D/d\ln a)(d\ln a/dt)D=dD/dt$, where $D$ is the linear theory
growth rate. Furthermore, either $\nabla\phi_{\mathrm{cm}}$ or
$\bm{r}_{\mathrm{cm}}$ in equation (1) may be omitted, since the mass averages
of both $\bm{r}-\bm{r}_{\mathrm{cm}}$ and of
$\nabla\phi-\nabla\phi_{\mathrm{cm}}$ vanish by definition.

For a homogeneous sphere, $R_{I}$ coincides with the sphere’s geometrical
radius $R$ (defined by $V=4\pi R^{3}/3$), and $\epsilon$ with its mean matter
overdensity $\delta_{R}$. In general, however, they are different. For spheres
centered on protohalos, which tend to be more centrally overdense, usually
$R_{I}<R$ (with the difference starting at first order in perturbations) and
$\epsilon>\delta_{R}$. At early times, the potential energy overdensity is
related to the total energy $E$ (kinetic plus potential) within $V$ as

$E=-\frac{4\pi
G}{3}\bar{\rho}\,MR_{I}^{2}\,\epsilon=-\frac{GM^{2}}{R}\,\frac{R_{I}^{2}}{R^{2}}\,\epsilon\,,$
(3)

This is $(5\epsilon/3)(R_{I}/R)^{2}$ times the potential energy of a
homogeneous sphere containing the same mass as $V$, that is
$-(3/5)\,(GM^{2}/R)$.

The dynamical role of $\epsilon$ is twofold. First, it represents for $R_{I}$
what $\delta_{R}$ represents for $R$, that is, it determines the time at which
$R_{I}$ turns around, according to the spherical collapse model: $R_{I}$ of a
patch with larger $\epsilon$ shrinks faster, and that patch collapses sooner.
Thus, the inertial radius of a volume $V$ that is a local maximum in
$\epsilon$ has a smaller collapse time than all neighbouring volumes of the
same mass. Secondly, for a sphere, the spatial gradient $\nabla\epsilon$ is
(proportional to) the sphere’s dipole moment. If this gradient vanishes, then
so does the dipole term of the gravitational potential at the surface. Hence,
if higher order multipoles are neglected, the initial velocities of the
particles at the surface all approximately converge to the center. For these
reasons, local maxima of the smoothed energy overdensity field are excellent
candidates for the centers of mass of initial patches that will evolve
dynamically into high density regions: i.e., protohaloes.

So far, like most of the literature to date, we have implicitly assumed $V$ to
be a sphere. However, different shapes of the same $V$ will generically have
different values of $\epsilon$. In particular, if one deforms a sphere’s
boundary so as to include larger values of
$(\bm{r}-\bm{r}_{\mathrm{cm}})\cdot\nabla\phi$ and exclude smaller ones, while
keeping the total mass fixed, then the new value of $\epsilon$ will be larger.
Thus, if one allows for such non-spherical deformations around an initially
spherical peak, one can always find a value of $\epsilon$ that exceeds that of
the spherical peak. However, any such deformation also changes $R_{I}^{2}$;
sufficiently large deformations may increase the denominator of $\epsilon$
more than the numerator. There is therefore a well defined, non-spherical
surface that maximises the value of $\epsilon$. As the inertial radius $R_{I}$
of the region within this surface shrinks even faster than the sphere’s, it is
tempting to take this (non-spherical) region as the best candidate for the
protohalo patch.

Figure 1: Isosurfaces of the infall potential $\mathcal{V}/R_{I}^{2}$, defined
in equation (6), centered on a protohalo patch. The values of $\epsilon$,
$\bm{r}_{\mathrm{cm}}$ and $\nabla\phi_{\mathrm{cm}}$ are the same (the ones
of the protohalo) for all surfaces.

###  2.2 Equipotential surfaces

We now want to characterize the surface that maximises $\epsilon$. Since we
just argued that there is always such a surface surrounding a spherical peak,
it should be sufficient to find that surface for which the infinitesimal
variation $\delta\epsilon$ vanishes for _any_ further infinitesimal
deformation that preserves the mass.

The variation of the mass $M$ under infinitesimal deformations of the volume
$V$ is given by the surface integral

$\delta M=\int_{S}\mathrm{d}\bm{S}\cdot\delta\bm{\lambda}\,\rho(\bm{r})$ (4)

where $S$ is the boundary surface of $V$,
$\mathrm{d}\bm{S}=\mathrm{d}S\,\bm{\hat{n}}$ is the normal to the surface
element, and $\delta\bm{\lambda}$ is the infinitesimal deformation vector
mapping each point on $S$ to its image. One must have $\delta M=0$ for the
deformation to preserve the mass.

Figure 2:  Comparison of protohalo patches (blue) with their equipotential
regions (orange; i.e., the region enclosed by isocontours of $\mathcal{V}$).
In all cases, particles are shown in the $x$-$y$ plane, projected in the $z$
direction. Top right corner of each panel gives the protohalo mass (in
$h^{-1}M_{\odot}$).

To compute the total variation of $\epsilon$, one must not only vary the
integral over $V$ appearing in equation (1) explicitly, but also the one in
$R_{I}^{2}$ (varying those in $\bm{r}_{\mathrm{cm}}$ and
$\nabla\phi_{\mathrm{cm}}$ gives a null net contribution). Using Leibnitz’s
rule, one gets

$\delta\epsilon=-\frac{5}{MR_{I}^{2}}\int_{S}\mathrm{d}\bm{S}\cdot\delta\bm{\lambda}\,\rho(\bm{r})\mathcal{V}(\bm{r})\,,$
(5)

where

$\mathcal{V}(\bm{r})\equiv(\bm{r}-\bm{r}_{\mathrm{cm}})\cdot\bigg{[}-(\nabla\phi-\nabla\phi_{\mathrm{cm}})+\frac{\epsilon}{3}(\bm{r}-\bm{r}_{\mathrm{cm}})\bigg{]}$
(6)

is the radial component of the acceleration relative to the center of mass,
minus its mass weighed average over $V$. It has the dimensions of, and for the
monopole term it actually equals, the potential relative to the center of mass
minus its average111Near the protohalo boundary
$\nabla\phi-\nabla\phi_{\mathrm{cm}}\sim\delta_{R}(\bm{r}-\bm{r}_{\mathrm{cm}})/3$
and $|\bm{r}-\bm{r}_{\mathrm{cm}}|\sim R\sim R_{I}$, so we expect
$\mathcal{V}/R_{I}^{2}\sim\epsilon-\delta_{R}$ to be of order unity (e.g Musso
& Sheth, 2021).. With a slightly stretched terminology, we dub
$\mathcal{V}(\bm{r})$ the _infall potential_.

To make $\epsilon$ stationary as per our original aim, we need to find the
surface $S$ such that $\delta\epsilon=0$ under _any_ deformation that
preserves the mass. By comparing equations (4) and (5), it follows immediately
that this must be the surface over which $\mathcal{V}(\bm{r})$ is constant,
since in this case one has $\delta\epsilon\propto\delta M=0$. That is, $S$ is
an iso-surface of $\mathcal{V}$; and, to the extent that $\mathcal{V}$ can be
related to the gravitational potential, an equipotential surface.

There will be of course many of these surfaces, usually nested in one another,
with different values of $\mathcal{V}$. Only one of them will contain exactly
the mass $M$. Increasing (or decreasing) $M$ will select a different surface,
usually enclosing (or enclosed by) the previous one. Each surface will also
have its own value of $\epsilon$, which is maximal at fixed $M$ but not
necessarily as $M$ varies. Just like in plain excursion sets, this value can
be associated to the time at which the mass $M$ is assembled. The set of these
nested iso-surfaces will thus provide a description of the halo’s mass
accretion history.

If $(\nabla\phi-\nabla\phi_{\mathrm{cm}})$ had only the monopole term, then
$\mathcal{V}(\bm{r})$ would be spherically symmetric, and all its iso-surfaces
would be spheres centered on $\bm{r}_{\mathrm{cm}}$. In general, the presence
of higher multipoles deforms them. Since
$-(\bm{r}-\bm{r}_{\mathrm{cm}})\cdot(\nabla\phi-\nabla\phi_{\mathrm{cm}})$ is
proportional to the infall velocity, and normally becomes more negative as the
distance grows, the boundary is pushed farther where the infall velocity is
larger than its spherically symmetric value. Particles that are infalling
faster will make it into the halo from further away. Thus, the minimum energy
principle describes the asphericity of protohalos as a response to the
anisotropy of the gravitational infall.

##  3 Measurements in simulations

We now test our minimum energy _ansatz_ using protohaloes from the Flora
simulation, the largest box in the SBARBINE suite (Despali et al., 2016). The
simulation evolved $1024^{3}$ dark matter particles each of mass $6.35\times
10^{11}h^{-1}M_{\odot}$ in a periodic cube of side
$L_{\mathrm{box}}=2h^{-1}$Gpc with a Planck13 background cosmology:
$\Omega_{m}=0.307$, $\Omega_{\Lambda}=0.693$, $\sigma_{8}=0.829$ and
$h=0.677$.

Our halo set contains 5378 haloes identified at $z=0$ using a Spherical
Overdensity (SO) halo finder with threshold of $319\times$ the background
density. The set includes all the 1378 haloes more massive than
$10^{15}h^{-1}M_{\odot}$, 2000 randomly chosen haloes with masses between
$10^{14}$ and $10^{15}h^{-1}M_{\odot}$, and 2000 randomly chosen haloes with
masses between $4\times 10^{13}$ and $10^{14}h^{-1}M_{\odot}$. We refer to the
patch defined by each halo’s particles in the initial conditions as the
‘protohalo’. For a subset of haloes (the 1387 most massive halos, and 1372 in
the intermediate mass bin) we also have counterparts identified with an
Ellipsoidal Overdensity (EO) halo finder (typically slightly less round, and
about 10 percent more massive). We used this control set to check that our
results do not depend strongly on the halo finder. To emphasize the connection
to energy peaks, we present results using $\sigma_{02}^{2}(R)\equiv\int
dk\,k^{2}P_{\rm Lin}(k)\,W_{2}^{2}(kR)$, where
$R\equiv(3M/4\pi\bar{\rho})^{1/3}$ and $W_{2}(x)\equiv 15\,j_{2}(x)/x^{2}$;
increasing $M$ decreases $\sigma_{02}$.

Figure 3:  _Top:_ Fraction of actual protohalo particles included in the
equipotential region. _Bottom:_ Alignment between the inertia tensor of the
protohalo and of the equipotential region (equation (9)). Both are shown as a
function of $\sigma_{02}$ (largest masses on the left)

We measure each protohalo’s center of mass position and velocity,
$\bm{r}_{\mathrm{cm}}$ and $\bm{v}_{\mathrm{cm}}$, by averaging over all its
particles. We then estimate its potential energy overdensity _tensor_ as

$\hat{\epsilon}_{ij}\equiv-3\,\frac{\sum_{n}[(\bm{r}-\bm{r}_{\mathrm{cm}})_{i}(\bm{v}-\bm{v}_{\mathrm{cm}})_{j}/fDH]_{n}}{\sum_{n}[(\bm{r}-\bm{r}_{\mathrm{cm}})\cdot(\bm{r}-\bm{r}_{\mathrm{cm}})]_{n}}\,,$
(7)

where $D$ is the $\Lambda$CDM density perturbation growth factor,
$f=\mathrm{d}\ln D/\mathrm{d}\ln a$ so that $\dot{D}=fDH$, and $n$ runs over
all $N$ particles in the protohalo. Next, with these values of
$\bm{r}_{\mathrm{cm}}$, $\bm{v}_{\mathrm{cm}}$ and
$\hat{\epsilon}=\mathrm{tr}(\bm{\mathsf{\hat{\epsilon}}})$, we construct

$\hat{\mathcal{V}}(\bm{r})\equiv(\bm{r}-\bm{r}_{\mathrm{cm}})\cdot\bigg{[}\frac{(\bm{v}-\bm{v}_{\mathrm{cm}})}{fDH}+\frac{\hat{\epsilon}}{3}(\bm{r}-\bm{r}_{\mathrm{cm}})\bigg{]}$
(8)

for each particle in a cube of side $3R$ centered on $\bm{r}_{\mathrm{cm}}$,
and select the $N$ particles having the lowest values. This singles out a
region that has the same mass as the protohalo, and is bound by an iso-surface
of $\mathcal{V}$, which we call the protohalo’s _equipotential region_.

Fig. 1 shows, for illustration, a few iso-surfaces of $\mathcal{V}$ (scaled by
$R_{I}^{2}$ to make it dimensionless) for a protohalo in the largest mass bin,
all relative to the _same_ values of $\bm{r}_{\mathrm{cm}}$,
$\bm{v}_{\mathrm{cm}}$ and $\hat{\epsilon}$ (the ones of the protohalo). We
refer to the one that encloses the _same_ number of particles (hence, the same
mass) as the protohalo as the ‘boundary’ of its equipotential region. As
argued in the previous section, this region has $\mathcal{V}/R_{I}^{2}\sim 1$
and the maximal value of $\epsilon$ for the given mass, and is the most
natural prediction for the protohalo boundary. Deep inside the region, the
equipotential surfaces break up into disconnected regions (not shown), which
we believe encodes the assembly history. To make this connection stronger,
however, one should not use the same $\bm{r}_{\mathrm{cm}}$,
$\bm{v}_{\mathrm{cm}}$ and $\hat{\epsilon}$ for the whole set, but recompute
them for each surface at every step.

Figure 2 compares the protohalo and equipotential regions for a few randomly
selected haloes in the three mass bins. Smaller blue dots represent protohalo
particles, and larger orange circles the ones of the equipotential region,
plotted in the $x$-$y$ plane and projected along the $z$ direction. By visual
inspection one can appreciate the excellent agreement between the two. It is
remarkable how our approach has captured strongly non-spherical features like
the lobes and protruding arms of most protohaloes.

Figure 4: _Top panel:_ ellipticity of protohaloes (blue dots) vs.
equipotential regions (orange). _Bottom:_ ratio of the two ellipticities. The
agreement is excellent, although equipotential surfaces are slightly more
elliptical. Figure 5: Reduced torque magnitudes (equation (11)) of actual
protohaloes vs. equipotential surfaces as a function of $\sigma_{02}$ (largest
masses on the left). _Top panel:_ torque of protohaloes (blue dots) and of
equiptotential surfaces (orange dots). _Bottom:_ Logarithm of the ratio of the
two magnitudes.

For a more quantitative comparison, we quantify the overlap of the two regions
in four different ways:
(_i_) the fraction of protohalo particles in the equipotential region (shown
in the top panel of Figure 3);
(_ii_) the matrix cosine between the inertia tensors
$\bm{\mathsf{T}}_{\mathrm{H}}$ and $\bm{\mathsf{T}}_{\mathrm{EQ}}$ of each
protohalo and equipotential region,

$\cos(\psi)\equiv\mathrm{Tr}(\bm{\mathsf{T}}_{\mathrm{H}}\cdot\bm{\mathsf{T}}_{\mathrm{EQ}})/\big{[}\mathrm{Tr}\big{(}\bm{\mathsf{T}}_{\mathrm{H}}^{2}\big{)}\mathrm{Tr}\big{(}\bm{\mathsf{T}}_{\mathrm{EQ}}^{2}\big{)}\big{]}^{1/2}\,$
(9)

(shown in Figure 3, bottom), where

$T_{ij}=\sum_{n}[(\bm{r}-\bm{r}_{\mathrm{cm}})_{i}(\bm{r}-\bm{r}_{\mathrm{cm}})_{j}]_{n}\,;$
(10)

(_iii_) comparison of ellipticities $e_{\rm H}$ and $e_{\rm EQ}$, where
$e\equiv\sqrt{1-a_{3}/a_{1}}$ and $a_{3}$ and $a_{1}$ are the smallest and
largest eigenvalues of each inertia tensor (Figure 4); and
(_iv_) comparison of the _reduced_ torques $\bm{\tilde{\tau}}_{\mathrm{H}}$
and $\bm{\tilde{\tau}}_{\mathrm{EQ}}$, defined as

$\tilde{\tau}_{k}\equiv\frac{\tau_{k}}{\mathrm{tr}(\bm{\mathsf{T}})}=-\frac{1}{3}\,\varepsilon_{kij}\,\hat{\epsilon}_{ij}$
(11)

with $\hat{\epsilon}_{ij}$ given in equation (7) (Figure 5). For _(i-iii)_ ,
the differences between the SO and EO halo sets were hardly noticeable, so we
only show the (larger) SO sample. However, the net torque on a sphere
vanishes, so the estimated torque is quite sensitive to small changes in the
halo (and hence protohalo) boundary. So, for _(iv)_ , we show results for the
EO sample.

Overall, protohalos tend to be slightly less elliptical than equipotential
surfaces ($e_{\rm H}<e_{\rm EQ}$ in bottom panel of Figure 4) but there is
otherwise excellent agreement between the two. As a result, equipotential
regions are more torqued (ratios in bottom panel of Figure 5 are slightly
greater than unity; these ratios are slightly less than 1 for SO haloes). This
small discrepancy is nearly mass independent.

Figure 6:  Same as Figure 2, but for low mass objects that are outliers ($\leq
40\%$ of overlap) in the top panel of Figure 3. Clearly, equipotential regions
(orange symbols) do not identify the actual protohalo patches (blue).

To explore this further, Figure 6 shows a few of the objects with the smallest
values (less than 0.4) in the top panel of Figure 3. Whereas protohalo
particles (blue) define a single blob, the equipotential regions are sometimes
split into two or more components. Clearly, our approach fails to describe
these objects, suggesting that corrections beyond linear order in
perturbations, which our approach ignores, can be important. While it is
certainly interesting to identify these corrections, they really only matter
for a small fraction of the full halo sample, so we believe it is fair to
conclude that, for the vast majority of halos, our approach identifies not
just the protohalo center of mass but also its shape.

##  4 Conclusions

We have demonstrated that protohaloes are very well described as regions that
maximize the energy overdensity $\epsilon$, bounded by surfaces of a suitably
defined (equation 6) infall potential (Figure 1). This approach captures both
the centers of mass of protohaloes and their boundaries (Figures 2–5).

By connecting protohalo shapes to the infall pattern, our approach provides a
framework not just for modeling halo shapes, but the anisotropy of infall
velocity around clusters and the cross-correlation between dynamical and weak
gravitational lensing mass estimates for clusters as well (e.g. White et al.,
2010), as well as how these correlate with the larger scale environment.

Going beyond a spherical model is crucial for predicting the initial torque.
Our minimum energy principle is therefore also suitable for making analytical
predictions of the angular momentum of a protohalo patch. This prediction uses
the energy tensor $\epsilon_{ij}$ (equations 7 and 11). In future work, we
will compare and contrast this with the predictions of Tidal Torque Theory
(Doroshkevich, 1970; White, 1984; Lee & Pen, 2001; Porciani et al., 2002;
Cadiou et al., 2021).

Although our approach proved to be robust to changes in how exactly halos are
identified in the evolved field (e.g. SO vs EO halo finders), we believe that
it can be most naturally compared with the objects found by the ‘boosted
potential’ halo finder recently proposed by Stücker et al. (2021). The two
approaches are in fact very similar, with the important difference that ours
(_i_) derives equipotential surfaces from the energy minimisation principle,
and (_ii_) holds in the initial conditions, and therefore provides a natural
framework for making analytical _predictions_ of halo statistics. Moreover, as
we noted, our equipotential approach encodes information about the assembly
history of each object (Figure 1 and associated discussion); since it is
similar in spirit to how halo substructures are identified in the boosted
potential approach, a comparison of our predictions with such measurements
should lead to interesting results, to be compared for instance with Cadiou et
al. (2020).

## Acknowledgements

Thanks to Giulia Despali for sharing (and helping with) the Flora simulation,
to Corentin Cadiou, Oliver Hahn, Dmitri Pogosyan for helpful comments, and to
the organisers and other participants of the KITP Cosmic Web workshop. Our
visit to KITP was supported by the National Science Foundation under
PHY-1748958.

## Data Availability

The data of the Flora simulation, and the post-processing quantities used in
this work, can be shared on reasonable request to the authors.

## References

  * Appel & Jones (1990) Appel L., Jones B., 1990, MNRAS, 245, 522
  * Bardeen et al. (1986) Bardeen J. M., Bond J. R., Kaiser N., Szalay A. S., 1986, Astrophys. J., 304, 15
  * Betancort-Rijo & López-Corredoira (2002) Betancort-Rijo J., López-Corredoira M., 2002, ApJ, 566, 623
  * Bonamigo et al. (2015) Bonamigo M., Despali G., Limousin M., Angulo R., Giocoli C., Soucail G., 2015, MNRAS, 449, 3171
  * Bond (1989) Bond J. R., 1989, in Lake Louise Winter Institute: Frontiers in Physics - From Colliders to Cosmology Lake Louise, Alberta, Canada, February 19-25, 1989\. pp 182–235
  * Bond & Myers (1996) Bond J. R., Myers S. T., 1996, Astrophys.J.Supp., 103, 1
  * Bond et al. (1991) Bond J., Cole S., Efstathiou G., Kaiser N., 1991, Astrophys.J., 379, 440
  * Borzyszkowski et al. (2017) Borzyszkowski M., Porciani C., Romano-Diaz E., Garaldi E., 2017, Mon. Not. Roy. Astron. Soc., 469, 594
  * Cadiou et al. (2020) Cadiou C., Pichon C., Codis S., Musso M., Pogosyan D., Dubois Y., Cardoso J. F., Prunet S., 2020, Mon. Not. Roy. Astron. Soc., 496, 4787
  * Cadiou et al. (2021) Cadiou C., Pontzen A., Peiris H. V., 2021, Mon. Not. Roy. Astron. Soc., 502, 5480
  * Chiueh & Lee (2001) Chiueh T., Lee J., 2001, ApJ, 555, 83
  * Despali et al. (2013) Despali G., Tormen G., Sheth R. K., 2013, MNRAS, 431, 1
  * Despali et al. (2016) Despali G., Giocoli C., Angulo R. E., Tormen G., Sheth R. K., Baso G., Moscardini L., 2016, MNRAS, 456, 2486
  * Doroshkevich (1970) Doroshkevich A. G., 1970, Astrofizika, 6, 581
  * Gunn & Gott (1972) Gunn J. E., Gott J. Richard I., 1972, Astrophys.J., 176, 1
  * Lam & Sheth (2008) Lam T. Y., Sheth R. K., 2008, MNRAS, 386, 407
  * Lee & Pen (2001) Lee J., Pen U.-L., 2001, ApJ, 555, 106
  * Ludlow et al. (2014) Ludlow A. D., Borzyszkowski M., Porciani C., 2014, Mon. Not. Roy. Astron. Soc., 445, 4110
  * Monaco (1997) Monaco P., 1997, MNRAS, 287, 753
  * Musso & Sheth (2021) Musso M., Sheth R. K., 2021, Mon. Not. Roy. Astron. Soc., 508, 3634
  * Nikakhtar et al. (2022) Nikakhtar F., Sheth R. K., Lévy B., Mohayaee R., 2022, Phys. Rev. Lett., 129, 251101
  * Paranjape & Sheth (2012) Paranjape A., Sheth R. K., 2012, MNRAS, 426, 2789
  * Paranjape et al. (2013) Paranjape A., Sheth R. K., Desjacques V., 2013, MNRAS, 431, 1503
  * Porciani et al. (2002) Porciani C., Dekel A., Hoffman Y., 2002, MNRAS, 332, 339
  * Press & Schechter (1974) Press W. H., Schechter P., 1974, Astrophys. J., 187, 425
  * Rossi et al. (2011) Rossi G., Sheth R. K., Tormen G., 2011, MNRAS, 416, 248
  * Salvador-Solé & Manrique (2021) Salvador-Solé E., Manrique A., 2021, ApJ, 914, 141
  * Sheth et al. (2001) Sheth R. K., Mo H., Tormen G., 2001, MNRAS, 323, 1
  * Stücker et al. (2021) Stücker J., Angulo R. E., Busch P., 2021, Mon. Not. Roy. Astron. Soc., 508, 5196
  * White (1984) White S. D. M., 1984, ApJ, 286, 38
  * White et al. (2010) White M., Cohn J. D., Smit R., 2010, MNRAS, 408, 1818