# On the probability of generating a primitive matrix††thanks: Jingwei CHEN.
Chongqing Key Laboratory of Automated Reasoning and Cognition, Chongqing
Institute of Green and Intelligent Technology, Chinese Academy of Sciences,
Chongqing $400714$, China; Chongqing College, University of Chinese Academy of
Sciences, Chongqing $400714$, China; Key Laboratory of Advanced Manufacturing
Technology of Ministry of Education, Guizhou Universityy, Guiyang $550025$,
China. Email<EMAIL_ADDRESS>
Yong FENG, Wenyuan WU. Chongqing Key Laboratory of Automated Reasoning and
Cognition, Chongqing Institute of Green and Intelligent Technology, Chinese
Academy of Sciences, Chongqing $400714$, China; Chongqing College, University
of Chinese Academy of Sciences, Chongqing $400714$, China. Email: {yongfeng,
<EMAIL_ADDRESS>
Yang LIU (Corresponding author). Information Science and Engineering,
Chongqing Jiaotong University, Chongqing 400074, China. Email:
<EMAIL_ADDRESS>
This research was supported by National Key Research and Development Project
(2020YFA0712303), NSFC (61903053), Youth Innovation Promotion Association of
CAS, Guizhou Science and Technology Program [2020]4Y056 and Chongqing Science
and Technology Program (cstc2021jcyj-msxmX0821, cstc2020yszx-jcyjX0005,
cstc2021yszx-jcyjX0004, 2022YSZX-JCX0011CSTB).

Jingwei CHEN  Yong FENG  Yang LIU  Wenyuan WU

###### Abstract

Given a $k\times n$ integer primitive matrix $\bm{A}$ (i.e., a matrix can be
extended to an $n\times n$ unimodular matrix over the integers) with the
maximal absolute value of entries $\|\bm{A}\|$ bounded by an integer $\lambda$
from above, we study the probability that the $m\times n$ matrix extended from
$\bm{A}$ by appending other $m-k$ row vectors of dimension $n$ with entries
chosen randomly and independently from the uniform distribution over
$\\{0\mbox{,\ }1\mbox{,\ }\ldots\mbox{,\ }\lambda-1\\}$ is still primitive. We
present a complete and rigorous proof of a lower bound on the probability,
which is at least a constant for fixed $m$ in the range $[k+1,n-4]$. As an
application, we prove that there exists a fast Las Vegas algorithm that
completes a $k\times n$ primitive matrix $\bm{A}$ to an $n\times n$ unimodular
matrix within expected $\widetilde{O}(n^{\omega}\log\|\bm{A}\|)$ bit
operations, where $\widetilde{O}$ is big-$O$ but without log factors, $\omega$
is the exponent on the arithmetic operations of matrix multiplication.

Keywords: Integer matrix, unimodular matrix, matrix completion, probabilistic
algorithm.

##  1 Introduction

A vector $\bm{x}\in\mathbb{Z}^{n}$ is called primitive if $\bm{x}=d\bm{y}$ for
$\bm{y}\in\mathbb{Z}^{n}$ and $d\in\mathbb{Z}$ implies $d=\pm 1$. More
generally, a matrix $\bm{A}\in\mathbb{Z}^{k\times n}$ with $k\leq n$ is called
primitive if $\bm{x}=\bm{yA}\in\mathbb{Z}^{n}$ for $\bm{y}\in\mathbb{Q}^{k}$
implies $\bm{y}\in\mathbb{Z}^{k}$; in this case we also say the $k$ rows of
$\bm{A}$ are primitive in $\mathbb{Z}^{n}$. In particular, an $n\times n$
primitive matrix over $\mathbb{Z}$ is also called _unimodular_ , i.e., an
integer square matrix with determinant $\pm 1$. It can be proved that a
$k\times n$ primitive matrix can always be extended to an $n\times n$
unimodular matrix over $\mathbb{Z}$; see, e.g., [MRW11].

Given a primitive matrix $\bm{A}\in\mathbb{Z}^{k\times n}$ with
$\|\bm{A}\|:=\max_{i,j}\lvert a_{i,j}\rvert$ bounded by an integer $\lambda$
from above, our focus in this paper will be on the probability of that the
$m\times n$ matrix extended from $\bm{A}$ by appending other $m-k$ vectors of
dimension $n$ with entries chosen randomly and independently from the uniform
distribution over $\Lambda:=\mathbb{Z}\cap[0\mbox{,\ }\lambda)$ is still
primitive.

###  1.1 Main results

In particular, we prove the following

######  Theorem 1.

Given a primitive matrix $\bm{A}\in\mathbb{Z}^{k\times n}$ with $\|\bm{A}\|$
bounded by an integer $\lambda$ from above and an integer $s$ with $0\leq
s\leq n-k-2$, let $\bm{B}\in\mathbb{Z}^{(n-s-1)\times n}$ be a matrix with
first $k$ rows copied from $\bm{A}$ and entries of the other rows chosen
randomly and independently from the uniform distribution over $\Lambda$. Then
the probability of the event that $\bm{B}$ is primitive is at least

$1-4\left(\frac{2}{3}\right)^{s+1}\left(1-\left(\frac{2}{3}\right)^{n-k-s-1}\right)-\dfrac{2(n-s)^{2}}{\lambda^{s+2}}\left(1-\frac{1}{\lambda^{n-k-s-1}}\right).$
(1)

Note that Theorem 1 holds for $k=0$ as well, which is the case of directly
choosing $n-s-1$ vectors with entries from the uniform distribution over
$\Lambda$. Roughly speaking, for this case, it was shown by Maze et al.
[MRW11] that when $\lambda\rightarrow\infty$, the limit probability of that an
$(n-s-1)\times n$ integer matrix with entries random chosen from the uniform
distribution over $\Lambda$ is primitive is

$\prod_{j=s+2}^{n}\dfrac{1}{\zeta(j)},$ (2)

where $\zeta(\cdot)$ is the Riemann’s zeta function. In this sense, Theorem 1
gives an effective lower bound on the probability for finite $\lambda$, and
hence will be useful in practice, especially in computer science.

From Eq. (1), the parameter $k$ plays a very limited role for the result. In
fact, one may easily obtain a simpler but worse bound:

$\displaystyle
1-4\left(\frac{2}{3}\right)^{s+1}-\frac{2(n-s)^{2}}{\lambda^{s+2}}.$ (3)

Note that this bound is independent of the parameter $k$. For example, if
$s=3$ is fixed and $\lambda\geq 3(n-3)^{2/5}$, then the bound (3) implies that
the resulting $(n-4)\times n$ matrix will be primitive with a probability at
least $0.2$.

Moreover, if $\lambda$ is large enough (with respect to $n$), this bound can
be further simplified as

$\displaystyle 1-(4+\delta)\left(\frac{2}{3}\right)^{s+1}$ (4)

for some $0<\delta<1$. Surprisingly, this oversimplified bound only depends on
$s$.

For given $k\mbox{,\ }n$ and $\lambda$, one can decide the smallest integer
$s\in[0,n-k-2]$ such that the lower bound given in Eq. (1) is a usable bound,
i.e., between $0$ and $1$. For instance, $s$ should be at least $3$ for
$4\left(\frac{2}{3}\right)^{s+1}<1$.

We remark that when $s=n-k-2$ the probability bound given in Theorem 1 matches
the empirical probability well according to our experiments in Section 4.

In addition, one may not further expect a constant probability for the case of
$s=-1$ (that corresponds to the resulting matrix is an $n\times n$ unimodular
matrix), since the natural density of random $n\times n$ unimodular matrices
is $0$; see [MRW11, Lemma 5].

###  1.2 Implications

As an application of Theorem 1, we present a fast Las Vegas algorithm
(Algorithm 2) that efficiently completes a primitive matrix
$\bm{A}\in\mathbb{Z}^{k\times n}$ to an $n\times n$ unimodular matrix $\bm{U}$
such that $\|\bm{U}\|\leq n^{O(1)}\|\bm{A}\|$. More specifically, we prove the
following

######  Theorem 2.

Given a primitive matrix $\bm{A}\in\mathbb{Z}^{k\times n}$, there exists a Las
Vegas algorithm that completes $\bm{A}$ to an $n\times n$ unimodular matrix
$\bm{U}$ such that $\|\bm{U}\|\leq n^{8}\|\bm{A}\|$ in an expected number of
$O(n^{\omega+\varepsilon}\log^{1+\varepsilon}\|\bm{A}\|)$ bit operations.

###  1.3 Techniques

The essential ingredient of our proof for Theorem 1 is adapted from [EGV00,
Section 6], which was used to analyze the expected number of nontrivial
invariant factors of a random integer matrix. The main idea is to give an
upper bound on the probability that the resulting $(n-s-1)\times n$ matrix is
not primitive. Based on the following lemma (whose proof can be obtained from,
e.g., [CS05, Fact 3]), the bound can be analyzed by localizing at each prime
$p$.

######  Lemma 1.

A $k\times n$ integer matrix $\bm{A}$ is not primitive if and only if there
exists at least one prime number $p$ such that the resulting matrix is not
full rank over the finite field $\mathbb{Z}_{p}:=\mathbb{Z}/p\mathbb{Z}$.

In addition, the algorithm for unimodular matrix completion is based on the
determinant reduction technique, which was originally introduced by Storjohann
in [Sto03, Section 15] for computing the determinant of a polynomial matrix,
with a worked example for integer matrix followed. More details about
determinant reduction for integer matrices and an iterated usage of this
technique are given in [Sto05, Section 13.2]. We give a full description of
the algorithm for the integer matrix case and present a detailed analysis in
Section 3.

###  1.4 Related work

Primitive and unimodular matrices have many applications in different areas.
For example, unimodular matrices can be applied to signal compression [PL02];
the lattice reduction algorithms [LLL82, Sch87] essentially produce a series
of unimodular matrices (linear transformations) to improve the basis quality
of the given lattice.

In particular, generating a primitive or unimodular matrix with given rows or
columns happens quite often in practice. For instance, one may need to
generate unimodular matrices with at least one column of all ones in linear
programming for simplex pivoting. In the literature, there exist many results
on what conditions should be satisfied for that a partial integral matrix can
be completed to a unimodular matrix. In 1956 Reiner [Rei56] proved that a row
vector can be completed to a unimodular matrix if and only if it is primitive.
Zhan [Zha06] proved that if $n$ entries of an $n\times n$ partial integral
matrix are prescribed and these $n$ entries do not constitute a row or a
column, then this matrix can be completed to a unimodular matrix. Fang [Fan07]
improved Zhan’s result by proving that if an $n\times n$ partial integral
matrix has $2n-3$ prescribed entries and any $n$ entries of these do not
constitute a row or a column, then it can be completed to a unimodular matrix.
Duffner and Silva [DS17] gave necessary and sufficient conditions for the
existence of unimodular matrices with a prescribed submatrix over a ring that
either is Hermite and Dedekind finite or has stable range one.

For the probability analysis, Maze et al. [MRW11] analyzed the natural density
of $k\times n$ primitive matrices. Guo et al. [GHL16] extended Maze et al.’s
result to a more general setting, where the natural density of $k\times n$
primitive matrices over all $m\times m$ ($m\geq\max\\{k\mbox{,\ }n\\}$)
integer matrices was considered. Note that the natural density can be
interpreted as a limit probability, where each matrix entry is randomly and
independently chosen from a set with an upper bound but the bound tends to
infinity. However, a finite version is usually preferable in practice, e.g.,
algorithm analysis. Therefore, the result in Theorem 1 will be useful.

Additionally, a somewhat “dual” case is considered in [FW14], where the
probability of that $m$ integeral vectors with bounded entries generate a same
lattice of rank $n$ was studied. In [FW14], the ideal choice is $m=n+1$, but a
theoretical lower bound on the probability was only proven for $m\geq 2n+1$.
Aggarwal and Regev in [AR16] and Kirshanova et al. in [KNSW20] considered a
closely related problem but the entries are randomly chosen from the discrete
gaussian distribution over $\mathbb{Z}^{n}$, where $m$ is even larger.

For algorithms, Randall [Ran91] presented an algorithm for generating random
matrices over a finite field with a given determinant, which naturally can be
used to generate unimodular matrices over finite fields. Kalaimani et al.
[KBS13] and Zhou and Labahn [ZL14] discussed algorithms for unimodular
completion of polynomial matrices. The unimodular matrix completion algorithm
discussed in this paper works for integer matrices, which is more efficient
than the standard method for this problem (see Remark 3 and 4).

Roadmap. We prove Theorem 1 in Section 2. Combining the determinant reduction
technique for integer matrices with Theorem 1, we give a fast algorithm for
the problem of unimodular matrix completion in Section 3. In Section 4, we
present an extensive experimental study on the probability that the resulting
$(n-s-1)\times n$ matrix is primitive and discuss some interesting problems
for further study.

##  2 Proof of Theorem 1

Given a primitive matrix $\bm{A}\in\mathbb{Z}^{k\times n}$ with $k<n$ and
$\|\bm{A}\|\leq\lambda$, we now consider to extend $\bm{A}$ to an
$(n-s-1)\times n$ matrix by choosing other $n-k-s-1$ vectors with entries
randomly and independently chosen from the uniform distribution over
$\Lambda=\mathbb{Z}\cap[0\mbox{,\ }\lambda)$, where the integer $s$ satisfies
$0\leq s\leq n-k-2$. Denote by $\bm{a}_{i}=(a_{i,j})_{1\leq j\leq n}$ the
$i$-th row of $\bm{A}$. Then $\|\bm{a}_{i}\|_{\infty}:=\max_{j}\\{\lvert
a_{i,j}\rvert\\}\leq\lambda$. We always assume that $\lambda\geq 2$ for
excluding the case of $\Lambda=\\{0\\}$. For convenience, we still use
$\bm{a}_{k+1}\mbox{,\ }\ldots\mbox{,\ }\bm{a}_{n-s-1}\in\mathbb{Z}^{n}$ to
denote the random vectors with each entry chosen randomly and independently
from the uniform distribution over $\Lambda$, and denote

$\displaystyle\bm{A}_{i}=\begin{pmatrix}\bm{a}_{1}\\\ \bm{a}_{2}\\\ \vdots\\\
\bm{a}_{i}\end{pmatrix}=\begin{pmatrix}a_{1,1}&a_{1,2}&\cdots&a_{1,n}\\\
a_{2,1}&a_{2,2}&\cdots&a_{2,n}\\\ \vdots&\vdots&&\vdots\\\
a_{i,1}&a_{i,2}&\cdots&a_{i,n}\end{pmatrix}\mbox{,\ }\quad i=k\mbox{,\
}k+1\mbox{,\ }\ldots\mbox{,\ }n-s-1.$

To prove Theorem 1, we firstly need to bound from above the probability $P$ of
the event that the matrix $\bm{A}_{n-s-1}$ is not primitive under the
assumption that the given matrix $\bm{A}$ is primitive.

######  Lemma 2.

Let all notations be as above. Then

$\displaystyle P\leq
4\left(\frac{2}{3}\right)^{s+1}\left(1-\left(\frac{2}{3}\right)^{n-k-s-1}\right)+\dfrac{2(n-s)^{2}}{\lambda^{s+2}}\left(1-\frac{1}{\lambda^{n-k-s-1}}\right).$

We prove Lemma 2 following the approach of Eberly et al. [EGV00, Section 6],
whose original goal was to bound the expected number of invariant factors for
random integer matrices. For $k\leq i\leq n-s-1$, we define the event

  1. $\bullet$

$\mathsf{MDep}_{i}$: There exists at least one prime number $p$ such that
$\mathrm{rank}(\bm{A}_{i})\leq i-1$ over $\mathbb{Z}_{p}$.

######  Remark 1.

The definition of $\mathsf{MDep}_{i}$ here is different from that in [EGV00,
Section 6], where $\mathsf{MDep}_{i}$ denotes the event that there exists at
least one prime $p$ such that $\mathrm{rank}(\bm{a}_{1}\mbox{,\
}\ldots\mbox{,\ }\bm{a}_{i})\leq i-2$ over $\mathbb{Z}_{p}$.

So the assumption that $\bm{A}$ is primitive is equivalent to that the event
$\neg\mathsf{MDep}_{k}$ happens. Under the assumption, we have

$\displaystyle P$
$\displaystyle=\Pr[\mathsf{MDep}_{n-s-1}]\leq\Pr[\mathsf{MDep}_{k+1}\vee\mathsf{MDep}_{k+2}\vee\cdots\vee\mathsf{MDep}_{n-s-1}]$
(5)
$\displaystyle=\Pr[\mathsf{MDep}_{k+1}\vee(\mathsf{MDep}_{k+2}\wedge\neg\mathsf{MDep}_{k+1})\vee\cdots\vee(\mathsf{MDep}_{n-s-1}\wedge\neg\mathsf{MDep}_{n-s-2})]$
$\displaystyle=\Pr[(\mathsf{MDep}_{k+1}\wedge\neg\mathsf{MDep}_{k})\vee(\mathsf{MDep}_{k+2}\wedge\neg\mathsf{MDep}_{k+1})\vee\cdots\vee(\mathsf{MDep}_{n-s-1}\wedge\neg\mathsf{MDep}_{n-s-2})]$
$\displaystyle\leq\sum_{i=k+1}^{n-s-1}\Pr[\mathsf{MDep}_{i}\wedge\neg\mathsf{MDep}_{i-1}]\leq\sum_{i=k+1}^{n-s-1}\Pr[\mathsf{MDep}_{i}\,|\,\neg\mathsf{MDep}_{i-1}].$

In order to bound $\Pr[\mathsf{MDep}_{i}\,|\,\neg\mathsf{MDep}_{i-1}]$ for
$k+1\leq i\leq n-s-1$, we need to introduce another useful event:

  1. $\bullet$

$\mathsf{Dep}_{i}$: The rows of $\bm{A}_{i}$ are linearly dependent over
$\mathbb{Q}$, i.e., $\mathrm{rank}(\bm{A}_{i})\leq i-1$ over $\mathbb{Q}$.

Now we have

$\displaystyle\Pr[\mathsf{MDep}_{i}\,|\,\neg\mathsf{MDep}_{i-1}]=$
$\displaystyle\Pr[\mathsf{MDep}_{i}\wedge(\mathsf{Dep}_{i}\vee\neg\mathsf{Dep}_{i})\,|\,\neg\mathsf{MDep}_{i-1}]$
(6) $\displaystyle=$
$\displaystyle\Pr[(\mathsf{MDep}_{i}\wedge\mathsf{Dep}_{i})\vee(\mathsf{MDep}_{i}\wedge\neg\mathsf{Dep}_{i})\,|\,\neg\mathsf{MDep}_{i-1}]$
$\displaystyle\leq$
$\displaystyle\Pr[(\mathsf{MDep}_{i}\wedge\mathsf{Dep}_{i})\,|\,\neg\mathsf{MDep}_{i-1}]+\Pr[(\mathsf{MDep}_{i}\wedge\neg\mathsf{Dep}_{i})\,|\,\neg\mathsf{MDep}_{i-1}].$

We first bound

$\Pr[(\mathsf{MDep}_{i}\wedge\mathsf{Dep}_{i})\,|\,\neg\mathsf{MDep}_{i-1}]\leq\Pr[\mathsf{Dep}_{i}\,|\,\neg\mathsf{MDep}_{i-1}].$

The latter is the probability of that $\mathrm{rank}(\bm{A}_{i})=i-1$ over
$\mathbb{Q}$ on the condition that $\bm{A}_{i-1}$ is primitive. From
$\bm{A}_{i-1}$ is primitive, it follows that there must exist $i-1$ columns of
$\bm{A}_{i-1}$ such that the submatrix consisting of the first $i-1$ rows of
these $i-1$ columns has rank $i-1$ over $\mathbb{Q}$. Denote by the set of
indices of these columns $C_{i-1}$. Now $\mathrm{rank}(\bm{A}_{i})=i-1$ over
$\mathbb{Q}$ implies that for all $j\notin C_{i-1}$, the entry $a_{i,j}$ must
be a linear combination of $a_{\ell,j}$’s for $\ell=1,\ldots,i-1$ with the
same rational coefficients that determined by those $a_{i,j}$’s with $j\in
C_{i-1}$. However, each entry of $\bm{a}_{i}$ are chosen randomly and
independently from the uniform distribution over $\Lambda$. Thus, for each
$j\notin C_{i-1}$, the likelihood that $a_{i,j}$ is equal to such a rational
linear combination is either $0$ or $\frac{1}{\lambda}$. Therefore

$\displaystyle\Pr[(\mathsf{MDep}_{i}\wedge\mathsf{Dep}_{i})\,|\,\neg\mathsf{MDep}_{i-1}]\leq\Pr[\mathsf{Dep}_{i}\,|\,\neg\mathsf{MDep}_{i-1}]\leq\left(\frac{1}{\lambda}\right)^{n-i+1}.$
(7)

To bound
$\Pr[(\mathsf{MDep}_{i}\wedge\neg\mathsf{Dep}_{i})\,|\,\neg\mathsf{MDep}_{i-1}]$,
let us consider primes $p<\lambda$ and primes $p\geq\lambda$, respectively.
For that, we define the following events:

  1. $\bullet$

$\mathsf{MDep}_{i}^{(p)}$: For prime $p$, $\mathrm{rank}(\bm{A}_{i})\leq i-1$
over $\mathbb{Z}_{p}$.

  2. $\bullet$

$\mathsf{MDep}_{i}^{(p<\lambda)}$: there exists a prime $p<\lambda$ such that
$\mathrm{rank}(\bm{A}_{i})\leq i-1$ over $\mathbb{Z}_{p}$.

###  2.1 The case of $p<\lambda$

We bound
$\Pr[(\mathsf{MDep}_{i}^{(p<\lambda)}\wedge\mathsf{Dep}_{i})\,|\,\neg\mathsf{MDep}_{i-1}]$
by $\Pr[\mathsf{MDep}_{i}^{(p<\lambda)}\,|\,\neg\mathsf{MDep}_{i-1}]$. The
latter one is the probability that there exists a prime $p<\lambda$ such that
$\mathrm{rank}(\bm{A}_{i})=i-1$ over $\mathbb{Z}_{p}$ under the condition that
$\bm{A}_{i-1}$ is primitive. Now we bound the latter one case by case.

If $\lambda=2$, then no such prime $p$ exists.

If $\lambda=3$, then $p=2$. Furthermore, the event
$(\mathsf{MDep}_{i}^{(p<3)}\,|\,\neg\mathsf{MDep}_{i-1})$ is that
$\mathrm{rank}(\bm{A}_{i})=i-1$ over $\mathbb{Z}_{2}$ assuming
$\mathrm{rank}(\bm{A}_{i-1})=i-1$ over $\mathbb{Z}_{p}$ for any prime $p$. It
follows from the assumption that there exist $i-1$ columns of $\bm{A}_{i}$
such that the submatrix consisting of the first $i-1$ rows of these $i-1$
columns has rank $i-1$ over $\mathbb{Z}_{2}$. Denote by the set of indices of
these columns $C_{i-1}$. Now, $\mathrm{rank}(\bm{A}_{i})=i-1$ for $p=2$ means
that for each $j\notin C_{i-1}$, the entry $a_{i,j}$ must be a linear
combination of $a_{\ell,j}^{\prime}s$ over $\mathbb{Z}_{2}$ for
$\ell=1,\ldots,i-1$. However, when $\lambda=3$ and $p=2$, we have

$\displaystyle\Pr_{x\leftarrow\Lambda}\left[x\equiv 0\mod
p\right]=\frac{2}{3}$

and

$\displaystyle\Pr_{x\leftarrow\Lambda}\left[x\equiv 1\mod
p\right]=\frac{1}{3},$

so for $\lambda=3$ it follows that

$\displaystyle\Pr\left[\mathsf{MDep}_{i}^{(p<3)}\,|\,\neg\mathsf{MDep}_{i-1}\right]\leq\left(\frac{2}{3}\right)^{n-i+1}.$

If $\lambda=4$, then $p=2$ or $p=3$. Similarly, we have

$\displaystyle\Pr\left[\mathsf{MDep}_{i}^{(p<4)}\,|\,\neg\mathsf{MDep}_{i-1}\right]\leq
2\left(\frac{1}{2}\right)^{n-(i-1)}\leq\left(\frac{2}{3}\right)^{n-i+1}$

for $i<n-1$.

If $\lambda=5$, then $p=2$ or $p=3$, and further we obtain

$\displaystyle\Pr\left[\mathsf{MDep}_{i}^{(p<5)}\,|\,\neg\mathsf{MDep}_{i-1}\right]\leq\left(\frac{3}{5}\right)^{n-(i-1)}+\left(\frac{2}{5}\right)^{n-(i-1)}\leq\left(\frac{2}{3}\right)^{n-i+1}$

for $i<n-1$.

If $\lambda=6$, then $p=2$ or $p=3$ or $p=5$ and we have

$\displaystyle\Pr\left[\mathsf{MDep}_{i}^{(p<6)}\,|\,\neg\mathsf{MDep}_{i-1}\right]\leq\left(\frac{1}{2}\right)^{n-i+1}+2\left(\frac{1}{3}\right)^{n-(i-1)}\leq\left(\frac{2}{3}\right)^{n-i+1}$

for $i<n-1$.

If $\lambda=7$, then $p=2$ or $p=3$ or $p=5$. It follows that

$\displaystyle\Pr\left[\mathsf{MDep}_{i}^{(p<7)}\,|\,\neg\mathsf{MDep}_{i-1}\right]\leq\left(\frac{4}{7}\right)^{n-i+1}+\left(\frac{3}{7}\right)^{n-(i-1)}+\left(\frac{2}{7}\right)^{n-(i-1)}\leq\left(\frac{2}{3}\right)^{n-i+1}$

for $i<n-1$.

If $\lambda\geq 8$, then $p=2$ or $p=3$ or $p=5$ or $p=7$, etc. Then for
$i<n-1$, we have

$\displaystyle\Pr\left[\mathsf{MDep}_{i}^{(p<\lambda)}\,|\,\neg\mathsf{MDep}_{i-1}\right]$
$\displaystyle\leq\left(\frac{1}{2}\right)^{n-i+1}+\left(\frac{3}{8}\right)^{n-(i-1)}+\left(\frac{1}{4}\right)^{n-(i-1)}+\sum_{p\geq
7}\left(\frac{2}{p-1}\right)^{n-i+1}$
$\displaystyle\leq\left(\frac{2}{3}\right)^{n-i+1}+\left(\frac{1}{3}\right)^{n-i+1}\cdot\sum_{p\geq
7}\frac{4}{(p-1)^{2}}$
$\displaystyle\leq\left(\frac{2}{3}\right)^{n-i+1}+4\cdot\left(\frac{1}{3}\right)^{n-i+1}\cdot\sum_{j=6}^{\infty}\frac{1}{j^{2}}$
$\displaystyle\leq\left(\frac{2}{3}\right)^{n-i+1}+4\cdot\left(\frac{1}{3}\right)^{n-i+1}\cdot\left(\zeta(2)-1-\frac{1}{4}-\frac{1}{9}-\frac{1}{16}-\frac{1}{25}\right)$
$\displaystyle\leq\left(\frac{2}{3}\right)^{n-i+1}+\frac{3}{4}\left(\frac{1}{3}\right)^{n-i+1},$

where $\zeta(\cdot)$ is Riemman’s zeta function and the fact that
$\left\lceil\lambda/p\right\rceil/\lambda$ does not increase with respect to
$\lambda$ and that $\left\lceil\lambda/p\right\rceil/\lambda\leq 2/(p-1)$ for
$p<\lambda$ were used. Therefore, for the case of $p<\lambda$ we proved the
following

######  Proposition 1.

Let $\lambda\geq 2$ be an integer and $k+1\leq i\leq n-3$, and suppose that
the event $\neg\mathsf{MDep}_{i-1}$ happens. The probability that there exists
any prime $p<\lambda$ such that $\mathrm{rank}(\bm{A}_{i})\leq i-1$ over
$\mathbb{Z}_{p}$ is at most
$\left(\frac{2}{3}\right)^{n-i+1}+\frac{3}{4}\left(\frac{1}{3}\right)^{n-i+1}$.
In particular,

$\displaystyle\Pr\left[(\mathsf{MDep}_{i}^{(p<\lambda)}\wedge\neg\mathsf{Dep}_{i})\,|\,\neg\mathsf{MDep}_{i-1}\right]\leq\left(\frac{2}{3}\right)^{n-i+1}+\frac{3}{4}\left(\frac{1}{3}\right)^{n-i+1}.$
(8)

###  2.2 The case of $p\geq\lambda$

Now our goal is to bound

$\displaystyle\begin{aligned}
\Pr\left[\left(\mathsf{MDep}_{i}^{(p\geq\lambda)}\wedge\neg\mathsf{Dep}_{i}\right)\,\bigg{|}\,\neg\mathsf{MDep}_{i-1}\right]&=\Pr\left[\left(\bigvee_{p\geq\lambda}\mathsf{MDep}_{i}^{(p)}\right)\wedge\neg\mathsf{Dep}_{i}\,\bigg{|}\,\neg\mathsf{MDep}_{i-1}\right]\\\
&=\Pr\left[\bigvee_{p\geq\lambda}\left(\mathsf{MDep}_{i}^{(p)}\wedge\neg\mathsf{Dep}_{i}\right)\,\bigg{|}\,\neg\mathsf{MDep}_{i-1}\right]\\\
&\leq\sum_{p\geq\lambda}\Pr\left[\left(\mathsf{MDep}_{i}^{(p)}\wedge\neg\mathsf{Dep}_{i}\right)\,\bigg{|}\,\neg\mathsf{MDep}_{i-1}\right]\end{aligned}$
(9)

where $p$ ranges all primes at least $\lambda$.

If $p\geq\lambda$ is a fixed prime, then the probability that
$\mathrm{rank}(\bm{A}_{i})\leq i-1$ over $\mathbb{Z}_{p}$ under the condition
$\mathrm{rank}(\bm{A}_{i-1})=i-1$ over $\mathbb{Q}$ is at most
$\left(\frac{1}{\lambda}\right)^{n-i+1}$, since the probability that a value
chosen randomly and independently from the uniform distribution on
$\Lambda=\mathbb{Z}\cap[0,\lambda)$ equals a given value in $\mathbb{Z}_{p}$
is either $0$ or $\frac{1}{\lambda}$, and hence at most $\frac{1}{\lambda}$.
So we have

$\displaystyle\Pr\left[(\mathsf{MDep}_{i}^{(p)}\wedge\neg\mathsf{Dep}_{i})\,|\,\neg\mathsf{MDep}_{i-1}\right]\leq\Pr\left[\mathsf{MDep}_{i}^{(p)}\,|\,\neg\mathsf{MDep}_{i-1}\right]\leq\left(\frac{1}{\lambda}\right)^{n-i+1}.$
(10)

Furthermore, we claim that for all primes $p>(i\cdot\lambda)^{i}$

$\displaystyle\Pr\left[\left(\mathsf{MDep}_{i}^{(p)}\wedge\neg\mathsf{Dep}_{i}\right)\,\bigg{|}\,\neg\mathsf{MDep}_{i-1}\right]=0.$
(11)

In fact, if there exists a prime $p>(i\cdot\lambda)^{i}$ that makes the event
$(\mathsf{MDep}_{i}^{(p)}\wedge\neg\mathsf{Dep}_{i})$ happen assuming that
$\bm{A}_{i-1}$ is primitive, then $p$ must divide all $i\times i$ minors of
$\bm{A}_{i}$. However, although $a_{i,j}$’s $(j=1,2,\ldots,n)$ are chosen
randomly and independently from the uniform distribution on $\Lambda$, once
$\neg\mathsf{Dep}_{i}$ happens, there exists at least one nonsingular $i\times
i$ submatrix in $\bm{A}_{i}$ over $\mathbb{Q}$. The absolute value of the
determinant of such a submatrix is at most
$i!\cdot\lambda^{i}\leq(i\cdot\lambda)^{i}$. Combining
$p>(i\cdot\lambda)^{i}$, $p$ divides the determinant, and the absolute value
of the determinant is bounded by $(i\cdot\lambda)^{i}$, we have the determiant
must be zero, which contradicts with non-singularity of the submatrix. Also,
note that the bound $(i\cdot\lambda)^{i}$ is independent of the choice of
$a_{i,j}$’s when such $a_{i,j}$’s make $\neg\mathsf{Dep}_{i}$ happen.
(Instead, when such $a_{i,j}$’s do not make $\neg\mathsf{Dep}_{i}$ happen, the
probability
$\Pr[(\mathsf{MDep}_{i}\wedge\mathsf{Dep}_{i})\,|\,\neg\mathsf{MDep}_{i-1}]$
is already discussed in (7).)

Now, the number of possible primes $p\geq\lambda$ that divides an integer
$d\leq(i\cdot\lambda)^{i}$ is at most $\log_{\lambda}(i\cdot\lambda)^{i}\leq
i(1+\log_{\lambda}i)$, where $d$ is in fact the greatest common divisor of all
nonsingular $i\times i$ submatrix of $\bm{A}_{i}$ over $\mathbb{Q}$. Although
the integer $d$ may vary as $a_{i,j}$’s are chosen randomly and independently
from the uniform distribution on $\Lambda$, it is always bounded by
$(i\cdot\lambda)^{i}$. Thus, the number of possible primes that we need to
consider is at most $i(1+\log_{\lambda}i)$. Thus, combining Eq. (9-11), we
have

$\displaystyle\begin{aligned}
\Pr\left[(\mathsf{MDep}_{i}^{(p\geq\lambda)}\wedge\neg\mathsf{Dep}_{i})\,|\,\neg\mathsf{MDep}_{i-1}\right]&\leq\sum_{p\geq\lambda}\Pr\left[\left(\mathsf{MDep}_{i}^{(p)}\wedge\neg\mathsf{Dep}_{i}\right)\,\bigg{|}\,\neg\mathsf{MDep}_{i-1}\right]\leq\left(i(1+\log_{\lambda}i)\right)\cdot\left(\frac{1}{\lambda}\right)^{n-i+1}.\end{aligned}$
(12)

Therefore, it follows from Eq. (8) and (12) that

$\displaystyle\begin{aligned}
&\Pr\left[(\mathsf{MDep}_{i}\wedge\neg\mathsf{Dep}_{i})\,|\,\neg\mathsf{MDep}_{i-1}\right]\\\
=&\Pr\left[(\mathsf{MDep}_{i}^{(p<\lambda)}\wedge\neg\mathsf{Dep}_{i})\,|\,\neg\mathsf{MDep}_{i-1}\right]+\Pr\left[(\mathsf{MDep}_{i}^{(p\geq\lambda)}\wedge\neg\mathsf{Dep}_{i})\,|\,\neg\mathsf{MDep}_{i-1}\right]\\\
\leq&\left(\frac{2}{3}\right)^{n-i+1}+\frac{3}{4}\left(\frac{1}{3}\right)^{n-i+1}+\left(i(1+\log_{\lambda}i)\right)\cdot\left(\frac{1}{\lambda}\right)^{n-i+1}.\end{aligned}$
(13)

#### Proof of Lemma 2

It follows from Eq. (6), (7), and (13) that for integer $\lambda\geq 2$,

$\displaystyle\Pr[\mathsf{MDep}_{i}|\neg\mathsf{MDep}_{i-1}]\leq\ $
$\displaystyle\Pr[(\mathsf{MDep}_{i}\wedge\mathsf{Dep}_{i})\,|\,\neg\mathsf{MDep}_{i-1}]+\Pr[(\mathsf{MDep}_{i}\wedge\neg\mathsf{Dep}_{i})\,|\,\neg\mathsf{MDep}_{i-1}]$
$\displaystyle\leq\ $
$\displaystyle\left(\frac{1}{\lambda}\right)^{n-i+1}+\left(\frac{2}{3}\right)^{n-i+1}+\frac{3}{4}\left(\frac{1}{3}\right)^{n-i+1}+\left(i(1+\log_{\lambda}i)\right)\cdot\left(\frac{1}{\lambda}\right)^{n-i+1}$
$\displaystyle\leq\ $
$\displaystyle\left(\frac{2}{3}\right)^{n-i+1}+\frac{3}{4}\left(\frac{1}{3}\right)^{n-i+1}+\left(i(1+\log_{\lambda}i)+1\right)\cdot\left(\frac{1}{\lambda}\right)^{n-i+1}.$

Therefore, under the assumption that $\neg\mathsf{MDep}_{k}$ happens, it
follows from Eq. (5) that

$\displaystyle P$
$\displaystyle=\Pr[\mathsf{MDep}_{n-s-1}]\leq\sum_{i=k+1}^{n-s-1}\Pr[\mathsf{MDep}_{i}|\neg\mathsf{MDep}_{i-1}]$
$\displaystyle\leq\sum_{i=k+1}^{n-s-1}\left(\left(\frac{2}{3}\right)^{n-i+1}+\frac{3}{4}\left(\frac{1}{3}\right)^{n-i+1}+\left(i(1+\log_{\lambda}i)+1\right)\cdot\left(\frac{1}{\lambda}\right)^{n-i+1}\right)$
$\displaystyle\leq
2\cdot\sum_{i=k+1}^{n-s-1}\left(\frac{2}{3}\right)^{n-i+1}+(n-s)(1+\log_{\lambda}(n-s-1))\sum_{i=k+1}^{n-s-1}\frac{1}{\lambda^{n-i+1}}$
$\displaystyle\leq
4\left(\frac{2}{3}\right)^{s+1}\left(1-\left(\frac{2}{3}\right)^{n-k-s-1}\right)+\dfrac{(n-s)^{2}}{(\lambda-1)\lambda^{s+1}}\left(1-\frac{1}{\lambda^{n-k-s-1}}\right)$
$\displaystyle\leq
4\left(\frac{2}{3}\right)^{s+1}\left(1-\left(\frac{2}{3}\right)^{n-k-s-1}\right)+\dfrac{2(n-s)^{2}}{\lambda^{s+2}}\left(1-\frac{1}{\lambda^{n-k-s-1}}\right).\qed$

Now, Theorem 1 is a direct consequence of Lemma 2.

######  Remark 2.

According to the above proof, we use localization at each prime $p$ and count
the number of possibilities over $\mathbb{Z}_{p}$. So the set $\Lambda$ can be
any set containing $\lambda$ contiguous integers, as in [EGV00, Section 6].

##  3 Unimodular matrix completion

When completing a given $k\times n$ primitive matrix to an $n\times n$
unimodular matrix, we first complete the input matrix to an $n\times n$ matrix
with $n-k$ vectors whose entries are chosen randomly and uniformly from the
uniform distribution over $\Lambda$, and then rectify the last four vectors
via the determinant reduction technique [Sto03, Section 15]. Repeat the above
process until the resulting matrix is unimodular.

Eq. (3) (a consequence of Theorem 1) with $s=3$ and $\lambda\geq 3(n-3)^{2/5}$
guarantees that the completed $(n-4)\times n$ matrix is still primitive with a
probability at least a constant, say $0.2$. Then after a constant number of
repeat, it is expected that we will obtain at least one $(n-4)\times 4$
primitive matrix, and the determinant reduction technique applies to produce a
unimodular matrix.

###  3.1 Hermite normal form

A matrix $\bm{A}\in\mathbb{Z}^{m\times n}$ of rank $r$ has a (row) Hermite
normal form (HNF) if there exists a square unimodular matrix $\bm{U}$ such
that $\bm{H}=\bm{UA}$ satisfies the following: there exist indices $1\leq
i_{1}<i_{2}<\cdots<i_{r}\leq n$ such that for $j=1,\ldots,r$, $h_{j,i_{j}}>0$,
$h_{j,k}=0$ if $k<i_{j}$ and $0\leq h_{\ell,i_{j}}<h_{j,i_{j}}$ if $\ell<j$;
the bottom $m-r$ rows of $\bm{H}$ are zero.

The matrix $\bm{H}$ is unique, denoted by $\mathrm{HNF}(\bm{A})$, although the
unimodular tranformations are usually not unique. With the definition of HNF,
we can character the primitive matrix as the following:

######  Lemma 3.

A $k\times n$ integer matrix $\bm{A}$ with $k<n$ is primitive iff
$\mathrm{HNF}(\bm{A}^{\mathrm{T}})=\binom{\bm{I}_{k}}{\bm{0}}$.

One can find a proof of this lemma in, e.g., [MRW11].

###  3.2 Determinant reduction

Given a nonsingular $\bm{A}\in\mathbb{Z}^{n\times n}$, the determinant
reduction introduced in [Sto03, Section 15] computes a matrix
$\bm{B}\in\mathbb{Z}^{n\times n}$, obtained from $\bm{A}$ by replacing the
last column, such that the last diagonal entry in the Hermite normal form of
$\bm{B}$ is one. The determinant reduction was originally presented for
integral polynomial matrix in [Sto03, Section 15], with a worked example for
integer matrix. Here, we give full details and analyses for the integer matrix
case.

Algorithm 1 (Determinant reduction)

0: A nonsingular integer matrix $\bm{A}\in\mathbb{Z}^{n\times n}$.

0: A matrix $\bm{B}\in\mathbb{Z}^{n\times n}$, with $\bm{B}$ equal to $\bm{A}$
except for possibly the last column, $\|\bm{B}\|\leq n^{2}\|\bm{A}\|$, and the
last diagonal entry of $\mathrm{HNF}(\bm{B})$ equal to one.

1: Set $\bm{C}_{n-1}$ to be the matrix consisting of the first $n-1$ columns
of $\bm{A}$. Compute a primitive vector $\bm{u}\in\mathbb{Z}^{n}$ such that
$\bm{u}\bm{C}_{n-1}=\bm{0}$.

2: Call an extended gcd algorithm to compute $\bm{b}\in\mathbb{Z}^{n}$ such
that $\bm{u}\bm{b}^{\mathrm{T}}=1$.

3: Set $\overline{\bm{A}}$ to be the $(n-1)\times(n-1)$ principal submatrix of
$\bm{A}$, and $\overline{\bm{b}}$ the vector consisting of the first $n-1$
entries of $\bm{b}$.

4: Compute
$\overline{\bm{q}}^{\mathrm{T}}:=\left\lceil\overline{\bm{A}}^{-1}\overline{\bm{b}}^{\mathrm{T}}\right\rfloor$.
//For a vector $\bm{v}$, $\left\lceil\bm{v}\right\rfloor$ means each entry
rounded.

5: Set $\bm{q}:=(\overline{\bm{q}}\mbox{,\ }0)$ and set $\bm{B}$ to be
$\bm{A}$ except replacing the last column by
$\bm{b}^{\mathrm{T}}-\bm{A}\bm{q}^{\mathrm{T}}$.

6: return $\bm{B}$.

######  Proposition 2.

Given a nonsingular $n\times n$ integer matrix $\bm{A}$, Algorithm 1 correctly
computes an $n\times n$ integer matrix $\bm{B}$ within
$O(n^{\omega+\varepsilon}\log^{1+\varepsilon}\|\bm{A}\|)$ bit operations,
where $\bm{B}$ satisfies the following:

  * •

$\bm{B}$ equals $\bm{A}$ except for possibly the last column,

  * •

the last diagonal entry of $\mathrm{HNF}(\bm{B})$ equals one,

  * •

$\|\bm{B}\|=O(n^{2}\|\bm{A}\|)$.

###### Proof.

The singularity of $\bm{A}$ implies that the vector $\bm{u}$ produced in Step
1 is the unique primitive vector in
$\\{\bm{x}\in\mathbb{Z}^{n}:\bm{xC}_{n-1}=\bm{0}\\}$. Let $\bm{U}$ be an
arbitrary unimodular matrix such that $\bm{H}:=\mathrm{HNF}(\bm{A})=\bm{UA}$.
Then the uniqueness of $\bm{u}$ implies that $\bm{u}$ must be the last row of
$\bm{U}$. By construction of $\bm{b}$ in Step 2, the matrix obtained from
$\bm{A}$ by replacing the last column with $\bm{b}^{\mathrm{T}}$ will have
$\mathrm{HNF}$ with the last diagonal entry one. The whole algorithm can be
expressed as the following equation

$\displaystyle\bm{B}=\begin{pmatrix}\overline{\bm{A}}&\overline{\bm{b}}{{}^{\mathrm{T}}}\\\
\overline{\bm{a}}&b_{n}\end{pmatrix}\cdot\begin{pmatrix}\bm{I}_{n-1}&-\overline{\bm{q}}{{}^{\mathrm{T}}}\\\
&1\end{pmatrix}=\begin{pmatrix}\overline{\bm{A}}&\overline{\bm{b}}{{}^{\mathrm{T}}}-\overline{\bm{A}}\overline{\bm{q}}{{}^{\mathrm{T}}}\\\
\overline{\bm{a}}&b_{n}-\overline{\bm{a}}\,\overline{\bm{q}}{{}^{\mathrm{T}}}\end{pmatrix},$

where $\overline{\bm{A}}$ is the $(n-1)\times(n-1)$ principal submatrix of
$\bm{A}$, $\overline{\bm{b}}$ is the vector consisting of the first $n-1$
entries of $\bm{b}$, and $\bar{\bm{a}}$ is the vector consisting of the first
$n-1$ entries of the last row of $\bm{A}$. Step 3-5 are to reduce the size of
$\overline{\bm{b}}$ by the columns of $\overline{\bm{A}}$. It follows from
$\bm{u}\bm{C}_{n-1}=\bm{0}$ that all entries of $\bm{u}\bm{A}$ are zero except
for possibly the last one. However, the last entry of $\bm{q}$ equals zero, so
we have $\bm{u}\bm{Aq}^{\mathrm{T}}=0$, and hence

$\displaystyle\bm{u}(\bm{b}^{\mathrm{T}}-\bm{Aq}^{\mathrm{T}})=\bm{u}\bm{b}^{\mathrm{T}}=1,$
(14)

which implies that $\bm{UB}$ is exactly $\bm{H}$ except for the last column
replaced by $(*,\ldots,*,1)^{\mathrm{T}}$, and hence the last diagonal of
$\mathrm{HNF}(\bm{B})$ must be one.

We now consider the size of $\bm{b}^{\mathrm{T}}-\bm{Aq}^{\mathrm{T}}$. First,

$\displaystyle\begin{aligned}
\left\|\overline{\bm{b}}^{\mathrm{T}}-\overline{\bm{A}}\overline{\bm{q}}^{\mathrm{T}}\right\|_{\infty}&=\left\|\overline{\bm{b}}^{\mathrm{T}}-\overline{\bm{A}}\left\lceil\overline{\bm{A}}^{-1}\overline{\bm{b}}^{\mathrm{T}}\right\rfloor\right\|_{\infty}\\\
&=\left\|\overline{\bm{b}}^{\mathrm{T}}-\overline{\bm{A}}\left(\overline{\bm{A}}^{-1}\overline{\bm{b}}^{\mathrm{T}}+\bm{\varepsilon}^{\mathrm{T}}\right)\right\|_{\infty}\leq\left\|\overline{\bm{A}}\bm{\varepsilon}^{\mathrm{T}}\right\|_{\infty}\leq\frac{n-1}{2}\|\overline{\bm{A}}\|,\end{aligned}$
(15)

where $\|\cdot\|_{\infty}$ is the $\ell_{\infty}$-norm of a vector and
$\|\bm{\varepsilon}^{\mathrm{T}}\|_{\infty}\leq\frac{1}{2}$ is used. In fact,
Step 4 is essentially the same as Babai’s rounding algorithm [Bab85]. Denote
$\bm{u}=(\overline{\bm{u}},u_{n})$. Without loss of generality, we can assume
that

$\displaystyle\|\bm{u}\|_{\infty}=u_{n}.$ (16)

Otherwise there must exist a permutation matrix $\bm{P}$ such that
$\bm{u}\bm{P}\bm{P}^{-1}\bm{C}_{n-1}=\bm{0}$, $\bm{u}\bm{P}$ satisfies Eq.
(16), and $\bm{P}^{-1}\bm{C}_{n-1}$ still corresponds to the first $n-1$
columns of $\bm{A}$ but with a certain column permutation. From Eq. (14), we
have
$u_{n}(b_{n}-\overline{\bm{a}}\,\overline{\bm{q}}^{\mathrm{T}})=1-\overline{\bm{u}}(\overline{\bm{b}}^{\mathrm{T}}-\overline{\bm{A}}\overline{\bm{q}}^{\mathrm{T}})$,
combining Eq. (16), which gives

$\displaystyle\left\lvert
b_{n}-\overline{\bm{a}}\,\overline{\bm{q}}^{\mathrm{T}}\right\rvert=\frac{1}{\left\lvert
u_{n}\right\rvert}\cdot\left\lvert
1-\overline{\bm{u}}(\overline{\bm{b}}^{\mathrm{T}}-\overline{\bm{A}}\overline{\bm{q}}^{\mathrm{T}})\right\rvert\leq
1+\frac{\|\overline{\bm{u}}\|_{\infty}}{\lvert
u_{n}\rvert}\cdot\|\overline{\bm{b}}^{\mathrm{T}}-\overline{\bm{A}}\overline{\bm{q}}^{\mathrm{T}}\|_{1}\leq
1+\frac{(n-1)^{2}}{2}\|\overline{\bm{A}}\|,$

where $\|\cdot\|_{1}$ is the $\ell_{1}$-norm of a vector, the first inequality
follows from that
$\left\lvert\bm{x}\bm{y}^{\mathrm{T}}\right\rvert\leq\left\|\bm{x}\right\|_{\infty}\cdot\left\|\bm{y}\right\|_{1}$
holds for any two vectors $\bm{x}$ and $\bm{y}$ of the same dimension, and the
second inequality follows from Eq. (15), (16), and $\|\bm{x}\|_{1}\leq
n\|\bm{x}\|_{\infty}$ for all $n$-dimensional vector $\bm{x}$. Therefore, the
resulting matrix $\bm{B}$ satisfies $\|\bm{B}\|\leq n^{2}\|\bm{A}\|$.

The cost of Algorithm 1 consists in nonsingular rational linear system solving
(Step 1 and 4) that can be finished by a Las Vegas algorithm in an expected
number of $O(n^{\omega+\varepsilon}\log^{1+\varepsilon}\|\bm{A}\|)$ bit
operations [MS04] and an extended gcd computation (Step 2) that can be
accomplished within $O(n^{2+\varepsilon}\log^{1+\varepsilon}\|\bm{A}\|)$ bit
operations [Sto05, Section 13.2]. Totally, Algorithm 1 costs at most
$O(n^{\omega+\varepsilon}\log^{1+\varepsilon}\|\bm{A}\|)$ bit operations. ∎

######  Corollary 1.

Given an $(n-1)\times n$ primitive matrix $\bm{A}$, there exists a Las Vegas
algorithm which completes $\bm{A}$ to an $n\times n$ unimodular matrix in an
expected number of $O(n^{\omega+\varepsilon}\log^{1+\varepsilon}\|\bm{A}\|)$
bit operations.

###### Proof.

Since $\bm{A}$ is primitive, Lemma 3 implies that
$\mathrm{HNF}(\bm{A}^{\mathrm{T}})$ has the form
$\binom{\bm{I}_{n-1}}{\bm{0}}$. Now, one can use the matrix
$\bm{B}=\binom{\bm{A}}{\bm{a}_{n}}^{\mathrm{T}}$ as the input for Algorithm 1,
where $\bm{a}_{n}$ can be an arbitrary integer vector such that the resulting
square matrix is nonsingular. Then, $\mathrm{HNF}(\bm{B})$ has the same form
as the identity matrix of order $n$ except for the last column. Thus, the
transpose of the output matrix of Algorithm 1 will be a unimodular completion
of $\bm{A}$. ∎

######  Remark 3.

To complete an $(n-1)\times n$ primtive matrix to an $n\times n$ unimodular
matrix, a standard method is the following: firstly compute $n$ determinants
of all $(n-1)\times(n-1)$ submatrix of the input matrix and then invoking an
extended euclidean algorithm will give the information of the last row.
However, this standard method can be finished in an expected number of
$O(n^{\omega+1+\varepsilon}\log^{1+\varepsilon}\|\bm{B}\|)$ bit operations,
even using the fast Las Vegas algorithm in [Sto05] for computing the
determinant of an integer matrix in an expected number of
$O(n^{\omega+\varepsilon}\log^{1+\varepsilon}\|\bm{B}\|)$ bit operations.

###  3.3 Proof of Theorem 2

We prove Theorem 2 by Algorithm 2, which uses the Iterated Determinant
Reduction presented in [Sto05, Section 13.2]. Given a nonsingular matrix
$\bm{A}\in\mathbb{Z}^{n\times n}$, the Iterated Determinant Reduction
technique computes a matrix $\bm{B}\in\mathbb{Z}^{n\times n}$, obtained from
$\bm{A}$ by replacing the last $d$ column, such that the last $d$ diagonal
entries in the Hermite normal form of $\bm{B}$ is one. The algorithm consists
of $d$ times calling of Algorithm 1, each followed with multiplying by a
permutation matrix $\bm{P}=\footnotesize\begin{pmatrix}\bm{0}&\bm{I}_{n-1}\\\
1&0\end{pmatrix}$ from right. Now we are ready to present our algorithm for
unimodular matrix completion (Algorithm 2).

Algorithm 2 (Unimodular matrix completion)

0: A primitive matrix $\bm{A}\in\mathbb{Z}^{k\times n}$ with $k<n$ and $n\geq
5$.

0: A unimodular matrix $\bm{U}\in\mathbb{Z}^{n\times n}$, with $\bm{U}$ equal
to $\bm{A}$ except for the last $n-k$ rows and $\|\bm{U}\|\leq
n^{8}\|\bm{A}\|$.

1: Repeat

2: Set $\lambda:=\max\\{\|\bm{A}\|,\lceil 3(n-3)^{2/5}\rceil\\}$ and
$\bm{B}:=\bm{A}^{\mathrm{T}}$.

3: Complete $\bm{B}$ as an $n\times n$ matrix with entries chosen randomly and
independently from the uniform distribution over $\\{0,1\mbox{,\
}\ldots\mbox{,\ }\lambda-1\\}$.

4: Repeat the following four times:

5: Set $\bm{B}$ as the output of Algorithm 1 with input $\bm{B}$ and then set
$\bm{B}:=\bm{BP}$.

6: Until $\det(\bm{B})=\pm 1$

7: return $\bm{U}:=(\bm{B}(\bm{P}^{-1})^{4})^{\mathrm{T}}$.

Denote $\bm{B}_{i}$ by the submatrix of $\bm{B}$ consisting of the first $i$
columns. The goal of the loop in Step 4 is to make the last four diagonal
entries of the Hermite normal form of $\bm{B}$ be one. First, we assume that
$\bm{B}_{n-4}^{\mathrm{T}}$ is primitive at the end of Step 3. Then Lemma 3
implies that

$\mathrm{HNF}(\bm{B})=\begin{pmatrix}\bm{I_{n-4}}&\ast&\ast&\ast&\ast\\\
&\ast&\ast&\ast&\ast\\\ &&\ast&\ast&\ast\\\ &&&\ast&\ast\\\ &&&&\ast\\\
\end{pmatrix}.$

After running Algorithm 1 once, the HNF of the resulting $\bm{B}$ has the
following form

$\mathrm{HNF}(\bm{B})=\begin{pmatrix}\bm{I_{n-4}}&\ast&\ast&\ast&\bm{0}\\\
&\ast&\ast&\ast&0\\\ &&\ast&\ast&0\\\ &&&\ast&0\\\ &&&&1\\\ \end{pmatrix}$

and

$\mathrm{HNF}(\bm{BP})=\begin{pmatrix}\bm{I_{n-3}}&\ast&\ast&\ast\\\
&\ast&\ast&\ast\\\ &&\ast&\ast\\\ &&&\ast\\\ \end{pmatrix}.$

Therefore, at the end of execution of Step 4, we have
$\mathrm{HNF}(\bm{B})=\bm{I_{n}}$ for the resluting $\bm{B}$, which implies
that $\bm{B}$ is unimodular.

Thanks to Eq. (3) with $s=3$ and $\lambda\geq 3(n-3)^{2/5}$, after each
execution of Step 3, $\bm{B}_{n-4}^{\mathrm{T}}$ is primitive with probability
at least $0.2$. Therefore, after $c$ (a constant) executions of Step 3, the
probability that there is at least one $\bm{B}_{n-4}^{\mathrm{T}}$ is
primitive is at least $1-0.8^{c}$. So after a constant number of executions of
the loop Step 1–6, it is expected that the output of Algorithm 2 is a
unimodular completion of $\bm{A}$. After each execution of Step 5, it follows
from Proposition 2 that the magnitude of the new column is bounded by
$O(n^{2}\|\bm{B}\|)$. Therefore, after four executions of Step 5, the
magnitude of the new column is bounded by $O(n^{8}\|\bm{A}\|)$. So the
correctness of Algorithm 2 is proved.

The main cost part of Algorithm 2 is a constant number of executions of the
loop Step 1–6. According to Proposition 2 and the algorithm for computing
determinant in [Sto05], each execution of the loop can be finished in an
expected number of $O(n^{\omega+\varepsilon}\log^{1+\varepsilon}\|\bm{A}\|)$
bit operations. As a consequence, Theorem 2 is proved. ∎

######  Remark 4.

Without the help of Theorem 1, one may use the Iterated Determinant Reduction
algorithm $n-k$ times for unimodular completion. This results in an algorithm
that completes $\bm{A}$ to an $n\times n$ unimodular matrix $\bm{U}$ with
$\|\bm{U}\|\leq n^{2(n-k)}\|\bm{A}\|$ in an expected number of
$O((n-k)n^{\omega+\varepsilon}\log^{1+\varepsilon}\|\bm{A}\|)$ bit operations.

##  4 Experiments and discussions

In this section, we present an experimental study on the empirical probability
for the case that is included in Theorem 1 ($s\geq 3$), and also for the case
that is not included in Theorem 1, e.g., the case of $s<3$. In all the
following experimental data, each empirical probability (labeled Exp.) are
obtained by running $10,000$ random tests on computer algebra system Maple and
counting the success rate, i.e., the proportion of primitive matrices among
all resulting matrices.

###  4.1 The case of $s\geq 3$

In Table 1 and 2 we study the for the case of $k=0$, so that we can compare
the empirical probability (column labeled Exp.), the lower bound given in
Theorem 1 (column labeled Th. 1), and the limit probability given in Eq. (2)
(column labeled Limit probability). We use three different bounds for
$\lambda$, namely, $\lambda=10^{5}$, $\lambda=10^{10}$ and $\lambda=10^{20}$.

Table 1: Average empirical probability vs the probability in Theorem 1 ( $k=0$ and $s=3$) | $\lambda=10^{5}$ | $\lambda=10^{10}$ | $\lambda=10^{20}$ |
---|---|---|---|---
$n$ | Exp. |  Th. 1 | Exp. |  Th. 1 | Exp. |  Th. 1 | Limit probability
$5$ | $0.9652$ | $0.7366$ | $0.9662$ | $0.7366$ | $0.9639$ | $0.7366$ | $0.9643$
$10$ | $0.9335$ | $0.2792$ | $0.9291$ | $0.2792$ | $0.9306$ | $0.2792$ | $0.9334$
$15$ | $0.9292$ | $0.2190$ | $0.9278$ | $0.2190$ | $0.9312$ | $0.2190$ | $0.9325$
$20$ | $0.9338$ | $0.2110$ | $0.9349$ | $0.2110$ | $0.9345$ | $0.2110$ | $0.9325$

Both Table 1 and 2 shows that the empirical probability is relatively near to
the limit probability, both of which are much better than our theoretical
lower bound given in Theorem 1. As indicated previously, $s=3$ is the smallest
$s$ such that the lower bound given in Theorem 1 is between $0$ and $1$ for
these experiments, while $s=n-k-2$ is the largest possible value.

Comparing Table 1 with 2 shows that larger $s$ implies that both higher
experimental success rate and better theoretical bound, however, smaller $s$
implies the resulting matrix is closer to a unimodular matrix. In particular,
for the case of $s=n-k-2$, the probability bound given in Theorem 1 matches
very well with both empirical probability and the limit probability.

Table 2: Average empirical probability vs the probability in Theorem 1 ($k=0$ and $s=n-k-2$) | $\lambda=10^{5}$ | $\lambda=10^{10}$ | $\lambda=10^{20}$ |
---|---|---|---|---
$n$ | Exp. |  Th. 1 | Exp. |  Th. 1 | Exp. |  Th. 1 | Limit probability
$5$ | $0.9652$ | $0.7366$ | $0.9662$ | $0.7366$ | $0.9639$ | $0.7366$ | $0.9643$
$10$ | $0.9995$ | $0.9653$ | $0.9990$ | $0.9653$ | $0.9990$ | $0.9653$ | $0.9990$
$15$ | $0.9999$ | $0.9954$ | $1.0000$ | $0.9954$ | $1.0000$ | $0.9954$ | $0.9999$
$20$ | $1.0000$ | $0.9993$ | $1.0000$ | $0.9993$ | $1.0000$ | $0.9993$ | $0.9999$

In addition, Table 1 and 2 also show that for different $\lambda$ with same
$n$, the theoretical bounds in column Th. 1 are almost the same and the data
of empirical probability are almost the same as well. This is because that for
a large enough $\lambda$, the bound given in Theorem 1 is almost independent
of $\lambda$, as indicated by the oversimplified bound given in (4). For this
reason, we fix $\lambda=10^{5}$ for all other experiments (Table 3–10).

Tables 3–6 are for the case of $k>0$, for which there does not exist a known
limit probability. We generate the initial primitive matrix as follows: We
first generate a $k\times n$ matrix, whose entries are chosen randomly and
independently from the uniform distribution over $[-\lambda\mbox{,\
}\lambda]\cap\mathbb{Z}$. If the matrix is not primitive, we regenerate a new
matrix until it is eventually primitive. For each initial matrix, we complete
it with uniformly random entries from $\Lambda$ to an $(n-s-1)\times n$ matrix
$10,000$ times and count the success rate. From these experiments, we can
observe a similar phenomenon as the case of $k=0$.

Table 3: Average empirical probability vs the probability in Theorem 1 ($k=1$, $s=3$) $n$ | $10$ | $15$ | $20$ | $25$ | $30$
---|---|---|---|---|---
Exp. | $0.9321$ | $0.9283$ | $0.9291$ | $0.9324$ | $0.9310$
Th. 1 | $0.3139$ | $0.2235$ | $0.2116$ | $0.2101$ | $0.2099$
Table 4: Average empirical probability vs the probability in Theorem 1 ($k=1$, $s=n-k-2$) $n$ | $5$ | $10$ | $15$ | $20$ | $25$ | $30$
---|---|---|---|---|---|---
Exp. | $0.8919$ | $0.9969$ | $0.9999$ | $1.0000$ | $1.0000$ | $1.0000$
Th. 1 | $0.6049$ | $0.9479$ | $0.9931$ | $0.9990$ | $0.9998$ | $0.9999$
Table 5: Average empirical probability vs the probability in Theorem 1 ($k=n/2$, $s=3$) $n$ | $16$ | $20$ | $24$ | $28$ | $32$ | $36$
---|---|---|---|---|---|---
Exp. | $0.9349$ | $0.9338$ | $0.9340$ | $0.9312$ | $0.9352$ | $0.9333$
Th. 1 | $0.3659$ | $0.2792$ | $0.2407$ | $0.2235$ | $0.2159$ | $0.2125$
Table 6: Average empirical probability vs the probability in Theorem 1 ($k=n/2$, $s=n-k-2$) $n$ | $16$ | $20$ | $24$ | $28$ | $32$ | $36$
---|---|---|---|---|---|---
Exp. | $0.9919$ | $0.9980$ | $0.9996$ | $0.9999$ | $1.0000$ | $1.0000$
Th. 1 | $0.9219$ | $0.9653$ | $0.9845$ | $0.9931$ | $0.9969$ | $0.9986$

Totally, the probability bound presented in Theorem 1 is tight, especially for
the case of large $s$. However, how to improve the theoretical bound for small
$s$ is an intriguing problem.

###  4.2 The case of $s<3$

For $s<3$, the lower bound on that the resulting $(n-s-1)\times n$ matrix is
primitive given in Theorem 1 will be negative and hence useless. Therefore, it
would be very interesting and useful to obtain a lower bound for the case of
$s<3$. Here we present some experimental results.

Table 7: Average empirical probability for the case of $s=2$ and $\lambda=10^{5}$ $n$ | $16$ | $20$ | $24$ | $28$ | $32$ | $36$
---|---|---|---|---|---|---
$k=0$ | $0.8599$ | $0.8543$ | $0.8575$ | $0.8604$ | $0.8628$ | $0.8643$
$k=1$ | $0.8654$ | $0.8618$ | $0.8580$ | $0.8671$ | $0.8646$ | $0.8620$
$k=\frac{n}{2}$ | $0.8609$ | $0.8611$ | $0.8651$ | $0.8621$ | $0.8682$ | $0.8662$
Table 8: Average empirical probability for the case of $s=1$ and $\lambda=10^{5}$ $n$ | $16$ | $20$ | $24$ | $28$ | $32$ | $36$
---|---|---|---|---|---|---
$k=0$ | $0.7201$ | $0.7177$ | $0.7124$ | $0.7227$ | $0.7125$ | $0.7110$
$k=1$ | $0.7103$ | $0.7141$ | $0.7154$ | $0.7129$ | $0.7118$ | $0.7192$
$k=\frac{n}{2}$ | $0.7168$ | $0.7212$ | $0.7210$ | $0.7226$ | $0.7106$ | $0.7154$
Table 9: Average empirical probability for the case of $s=0$ and $\lambda=10^{5}$ $n$ | $16$ | $20$ | $24$ | $28$ | $32$ | $36$
---|---|---|---|---|---|---
$k=0$ | $0.4365$ | $0.4363$ | $0.4353$ | $0.4435$ | $0.4434$ | $0.4377$
$k=1$ | $0.4323$ | $0.4337$ | $0.4345$ | $0.4451$ | $0.4336$ | $0.4440$
$k=\frac{n}{2}$ | $0.4385$ | $0.4371$ | $0.4290$ | $0.4330$ | $0.4427$ | $0.4330$
Table 10: Average empirical probability for the case of $k=1$ and $\lambda=10^{5}$ with fixed $(1\mbox{,\ }1\mbox{,\ }\ldots\mbox{,\ }1)\in\mathbb{Z}^{n}$ as the $k\times n$ matrix to be completed. $n$ | $16$ | $20$ | $24$ | $28$ | $32$ | $36$
---|---|---|---|---|---|---
$s=0$ | $0.4330$ | $0.4354$ | $0.4342$ | $0.4381$ | $0.4284$ | $0.4407$
$s=1$ | $0.7163$ | $0.7120$ | $0.7209$ | $0.7198$ | $0.7218$ | $0.7159$
$s=2$ | $0.8629$ | $0.8641$ | $0.8649$ | $0.8673$ | $0.8616$ | $0.8568$

All test examples in Table 7–9 are generated with the same method described in
the last subsection, i.e., the input matrices are randomly chosen. In Table
10, all test examples are fixed to a $1\times n$ primitive row $(1\mbox{,\
}1\mbox{,\ }\ldots\mbox{,\ }1)$, but still with $\lambda=10^{5}$. All of these
tables show that for $0\leq s\leq 2$, the empirical probability of that the
resulting matrix is primitive is relatively high, similar with that of the
case $s\geq 3$ shown in the last section. However, an effective theoretical
lower bound on the probability for this case is left open.

#### Acknowledgments.

The authors would like to thank two anonymous referees for helpful comments
and suggestions that makes the presentation of this paper clearer.

## References

  * [AR16] Divesh Aggarwal and Oded Regev.  A note on discrete gaussian combinations of lattice vectors.  Chicago Journal of Theoretical Computer Science, 2016(Article 07):1–11, 2016.  https://doi.org/10.4086/cjtcs.2016.007.
  * [Bab85] László Babai.  On Lovász’ lattice reduction and the nearest lattice point problem.  In K. Mehlhorn, editor, Proc STACS ’85 (January 3–5, 1985, Saarbrücken, Germany), volume 182 of Lecture Notes in Computer Science, pages 13–20. Springer, Heidelberg, 1985.  https://doi.org/10.1007/BFb0023990.
  * [CS05] Zhuliang Chen and Arne Storjohann.  Lattice compression of integer matrices.  https://cs.uwaterloo.ca/~astorjoh/latcomp.ps, 2005.
  * [DS17] M. Graça Duffner and Fernando C. Silva.  On the existence of unimodular matrices with a prescribed submatrix.  Linear Algebra and its Applications, 515:321–330, 2017.  https://doi.org/10.1016/j.laa.2016.11.015.
  * [EGV00] Wayne Eberly, Mark Giesbrecht, and Gilles Villard.  Computing the determinant and Smith form of an integer matrix.  In Proceedings of the 41st Annual Symposium on Foundations of Computer Science (12–14 November 2000, Redondo Beach, USA), pages 675–685. IEEE Computer Society, Los Alamitos, 2000.  https://doi.org/10.1109/SFCS.2000.892335.
  * [Fan07] Maozhong Fang.  On the completion of a partial integral matrix to a unimodular matrix.  Linear Algebra and its Applications, 422(1):291–294, 2007.  https://doi.org/10.1016/j.laa.2006.10.012.
  * [FW14] Felix Fontein and Pawel Wocjan.  On the probability of generating a lattice.  Journal of Symbolic Computation, 64:3–15, 2014.  https://doi.org/10.1016/j.jsc.2013.12.002.
  * [GHL16] Xiangqian Guo, Fengdan Hou, and Xuewen Liu.  Natural density of integral matrices that can be extended to invertible integral matrices.  Linear and Multilinear Algebra, 64(9):1878–1886, 2016.  https://doi.org/10.1080/03081087.2015.1127316.
  * [KBS13] Rachel K. Kalaimani, Madhu N. Belur, and Sivaramakrishnan Sivasubramanian.  Generic pole assignability, structurally constrained controllers and unimodular completion.  Linear Algebra and its Applications, 439(12):4003–4022, 2013.  https://doi.org/10.1016/j.laa.2013.10.004.
  * [KNSW20] Elena Kirshanova, Huyen Nguyen, Damien Stehlé, and Alexandre Wallet.  On the smoothing parameter and last minimum of random orthogonal lattices.  Designs, Codes and Cryptography, 88(5):931–950, 2020.  https://doi.org/10.1007/s10623-020-00719-w.
  * [LLL82] Arjen K. Lenstra, Hendrik W. Lenstra, and László Lovász.  Factoring polynomials with rational coefficients.  Mathematische Annalen, 261(4):515–534, 1982.  https://doi.org/10.1007/BF01457454.
  * [MRW11] Gérard Maze, Joachim Rosenthal, and Urs Wagner.  Natural density of rectangular unimodular integer matrices.  Linear Algebra and its Applications, 434(5):1319–1324, 2011.  https://doi.org/10.1016/j.laa.2010.11.015.
  * [MS04] Thom Mulders and Arne Storjohann.  Certified dense linear system solving.  Journal of Symbolic Computation, 37(4):485–510, 2004.  https://doi.org/10.1016/j.jsc.2003.07.004.
  * [PL02] See-May Phoong and Yuan-Pei Lin.  Application of unimodular matrices to signal compression.  In Proceedings of 2002 IEEE International Symposium on Circuits and Systems (May 26–29, 2002, Phoenix-Scottsdale, USA), pages 837–840. IEEE, Piscataway, 2002.  https://doi.org/10.1109/ISCAS.2002.1009971.
  * [Ran91] Dana Randall.  Efficient generation of random nonsingular matrices.  Technical Report UCB/CSD-91-658, EECS Department, University of California, Berkeley, Nov 1991.  http://www2.eecs.berkeley.edu/Pubs/TechRpts/1991/6147.html.
  * [Rei56] Irving Reiner.  Unimodular complements.  The American Mathematical Monthly, 63(4):246–247, 1956.  https://doi.org/10.2307/2310351.
  * [Sch87] Claus-Peter Schnorr.  A hierarchy of polynomial time lattice basis reduction algorithms.  Theoretical Computer Science, 53(2–3):201–224, 1987.  https://doi.org/10.1016/0304-3975(87)90064-8.
  * [Sto03] Arne Storjohann.  High-order lifting and integrality certification.  Journal of Symbolic Computation, 36(3):613–648, 2003.  https://doi.org/10.1016/S0747-7171(03)00097-X.
  * [Sto05] Arne Storjohann.  The shifted number system for fast linear algebra on integer matrices.  Journal of Complexity, 21(4):609–650, 2005.  https://doi.org/10.1016/j.jco.2005.04.002.
  * [Zha06] Xingzhi Zhan.  Completion of a partial integral matrix to a unimodular matrix.  Linear Algebra and its Applications, 414(1):373–377, 2006.  https://doi.org/10.1016/j.laa.2005.10.013.
  * [ZL14] Wei Zhou and George Labahn.  Unimodular completion of polynomial matrices.  In Katsusuke Nabeshima, editor, Proceedings of the 39th International Symposium on Symbolic and Algebraic Computation (July 23–25, 2014, Kobe, Japan), pages 413–420. ACM, New York, 2014.  https://doi.org/10.1145/2608628.2608640.


# Beyond the Local Void: A comprehensive view on the origins of the Amaterasu
particle

Nadine Bourriche Max Planck Institute for Physics
Boltzmannstraße 8, 85748 Garching, Germany  Technical University of Munich
James-Franck-Straße 1, 85748 Garching  Francesca Capel Max Planck Institute
for Physics
Boltzmannstraße 8, 85748 Garching, Germany

###### Abstract

We use the reconstructed properties of the Amaterasu particle, the second-
highest energy cosmic ray ever detected, to map out three-dimensional
constraints on the location of its unknown source. We highlight possible
astrophysical sources that are compatible with these regions and requirements.
Among these, M82, a powerful starburst galaxy, stands out as a strong
candidate due to its position and proximity. To derive our constraints, we use
CRPropa 3 to model all relevant propagation effects, including deflections in
the Galactic and extra-Galactic magnetic fields. We consider key input
quantities such as source distance, position, energy, and the strength and
coherence length of the extra-Galactic magnetic field as free parameters. We
then infer constraints on these parameters by applying approximate Bayesian
computation. We present our results, demonstrating the impact of different
assumptions for the arrival mass of the Amaterasu particle and the systematic
uncertainties on the energy scale.

Particle astrophysics(96) — Cosmic rays (329) — Ultra-high-energy cosmic
radiation(1733) — Bayesian statistics(1900)

††software: astropy (Astropy Collaboration et al., 2013, 2018), numpy (Harris
et al., 2020), matplotlib (Hunter, 2007), scipy (Virtanen et al., 2020), h5py
(Collette, 2013; Collette et al., 2023), seaborn (Waskom, 2021), cartopy (Met
Office, 2024), CRPropa 3 (Batista et al., 2022)

##  1 Introduction

Ultra-high energy cosmic rays (UHECRs) are charged particles with energies
that exceed $E\geq 10^{18}$ eV. The origin of these particles is still unknown
despite numerous efforts to identify their sources. These searches are
challenged by the complex nature of UHECR propagation, including energy losses
and deflections by intervening magnetic fields. Although advancements in the
search for sources of UHECRs have been made, including the possible
association of UHECRs with starburst galaxies (SBGs) and active galactic
nuclei (AGN) catalogs, these models cannot explain the UHECR data at the
highest energies (Aab et al., 2018; Abbasi et al., 2018; Capel & Mortlock,
2019; Abdul Halim et al., 2024). While few in number, events at the highest
energies have the potential to be the most constraining in terms of their
origins. Their high energies suggest nearby horizons and trajectories are less
affected by magnetic fields, even when considering the evidence for the
increasing mass and charge of UHECRs as a function of energy (Batista et al.,
2019). As such, studying individual events at the highest energies is a
complementary and powerful way to search for the sources of UHECRs (Globus et
al., 2023; Bourriche & Capel, 2023).

In this work, we focus on the recent detection of the Amaterasu particle by
the Telescope Array Collaboration (TA, Abbasi et al. 2023). Amaterasu was
detected by the surface detector of the Telescope Array with an energy of
$\mathrm{E=244\pm 29(stat.)^{+51}_{-76}(syst.)}$ EeV and an arrival direction
of ${\mathrm{R.A.},\mathrm{Dec.}=(255.9\pm 0.6,16.1\pm 0.5)^{\circ}}$, making
it the second-highest energy particle ever detected and the highest energy
particle ever detected in an actively operating observatory. The high energy
scale makes Amaterasu a strong candidate for an individual event-based UHECR
source search. Previous works have investigated the possible source of
Amaterasu by studying its compatibility with models for UHECR production in
nearby galaxies (Kuznetsov, 2023) and by estimating its deflection and horizon
through backtracking and 1D simulations (Abbasi et al., 2023; Unger & Farrar,
2024a). The results suggest that Amaterasu’s detected direction does not
correspond to any known galaxy but seems to come from the Local Void, an
especially low-density region of the Universe (Tully et al., 2008). This
conclusion has led to the proposal of past astrophysical transient sources
(Farrar, 2024), ultraheavy cosmic rays (Zhang et al., 2024), magnetic
monopoles (Frampton & Kephart, 2024), Lorentz invariance violation (Lang,
2024), and superheavy dark matter (Sarmah et al., 2024) as possible
explanations for the Amaterasu observation.

We revisit this puzzle by leveraging Approximate Bayesian Computation (ABC,
Beaumont 2019) to enable the use of 3D simulations with CRPropa 3 (Batista et
al., 2022) and the inclusion of a comprehensive description of UHECR
propagation in our inference. We also include the statistical and systematic
uncertainties in the reconstructed properties of Amaterasu, as well as key
modelling uncertainties regarding the Galactic and extra-Galactic magnetic
fields and possible energy of Amaterasu at its source. The result of applying
ABC in this context is a 3D posterior distribution of possible source
locations that are consistent with the observed energy and direction of
Amaterasu, marginalising over the various sources of uncertainty. This paper
is organised as follows: in Section 2, we describe the CRPropa 3 simulation
setup, the free parameters considered, and our choice of priors, in Section 3
we detail our implementation of ABC to exploit the results of the simulations,
in Section 4 we present our results and in Section 5 we discuss their
implications in the framework of possible astrophysical sources.

##  2 Physical model

To model the propagation of UHECRs, we use CRPropa 3 to perform 3D
simulations. We include all relevant particle interactions: photo-pion
production, photo-disintegration, electron-pair production, as well as nuclear
decay and adiabatic losses. The extra-Galactic magnetic field (EGMF) is
modelled as a Gaussian random field with a Kolmogorov turbulence spectrum on a
regular 3D grid with a spacing of 25 kpc. The EGMF is parameterised by the
root mean square field strength, $B_{\mathrm{rms}}$, and the coherence length,
$L_{\mathrm{c}}$. For the Galactic magnetic field (GMF), we use the JF12
model111We note that at the time of this work, the model for the coherent GMF
described in Unger & Farrar (2024b) is not yet publicly available in the suite
of CRPropa 3 lens models. (Jansson & Farrar, 2012) as implemented in CRPropa
3, considering both the regular and turbulent components. In addition to the
free parameters $B_{\mathrm{rms}}$ and $L_{\mathrm{c}}$, we consider different
realisations of the turbulent fields for both the GMF and EGMF at each
iteration of our ABC algorithm (described further in Section 3) to account for
their stochastic nature in our final results.

We assume that Amaterasu leaves the source as an iron nucleus. This assumption
is in line with the findings of previous work and is motivated by the fact
that heavier source compositions are more easily accelerated to high energies.
Additionally, the observed UHECR spectrum at the highest energies is well-
described by models assuming heavy source compositions (see, e.g., Abdul Halim
et al. 2024; Unger & Farrar 2024a and references therein). Finally, the lack
of deflections expected from a lighter particle further complicates an
astrophysical interpretation (Abbasi et al., 2023; Kuznetsov, 2023). While we
fix the particle type at the source, we leave the energy at the source as a
free parameter, $E_{\mathrm{src}}$, to explore its impact on our results. In
principle, our framework can treat the particle type at the source as a free
parameter, but we choose to fix it here due to the weak constraints expected
from the data and for the clarity of the resulting interpretation.

To map out the possible volume of space consistent with the measured energy
and arrival direction of Amaterasu, we consider the source location as a set
of free parameters: the galactic longitude and latitude, $(l,b)$, and the
distance, $D_{\mathrm{src}}$. This gives us a total of six free parameters:
$l$, $b$, $D_{\mathrm{src}}$, $E_{\mathrm{src}}$, $B_{\mathrm{rms}}$ and
$L_{c}$. As we perform inference within a Bayesian framework, we define priors
for these free parameters. Our prior choices are made to be constraining
enough to include important physical information but wide enough to avoid
driving the resulting inference where possible.

The EGMF remains relatively poorly understood (Durrer & Neronov, 2013).
Outside of galaxy clusters, we expect $\left\langle
B_{\mathrm{rms}}^{2}||L_{\mathrm{c}}\right\rangle^{1/2}\lesssim
10^{-8}\,\text{G}\,\text{Mpc}^{1/2}$, where $L_{\mathrm{c}}\lesssim 1$ Mpc
(Kotera & Olinto, 2011). To reflect the large uncertainties while satisfying
this constraint, we choose uniform priors over the EGMF parameters such that
$B_{\mathrm{rms}}=[0,10]$ nG and $L_{\mathrm{c}}=[50,1000]$ kpc. The lower
bounds of these priors are chosen such that the EGMF deflections in this range
are expected to be negligible compared to the other included modelling
uncertainties.

For the energies of Amaterasu at the source, we choose a power law spectrum
prior of $E^{-\Gamma}$, where $\Gamma=1$. This is the best-fit spectral index
obtained by the combined composition and spectrum fit of the Pierre Auger data
(Aab et al., 2017), and also the same assumption as made in (Unger & Farrar,
2024a). While we could let $\Gamma$ be a free parameter, we do not expect to
be able to constrain the spectral shape with a single event and so choose to
fix it here. The minimum energy at the source is set to be 3$\sigma$ below the
reconstructed energy of Amaterasu, $E$, where $\sigma$ is the statistical
uncertainty of the energy as measured by TA. This choice of minimum energy is
driven by the fact that arrival energies more than $3\sigma$ lower than the
detected energy will be rejected by our ABC algorithm. The maximum energy is
chosen to be 3 ZeV to explore the impact of larger post-GZK horizon (Greisen,
1966) energies on our results.

The source location prior is defined in 2 steps. The prior for
$D_{\mathrm{src}}$ is set such that source positions are uniformly sampled
within a spherical volume. We also only consider distances in the range
$[1,D_{\mathrm{max}}]$ Mpc, with the choice of $D_{\mathrm{max}}$ detailed in
Section 3. The priors of $l$ and $b$ of the source are determined by taking
into consideration both the GMF and EGMF. We estimate the impact of the GMF by
backtracking (Thielheim & Langhoff, 2002). We model the galaxy as a sphere of
$20$ kpc, we place the Earth $8.5$ kpc from the center. Starting from the
detected arrival direction of Amaterasu, we backtrack multiple realisations of
an iron nucleus to the Galactic boundary with energies sampled from a Gaussian
distribution centred on 244 EeV with a standard deviation corresponding to the
reported statistical uncertainty of 29 EeV. By assuming Amaterasu arrives as
an iron nucleus at Earth, we consider the largest possible deflections. We
then estimate the average deflection angle caused by the EGMF using equation
2.6 from Harari et al. (2002). We use the arrival direction at the Galactic
boundary, along with the average deflections expected from the EGMF, to
construct the distribution from which we then sample prior source positions.
As the EGMF deflections depend on the other free parameters
$D_{\mathrm{src}}$, $E_{\mathrm{src}}$, $B_{\mathrm{rms}}$ and
$L_{\mathrm{c}}$, the position prior is slightly different at every iteration
of the ABC algorithm. We verify that the overall prior resulting from running
our algorithm covers much of the sky, and the resulting distributions shown in
Section 4 are more constraining in all cases.

##  3 Methods

In this work, we take a simulation-based approach to inference, using ABC to
estimate the posterior distribution of our free parameters. The advantage of
using ABC in this context, as opposed to traditional methods, is that we do
not need to define the likelihood function that connects the observed data to
the model parameters explicitly; instead, our choice of simulation defines it
implicitly, as defined in Section 2. ABC proceeds by sampling parameter values
from the prior distributions, running a simulation based on these parameters
to generate a simulated data set, and then comparing this simulated data to
the observed data. Proposed parameter values are then accepted or rejected
based on how closely the simulated data match the observations, thereby
approximating the posterior distribution.

We use the reported arrival direction and reconstructed energy of Amaterasu as
data. We define a distance metric for comparison with simulated data as
$\mathrm{|X_{sim}-X_{obs}|}$, where $\mathrm{X_{sim}}$ is the simulated data
and $\mathrm{X_{obs}}$ is the observed one, where $X$ is the galactic
latitude, longitude, or the energy. To compare the two, a tolerance
$\mathrm{\epsilon}=3\sigma$ is introduced, with $\sigma$ denoting the
reconstruction uncertainty on the data. A set of parameters is accepted if the
resulting simulation produces an event with its arrival energy and direction
within $\mathrm{3\sigma}$ of those observed for Amaterasu. Each accepted
parameter set is weighted based on how well the resulting events match the
observed data. The weights for each observable are defined as
$\mathrm{\omega_{X}=\frac{1}{\sigma_{X}\sqrt{2\pi}}\exp\left(-0.5\left(\frac{X_{obs}-X_{sim}}{\sigma_{X}}\right)^{2}\right)}$,
and then combined them into a total weight
$\mathrm{\omega=\omega_{l}\cdot\omega_{b}\cdot\omega_{E}}$.

For each proposed combination of parameters, we simulate as described in
Section 2, with $\mathrm{10^{6}}$ particles emitted isotropically from a
single source location. This number gives a good compromise concerning
computational efficiency and results in reasonable numbers of detected events.
In the rare cases that we accept more than a single event from one simulation
iteration, one is randomly chosen to define the resulting weights.

As discussed in Abbasi et al. (2023), the reconstructed energy reported for
the Amaterasu particle is subject to large systematic uncertainties. To take
these uncertainties into account in our analysis, we consider 2 different
cases for the detected energy: 1) the nominal case with
${E_{\mathrm{nom}}=244\pm 29}$ EeV, and 2) the lower end of the systematic
range with ${E_{\mathrm{low}}=168\pm 29}$ EeV. In both cases, the statistical
uncertainties are treated equivalently, and we do not consider the higher end
of the reported systematic range as these results will only be more
constraining than those in case 1). For case 1), we set the corresponding
${D_{\mathrm{max}}=12}$ Mpc, and for case 2) we use ${D_{\mathrm{max}}=15}$
Mpc. These values are based on the expected maximum horizon for an iron
nucleus at the relevant energy scale (see, e.g., Bourriche & Capel 2023).

##  4 Results

We apply our ABC approach to the nominal and low energy cases, obtaining a
total of 512 and 519 parameter sets, respectively222For all the proposed
parameter sets resulting from our prior, only a fraction of $\sim 10^{-3}$
result in simulations with accepted events. making the implementation of our
ABC algorithm computationally challenging.. We summarise the resulting
posterior distributions for $l$, $b$, $D_{\mathrm{src}}$, $E_{\mathrm{src}}$,
$B_{\mathrm{rms}}$ and $L_{\mathrm{c}}$ below. In the visualisation of our
results, we use Gaussian kernel density estimation, making sure to account for
hard parameter boundaries where relevant. We also verify that the resulting
distributions are robust to reasonable variations in the bandwidth.

Figure 1: Sky maps resulting from the nominal energy run (case 1) showing the
possible source positions of Amaterasu in galactic coordinates. The magenta
star marks the measured arrival direction of Amaterasu, and the dashed line
outlines the Local Void. The circular markers show galaxies within the
accepted $D_{\mathrm{src}}$ range, with larger markers indicating SBGs and AGN
and smaller markers showing quiescent galaxies. In Fig.1(a), the green
contours outline the 10%, 30%, 70%, and 90% regions of highest posterior
density, independent of the composition at arrival. Fig 1(b) shows the
composition-dependent posterior distribution conditional on the particle
arriving with a mass number, $A$, in one of three groups. The orange contours
show the case for light elements with $A<4$, pink for $4\leq A<28$, and blue
for $A\geq 28$. The contour levels are as in Fig. 1(a).

Fig. 1 shows the sky maps in Galactic coordinates for case 1) with
${E_{\mathrm{nom}}=244\pm 29}$ EeV. In Fig. 1(a), the total posterior
distribution is compared to known astrophysical sources, with the SBG source
list from Aab et al. (2018) and AGN from Baumgartner et al. (2013) and Ajello
et al. (2017) shown. We also include quiescent galaxies from the 2MASS survey
(Huchra et al., 2012). All objects are colour-coded according to their
distance from Earth. Fig. 1(b) shows the same plot for the case of the
posterior distribution assuming different arrival particle compositions. We
find that most of the accepted simulations result in events that are Helium or
fragments of it, but heavier arrival compositions are also possible. We divide
the accepted parameters based on the accepted events type into three groups
based on their mass number, $A$: light ($A<4$), medium ($4\leq A<28$), and
heavy ($A\geq 28$). These groups are chosen to contain a comparable number of
accepted parameter sets. We clearly see that UHECRs that are heavier on
arrival suffer larger deflections because of their stronger charge, shifting
the center of the posterior further from the detected arrival direction of
Amaterasu.

The corresponding sky maps for case 2) with ${E_{\mathrm{low}}=168\pm 29}$ EeV
are shown in Fig. 2. The resulting posterior distributions for the source
direction follow a similar trend to those shown in Fig. 1, but are more
extended since lower energies lead to larger $D_{\mathrm{src}}$ and thus
larger magnetic deflections.

Figure 2: Sky maps of the low energy run (case 2). The layout is as in Fig. 1,
with Fig. 2(a) showing the total posterior distribution and Fig. 2(b) showing
the composition-dependent distributions.

Fig. 3(c) and 3(d) show the marginal posterior distributions of
$D_{\mathrm{src}}$ for the case of $E_{\mathrm{nom}}$ and $E_{\mathrm{low}}$,
respectively. The distributions are shown as stacked histograms with the
contribution from assuming different arrival mass groups highlighted in
different colours, as in Figs. 1(b) and 2(b). For $E_{\mathrm{nom}}$, the
source distance is relatively constrained, with the most probable distance
around ${D_{\mathrm{src}}\sim 2}$ Mpc. We also see that a medium composition
dominates at closer distances while a heavier composition is preferred for
${D_{\mathrm{src}}>6}$ Mpc. The case of $E_{\mathrm{low}}$ is more complex,
with a most probable distance of ${D_{\mathrm{src}}\sim 4}$ Mpc, but a longer
tail out to over 14 Mpc and a mixed arrival composition across the distance
range.

Figure 3: The marginal posterior distribution of $D_{\mathrm{src}}$ is shown
for the nominal and low energy cases in panels (a) and (b), respectively. The
results are shown as a stacked histogram to highlight the relative
contributions of $D_{\mathrm{src}}$ values that lead to accepted events in
different arrival mass groups. Orange, pink, and blue bars indicate arrival
mass numbers of ${A<4}$, ${4\leq A<28}$, and ${A\geq 28}$, respectively.

We summarise the results for the remaining free parameters in Fig. 4. To do
so, we combine the EGMF parameters to reflect their impact on UHECR
deflections and convert $E_{\mathrm{src}}$ into the source rigidity,
${R=E_{\mathrm{src}}/eZ}$, where $eZ$ is the particle charge. The preferred
rigidity is $\mathrm{R\approx 10^{20.2}}$ V and $\mathrm{R\approx 10^{18.8}}$
V for the nominal energy and $\mathrm{R\approx 10^{20}}$ V for the low energy
case. The results on $\mathrm{B_{rms}\sqrt{Lc}}$ show a preference for smaller
values, but no values are excluded.

Figure 4: The joint marginal posterior distribution for the source rigidity,
$R$, and the EGMF parameters, $B_{\mathrm{rms}}\sqrt{L_{\mathrm{c}}}$ for the
nominal (a) and low energy (b) runs. The contours show the 10%, 30%, 70%, and
90% regions of highest posterior density.

##  5 Discussion

Our combination of 3D CRPropa 3 simulations with ABC allows us to include a
more realistic physical model for UHECR propagation and directly constrain the
parameters of this model and their correlations.

When considering all possible arrival masses and $E_{\mathrm{nom}}$, as shown
in figure 1(a), there is no significant overlap of the $(l,b)$ posterior
distribution with known objects and the preferred direction is largely
contained within the Local Void, consistent with the findings of Abbasi et al.
(2023). However, if we consider the case of heavier arrival compositions, as
shown in Fig. 1(b), our contours overlap with three astrophysical objects. Two
quiet galaxies are within the 70% contour, of which one, NGC 6789, is an
extremely close $\mathrm{D\leq 4}$ Mpc AGN candidate. AGN can potentially
accelerate UHECRs in relativistic jet flows or via internal trans relativistic
shocks (Rieger, 2022). However, NGC 6789 does not show signs of intensive past
activity like Cen A. The third source, NGC 6946, is an SBG that lies within
our 90% contour. NGC 6946 has also been found close to the region of possible
source positions in Unger & Farrar (2024a). However, it was disfavoured as a
convincing source candidate as it contributes only 3% to the summed 1.4 GHz
radio flux of SBGs within horizons similar to ours.

The lower energy case reveals many more objects consistent with our posterior
distribution. In this case, shown in figure 2(a), our 70% contour overlaps
with multiple astrophysical objects that are within ${D_{\mathrm{src}}\leq
14}$ Mpc. However, all of these are galaxies that appear relatively inactive
and do not belong to a source class that satisfies the Hillas criterion
(Hillas, 1984). The 90% contour contains three known SBGs: M81, M82, and NGC
6946. M82 is a powerful and nearby SBG at a distance of only 3.6 Mpc. It lies
a few degrees from the TA hotspot and is commonly invoked as a UHECR source
candidate (Abbasi et al., 2014; He et al., 2016). M82 is also the brightest
SBG within our horizon, accounting for over 18% of the total 1.4 GHz radio
flux (Aab et al., 2018), where the 1.4 GHz radio flux is often used as a proxy
of UHECR production (Aab et al., 2018; Abdul Halim et al., 2024).

The mass-dependent posterior contours for the low energy case continue to
reveal further astrophysical source candidates, particularly for the heavier
arrival compositions, as shown in Fig. 2(b). At the center of the
distribution, there are several quieter galaxies. Additionally, we see three
SGBs within the 30% region: M81, M82, and NGC 6946. Several more SGBs and AGN,
such as NGC 2403, NGC 6946, NGC 1569, IC142, and Maffei 2, are contained
within the 70% and 90% contours.

Our results offer alternatives to the origins of Amaterasu outside of the
Local Void when considering key modelling uncertainties, the possibility of a
heavier arrival mass, and/or the lower end of the systematic energy
uncertainty. While there are currently no strong constraints on the arrival
mass composition, it seems that the UHECR composition evolves with energy and
we can expect to see heavier elements at the highest energies (Mayotte et al.,
2023). Considering its characteristics and its proximity to the center of our
distribution, M82 appears to be a strong astrophysical candidate.

There are several possible phenomena that can drive UHECR acceleration in
SBGs, such as superwinds caused by the additive action of stellar winds and
supernovae (Anchordoqui, 2019). Another explanation for the association
between UHECRs and SBGs that also includes the action of AGN is offered by the
echo model (Bell & Matthews, 2022).

Alternatively, Amaterasu could have been produced by a past transient event.
The average propagation time of our accepted events for the nominal and the
low energy case throughout Galactic and extra-Galactic magnetic fields is
${\approx 48}$ Myrs and ${\approx 209}$ Myrs, respectively. Evidence from the
distribution of gas around the M81/M82 systems shows that the two galaxies
experienced a tidal encounter within the last $\mathrm{\approx 200}$ Myrs (Yun
et al., 1994). Other transient events that could be responsible for UHECR
acceleration are tidal disruption events, gamma-ray bursts, or binary neutron
star mergers (Biehl et al., 2018; Batista et al., 2019; Farrar, 2024).

##  6 Conclusions

Through a comprehensive approach that couples a full 3D UHECR propagation
model with simulation-based inference, we are able to constrain a volume of
possible source locations consistent with the observed Amaterasu particle. We
assume that the Amaterasu particle is an iron nucleus at the source and
marginalise over uncertainties in its source energy and the modelling of
Galactic and extra-Galactic magnetic fields. Considering the systematic
uncertainty on the reconstructed energy and the possibility of a heavier
arrival mass, $A\geq 28$, our results are consistent with regions of space
beyond the Local Void. Several astrophysical sources lie within the proposed
volume, including the starburst galaxy M82, a nearby SBG only a few Mpc away
from us that is also the brightest SBG in the 1.4 GHz radio flux within our
horizon. Alternatively, Amaterasu could have been produced by a transient
event in one of the many quiescent galaxies contained within our volume.

Our results demonstrate the importance of resolving the UHECR energy and
arrival mass or charge at the highest energies. We see a significant impact on
the possible volumes consistent with the Amaterasu particle for the different
cases considered here, and a more constrained measurement of the nominal
energy reported for Amaterasu or evidence for a lighter arrival composition at
these energies would leave little room for origins outside of the Local Void.
Planned upgrades to existing observatories such as TAx4 (Fujisue, 2023) and
AugerPrime (Castellina, 2019) will soon augment our knowledge, and we expect
further advances from next-generation facilities such as POEMMA and GRAND
(Coleman et al., 2023).

## acknowledgments

We thank M. Unger, T. K. Bister, A. Fedynitch, and A. van Vliet for their
valuable input during our discussions, as well as the Max Planck Computing and
Data Facility for the use of the Raven HPC system. N. Bourriche acknowledges
the financial support from the Excellence Cluster ORIGINS, which is funded by
the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under
Germany’s Excellence Strategy - EXC-2094-390783311.

## References

  * Aab et al. (2017) Aab, A., Abreu, P., Aglietta, M., et al. 2017, Journal of Cosmology and Astroparticle Physics, 2017, 038, doi: 10.1088/1475-7516/2017/04/038
  * Aab et al. (2018) —. 2018, The Astrophysical Journal Letters, 853, L29, doi: 10.3847/2041-8213/aaa66d
  * Abbasi et al. (2014) Abbasi, R. U., Abe, M., Abu-Zayyad, T., et al. 2014, The Astrophysical Journal Letters, 790, L21, doi: 10.1088/2041-8205/790/2/l21
  * Abbasi et al. (2018) —. 2018, The Astrophysical Journal Letters, 867, L27, doi: 10.3847/2041-8213/aaebf9
  * Abbasi et al. (2023) Abbasi, R. U., Allen, M. G., Arimura, R., et al. 2023, Science, 382, 903, doi: 10.1126/science.abo5095
  * Abdul Halim et al. (2024) Abdul Halim, A., Abreu, P., Aglietta, M., et al. 2024, Journal of Cosmology and Astroparticle Physics, 2024, 022, doi: 10.1088/1475-7516/2024/01/022
  * Ajello et al. (2017) Ajello, M., Atwood, W. B., Baldini, L., et al. 2017, The Astrophysical Journal Supplement Series, 232, 18, doi: 10.3847/1538-4365/aa8221
  * Anchordoqui (2019) Anchordoqui, L. A. 2019, Physics Reports, 801, 1, doi: 10.1016/j.physrep.2019.01.002
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Batista et al. (2019) Batista, R. A., Biteau, J., Bustamante, M., et al. 2019, Frontiers in Astronomy and Space Sciences, 6, 23, doi: 10.3389/fspas.2019.00023
  * Batista et al. (2022) Batista, R. A., Tjus, J. B., Dörner, J., et al. 2022, Journal of Cosmology and Astroparticle Physics, 2022, 035, doi: 10.1088/1475-7516/2022/09/035
  * Baumgartner et al. (2013) Baumgartner, W. H., Tueller, J., Markwardt, C. B., et al. 2013, The Astrophysical Journal Supplement Series, 207, 19, doi: 10.1088/0067-0049/207/2/19
  * Beaumont (2019) Beaumont, M. A. 2019, Annual Review of Statistics and Its Application, 6, 379, doi: https://doi.org/10.1146/annurev-statistics-030718-105212
  * Bell & Matthews (2022) Bell, A. R., & Matthews, J. H. 2022, Monthly Notices of the Royal Astronomical Society, 511, 448, doi: 10.1093/mnras/stac031
  * Biehl et al. (2018) Biehl, D., Boncioli, D., Lunardini, C., & Winter, W. 2018, Scientific Reports, 8, 10828, doi: 10.1038/s41598-018-29022-4
  * Bourriche & Capel (2023) Bourriche, N., & Capel, F. 2023, Proceedings of 38th International Cosmic Ray Conference — PoS(ICRC2023), 362, doi: 10.22323/1.444.0362
  * Capel & Mortlock (2019) Capel, F., & Mortlock, D. J. 2019, Monthly Notices of the Royal Astronomical Society, 484, 2324, doi: 10.1093/mnras/stz081
  * Castellina (2019) Castellina, A. 2019, EPJ Web of Conferences, 210, 06002, doi: 10.1051/epjconf/201921006002
  * Coleman et al. (2023) Coleman, A., Eser, J., Mayotte, E., et al. 2023, Astroparticle Physics, 149, 102819, doi: 10.1016/j.astropartphys.2023.102819
  * Collette (2013) Collette, A. 2013, Python and HDF5 (O’Reilly)
  * Collette et al. (2023) Collette, A., Kluyver, T., Caswell, T. A., et al. 2023, h5py/h5py: 3.8.0, 3.8.0, Zenodo, doi: 10.5281/zenodo.7560547
  * Durrer & Neronov (2013) Durrer, R., & Neronov, A. 2013, The Astronomy and Astrophysics Review, 21, 62, doi: 10.1007/s00159-013-0062-7
  * Farrar (2024) Farrar, G. R. 2024.  https://arxiv.org/abs/2405.12004
  * Frampton & Kephart (2024) Frampton, P. H., & Kephart, T. W. 2024, The Amaterasu Cosmic Ray as a Magnetic Monopole and Implications for Extensions of the Standard Model.  https://arxiv.org/abs/2403.12322
  * Fujisue (2023) Fujisue, K. 2023, Proceedings of 38th International Cosmic Ray Conference — PoS(ICRC2023), 308, doi: 10.22323/1.444.0308
  * Globus et al. (2023) Globus, N., Fedynitch, A., & Blandford, R. D. 2023, The Astrophysical Journal, 945, 12, doi: 10.3847/1538-4357/acaf5f
  * Greisen (1966) Greisen, K. 1966, Physical Review Letters, 16, 748, doi: 10.1103/physrevlett.16.748
  * Harari et al. (2002) Harari, D., Mollerach, S., Roulet, E., & Sánchez, F. 2002, Journal of High Energy Physics, 2002, 045, doi: 10.1088/1126-6708/2002/03/045
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * He et al. (2016) He, H.-N., Kusenko, A., Nagataki, S., et al. 2016, Phys. Rev. D, 93, 043011, doi: 10.1103/PhysRevD.93.043011
  * Hillas (1984) Hillas, A. M. 1984, Annual Review of Astronomy and Astrophysics, 22, 425, doi: 10.1146/annurev.aa.22.090184.002233
  * Huchra et al. (2012) Huchra, J. P., Macri, L. M., Masters, K. L., et al. 2012, The Astrophysical Journal Supplement Series, 199, 26, doi: 10.1088/0067-0049/199/2/26
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Jansson & Farrar (2012) Jansson, R., & Farrar, G. R. 2012, The Astrophysical Journal, 757, 14, doi: 10.1088/0004-637x/757/1/14
  * Kotera & Olinto (2011) Kotera, K., & Olinto, A. V. 2011, Annual Review of Astronomy and Astrophysics, 49, 119, doi: 10.1146/annurev-astro-081710-102620
  * Kuznetsov (2023) Kuznetsov, M. Y. 2023, arXiv, doi: 10.48550/arxiv.2311.14628
  * Lang (2024) Lang, R. G. 2024, New physics as a possible explanation for the Amaterasu particle.  https://arxiv.org/abs/2405.03528
  * Mayotte et al. (2023) Mayotte, E. W., Collaboration, T. P. A., Halim, A. A., et al. 2023, Proceedings of 38th International Cosmic Ray Conference — PoS(ICRC2023), 365, doi: 10.22323/1.444.0365
  * Met Office (2024) Met Office. 2024, Cartopy: a cartographic python library with a Matplotlib interface, Exeter, Devon.  https://scitools.org.uk/cartopy
  * Rieger (2022) Rieger, F. M. 2022, Universe, 8, 607, doi: 10.3390/universe8110607
  * Sarmah et al. (2024) Sarmah, P., Das, N., Borah, D., Chakraborty, S., & Mehta, P. 2024, The Amaterasu particle: constraining the superheavy dark matter origin of UHECRs.  https://arxiv.org/abs/2406.03174
  * Thielheim & Langhoff (2002) Thielheim, K. O., & Langhoff, W. 2002, Journal of Physics A: General Physics, 1, 694, doi: 10.1088/0305-4470/1/6/308
  * Tully et al. (2008) Tully, R. B., Shaya, E. J., Karachentsev, I. D., et al. 2008, The Astrophysical Journal, 676, 184, doi: 10.1086/527428
  * Unger & Farrar (2024a) Unger, M., & Farrar, G. R. 2024a, The Astrophysical Journal Letters, 962, L5, doi: 10.3847/2041-8213/ad1ced
  * Unger & Farrar (2024b) —. 2024b, The Coherent Magnetic Field of the Milky Way.  https://arxiv.org/abs/2311.12120
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Waskom (2021) Waskom, M. L. 2021, Journal of Open Source Software, 6, 3021, doi: 10.21105/joss.03021
  * Yun et al. (1994) Yun, M. S., Ho, P. T. P., & Lo, K. Y. 1994, Nature, 372, 530, doi: 10.1038/372530a0
  * Zhang et al. (2024) Zhang, B. T., Murase, K., Ekanger, N., Bhattacharya, M., & Horiuchi, S. 2024, Ultraheavy Ultrahigh-Energy Cosmic Rays.  https://arxiv.org/abs/2405.17409