# Branching random walks with regularly varying perturbations

Krzysztof Kowalski111University of Wrocław, Poland. Email:
<EMAIL_ADDRESS>Supported by the National Science Center,
Poland (OPUS, grant number 2019/33/B/ST1/00207)
MSC2020 subject classification: 60J80, 60F05, 60F15
keywords: branching random walk, perturbations, maximal position, central
limit theorem, strong limit, derivative martingale, additive martingale,
heavy-tailed distributions

###### Abstract

We consider a modification of classical branching random walk, where we add
i.i.d. perturbations to the positions of the particles in each generation. In
this model, which was introduced and studied by Bandyopadhyay and Ghosh
(2023), perturbations take form $\frac{1}{\theta}\log\frac{X}{E}$, where
$\theta$ is a positive parameter, $X$ has arbitrary distribution $\mu$ and $E$
is exponential with parameter 1, independent of $X$. Working under finite mean
assumption for $\mu$, they proved almost sure convergence of the rightmost
position to a constant limit, and identified the weak centered asymptotics
when $\theta$ does not exceed certain critical parameter $\theta_{0}$. This
paper complements their work by providing weak centered asymptotics for the
case when $\theta>\theta_{0}$ and extending the results to $\mu$ with
regularly varying tails. We prove almost sure convergence of the rightmost
position and identify the appropriate centering for the weak convergence,
which is of form $\alpha n+c\log n$, with constants $\alpha$, $c$ depending on
the ratio of $\theta$ and $\theta_{0}$. We describe the limiting distribution
and provide explicitly the constants appearing in the centering.

##  1 Introduction

Branching random walk on the real line, abbreviated as BRW, is constructed as
follows. The process starts with a single particle placed at 0. Given a point
process ${\mathcal{Z}}=\sum_{k=1}^{N}\delta_{\xi_{k}}$ on $\mathbb{R}$, where
$N$ is a random variable on $\mathbb{N}_{0}$, the original particle at time 1
dies and gives birth to $N$ particles positioned according to ${\mathcal{Z}}$.
These particles are called the first generation of the process. At time 2,
each of these particles reproduces independently and has offspring with
positions relative to their parent’s position given by an independent copy of
${\mathcal{Z}}$.

The process continues infinitely. As a result we obtain a marked tree
$\left(S,\mathbb{T}\right)$, where the tree $\mathbb{T}$ is the set of all
particles equipped in the natural tree structure, and $S_{v}$ is the position
of a given particle $v\in\mathbb{T}$. We write $|v|$ for the generation of
$v$. For a BRW with displacements given by $\xi$, let
$R_{n}=\sup_{|v|=n}S_{v}$ denote the position of the most right particle at
time $n$. There is a rich literature concerning the asymptotic behaviour of
$R_{n}$. In 1976 Biggins [5] proved under minimal assumptions the law of large
numbers for $R_{n}$, i.e. $\frac{R_{n}}{n}$ converges almost surely to a
constant $\alpha$. The corresponding limit theorem was proved by Aïdékon [1]
in 2013, who showed that $R_{n}-\alpha n+c\log n$ converges in distribution to
a random shift of the Gumbel distribution. We refer to Shi [12] for an
extensive description of recent results on branching random walks.

In this paper we consider a perturbed branching random walk $S^{*}$, which is
a modification of $S$, in which we add a random perturbation to the position
of every particle, i.e.

$S^{*}_{v}=S_{v}+X_{v}$

where $\\{X_{v}\\}_{v\in\mathbb{T}}$ are i.i.d. random variables independent
of $S$. Note that the perturbation added to the position of a vertex
$v\in\mathbb{T}$ does not influence the positions of its offspring, which
explains that the process is sometimes called last progeny modified branching
random walk. In this paper we study the model introduced by Bandyopadhyay and
Ghosh in [3], where the perturbations have the form

$X_{v}(\theta)=\frac{1}{\theta}\log\frac{Y_{v}}{E_{v}}$

for a given positive real number $\theta$, and $\\{Y_{v}\\}_{v\in\mathbb{T}}$
that are independent positive random variables with distribution $\mu$,
independent of $\\{E_{v}\\}_{v\in\mathbb{T}}$, which are i.i.d. with
distribution $\text{Exp}{\left(1\right)}$. This model was further studied in
the context of large deviations [10] and inhomogeneous time BRW [4]. Main
motivation for considering it comes from the connection between the supremum
of the perturbed BRW $R^{*}_{n}(\theta)=\sup_{|v|=n}S^{*}_{v}$ and random
weighted sums. More precisely, Theorem 3.6 in [3] states, that

$\displaystyle\theta R^{*}_{n}(\theta)\stackrel{{\scriptstyle d}}{{=}}\log
Y_{n}(\theta)-\log E$ (1.1)

where $Y_{n}(\theta)=\sum_{|v|=n}e^{\theta S_{v}}Y_{v}$ and $E$ is exponential
with parameter 1, independent of $Y_{n}(\theta)$. We will sometimes write
$R^{*}_{n}=R^{*}_{n}(\theta)$ and $Y_{n}=Y_{n}(\theta)$ if the parameter is
clear from the context. The asymptotics of $R^{*}_{n}$ will be related very
closely to the behaviour of $Y_{n}$. As it turns out, the sequence $Y_{n}$ is
quite well described in the literature, see e.g. [7] and [8]. Properties of
$R_{n}^{*}$ are obviously determined by the parameter $\theta$. It turns out
that one needs to control its position with respect to the critical parameter
$\theta_{0}$ defined as

$\theta_{0}=\inf\left\\{\theta>0:\nu(\theta)=\theta\nu^{\prime}(\theta)\right\\}$

where

$\nu(\theta)=\log\mathbb{E}\left[\sum_{i=1}^{N}e^{\theta\xi_{i}}\right]$

is the log-Laplace transform of $\mathcal{Z}$, and
$\nu^{\prime}(\theta)=e^{-\nu(\theta)}\mathbb{E}\left[\sum_{i=1}^{N}\xi_{i}e^{\theta\xi_{i}}\right]$.
Note that $\nu$ does not have to be differentiable at $\theta$ for this
quantity to exist, and that in general $\theta_{0}$ may be infinite.

In [3] branching random walks with such perturbations where studied in the
case when $\mu$ has finite mean. In particular, the authors proved that

$\displaystyle\frac{R_{n}^{*}(\theta)}{n}\xrightarrow[n\rightarrow\infty]{a.s.}\begin{cases}\frac{\nu\left(\theta\right)}{\theta}&\theta<\theta_{0}\\\
\frac{\nu(\theta_{0})}{\theta_{0}}&\theta\geq\theta_{0}\end{cases}$

and identified weak centered asymptotics for $\theta\leq\theta_{0}$. However,
the result for $\theta>\theta_{0}$ was only obtained for the degenerated
perturbations with $\mu=\delta_{1}$. The goal of this paper is to complete the
results from [3] by providing the missing weak centered asymptotics for the so
called above the boundary case, and to extend them to $\mu$ with infinite
mean, with special focus on distributions with regularly varying tails.

##  2 Main results

Let $\gamma\in\left(0,1\right)$. Our main assumption for $\mu$ is as follows:

$x^{\gamma}\left(1-F(x)\right)\xrightarrow[x\rightarrow+\infty]{}c_{+}>0$ (H)

where $F$ is the probability distribution function of $\mu$.

This assumption tells us that $\mu$ belongs to the domain of attraction of a
stable law with characteristic function

$\tilde{g}\left(t\right)=e^{-k|t|^{\gamma}\left(1-i\tan{\left(\frac{\pi\gamma}{2}\right)}\text{sign}t\right)}$
(2.1)

where

$k=\frac{\pi
c_{+}}{2\Gamma\left(\gamma\right)\sin\left(\pi\gamma/2\right)}>0.$

Furthermore, it means that if $Y$ has distribution $\mu$, then
$\mathbb{E}\left[Y^{\gamma}\right]=\infty$, but
$\mathbb{E}\left[Y^{r}\right]<\infty$ for any $r\in\left(0,\gamma\right)$.
$(\textbf{H})$ will be assumed in most of the paper, however the result for
the above the boundary case will be stated under more general assumption.

We assume for the rest of the paper that
$\text{supp}\left(\mu\right)\subset\mathbb{R}_{+}$, the system survives with
probability 1 (i.e. $\mathbb{P}(N=0)=0$), and $\nu$ is finite on some open
interval $I$, with $0\in I$. Since $\theta_{0}$ is finite, the last assumption
guarantees, by convexity of $\nu$, that $\nu$ is differentiable on
$\left(-s,\theta_{0}\right)$ for some $s>0$, and has the left derivative at
$\theta_{0}$. One can also characterize $\theta_{0}$ as the argument
minimizing $\frac{\nu(\theta)}{\theta}$ over $\theta>0$. Throughout this
paper, existence of finite $\theta_{0}$ will only be assumed when necessary.

As proved in [5], if $\theta_{0}$ is finite, then

$\displaystyle\frac{R_{n}}{n}\xrightarrow[n\rightarrow\infty]{}\frac{\nu(\theta_{0})}{\theta_{0}},\quad\mbox{a.s.}$
(2.2)

For $\theta$ such that $\nu(\theta)<\infty$, let

$\displaystyle W_{n}(\theta)=e^{-n\nu(\theta)}\sum_{|v|=n}e^{\theta S_{v}}.$
(2.3)

$W_{n}(\theta)$ is called the additive martingale associated with $S$. We
denote $W_{n}=W_{n}(\theta_{0})$. Note that as a positive martingale
$W_{n}(\theta)$ converges almost surely to some finite limit. If
$\nu^{\prime}(\theta)<\infty$, then Biggins martingale convergence theorem [6]
states, that the almost sure limit of $W_{n}(\theta)$ is non-degenerate if and
only if $\nu^{\prime}(\theta)<\nu(\theta)/\theta$ and

$\mathbb{E}\left[W_{1}(\theta)\log_{+}W_{1}(\theta)\right]<\infty.$ (2.4)

Furthermore, the limit is then positive almost surely. The first condition is
equivalent to $\theta<\theta_{0}$, thus

$\displaystyle
W_{n}(\theta)\xrightarrow[n\rightarrow\infty]{a.s.}\begin{cases}W^{\infty}_{\theta}&\mbox{if
$\theta<\theta_{0}$ and \eqref{eq:h2} is satisfied,}\\\
0&\mbox{otherwise.}\end{cases}$ (2.5)

where $W_{\theta}^{\infty}$ is finite and positive almost surely. We also
define the derivative martingale associated with $S$ as

$D_{n}=-\sum_{|v|=n}\left(\theta_{0}S_{v}-n\nu(\theta_{0})\right)e^{\theta_{0}S_{v}-n\nu(\theta_{0})}$

As seen in Proposition A.3 from [1], under assumptions

(L1)

$\theta_{0}<\infty$ and

$\displaystyle\mathbb{E}\left[\sum_{i=1}^{N}e^{\theta_{0}\xi_{i}}\xi_{i}^{2}\right]<\infty.$

(L2)

$\theta_{0}<\infty$, and for
$\tilde{X}=\sum_{i=1}^{N}e^{\theta_{0}\xi_{i}}\xi_{i}$,
$X=\sum_{i=1}^{N}e^{\theta_{0}\xi_{i}}$

$\displaystyle\mathbb{E}\left[\tilde{X}\log_{+}\tilde{X}\right]<\infty,$
$\displaystyle\mathbb{E}\left[X\log_{+}^{2}X\right]<\infty,$

where $\log_{+}x=\max\\{0,\log x\\}$ and
$\log^{2}_{+}x=\left(\log_{+}x\right)^{2}$, we have

$\displaystyle D_{n}\xrightarrow[n\rightarrow\infty]{a.s.}D_{\infty}$ (2.6)

where $D_{\infty}$ is finite and positive almost surely.

These two martingales are connected through Theorem 1.1 from [2], which states
that under $(\textbf{L1})$ and $(\textbf{L2})$

$\displaystyle
n^{\frac{1}{2}}W_{n}(\theta_{0})\xrightarrow[n\rightarrow\infty]{\mathbb{P}}c_{\infty}D_{\infty}$
(2.7)

where

$c_{\infty}=\left(\frac{2}{\pi\sigma^{2}}\right)^{\frac{1}{2}}\quad\mbox{and}\quad\sigma^{2}=\mathbb{E}\left[\sum_{i=1}^{N}\left(\theta_{0}\xi_{i}-\nu(\theta_{0})\right)^{2}e^{\theta_{0}\xi_{i}-\nu(\theta_{0})}\right].$

For more results on these martingales and their limits see for example Chapter
3 in [12].

We are now ready to present our results, starting with the almost sure
convergence.

######  Theorem 2.1.

Assume that condition $(\textbf{H})$ is satisfied and
$\mathbb{E}\left[W_{1}(\gamma\theta)\log_{+}W_{1}(\gamma\theta)\right]<\infty$
for $\theta<\frac{\theta_{0}}{\gamma}$. Then

$\displaystyle\frac{R_{n}^{*}}{n}\xrightarrow[n\rightarrow\infty]{a.s.}\frac{\nu\left(\gamma\theta\right)}{\gamma\theta}.$

######  Theorem 2.2.

If $\theta_{0}\leq\theta$ and $\mu$ has finite r-th moments for all
$r<\frac{\theta_{0}}{\theta}$, then

$\displaystyle\frac{R_{n}^{*}}{n}\xrightarrow[n\rightarrow\infty]{a.s.}\frac{\nu(\theta_{0})}{\theta_{0}}.$

In particular, the conditions of Theorem 2.2 hold if $\mu$ satisfies
$(\textbf{H})$ and $\theta\geq\frac{\theta_{0}}{\gamma}$.

The results concerning convergence in distribution we split into three cases.

######  Theorem 2.3 (Below the boundary case).

Assume that condition $(\textbf{H})$ is satisfied and
$\mathbb{E}\left[W_{1}(\gamma\theta)\log_{+}W_{1}(\gamma\theta)\right]<\infty$
for $\theta<\frac{\theta_{0}}{\gamma}$. Then

$R^{*}_{n}-n\frac{\nu\left(\gamma\theta\right)}{\gamma\theta}\xrightarrow[n\rightarrow\infty]{d}\frac{1}{\theta}\left(\log
H_{\theta}-\log E\right),$

where $H_{\theta}$ is finite and positive almost surely, and $E$ is
exponential with intensity 1, independent of $H_{\theta}$. Furthermore,
$H_{\theta}$ has the characteristic function
$\mathbb{E}\left[\tilde{g}\left(t\left(W^{\infty}_{\gamma\theta}\right)^{\gamma}\right)\right]$
where $W^{\infty}_{\gamma\theta}$ is the limit from (2.5) and $\tilde{g}$ is
defined in (2.1).

######  Theorem 2.4 (The boundary case).

Assume $(\textbf{L1})$ and $(\textbf{L2})$. If $\mu$ satisfies $(\textbf{H})$,
then for $\theta=\frac{\theta_{0}}{\gamma}$

$R^{*}_{n}-n\frac{\nu\left(\theta_{0}\right)}{\theta_{0}}+\frac{1}{2\theta_{0}}\log
n\xrightarrow[n\rightarrow\infty]{d}\frac{1}{\theta}\left(\log
H_{\theta_{0}}-\log E\right)$

where $H_{\theta_{0}}$ is finite and positive almost surely and $E$ is
exponential with intensity 1, independent of $H_{\theta_{0}}$. Furthermore,
$H_{\theta_{0}}$ has the characteristic function
$\mathbb{E}\left[\tilde{g}\left(t(c_{\infty}D_{\infty})^{\gamma}\right)\right]$,
where $\tilde{g}$ is defined in (2.1).

######  Theorem 2.5 (Above the boundary case).

Assume $(\textbf{L1})$, $(\textbf{L2})$ and that for all $s\in\mathbb{R}$,
$\mathbb{P}\left(\xi_{1},\xi_{2},\dots,\in s\mathbb{Z}\right)<1$. If
$\theta>\theta_{0}$ and $\mu$ has a finite $r$-th moment for some
$r>\frac{\theta_{0}}{\theta}$ and is not concentrated on a single point, then

$R_{n}^{*}-n\frac{\nu(\theta_{0})}{\theta_{0}}+\frac{3\log
n}{2\theta_{0}}\xrightarrow[n\rightarrow\infty]{d}\frac{1}{\theta}\left(\log
Z-\log E\right)$

where $Z$ is finite and positive almost surely and $E$ is exponential with
intensity 1, independent of $Z$. Furthermore, $Z\stackrel{{\scriptstyle
d}}{{=}}(D_{\infty})^{\frac{\theta}{\theta_{0}}}Z_{\frac{\theta_{0}}{\theta}}$,
where $D_{\infty}$ is the almost sure limit of the derivative martingale
defined in (2.6) and $Z_{\frac{\theta_{0}}{\theta}}$ is strictly
$\frac{\theta_{0}}{\theta}$-stable independent of $D_{\infty}$.

In particular, if $\mu$ has finite mean and $\theta>\theta_{0}$, or if $\mu$
satisfies $(\textbf{H})$ and $\theta>\frac{\theta_{0}}{\gamma}$, then the
assumption in the last theorem is satisfied. It is worth noting that the
logarithmic correction term in Theorem 2.5 is the same as in the case of a
classic BRW without perturbations (see [1]). This result is also more
extensive than Theorem 2.6 in [3], where the asymptotics were only given for
the case $\mu=\delta_{1}$.

##  3 Proofs of Theorems 2.3, 2.4, and 2.5

First we recall Lemma 4.1 from [7] (presented here in a slightly more
accessible form for our use), that will be useful to understand behaviour of
the asymptotics of $Y_{n}$.

######  Lemma 3.1.

Let $\\{Y_{v}\\}_{v\in\mathbb{T}}$ be i.i.d. random variables with
distribution $\mu$ satisfying $(\textbf{H})$, and
$\\{A_{v}\\}_{v\in\mathbb{T}}$ be a sequence of positive random variables,
independent of $\\{Y_{v}\\}_{v\in\mathbb{T}}$, such that

$\sum_{|v|=n}A_{v}^{\gamma}\xrightarrow[n\rightarrow\infty]{\mathbb{P}}A\quad\mbox{and}\quad\sup_{|v|=n}A_{v}\xrightarrow[n\rightarrow\infty]{\mathbb{P}}0$

for some positive random variable $A$. Then

$\sum_{|v|=n}A_{v}Y_{v}\xrightarrow[n\rightarrow\infty]{d}H$

where $H$ has the characteristic function
$\varphi_{H}\left(t\right)=\mathbb{E}\left[\tilde{g}(tA^{\frac{1}{\gamma}})\right]$.

###### Proof of Theorem 2.3.

First we will prove that

$Y_{n}(\theta)e^{-n\frac{\nu\left(\gamma\theta\right)}{\gamma}}\xrightarrow[n\rightarrow\infty]{d}H_{\theta}$
(3.1)

where $H_{\theta}$ has the characteristic function
$\mathbb{E}\tilde{g}\big{(}t(W^{\infty}_{\gamma\theta})^{\frac{1}{\gamma}}\big{)}$
and moreover $H_{\theta}$ is positive almost surely. For this purpose we will
use Lemma 3.1 for $A_{v}=e^{\theta
S_{v}-|v|\frac{\nu\left(\theta\gamma\right)}{\gamma}}$. To check its
hypotheses observe

$\sup_{|v|=n}A_{v}=\sup_{|v|=n}e^{n\theta\left(\frac{S_{v}}{n}-\frac{\nu\left(\theta\gamma\right)}{\theta\gamma}\right)}=e^{n\theta\left(\frac{R_{n}}{n}-\frac{\nu\left(\theta\gamma\right)}{\theta\gamma}\right)}.$

As
$\frac{R_{n}}{n}\xrightarrow[n\rightarrow\infty]{a.s.}\frac{\nu(\theta_{0})}{\theta_{0}}$
and $\theta_{0}$ is the argument minimizing $\frac{\nu\left(t\right)}{t}$,
then
$n\theta\left(\frac{R_{n}}{n}-\frac{\nu\left(\theta\gamma\right)}{\theta\gamma}\right)\xrightarrow[n\rightarrow\infty]{a.s.}-\infty$
and so $\sup_{|v|=n}A_{v}\xrightarrow[n\rightarrow\infty]{a.s.}0$.
Furthermore, in view of (2.5),

$\sum_{|v|=n}A_{v}^{\gamma}=W_{n}\left(\theta\gamma\right)$

converges to $W^{\infty}_{\gamma\theta}$, which is positive almost surely.
Summarizing, Lemma 3.1 entails (3.1)

To see that $H_{\theta}$ is positive almost surely, choose any
$\varepsilon>0$. Then using Exercise 3.3.2 from [9]

$\displaystyle\mathbb{P}\left(H_{\theta}=0\right)$
$\displaystyle=\lim_{T\rightarrow\infty}\frac{1}{2T}\int_{-T}^{T}\mathbb{E}\left[\tilde{g}\big{(}t(W^{\infty}_{\gamma\theta})^{\frac{1}{\gamma}}\big{)}\right]dt$
$\displaystyle\leq\lim_{T\rightarrow\infty}\frac{1}{2T}\int_{-T}^{T}\mathbb{E}\left[e^{-k|t|^{\gamma}W^{\infty}_{\gamma\theta}}\right]dt$
$\displaystyle=\lim_{T\rightarrow\infty}\frac{1}{2T}\int_{-T}^{T}\left(\mathbb{E}\left[e^{-k|t|^{\gamma}|W^{\infty}_{\gamma\theta}|};W^{\infty}_{\gamma\theta}\leq\varepsilon\right]+\mathbb{E}\left[e^{-k|t|^{\gamma}W^{\infty}_{\gamma\theta}};W^{\infty}_{\gamma\theta}>\varepsilon\right]\right)dt$
$\displaystyle\leq\mathbb{P}\left(W^{\infty}_{\gamma\theta}\leq\varepsilon\right)+\lim_{T\rightarrow\infty}\frac{1}{2T}\int_{-T}^{T}e^{-k|t|^{\gamma}\varepsilon}dt$

Now the function $t\rightarrow e^{-k\varepsilon|t|^{\gamma}}$ is integrable,
hence the limit of the second term is 0 for any $\varepsilon$. The first term
can be made arbitrarily small through the choice of $\varepsilon$. Since we
know that $W^{\infty}_{\gamma\theta}$ is positive almost surely, we conclude
positivity of $H_{\theta}$.

Next, recalling (1.1)

$R^{*}_{n}-n\frac{\nu\left(\gamma\theta\right)}{\gamma\theta}\stackrel{{\scriptstyle
d}}{{=}}\frac{1}{\theta}\left(\log Y_{n}(\theta)-\log
E\right)-n\frac{\nu\left(\gamma\theta\right)}{\gamma\theta}=\frac{1}{\theta}\left(\log
Y_{n}(\theta)e^{-n\frac{\nu\left(\gamma\theta\right)}{\gamma}}-\log E\right)$

with $E\sim Exp\left(1\right)$ independent of $Y_{n}$. Finally, by (3.1) and
the continuous mapping theorem

$\log\left(Y_{n}(\theta)e^{-n\frac{\nu\left(\gamma\theta\right)}{\gamma}}\right)\xrightarrow[n\rightarrow\infty]{d}\log
H_{\theta}$

where distribution of $H_{\theta}$ is as specified in the statement. ∎

###### Proof of Theorem 2.4(the boundary case).

Define $A_{v}=e^{\theta
S_{v}-n\frac{\nu\left(\theta_{0}\right)\theta}{\theta_{0}}+\frac{\theta}{2\theta_{0}}\log
n}$ then, by (2.7), we have

$\sum_{|v|=n}A^{\gamma}_{v}=n^{\frac{1}{2}}W_{n}\xrightarrow[n\rightarrow\infty]{\mathbb{P}}\left(\frac{2}{\pi\sigma^{2}}\right)^{\frac{1}{2}}D_{\infty}=Z_{\theta_{0}}$

The second condition of Lemma 3.1,
$\sup_{|v|=n}A_{v}\xrightarrow[n\rightarrow\infty]{\mathbb{P}}0$, follows by
applying Proposition A.3. in [11] to
$V_{u}=\theta_{0}S_{u}-|u|\nu(\theta_{0})$. Then once again by Lemma 3.1

$Y_{n}(\theta)n^{\frac{\theta}{2\theta_{0}}}e^{-n\frac{\theta\nu(\theta_{0})}{\theta_{0}}}=\sum_{|v|=n}A_{v}Y_{v}\xrightarrow[n\rightarrow\infty]{d}H_{\theta_{0}}.$
(3.2)

Next, by (1.1)

$\displaystyle
R^{*}_{n}-n\frac{\nu\left(\theta_{0}\right)}{\theta_{0}}+\frac{1}{2\theta_{0}}\log
n\stackrel{{\scriptstyle d}}{{=}}\frac{1}{\theta}\left(\log Y_{n}(\theta)-\log
E\right)-n\frac{\nu\left(\theta_{0}\right)}{\theta_{0}}+\frac{1}{2\theta_{0}}\log
n$
$\displaystyle=\frac{1}{\theta}\left(\log\big{(}n^{\frac{\theta}{2\theta_{0}}}Y_{n}(\theta)e^{-n\frac{\theta\nu(\theta_{0})}{{\theta_{0}}}}\big{)}-\log
E\right)$

and by (3.2) and the continuous mapping theorem

$\log\big{(}n^{\frac{\theta}{2\theta_{0}}}Y_{n}(\theta)e^{-n\frac{\theta\nu(\theta_{0})}{{\theta_{0}}}}\big{)}\xrightarrow[n\rightarrow\infty]{d}\log
H_{\theta_{0}}$

where distribution of $H_{\theta_{0}}$ is as specified in the statement. ∎

###### Proof of Theorem 2.5 (above the boundary case).

The proof relies on Proposition 3.2 in [8]. The assumptions for Theorem 2.5
with condition $\theta_{0}\nu^{\prime}(\theta_{0})=\nu(\theta_{0})$ for $S$
are equivalent to assumptions (A1) through (A3) from [8] for a BRW

$V_{u}=-\theta\left(S_{u}-|u|\frac{\nu(\theta_{0})}{\theta_{0}}\right)$

with critical parameter $\vartheta=\frac{\theta_{0}}{\theta}$. Proposition 3.2
in [8] entails

$n^{\frac{3}{2}\vartheta}\sum_{|u|=n}e^{-V_{u}}Y_{u}\xrightarrow[n\rightarrow\infty]{d}Z$

where Z is positive almost surely. Furthermore, by equation (1.13) in [8], we
have $Z\stackrel{{\scriptstyle d}}{{=}}D^{\frac{1}{\vartheta}}Z_{\vartheta}$,
where $D$ is the limit of the derivative martingale associated with
$-\vartheta V$, and $Z_{\vartheta}$ is strictly $\vartheta$-stable independent
of $D$. If we let $\psi$ be the log-Laplace transform of $-\vartheta V$, then
it satisfies the equation $\psi(1)=0=\psi^{\prime}(1)$, so the derivative
martingale associated with $-\vartheta V$ is

$\sum_{|u|=n}\vartheta V_{u}e^{-\vartheta
V_{u}}=-\sum_{|v|=n}\left(\theta_{0}S_{v}-n\nu(\theta_{0})\right)e^{\theta_{0}S_{v}-n\nu(\theta_{0})}=D_{n}$

so $D^{\frac{1}{\vartheta}}=(D_{\infty})^{\frac{\theta}{\theta_{0}}}$.
Therefore

$n^{\frac{3\theta}{2\theta_{0}}}Y_{n}(\theta)e^{n\frac{\theta\nu\left(\theta_{0}\right)}{\theta_{0}}}=n^{\frac{3}{2}\eta}\sum_{|u|=n}e^{V_{u}}Y_{u}\xrightarrow[n\rightarrow\infty]{d}Z$

where the distribution of $Z$ is as in the statement. By (1.1) we have

$\displaystyle R_{n}^{*}-n\frac{\nu(\theta_{0})}{\theta_{0}}+\frac{3\log
n}{2\theta_{0}}$ $\displaystyle\stackrel{{\scriptstyle
d}}{{=}}\frac{1}{\theta}\left(\log Y_{n}(\theta)-\log
E\right)-n\frac{\nu(\theta_{0})}{\theta_{0}}+\frac{3\log n}{2\theta_{0}}$
$\displaystyle=\frac{1}{\theta}\left(\log
n^{\frac{3\theta}{2\theta_{0}}}Y_{n}(\theta)e^{-n\frac{\theta\nu\left(\theta_{0}\right)}{\theta_{0}}}-\log
E\right)$

and by the continuous mapping theorem

$\log
n^{\frac{3\theta}{2\theta_{0}}}Y_{n}(\theta)e^{-n\frac{\theta\nu\left(\theta_{0}\right)}{\theta_{0}}}\xrightarrow[n\rightarrow\infty]{d}\log
Z$

which completes the proof. ∎

##  4 Proof of Theorems 2.1 and 2.2

We start with the following Lemma, which gives the convergence in probability.

######  Lemma 4.1.

(a) If conditions of Theorem 2.1 hold, then

$\displaystyle\frac{R_{n}^{*}}{n}\xrightarrow[n\rightarrow\infty]{\mathbb{P}}\frac{\nu\left(\gamma\theta\right)}{\gamma\theta}$

(b) If conditions of Theorem 2.2 hold, then

$\frac{R_{n}^{*}}{n}\xrightarrow[n\rightarrow\infty]{\mathbb{P}}\frac{\nu(\theta_{0})}{\theta_{0}}$

###### Proof.

Let $\beta=\gamma\theta$ in case (a) and $\beta=\theta_{0}$ in case (b). We
will prove first that

$\frac{1}{n\theta}\log\big{(}Y_{n}(\theta)e^{-n\frac{\theta\nu\left(\beta\right)}{\beta}}\big{)}\xrightarrow[n\rightarrow\infty]{\mathbb{P}}0$
(4.1)

We consider first the case (b). Fix an arbitrary $\varepsilon>0$ and choose
$\delta<\frac{\theta_{0}}{\theta}$ satisfying

$\frac{\nu(\theta_{0})}{\theta_{0}}-\frac{\nu\left(\delta\theta\right)}{\delta\theta}+\varepsilon>0.$

The Markov inequality yields

$\displaystyle\mathbb{P}\left(\frac{1}{n\theta}\log\big{(}Y_{n}(\theta)e^{-n\frac{\theta\nu\left(\theta_{0}\right)}{\theta_{0}}}\big{)}>\varepsilon\right)$
$\displaystyle=\mathbb{P}\left(\delta\log\big{(}Y_{n}(\theta)e^{-n\frac{\theta\nu\left(\theta_{0}\right)}{\theta_{0}}}\big{)}>n\theta\delta\varepsilon\right)$
$\displaystyle=\mathbb{P}\left(Y_{n}(\theta)^{\delta}e^{-n\delta\frac{\theta\nu\left(\theta_{0}\right)}{\theta_{0}}}>e^{n\theta\delta\varepsilon}\right)$
$\displaystyle\leq\mathbb{E}[Y_{n}\left(\theta\right)^{\delta}]e^{-\delta
n\theta\left(\frac{\nu(\theta_{0})}{\theta_{0}}+\varepsilon\right)}.$

Since $\delta<1$ and for any $v$ the random variable $Y_{v}$ is independent of
$S_{v}$, we obtain

$\mathbb{E}[Y_{n}(\theta)^{\delta}]=\mathbb{E}\bigg{[}\bigg{(}\sum_{|v|=n}e^{\theta
S_{v}}Y_{v}\bigg{)}^{\delta}\bigg{]}\leq\mathbb{E}\bigg{[}\sum_{|v|=n}e^{\theta\delta
S_{v}}Y_{v}^{\delta}\bigg{]}=e^{n\nu\left(\theta\delta\right)}\mathbb{E}[Y^{\delta}],$
(4.2)

where the last expectation is finite. Summarizing

$\displaystyle\mathbb{P}\left(\frac{1}{n\theta}\log\big{(}Y_{n}(\theta)e^{-n\frac{\theta\nu\left(\theta_{0}\right)}{\theta_{0}}}\big{)}>\varepsilon\right)\leq\mathbb{E}[Y^{\delta}]e^{-\delta
n\theta\left(\frac{\nu(\theta_{0})}{\theta_{0}}-\frac{\nu\left(\delta\theta\right)}{\delta\theta}+\varepsilon\right)}$

and thanks to our choice of $\delta$ the above expression converges to 0 as
$n$ tends to $+\infty$. To prove the remaining bound, denote
$v_{n}=\arg\max_{|v|=n}S_{v}$ . Since

$\displaystyle\frac{1}{n\theta}\log\big{(}Y_{n}(\theta)e^{-n\frac{\theta\nu(\theta_{0})}{\theta_{0}}}\big{)}\geq\frac{1}{n\theta}\log\big{(}e^{n\theta\left(\frac{S_{v_{n}}}{n}-\frac{\nu(\theta_{0})}{\theta_{0}}\right)}Y_{v_{n}}\big{)}=\frac{R_{n}}{n}-\frac{\nu(\theta_{0})}{\theta_{0}}+\frac{1}{n\theta}\log
Y_{v_{n}},$

for any parameters $0<\delta<\varepsilon$ we obtain

$\displaystyle\mathbb{P}\left(\frac{1}{n\theta}\log\big{(}Y_{n}(\theta)e^{-n\frac{\theta\nu\left(\theta_{0}\right)}{\theta_{0}}}\big{)}<-\varepsilon\right)$
$\displaystyle\leq\mathbb{P}\left(\frac{R_{n}}{n}-\frac{\nu(\theta_{0})}{\theta_{0}}+\frac{1}{n\theta}\log
Y_{v_{n}}<-\varepsilon\right)$
$\displaystyle=\mathbb{P}\left(e^{n\theta\left(\frac{R_{n}}{n}-\frac{\nu(\theta_{0})}{\theta_{0}}\right)}Y_{v_{n}}<e^{-\varepsilon
n\theta}\right)$
$\displaystyle=\mathbb{P}\left(e^{n\theta\left(\frac{R_{n}}{n}-\frac{\nu(\theta_{0})}{\theta_{0}}+\varepsilon\right)}Y_{v_{n}}<1\right)$
$\displaystyle\leq\mathbb{P}\left(e^{n\theta\delta}Y_{v_{n}}<1\right)+\mathbb{P}\left(\frac{R_{n}}{n}-\frac{\nu(\theta_{0})}{\theta_{0}}+\varepsilon<\delta\right)$

Now, since $\frac{R_{n}}{n}$ converges almost surely to
$\frac{\nu(\theta_{0})}{\theta_{0}}$ and $\delta<\varepsilon$, the second term
converges to 0. For the first term we have

$\mathbb{P}\left(e^{n\theta\delta}Y_{v_{n}}<1\right)=\mathbb{P}\left(Y<e^{-n\theta\delta}\right)\rightarrow
0.$

Thus, we conclude the proof of (4.1) for case (b).

For case (a), by Theorem 2.3, $\log
Y_{n}(\theta)e^{-n\frac{\nu\left(\gamma\theta\right)}{\gamma}}$ converges in
distribution to $\log H_{\theta}$ and this limit is finite almost surely.
Therefore
$\frac{1}{n\theta}\log\big{(}Y_{n}(\theta)e^{-n\frac{\nu\left(\gamma\theta\right)}{\gamma}})$
converges in distribution to 0, and hence the convergence holds in probability
as well. Thus, the proof of (4.1) is completed.

To prove Lemma 4.1 notice that using (1.1) we can write

$\displaystyle\frac{R^{*}_{n}(\theta)}{n}$
$\displaystyle\stackrel{{\scriptstyle d}}{{=}}\frac{\log
Y_{n}(\theta)}{n\theta}-\frac{\log
E}{n\theta}=\frac{1}{n\theta}\log\big{(}Y_{n}(\theta)e^{-n\frac{\theta\nu\left(\beta\right)}{\beta}}\big{)}+\frac{\nu\left(\beta\right)}{\beta}-\frac{\log
E}{n\theta}$

Now $\frac{\log E}{n\theta}$ converges to 0 almost surely and by (4.1),
$\frac{1}{n\theta}\log
Y_{n}(\theta)e^{-n\frac{\theta\nu\left(\beta\right)}{\beta}}$ converges to 0
in probability. That completes the proof of the Lemma. ∎

###### Proof of Theorems 2.1 and 2.2 (almost sure convergence).

To prove the almost sure convergence we adopt here the arguments given in the
proof of Theorem 2.1 in [3]. For the sake of completeness, we present a
complete proof. Again let $\beta=\gamma\theta$ if conditions of Theorem 2.1
are satisfied and $\beta=\theta_{0}$ if conditions of Theorem 2.2 are
satisfied. We start with the upper bound

$\limsup\limits_{n\rightarrow\infty}{\frac{R^{*}_{n}(\theta)}{n}}\leq\frac{\nu\left(\beta\right)}{\beta}\hskip
14.22636pt\text{ a.s.}$ (4.3)

Fix any $\varepsilon>0$. By (1.1) and the Markov inequality we get that for
any $\delta<\min(\frac{\theta_{0}}{\theta},1)$

$\displaystyle\mathbb{P}\left(\frac{R^{*}_{n}(\theta)}{n}-\frac{\nu\left(\beta\right)}{\beta}>\varepsilon\right)$
$\displaystyle=\mathbb{P}\left({\theta\delta
R^{*}_{n}(\theta)}-\frac{\theta\delta
n\nu\left(\beta\right)}{\beta}>n\delta\theta\varepsilon\right)$
$\displaystyle=\mathbb{P}\left(\log\frac{Y_{n}(\theta)^{\delta}}{E^{\delta}}-\frac{\theta\delta
n\nu\left(\beta\right)}{\beta}>n\delta\theta\varepsilon\right)$
$\displaystyle\leq e^{-\delta
n\theta\left(\frac{\nu\left(\beta\right)}{\beta}+\varepsilon\right)}\mathbb{E}\left[E^{-{\delta}}\right]\mathbb{E}\left[Y_{n}(\theta)^{\delta}\right]$
$\displaystyle\leq e^{-\delta
n\theta\left(\frac{\nu\left(\beta\right)}{\beta}-\frac{\nu\left(\theta\delta\right)}{\theta\delta}+\varepsilon\right)}\Gamma\left(1-\delta\right)\mathbb{E}\left[Y^{\delta}\right],$

where the last inequality follows from (4.2).

Since $\nu$ is continuous, we can choose $\delta$ so that

$\frac{\nu\left(\beta\right)}{\beta}-\frac{\nu\left(\delta\theta\right)}{\delta\theta}+\varepsilon>0,$

therefore the series

$\sum_{n=1}^{\infty}\mathbb{P}\left(\frac{R^{*}_{n}(\theta)}{n}-\frac{\nu\left(\beta\right)}{\beta}>\varepsilon\right)$

converges. The Borel-Cantelli lemma and arbitrariness of $\varepsilon$ entails
(4.3).

Finally our goal is to prove the lower bound

$\liminf_{n\rightarrow\infty}\frac{R_{n}^{*}\left(\theta\right)}{n}\geq\frac{\nu\left(\beta\right)}{\beta}\hskip
14.22636pt\text{ a.s.}$ (4.4)

For $u$ such that $|u|=m\leq n$, we define

$R_{n-m}^{*\left(u\right)}\left(\theta\right):=\max_{v>u,|v|=n}\left(S\left(v\right)+\frac{1}{\theta}\log\left(Y_{v}/E_{v}\right)\right)-S\left(u\right).$

Note that $\\{R_{n-m}^{*\left(u\right)}\left(\theta\right)\\}_{|u|=m}$ are
i.i.d. and have the same distribution as $R_{n-m}^{*}\left(\theta\right)$.
Now,

$\displaystyle R_{n}^{*}\left(\theta\right)$
$\displaystyle=\max_{|u|=m}\max_{v>u,|v|=n}\left(S\left(v\right)+\frac{1}{\theta}\log\left(Y_{v}/E_{v}\right)\right)$
$\displaystyle=\max_{|u|=m}\left(S\left(u\right)+R_{n-m}^{*\left(u\right)}\left(\theta\right)\right)$
$\displaystyle\geq
S\left(\tilde{u}_{m}\right)+\max_{|u|=m}\left(R_{n-m}^{*\left(u\right)}\left(\theta\right)\right),$

where

$\tilde{u}_{m}:=\arg\max_{|u|=m}\left(R_{n-m}^{*\left(u\right)}\left(\theta\right)\right).$

Now, for any $\varepsilon\in\left(0,1\right)$ and small $s$ such that
$\nu(-{s}/2)$ is finite,

$\displaystyle\mathbb{P}\left(\frac{R_{n}^{*}\left(\theta\right)}{n}-\frac{\nu\left(\beta\right)}{\beta}<-\varepsilon\right)$
$\displaystyle\leq\,\,\mathbb{P}\left(S\left(\tilde{u}_{[\sqrt{n}]}\right)+\max_{|u|=[\sqrt{n}]}\left(R_{n-[\sqrt{n}]}^{*\left(u\right)}\left(\theta\right)\right)<n\left(\frac{\nu\left(\beta\right)}{\beta}-\varepsilon\right)\right)$
$\displaystyle\leq\,\,\mathbb{P}\left(\max_{|u|=[\sqrt{n}]}\left(R_{n-[\sqrt{n}]}^{*\left(u\right)}\left(\theta\right)\right)<n\left(\frac{\nu\left(\beta\right)}{\beta}-\frac{\varepsilon}{2}\right)\right)+\mathbb{P}\left(S\left(\tilde{u}_{[\sqrt{n}]}\right)<-\frac{n\varepsilon}{2}\right)$
$\displaystyle\leq\,\,\mathbb{E}\left[\mathbb{P}\left(R_{n-[\sqrt{n}]}^{*}\left(\theta\right)<n\left(\frac{\nu\left(\beta\right)}{\beta}-\frac{\varepsilon}{2}\right)\right)^{N_{[\sqrt{n}]}}\right]+e^{-n\varepsilon
s/4}\cdot\mathbb{E}\left[e^{-sS\left(\tilde{u}_{[\sqrt{n}]}\right)/2}\right],$

where $N_{k}$ is the number of offspring in $k$-th generation. Recalling Lemma
4.1, for all large enough $n$,

$\mathbb{P}\left(R_{n-[\sqrt{n}]}^{*}\left(\theta\right)<n\left(\frac{\nu\left(\beta\right)}{\beta}-\frac{\varepsilon}{2}\right)\right)<\varepsilon.$

Observe that $N_{k}<n$ implies at least $k-\lceil\log_{2}n\rceil$ many
particles up to generation $k$ give birth to only one offspring. Indeed,
denote as $A_{k}$ the number of particles up to $k$-th generation that had
only one child. Then clearly

$N_{k}<n\implies k-\lceil\log_{2}N_{k}\rceil\geq k-\lceil\log_{2}n\rceil$

so it suffices to show $k-\lceil\log_{2}N_{k}\rceil\leq A_{k}$. This on the
other hand is equivalent to $2^{k}\leq
2^{\lceil\log_{2}N_{k}\rceil}2^{A_{k}}$, which is obviously true if $A_{k}\geq
k$. If $A_{k}<k$, then it means that there is at least $k-A_{k}$ generations
in which all the particles have at least 2 children, so $N_{k}\geq
2^{k-A_{k}}$, which again gives the desired result. Therefore

$\mathbb{P}\left(N_{[\sqrt{n}]}<n\right)\leq\left(\mathbb{P}\left(N=1\right)\right)^{[\sqrt{n}]-\lceil\log_{2}n\rceil}.$

so

$\mathbb{E}\left[\mathbb{P}\left(R_{n-[\sqrt{n}]}^{*}\left(\theta\right)<n\left(\frac{\nu\left(\beta\right)}{\beta}-\frac{\varepsilon}{2}\right)\right)^{N_{[\sqrt{n}]}}\right]\leq\mathbb{E}[\varepsilon^{N_{[\sqrt{n}]}}]\\\
\leq\mathbb{E}\left[\varepsilon^{n};N_{[\sqrt{n}]}\geq
n\right]+\mathbb{P}\left(N_{[\sqrt{n}]}<n\right)\leq\varepsilon^{n}+\left(\mathbb{P}\left(N=1\right)\right)^{[\sqrt{n}]-\lceil\log_{2}n\rceil}\leq\varepsilon_{1}^{\sqrt{n}}$

for some $\varepsilon_{1}<1$. To estimate the second term, we bound supremum
by the sum and we have

$\mathbb{E}\left[e^{-sS\left(\tilde{u}_{[\sqrt{n}]}\right)/2}\right]\leq\mathbb{E}\bigg{[}\sum_{|v|=\lceil\sqrt{n}\rceil}e^{-\frac{s}{2}S_{v}}\bigg{]}=e^{[\sqrt{n}]\nu\left(-s/2\right)}.$

Therefore we have for all large enough $n$,

$\mathbb{P}\left(\frac{R_{n}^{*}\left(\theta\right)}{n}-\frac{\nu(\beta)}{\beta}<-\varepsilon\right)\leq\,\varepsilon_{1}^{\sqrt{n}}+e^{-n\varepsilon
s/4+[\sqrt{n}]\nu\left(-s/2\right)}.$

Since for every $\varepsilon\in\left(0,1\right)$,

$\sum_{n=1}^{\infty}\mathbb{P}\left(\frac{R_{n}^{*}\left(\theta\right)}{n}-\frac{\nu(\beta)}{\beta}<-\varepsilon\right)<\infty,$

using Borel-Cantelli Lemma once again we deduce (4.4), completing the proof. ∎

## References

  * [1] E. Aïdékon (2013). _Convergence in law of the minimum of a branching random walk._ Annals of Probability 2013, Vol. 41, No. 3A, 1362-1426 MR3098680
  * [2] E. Aïdékon, Z. Shi (2014). _The Seneta-Heyde scaling for the branching random walk._ Annals of Probability 2014, Vol. 42, No. 3, 959–993 MR3189063
  * [3] A. Bandyopadhyay, P. P. Ghosh (2023). _Right-Most Position of a Last Progeny Modified Branching Random Walk_ , arXiv:2106.02880v2
  * [4] A. Bandyopadhyay, P. P. Ghosh (2023). _Right-most position of a last progeny modified time inhomogeneous branching random walk_ , Statistics & Probability Letters Vol. 193, February 2023, 109697 MR4499430
  * [5] J. D. Biggins (1976). _The First- and Last-Birth Problems for a Multitype Age-Dependent Branching Process_ Advances in Applied Probability Vol. 8, No. 3 (Sep., 1976), pp. 446-459 MR0420890
  * [6] J. D. Biggins(1977). _Martingale convergence in the branching random walk._ J. Appl. Probab. 14, 25–37. MR0433619
  * [7] K. Bogus, D. Buraczewski and A. Marynych (2020). _Self-similar solutions of kinetic-type equations: The boundary case_ , Stochastic Processes and their Applications 130 (2020) 677–693 MR4046515
  * [8] D. Buraczewski, K. Kolesko, M. Meiners (2021). _Self-similar solutions to kinetic-type evolution equations: beyond the boundary case._ Electron. J. Probab. 26 1 - 18, 2021. MR4216515
  * [9] R. Durrett (2019). _Probability: Theory and Examples_ Cambridge University Press, 2019 MR3930614
  * [10] P. P. Ghosh (2022) _Large deviations for the right-most position of a last progeny modified branching random walk._ Electron. Commun. Probab. 27 1 - 13, 2022. MR4375913
  * [11] A. Iksanov, K. Kolesko, M. Meiners (2018). _Fluctuations of Biggins’ martingales at complex parameters_ , Annales de l Institut Henri Poincaré Probabilités et Statistiques 56(4):2445-2479 MR4164844
  * [12] Z. Shi, _Branching random walks_ (2015). École d’Été de Probabilités de Saint-Flour XLII – 2012, Springer MR3444654


Void-Finding Systematics using Crossing Numbers

Dahlia Veyrat

0000-0001-8101-2836]Dahlia Veyrat
Department of Physics & Astronomy, University of Rochester, 500 Wilson Blvd., Rochester, NY  14627

0000-0002-9540-546X]Kelly A. Douglass
Department of Physics & Astronomy, University of Rochester, 500 Wilson Blvd., Rochester, NY  14627

0000-0001-5537-4710]Segev BenZvi
Department of Physics & Astronomy, University of Rochester, 500 Wilson Blvd., Rochester, NY  14627

We study how well void-finding algorithms identify cosmic void regions and whether we can quantitatively and qualitatively describe their biases by comparing the voids they find with dynamical information from the underlying matter distribution. Using the ORIGAMI algorithm to determine the number of dimensions along which dark matter particles have undergone shell-crossing (crossing number) in $N$-body simulations from the AbacusSummit simulation suite, we identify dark matter particles which have undergone no shell crossing as belonging to voids. We then find voids in the corresponding halo distribution using two different void-finding algorithms: VoidFinder and , a ZOBOV-based algorithm. The resulting void catalogs are compared to the distribution of dark matter particles to examine how their crossing numbers depend on void proximity. While both algorithms' voids have a similar distribution of crossing numbers near their centers, we find that beyond 0.25 times the effective void radius, voids found by VoidFinder exhibit a stronger preference for particles with low crossing numbers than those found by . We examine two possible methods of mitigating this difference in efficacy between the algorithms. While we are able to partially mitigate the ineffectiveness of by using distance from the void edge as a measure of centrality, we conclude that VoidFinder more reliably identifies dynamically-distinct regions of low crossing number.

§ INTRODUCTION

On cosmic scales, the structure of matter in the observable universe exhibits a complex, web-like distribution [Bond et al., 1996], with distinct filaments stretching between dense clusters of galaxies. In the space between these superstructures, galaxy redshift surveys have found vast, relatively empty regions containing very few galaxies [Jõeveer et al., 1978, Gregory & Thompson, 1978, Kirshner et al., 1981]. These cosmic voids, which occupy most of the volume of the universe [de Lapparent et al., 1986, Geller & Huchra, 1989], are the result of the gravitational instability in primordial underdense regions, similar to the formation of clusters from gravitational collapse of primordial overdensities [Zeldovich, 1970, van de Weygaert & Platen, 2011].

The emptiness of voids provides a unique environment of gravitational evolution for cosmological and astrophysical studies. Cosmologically, the lack of large-scale gravitational collapse causes dynamics within voids to remain in the linear regime for a relatively longer time [Goldberg & Vogeley, 2004]. Their uniqueness as cosmic structures has made voids well-suited to studies of the Alcock-Paczyński effect <cit.>, dark energy <cit.>, baryon acoustic oscillations <cit.>, and weak lensing <cit.>. Additionally, the uniqueness of voids as an intergalactic environment results in measurably different evolution and properties of galaxies within them <cit.>.

A number of different types of algorithms have been used to detect voids in the cosmic web. The VoidFinder algorithm [El-Ad & Piran, 1997, Hoyle & Vogeley, 2002] finds relatively empty spheres in the distribution of galaxies and combines them into individual voids. Another common strategy is to link low-density regions together using a watershed algorithm. Several implementations using a Voronoi tessellation to approximate local density exist <cit.>.

While the aforementioned algorithms find voids geometrically, voids are expected to be dynamically unique structures [Sheth & van de Weygaert, 2004]. With simulations, it is possible to define voids relative to the dynamics of the matter within them. Cosmological simulations have been used to predict properties of voids [Ricciardelli et al., 2014, Hamaus et al., 2014] and forecast the cosmological constraining power of voids [Pisani et al., 2015].

One method for identifying dynamical voids in simulations is the computation of the number of dimensions along which matter has gravitationally collapsed [Falck et al., 2012]. Voids are defined as regions undergoing no shell-crossing.

In this work, we investigate the accuracy of void-finding algorithms in detecting dynamical void regions in the matter distribution from dark matter tracer information. We use an $N$-body simulation from the AbacusSummit simulation suite [Maksimova et al., 2021], on which the ORIGAMI algorithm [Falck et al., 2012] is run to determine their crossing numbers — the number of dimensions along which they have undergone shell-crossing (that is to say, the number of dimensions of gravitational collapse). Then, using the crossing numbers to classify particles as belonging to voids, walls, filaments, or clusters, we compare this classification with the voids identified by two different void-finding algorithms in the Void Analysis Software Toolkit <cit.>: VoidFinder and . This allows us to quantify the relative accuracy of these algorithms in detecting dynamically distinct regions dominated by low crossing numbers.

The paper is organized as follows. In Section <ref>, we describe the theory of void evolution relevant to our study, and in Section <ref>, we describe the AbacusSummit simulation suite and the properties of the simulation used in the crossing number analysis. Sections <ref> and <ref> present the algorithms used to compute crossing numbers of dark matter particles and to find voids in the distribution of halos, respectively. Results are discussed in Section <ref>, examining the relationships between void regions defined by void-finding algorithms and the void particles identified by ORIGAMI. In Section <ref>, we explore ways to mitigate the relatively poor classification of the watershed algorithm, .

§ EXCURSION FORMALISM

The excursion-set formalism, first proposed by [Press & Schechter, 1974], provides an analytical model of the gravitational collapse and virialization of dark matter halos. [Bond et al., 1991] developed the excursion model to describe the halo mass function, but the model was not applied to a theoretical description of voids until the pioneering work of [Sheth & van de Weygaert, 2004].

The excursion-set model of voids describes the evolution of an initial underdensity in the matter distribution. [Sheth & van de Weygaert, 2004] find that, as matter is attracted to the overdense surroundings, the typical shape of underdensities becomes more spherical and can be effectively described by a series of spherical shells. More central shells experience a stronger “repulsion” from void centers, leading to a distinct “wall” feature — a build up of expanding matter near the void edge as shells approach and cross each other. We expect, then, to find a high but narrow preference for high dark matter particle crossing numbers (the number of dimensions along which they have undergone shell-crossing) near void edges, along with the preference for low crossing numbers in voids' more central regions.

§ SIMULATIONS

The AbacusSummit simulation suite [Maksimova et al., 2021] is a set of over 100 periodic $N$-body simulations spanning several box sizes up to $(7.5\text{ \hGpc})^3$, particle mass resolutions down to $\sim3\times10^8 M_\odot$, and implementing several cosmologies. AbacusSummit was produced with the Abacus $N$-body algorithm [Garrison et al., 2021] at the Oak Ridge Leadership Computing Facility's Summit supercomputer and is the largest high-accuracy $N$-body data set produced to date.

We use an AbacusSummit Hugebase simulation, which evolves $2304^3$ dark matter particles of mass $\sim5\times10^{10}M_\odot$ within a periodic box of width 2 . This simulation uses Planck 2018 $\Lambda$CDM cosmology: $\Omega_m=0.315$, $h=0.674$, $\sigma_8=0.811$, $n_s=0.965$ [Planck Collaboration et al., 2020]. In addition to full particle timeslices, the AbacusSummit data products include halo catalogs produced using the CompaSO halo-finder [Hadzhiyska et al., 2021]. We run all void-finding algorithms on the corresponding halo catalog of the Hugebase simulation, which contains $\sim1.8\times10^7$ halos.

§ CROSSING NUMBERS USING ORIGAMI

The gravitational collapse of the matter distribution in these $N$-body simulations results in the crossing of dark matter particle positions as their dynamics become nonlinear. The number of dimensions of shell-crossing, and thus of nonlinearity, determines the classification of dark matter particles into various types of large-scale structures: clusters from three-dimensional collapse, filaments from two-dimensional collapse, and walls from one-dimensional collapse. Particles which have undergone no shell-crossing are defined as belonging to cosmic voids. The ORIGAMI algorithm [Falck et al., 2012] computes the number of shell-crossing dimensions, referred to as the crossing number, or CN, by comparing the final relative positions of pairs of particles with their relative positions on the initial grid.

Running the base ORIGAMI algorithm on simulations with more than $512^3$ particles requires a prohibitively large amount of memory. We introduce modifications that compute crossing numbers for subvolumes of periodic simulations by implementing a buffer zone and subsequently merge the subvolume results, allowing ORIGAMI to be run on larger simulations. Crossing numbers in the AbacusSummit simulations were computed using this modified version of ORIGAMI[The modified ORIGAMI is available for download at <https://github.com/dveyrat/origami/tree/subdivide>.].

§ VAST

A 100 $\times$100 subsection of a 10 thick slice of the AbacusSummit Hugebase simulation. The locations of particles with crossing number 0 are shown in blue, 1 in orange, 2 in red, and 3 in violet. The intersections of voids with the center of the slice are shown in blue for VoidFinder (top), /VIDE (middle), and /REVOLVER (bottom).

We identify voids in the distribution of halos using the Void Analysis Software Toolkit[VAST is available for download at <https://github.com/desi-ur/vast>.] <cit.>. VAST includes an implementation of the VoidFinder algorithm, which grows and merges nearly empty spheres, and , an implementation of the ZOBOV algorithm [Neyrinck, 2008], which uses a watershed method to identify voids.  The voids found by each algorithm are shown in Figure <ref>, overlaying the dark matter particles within a 10 -thick slice of the AbacusSummit Hugebase simulation.

§.§ VoidFinder

VoidFinder, which was originally described by [El-Ad & Piran, 1997], begins by removing isolated tracers from the catalog, defined as those with a distance to their third-nearest neighbor $d_\text{3NN}>\overline{d_\text{3NN}}+1.5\sigma_{d_\text{3NN}}$. The remaining tracers are placed on a three-dimensional grid, and a sphere is grown from each empty cell until its surface is bounded by four tracers. Next, maximal spheres, defined as spheres with radii $r\geq10$that do not overlap a larger maximal sphere by more than 10% of their volume, are selected from the population of spheres. Each maximal sphere is identified as belonging to a void, and that void is built from the union of all spheres which intersect exactly one maximal sphere by at least 50% of their volume. See [Hoyle & Vogeley, 2002] for a more detailed description of VoidFinder.

§.§ Voronoi Voids ()

Voronoi Voids, or , is a Python implementation of the ZOBOV algorithm, first described in [Neyrinck, 2008]. begins by creating a Voronoi tesselation of the tracer distribution, whose cell volumes are used as an estimate of local density. A watershed algorithm is then used to combine groups of adjacent Voronoi cells into “zones,” where each cell is put into the same zone as its least dense neighboring cell, and any cell which is less dense than any of its neighbors is identified as a local density minimum that serves as its zone's central cell.

Adjacent zones are then linked together by the least-dense pair of adjacent cells between them, whose density is referred to as the “linking density.” The collections of linked zones (including single zones) are identified as voids, creating a hierarchy of voids up to one containing the entire catalog. Several methods exist to prune the full hierarchy and extract meaningful voids from it <cit.>.  We use two of the more popular pruning methods: VIDE pruning <cit.>, which sets both a maximum linking density for the zone-linking step and a maximum central density for each void, and REVOLVER pruning (/REVOLVER, [Nadathur et al., 2019]), which labels each zone as a unique void and removes any voids with effective radius less than the median void effective radius.

§ RESULTS

Crossing number distributions
particles in    particles not in     particles in     particles not in     particles in     particles not in
crossing     all DM     VoidFinder     VoidFinder     /VIDE     /VIDE     /REVOLVER     /REVOLVER
number     particles     voids     voids     voids     voids     voids     voids
0     30.8%     46.9%     16.6%     30.6%     31.3%     31.3%     28.7%

1     19.4%     22.4%     16.6%     19.4%     19.4%     19.4%     19.1%

2     18.2%     15.6%     20.5%     18.2%     18.1%     18.1%     18.5%

3     31.6%     15.0%     46.3%     31.8%     31.2%     31.2%     33.6%

Distribution of crossing numbers within and outside of voids found by VoidFinder and , as well as the total distribution in the AbacusSummit Hugebase simulation. All uncertainties are below 0.1%, estimated using several AbacusSummit Hugebase realizations.

410,852 voids were found in the Abacus Hugebase simulation by VoidFinder, 79,721 by /VIDE, and 43,430 by /REVOLVER. The fraction of dark matter particles with each crossing number in the AbacusSummit Hugebase simulation, and of those in voids defined by each of the void-finding algorithms studied, can be found in Table <ref>. Throughout the entire simulation, we find that $\sim$31% of the dark matter particles have undergone no shell-crossing, $\sim$19% have undergone shell-crossing along one dimension, $\sim$18% along two dimensions, and $\sim$32% along three dimensions. These results are in good agreement with the crossing number distribution trends in [Falck & Neyrinck, 2015] for similar particle resolutions. Our results, using a simulation with an initial grid spacing $L/N=(1000$ )$/1152$, are similar to those in [Falck & Neyrinck, 2015] for initial grid spacing $L/N=(100$ $)/128$ but with a shift towards more void (CN=0) particles and fewer cluster (CN=3) particles, agreeing with the trend for greater initial grid spacing.

§.§ Crossing Number Distribution within Voids

Normalized void density profiles by crossing number for VoidFinder (top), /VIDE (middle), and /REVOLVER (bottom) voids. All uncertainties are negligibly small, estimated using several AbacusSummit Hugebase realizations.

The distributions of crossing numbers by normalized distance from void centers found in real space by VoidFinder (top), /VIDE (middle), and /REVOLVER (bottom).

We expect the distribution of crossing numbers in void regions to contain an excess of CN=0 particles, because any dark matter particles which are identified by ORIGAMI to have undergone shell-crossing are expected to be part of a gravitationally bound structure: either wall, filament, or cluster. To determine the effectiveness of VoidFinder and at detecting regions with low crossing number, we examine the distribution of crossing numbers inside and outside voids. We observe that the distributions of dark matter particles with CN=0 (void particles) and CN=3 (cluster particles) differ significantly when the particles are categorized by whether they are inside or outside a void (Table <ref>). The distributions of CN=1 and CN=2 particles differ relatively little between void/non-void regions. While particles in voids are more likely to have CN=0 and less likely to have CN=3 than the background distribution, this preference is significantly stronger in VoidFinder voids, where $\sim$47% of particles have CN=0 (compared to $\sim$31% in /VIDE and /REVOLVER voids, and $\sim$31% overall) and only $\sim$15% have CN=3 (compared to $\sim$32% in /VIDE voids, $\sim$31% in /REVOLVER voids, and $\sim$32% overall).

Figure <ref> shows the crossing number density profiles of dark matter particles around voids, and Figure <ref> shows how the distribution of crossing numbers is related to the normalized distance to void centers, $r/R_\text{eff}$, where the effective radius $R_\text{eff}$ is defined as the radius of a sphere with the same volume as the void. We find that the particles closest to void centers are significantly more likely to have crossing number 0 than 3 for both VoidFinder and . However, while this preference for low crossing numbers extends to nearly $r=R_\text{eff}$ for VoidFinder voids, the distribution around voids begins to gradually shift toward the background crossing number distribution at $r\approx0.25R_\text{eff}$ for /VIDE voids, and $r\approx0R_\text{eff}$ for /REVOLVER voids. At a normalized distance of $r\approx0.8R_\text{eff}$, the distribution of crossing numbers in voids is nearly identical to the background.

The expected wall feature of voids is visually evident in both Figures <ref> and <ref>, but it is more prominent around voids found by VoidFinder. There is a clear, sharp increase in preference for higher crossing numbers near $r=R_\text{eff}$, in good agreement with the theory developed in [Sheth & van de Weygaert, 2004], which predicts a build-up and crossing of mass shells as they expand radially outward from void centers (see Section <ref> for details). The wall feature around voids is both less pronounced and less localized around $r=R_\text{eff}$.

We also observe a sharp feature at $r=0.25R_\text{eff}$ in the distributions of crossing numbers in /VIDE voids. This is a result of the central density cut present in the algorithm, which removes all voids with density above some threshold within $r<0.25R_\text{eff}$. As seen in the bottom panel of Figures <ref> and <ref>, The distribution of crossing numbers around voids found by /REVOLVER, which does not include a central density cut, lacks this feature.

§.§ CN=0 Particles not in Voids

While the majority of particles with crossing number 0 are located within VoidFinder voids, $\sim$28% do not fall within a void. Further, while particles with lower crossing numbers are more likely to be found within voids than those with higher crossing numbers, $\sim$36% of CN=0 particles do not fall within a /VIDE void and $\sim$15% do not fall within a /REVOLVER void. To determine whether these particles belong to regions we expect to be classified as void, we investigate their local environments.

For each dark matter particle, we investigate its local environment by identifying all particles within 2of it, and counting these particles by crossing number. In a void-like environment, we expect fewer neighboring particles with higher crossing numbers, as well as a lower density of particles overall, leading to a lower particle count for CN0 particles relative to CN=0 and a lower total particle count.

Local environment of CN=0 particles
2cVoidFinder     2c/VIDE
crossing number     all     inside voids     outside voids     inside voids     outside voids     inside voids     outside voids
0     19 (11,28)     18 (11,27)     21 (13,30)     19 (11,28)     19 (11,28)     19 (11,28)     19 (11,28)

1-3     11 (0,61)     7 (0,41)     27 (3,121)     11 (0,61)     10 (0,59)     10 (0,60)     12 (0,65)

total     32 (13,86)     28 (12,67)     51 (20,146)     33 (14,87)     32 (13,85)     32 (13,86)     34 (14,90)

Median number of particles within 2 of particles with CN=0 that fall either inside or outside a void. Values in parentheses indicate the boundaries of the central 68% of the distribution.

The crossing numbers within 2 of CN=0 particles inside and outside of voids are shown in Table <ref> for both VoidFinder and . Overall, CN=0 particles are found to be located in void-like environments with little clustering (CN0 particles). CN=0 particles outside of VoidFinder voids, however, tend to be located in environments with more neighboring particles overall and several times more CN0 particles than in the local environments of CN=0 particles in VoidFinder voids. This suggests that, while dynamical information alone classifies all CN=0 particles as belonging to void regions, many are located in significantly denser environments and subsequently are not expected to belong to cosmic voids. These are therefore expected to not be located in VoidFinder voids.

The local environments of CN=0 particles in voids, however, are much more similar to those of their undetected counterparts and of CN=0 particles as a whole. Like the overall crossing number distributions in Table <ref>, this indicates worse classification of void regions by , with no evident environmental distinction between the
CN=0 particles inside voids and those outside.

§.§ Void Finding in Redshift Space

As an additional check, both void-finding algorithms were run on a version of the halo catalog modified using peculiar velocity information to simulate the effect of peculiar velocities on apparent positions, an effect known as redshift-space distortions (RSDs). Because peculiar velocities within voids tend to be directed outwards, voids appear elongated in the line-of-sight direction in redshift space, and we expect a smoothing of features in the radial distribution of crossing numbers around voids.

Figure <ref> shows how the distribution of crossing numbers is related to normalized distance to void centers found in redshift space. While the distributions are similar to those shown in Figure <ref> at low and high $r/R_\textnormal{eff}$, the distinct features at $r/R_\textnormal{eff}=1$ for VoidFinder and $r/R_\textnormal{eff}=0.25$ for /VIDE appear to be smoothed out, as expected due to the distortion of the voids along the line of sight. This suggests that performing reconstruction of real-space positions in a redshift catalog before running a void-finding algorithm may improve the voids found.

The distributions of crossing numbers by normalized distance from void centers found in redshift space by VoidFinder (top), /VIDE (middle), and /REVOLVER (bottom).

§.§ Summary of Results

We expect void regions to contain an excess of particles with crossing number 0 relative to the background distribution, and both voids found by VoidFinder and voids found by contain such an excess. In voids, however, there is a strong excess only in the most central ($r<0.25R_\textnormal{eff}$) regions, while in VoidFinder voids it extends nearly to $R_\textnormal{eff}$. Further, there is an excess of CN=3 (cluster) particles near the edge of VoidFinder voids, in agreement with dynamical theories of void evolution. We conclude that VoidFinder identifies void regions more accurately than .

§ MITIGATING POOR CLASSIFICATION OF VOIDS

§.§ Linking Density

While the algorithm is mostly parameter-free, there is an input parameter in /VIDE: the maximum zone-linking density, which limits the density of Voronoi cells which can link adjacent zones into voids. The maximum linking density (see Section <ref> for details) is defined relative to the mean number density of the catalog, $\overline{n}$, and is equal to $0.2\overline{n}$ by default, but can be as low as $0$ (allowing no linking of zones, similar to /REVOLVER) or be undefined (leading to all zones being linked together into a single void).

While a lower maximum linking density is expected to limit the merging of voids across denser regions, the growth of voids up to density maxima occurs during the zone creation step, rather than in the subsequent zone-merging step. Consequently, varying the amount of zone merging that occurs is not expected to mitigate the inclusion of the dense shell-crossing region in the void volume of /VIDE voids.

We examine the effect of the maximum linking density on the distribution of crossing numbers within /VIDE voids. These results are shown in Table <ref>. The percent of particles in voids with a given crossing number does not vary by more than $\sim$1% across different linking densities. This is a minor shift compared to the difference between VoidFinder voids and /VIDE voids ($\sim$12% for both CN=0 and CN=3; see Table <ref>), and still leaves the /VIDE distribution relatively similar to the distribution of all crossing numbers.

Crossing number distributions
4cmaximum linking density
crossing number     $0.1\overline{n}$     $0.2\overline{n}$     $0.5\overline{n}$     $1.0\overline{n}$
0     32.2%     33.0%     33.2%     33.2%

1     20.9%     20.9%     20.8%     20.8%

2     18.5%     18.4%     18.2%     18.2%

3     28.4%     27.7%     27.8%     27.8%

Effect of maximum linking density on the distribution of crossing numbers within /VIDE voids. Void-finding was done using a (128 )$^3$ subvolume of the AbacusSummit Hugebase simulation.

§.§ Depth-in-Void

Distribution of crossing numbers by normalized depth from void edges within VoidFinder (top), /VIDE (middle), and /REVOLVER (bottom) voids. Depths were computed only for a (1 )$^3$ subvolume of the AbacusSummit Hugebase simulation.

While the centers of /VIDE voids have similar crossing number distributions to those of VoidFinder voids, this is only out to roughly $0.25R_\text{eff}$ (see the middle panel of Figure <ref>). Restricting the analysis to central regions using distance from the void center is therefore not an effective modification of , as a significant amount of the void volume is omitted, restricting analysis to a small subset of void regions.

Because voids generally have irregular shapes, distance from the edge (depth-in-void) may serve as a better measure of centrality than distance from the center, which is more effective for the more spherical voids found by VoidFinder. [Zaidouni et al., 2023] show that void galaxies with a normalized distance from the /VIDE void edge of at least $0.4d_\text{max}$ have a distribution of astrophysical properties similar to those of galaxies within VoidFinder voids.

We examine the distribution of crossing numbers at different depths within voids to determine whether depth is a better measure of centrality than distance from the center. The results are shown in Figure <ref> for both pruning methods as well as VoidFinder. The previously observed trend of lower crossing numbers for more centrally located void particles is apparent, but just as in the distribution based on distance from the void center (Figure <ref>), the voids do not exhibit the sharp wall feature observed in VoidFinder voids. Further, the central regions of voids defined using depths ($1-d/d_\text{max}<0.25$) do not favor low crossing numbers as strongly as central regions defined using distance from void centers ($r/R_\text{eff}<0.25$). These results indicate there is no clear method for extracting the dynamically-distinct parts of the void regions from voids. The shift back to lower crossing numbers at the edge of voids can be attributed to the fact that these void edges are made up of boundaries between halos' Voronoi cells, an inherently low-density region.

§ CONCLUSIONS

We study the relative accuracy of two different void-finding algorithms in detecting dynamical void regions by examining the positions of their voids relative to non-clustering regions in the dark matter distribution. Using a (2 )$^3$ $N$-body simulation from the AbacusSummit simulation suite [Maksimova et al., 2021], the evolutionary history of each dark matter particle was identified using the ORIGAMI algorithm [Falck et al., 2012], which counts the number of dimensions along which dark matter particles have undergone shell-crossing. Identifying particles with crossing number 0 as belonging to voids, we compare their positions with voids found in the corresponding halo distribution by VoidFinder and , two void-finding algorithms implemented in VAST [Douglass et al., 2022].

We find that while both void-finding algorithms produce voids with a central bias towards lower crossing numbers, this bias gradually diminishes at greater distances from void centers, while the preference remains strong up to the edge of VoidFinder voids. This suggests that includes the shell-crossing region surrounding voids as part of the void volume. Further, at the void edge defined as $r/R_\text{eff}\approx1$, voids found by VoidFinder have a preference for higher crossing numbers that is even stronger than the background distribution at larger radii, in good agreement with the shell-crossing predictions of the excursion-set formalism by [Sheth & van de Weygaert, 2004]. This feature is absent in the distribution of crossing numbers around voids found by .

Given the relative inability of to identify dynamically-distinct voids, we attempt several methods to improve the classification of void regions by . The /VIDE algorithm has one user-set parameter, the maximum linking density. We find that varying the maximum linking density has little effect on the distribution of crossing numbers within voids. We also examine the crossing number distribution as a function of distance from the void edge (depth) rather than distance from the void center. While void depth has been used to provide a definition of a central region with a stronger preference for void-like galaxy properties [Zaidouni et al., 2023], this did not improve the crossing number profiles of /VIDE voids. We conclude that VoidFinder more effectively identifies the dynamically-distinct regions primarily occupied by dark matter particles with low crossing number.

The AbacusSummit Hugebase simulation analyzed in this work is one of the largest to be used in a crossing number analysis. We modified the ORIGAMI algorithm (available for downloads at <https://github.com/dveyrat/origami/tree/subdivide>) to allow it to be run on this $2304^3$ particle set. These modifications enable similar analyses using other AbacusSummit simulations, such as studies of the effects of different void-finding algorithms on cosmological constraints.

§ ACKNOWLEDGEMENTS

The authors would like to thank Stephen W. O'Neill, Jr. for his help improving the VoidFinder algorithm in VAST, and Michael Vogeley and Fiona Hoyle for the original VoidFinder code and void-finding expertise. D.V. and S.B. acknowledge support from the U.S. Department of Energy Office of High Energy Physics under the grant DE-SC0008475.  K.D. and D.V. acknowledge support from grant 62177 from the John Templeton Foundation.

This research used resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility located at Lawrence Berkeley National Laboratory.

[Bond et al., 1991]
Bond, J. R., Cole, S., Efstathiou, G., & Kaiser, N. 1991, , 379,
440, 10.1086/170520

[Bond et al., 1996]
Bond, J. R., Kofman, L., & Pogosyan, D. 1996, , 380, 603,

[Chantavat et al., 2017]
Chantavat, T., Sawangwit, U., & Wandelt, B. D. 2017, , 836, 156,

[de Lapparent et al., 1986]
de Lapparent, V., Geller, M. J., & Huchra, J. P. 1986, , 302, L1,

[Douglass et al., 2019]
Douglass, K. A., Smith, J. A., & Demina, R. 2019, , 886, 153,

[Douglass et al., 2022]
Douglass, K. A., Veyrat, D., O'Neill, S. W., et al. 2022, Journal of Open
Source Software, 7, 4033, 10.21105/joss.04033

[El-Ad & Piran, 1997]
El-Ad, H., & Piran, T. 1997, , 491, 421, 10.1086/304973

[Falck & Neyrinck, 2015]
Falck, B., & Neyrinck, M. C. 2015, , 450, 3239,

[Falck et al., 2012]
Falck, B. L., Neyrinck, M. C., & Szalay, A. S. 2012, , 754, 126,

[Garrison et al., 2021]
Garrison, L. H., Eisenstein, D. J., Ferrer, D., Maksimova, N. A., & Pinto,
P. A. 2021, Monthly Notices of the Royal Astronomical Society, 508, 575,

[Geller & Huchra, 1989]
Geller, M. J., & Huchra, J. P. 1989, Science, 246, 897,

[Goldberg & Vogeley, 2004]
Goldberg, D. M., & Vogeley, M. S. 2004, , 605, 1,

[Gregory & Thompson, 1978]
Gregory, S. A., & Thompson, L. A. 1978, , 222, 784,

[Habouzit et al., 2020]
Habouzit, M., Pisani, A., Goulding, A., et al. 2020, Monthly Notices of the
Royal Astronomical Society, 493, 899, 10.1093/mnras/staa219

[Hadzhiyska et al., 2021]
Hadzhiyska, B., Eisenstein, D., Bose, S., Garrison, L. H., & Maksimova, N.
2021, Monthly Notices of the Royal Astronomical Society, 509, 501,

[Hamaus et al., 2016]
Hamaus, N., Pisani, A., Sutter, P. M., et al. 2016, , 117, 091302,

[Hamaus et al., 2014]
Hamaus, N., Wandelt, B. D., Sutter, P. M., Lavaux, G., & Warren,
M. S. 2014, , 112, 041304, 10.1103/PhysRevLett.112.041304

[Hoyle et al., 2005]
Hoyle, F., Rojas, R. R., Vogeley, M. S., & Brinkmann, J. 2005, The
Astrophysical Journal, 620, 618, 10.1086/427176

[Hoyle & Vogeley, 2002]
Hoyle, F., & Vogeley, M. S. 2002, , 566, 641, 10.1086/338340

[Jõeveer et al., 1978]
Jõeveer, M., Einasto, J., & Tago, E. 1978, , 185, 357,

[Kirshner et al., 1981]
Kirshner, R. P., Oemler, A., J., Schechter, P. L., & Shectman, S. A.
1981, , 248, L57, 10.1086/183623

[Lavaux & Wandelt, 2012]
Lavaux, G., & Wandelt, B. D. 2012, , 754, 109,

[Maksimova et al., 2021]
Maksimova, N. A., Garrison, L. H., Eisenstein, D. J., et al. 2021, Monthly
Notices of the Royal Astronomical Society, 508, 4017,

[Mao et al., 2017]
Mao, Q., Berlind, A. A., Scherrer, R. J., et al. 2017, , 835, 160,

[Melchior et al., 2014]
Melchior, P., Sutter, P. M., Sheldon, E. S., Krause, E., & Wandelt,
B. D. 2014, , 440, 2922, 10.1093/mnras/stu456

[Nadathur et al., 2019]
Nadathur, S., Carter, P. M., Percival, W. J., Winther, H. A., &
Bautista, J. E. 2019, , 100, 023504,

[Neyrinck, 2008]
Neyrinck, M. C. 2008, , 386, 2101,

[Patiri et al., 2006]
Patiri, S. G., Prada, F., Holtzman, J., Klypin, A., & Betancort-Rijo, J. 2006,
Monthly Notices of the Royal Astronomical Society, 372, 1710,

[Pisani et al., 2015]
Pisani, A., Sutter, P. M., Hamaus, N., et al. 2015, , 92, 083531,

[Planck Collaboration et al., 2020]
Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020, , 641,
A6, 10.1051/0004-6361/201833910

[Press & Schechter, 1974]
Press, W. H., & Schechter, P. 1974, , 187, 425, 10.1086/152650

[Ricciardelli et al., 2014]
Ricciardelli, E., Quilis, V., & Varela, J. 2014, , 440, 601,

[Rojas et al., 2005]
Rojas, R. R., Vogeley, M. S., Hoyle, F., & Brinkmann, J. 2005, ,
624, 571, 10.1086/428476

[Sheth & van de Weygaert, 2004]
Sheth, R. K., & van de Weygaert, R. 2004, Monthly Notices of the Royal
Astronomical Society, 350, 517, 10.1111/j.1365-2966.2004.07661.x

[Sutter et al., 2012]
Sutter, P. M., Lavaux, G., Wandelt, B. D., & Weinberg, D. H. 2012,
, 761, 187, 10.1088/0004-637X/761/2/187

[Sutter et al., 2014]
Sutter, P. M., Pisani, A., Wandelt, B. D., & Weinberg, D. H. 2014,
, 443, 2983, 10.1093/mnras/stu1392

[Sutter et al., 2015]
Sutter, P. M., Lavaux, G., Hamaus, N., et al. 2015, Astronomy and
Computing, 9, 1, 10.1016/j.ascom.2014.10.002

[van de Weygaert & Platen, 2011]
van de Weygaert, R., & Platen, E. 2011, in International Journal of Modern
Physics Conference Series, Vol. 1, International Journal of Modern Physics
Conference Series, 41–66, 10.1142/S2010194511000092

[Verza et al., 2019]
Verza, G., Pisani, A., Carbone, C., Hamaus, N., & Guzzo, L. 2019,
, 2019, 040, 10.1088/1475-7516/2019/12/040

[Zaidouni et al., 2023]
Zaidouni, F., Veyrat, D., Douglass, K. A., & BenZvi, S. 2023, In

[Zeldovich, 1970]
Zeldovich, Y. B. 1970, Astron. Astrophys. 5: 84-9(Mar 1970).

[Zhao et al., 2020]
Zhao, C., Chuang, C.-H., Kitaura, F.-S., et al. 2020, , 491,
4554, 10.1093/mnras/stz3339

[Zhao et al., 2021]
Zhao, C., Variu, A., He, M., et al. 2021, arXiv e-prints,