# Can ordered precursors promote the nucleation of solid solutions?

Caroline Desgranges  Department of Chemistry, New York University, New York,
New York 10003, United States  Department of Chemistry, University of North
Dakota, Grand Forks North Dakota 58202, United States  Jerome Delhommelle
Department of Chemistry, New York University, New York, New York 10003, United
States  Department of Chemistry, University of North Dakota, Grand Forks North
Dakota 58202, United States

###### Abstract

Crystallization often proceeds through successive stages that lead to a
gradual increase in organization. Using molecular simulation, we determine the
nucleation pathway for solid solutions of copper and gold. We identify a new
nucleation mechanism (liquid$\to$$L1_{2}$ precursor$\to$solid solution),
involving a chemically ordered intermediate that is more organized than the
end product. This nucleation pathway arises from the low formation energy of
$L1_{2}$ clusters which, in turn, promote crystal nucleation. We also show
that this mechanism is composition-dependent since the high formation energy
of other ordered phases precludes them from acting as precursors.

The emergence of order from disorder has long fascinated scientists.Bernstein
and Bernstein (2002); Sosso et al. (2016) For instance, crystal nucleation, i.
e. formation of a crystal nucleus from a disordered liquid, has attracted
considerable attention since the end of the $19^{th}$ century and remains to
this day an extremely challenging phenomenon. As acknowledged in Ostwald’s
famous step rule,Ostwald (1897) crystallization can be a multi-step process,
resulting in a gradual onset of order until the stable crystal polymorph is
obtained. There has been tremendous progress in our understanding of this
process in recent years. For instance, density functional theory and computer
simulations have shed light on the gradual progress of the disordered liquid
phase towards the stable crystal phase. In such systems, density fluctuations
in the liquid give rise to a crystal nucleus of a metastable form, whose
features (density and order) are intermediate between the liquid and the
stable polymorph, followed by the formation of the stable polymorph.Oxtoby
(1998); ten Wolde et al. (1995); Desgranges and Delhommelle (2007a) Examples
of this progressive onset of order include the pre-organization of liquid
pockets, with a structural order intermediate between liquid and
crystal,Lechner et al. (2011); Desgranges and Delhommelle (2011) before the
formation of a crystal nucleus, as well as the occurrence of phase separation
and demixing prior to crystallization, as observed for polymer mixturesTanaka
and Nishi (1985) and metal alloys.Desgranges and Delhommelle (2014) Here we
focus on the crystal nucleation of solid solutions of metals and identify a
new nucleation pathway that does not involve the expected gradual increase of
order during crystal nucleation. On the contrary, we find that the nucleation
pathway can involve a more ordered intermediate structure, since it exhibits
both structural and chemical order, than the solid solution, which only has
structural order.

The formation and properties of solid solutions have attracted considerable
attention in recent years.Zhang et al. (2017) For instance, solid solutions of
semiconductors, such as $GaN$ and $ZnO$, are remarkably stable and efficient
photocatalysts during water splitting.Maeda et al. (2005) Similarly, solid
solutions of $MoS_{2}$ and $P$ are highly efficient and cost-effective
materials for hydrogen production.Ye et al. (2016) Solid solutions of metals
also play a key role in the development of multi-principal elements alloys,
including high entropy alloys with a superior mechanical performance.Senkov et
al. (2015); Wu et al. (2016a) At the other end of the spectrum, metal alloys
can also form highly ordered structures with both structural and chemical
orders as in the $L1_{2}$ and $L1_{0}$ phases, with unique magnetic properties
as, for instance, for $FePt$ nanoparticles.Sun et al. (2000) Understanding the
interplay between structural and chemical order is thus key for a wide range
of applications.

We focus here on the nucleation process for different compositions of the
copper-gold alloy. These systems are characterized by a negative enthalpy of
mixing and exhibit solid solutions and ordered phases over a wide composition
range, with an order-disorder transition taking place beyond a specific
temperature, as shown in the pioneering work of Shockley Shockley (1938).
Specifically, we examine this process for a mole fraction of $x_{Cu}=0.25$
($CuAu_{3}$ alloy), $x_{Cu}=0.5$ ($CuAu$ alloy) and $x_{Cu}=0.75$ ($Cu_{3}Au$
alloy). $CuAu$ has an ordered $L1_{0}$ phase, stable below $683$K, and
undergoes an order-disorder transition into the disordered close-packed (CP)
structure above this temperature. Similarly, $CuAu_{3}$ and $Cu_{3}Au$ have an
ordered $L1_{2}$ phase and the disordered CP structure becomes stable above
$663$K and $500$K, respectively.Ozoliņš et al. (1998)

We use molecular simulation to study the nucleation process in $CuAu$ alloys.
To model the interactions between metal atoms, we employ the many-body quantum
corrected Sutton-Chen embedded atoms model (qSC-EAM)Luo et al. (2003). The
qSC-EAM force field performs very well when compared to the experiment for the
lattice parameters, cohesive energies, surface energies as well as melting
points,Kart et al. (2005); Sup  which are key parameters for the
crystallization process. Simulations of the nucleation process are carried out
in the isothermal-isobaric ensemble within a hybrid Monte Carlo-Molecular
Dynamics framework as in previous work on bimetallic systemsDesgranges and
Delhommelle (2018, 2014, 2016) at $1$ bar and a temperature 45% below the
melting point obtained for each composition of the alloy.Kart et al. (2005)
The temperature is set to $754$K for $CuAu_{3}$, $748$K for $CuAu$ and $737$K
for $Cu_{3}Au$, i.e. above the order-disorder transition temperatures for all
systems. During crystal nucleation, the system has to overcome a large free
energy barrier, that corresponds to the formation of a nucleus of a critical
size.ten Wolde et al. (1995) Non-Boltzmann sampling methods thus need to be
employed to sample the entire nucleation pathway and gain access to the free
energy barrier of nucleation. Here we use the umbrella sampling (US)
technique,Torrie and Valleau (1974) in which the sampling of configurations
with low Boltzmann weights is bolstered by a harmonic function of the order
parameterSteinhardt et al. (1983) $Q_{6}$. The $Q_{6}$ order parameter
measures the rate of crystallinity as shown in prior experimentalGasser et al.
(2001) and simulation studies.ten Wolde et al. (1995); Desgranges and
Delhommelle (2007a) This approach has been shown to be very efficient for
nucleation in single-component systemsten Wolde et al. (1995); Auer and
Frenkel (2001); Desgranges and Delhommelle (2007b); Beckham et al. (2008), as
well as in binary systems and bimetallic alloys.Valeriani et al. (2005);
Desgranges and Delhommelle (2018, 2014) We add that Frenkel, Dijkstra and co-
workers have carried out thorough comparisons of several simulations methods
for crystal nucleation.Valeriani et al. (2005); Filion et al. (2010)
Specifically, Dijkstra and co-workers performed US simulations to determine
the free energy barrier, nucleation rate, size and structure of the crystal
nucleus. They showed that the US results were in excellent agreement with
those obtained using other approaches, such as forward-flux sampling and
molecular dynamics simulations.Filion et al. (2010) Similar conclusions were
drawn by Frenkel and co-workers on crystal nucleation in molten sodium
chloride.Valeriani et al. (2005)

We start by presenting the results for the free energy barrier of nucleation.
Fig. 1(a) shows the free energy profiles as a function of the order parameter.
The profiles are qualitatively the same for the three mole fractions. First,
the free energy profile starts from a minimum, corresponding to the metastable
liquid phase, which is disordered and associated with a low value of the order
parameter. Then, as the order parameter increases and the crystal nucleus
forms, the free energy increases until it reaches a maximum, corresponding to
the formation of a critical nucleus. There is, however, a striking difference
between the results obtained for $x_{Cu}=0.5$ and the other two mole
fractions. The height of the free energy barrier for $CuAu$ is $21.5\pm
2~{}k_{B}T$, almost twice that for the other two systems. Moreover, the free
energy barriers for $CuAu_{3}$ and $Cu_{3}Au$ are roughly the same ($12.1\pm
1~{}k_{B}T$ and $11.8\pm 1~{}k_{B}T$, respectively). This suggests that the
nucleation mechanism in $CuAu$ is very different from the other two systems,
and that nucleation proceeds in a similar fashion for $Cu_{3}Au$ and
$CuAu_{3}$.

(a) (b)

Figure 1: Crystal nucleation in $CuAu$ alloys. (a) Free energy barrier of
nucleation. (b) (top) Size and structural (inset plot) features for the
critical nucleus and (bottom) fraction of ordered structure in the critical
nucleus ($L1_{0}$ for $x_{Cu}=0.5$ or $L1_{2}$ for $x_{Cu}=0.25$ and
$x_{Cu}=0.75$).

To identify the reason for these different behaviors, we examine the features
of the system at the end of the nucleation step. We carry out a detailed
analysis of the size and structure of the critical nucleus as in prior
work.Auer and Frenkel (2001); Desgranges and Delhommelle (2018, 2014) The size
of the critical nucleus is plotted, for each composition, in Fig. 1(b). This
graph shows that the size of the critical nucleus increases steadily with the
mole fraction in $Cu$. This increase in size is correlated with the steady
decrease in the free energy difference between the solid and liquid
$\Delta\mu$ that drives the nucleation process (numerical data are provided in
the supplemental information). This is, however, a monotonic trend that cannot
account for the unexpectedly high free energy for $x_{Cu}=0.5$, or for the
similar results obtained for $x_{Cu}=0.25$ and $x_{Cu}=0.75$. We also show the
structure of the critical nucleus in Fig. 1(b) as a function of $x_{Cu}$, with
a steady decrease in the fraction of atoms with a close-packed-like (CP-like)
environment with $x_{Cu}$. Again, we find a monotonic trend as a function of
$x_{Cu}$ that does not provide a rationale for the free energy results. We
finally recall that crystal nucleation is carried out here above the order-
disorder transition for all compositions. As shown at the bottom of Fig. 1(b),
the fraction of ordered structures, be it $L1_{2}$ or $L1_{0}$, is very small
(less than $10$%) in the critical nucleus for all $x_{Cu}$. However, we do
observe a non-monotonic dependence of this fraction as a function of $x_{Cu}$,
which is a first clue that suggests that ordered structures play a role in the
behavior observed for the free energy.

A closer inspection of the free energy profiles of Fig. 1(a) reveals that,
even for low values of $Q_{6}$, the slope of the free energy profile is much
steeper for $x_{Cu}=0.5$. This suggests a different nucleation mechanism in
the early stages of the process. It has been known since the pioneering work
of FranckFranck (1952) that locally favored structures of low energy, such as
icosahedra, form in the supercooled liquid and can impact nucleation. Kawasaki
and TanakaKawasaki and Tanaka (2010) have recently confirmed that the
supercooled liquid is spanned by transient structures, including medium-range
crystal-like domains, that can promote crystal nucleation. Can the presence of
different types of low energy precursors in the supercooled liquids account
for the different nucleation mechanisms?

(a) (b)

Figure 2: Supercooled liquids for the copper-gold alloy. (a) $Cu-Cu$ pair
correlation functions. (b) Number density distributions of clusters formed by
connected ISROs via volume sharing, with (inset plot) ISRO formation energy as
a function of the node degree $k$

We start with the characterization of icosahedral order in the supercooled
liquids. Icosahedral order is often found in supercooled metallic
liquids,Liang et al. (2014) most notably when there is a significant size
difference as here between Au and Cu with $r_{Au}:r_{Cu}=1.13$. The radial
distribution functions in Fig. 2(a) reveal the onset of splitting for the
second peak for the three supercooled liquids. Such a behavior has been linked
to the existence of icosahedral order.Liang et al. (2014) Since icosahedral
order is not structurally compatible with crystalline order, the presence of
icosahedral clusters (icosahedral short-range order or ISRO) and of connected
icosahedra (medium-range order or MRO) is an inhibitor of crystal
nucleationDesgranges and Delhommelle (2018) and often leads to the onset of
glassy behavior.Leocmach and Tanaka (2012) We therefore test the following
hypothesis: is the increase in height of the free energy barrier for
$x_{Cu}=0.5$ caused by icosahedral order? We determine the number density of
clusters of connected icosahedra as a function of the node degree $k$
following the analysis carried out in prior work.Wu et al. (2016b); Desgranges
and Delhommelle (2018, 2007c) We also calculate the formation energy of these
clusters as a function of $k$, choosing as the reference for the energy of
each atom the value for crystals of pure $Cu$ and $Au$ in line with Wu et
al.Wu et al. (2016b) Results are shown in Fig. 2(b) and confirm the presence
of icoasahedral MRO. We also observe that the $x_{Cu}=0.5$ supercooled liquid
exhibits the same qualitative features as for the other two mole fractions.
The number density plot for connected icosahedra are very similar for
$x_{Cu}=0.25$ and $x_{Cu}=0.5$, while the number density is slightly greater
for $x_{Cu}=0.75$ as a result of the lower formation energy especially at low
$k$ values. However, the formation energy varies very slowly as a function of
$k$ for all three systems. For instance, for $x_{Cu}=0.5$, it decreases from
$0.12eV$ ($k=0$) to $0.08eV$ ($k=5$) and never reaches a level that can
actually impact nucleation. Previous work on supercooled liquids of $CuZr$ and
$Ag_{6}Cu_{4}$ have reported 2- to 4-fold decreases in the formation energy Wu
et al. (2016b) over the same range of $k$, and values of the order of $0.05eV$
for $k=5$ for the development of MRO structures large enough to impact the
height of the free energy barrier of nucleation. Desgranges and Delhommelle
(2018) We thus conclude that the presence of icosahedral SRO and MRO does not
account for the different nucleation mechanisms observed here.

Figure 3: Formation energy (left y-axis) of ordered $L1_{2}$ and $L1_{0}$
clusters, together with the number of ordered clusters (right y-axis) per unit
volume and unit time.

Given the ability of copper-gold alloys to form ordered structures, we now
look for the possible presence of clusters with a $L1_{0}$ or a $L1_{2}$
structure in the supercooled liquids using a map based on Steinhardt order
parameters $(\hat{q_{4}},\hat{w_{4}})$ (details are provided in the
supplemental material). The formation energies are given in Fig. 3. Very
interestingly, the results show that the formation energies for $L1_{2}$
clusters are much lower than for the $L1_{0}$ cluster. This provides a first
insight into the different behavior that takes place in the supercooled liquid
for $x_{Cu}=0.25$ and $x_{Cu}=0.75$, which can give rise to $L1_{2}$ clusters,
and for $x_{Cu}=0.5$, in which $L1_{0}$ is the ordered structure. Furthermore,
the formation energy for $L1_{2}$ clusters is below $0.1eV$ for both
$x_{Cu}=0.25$ and $x_{Cu}=0.75$, which is lower than the formation energy of
ISRO ($k=0$) and of the connected icosahedra MRO with the highest number
density ($k\leq 3$). This means that the formation of $L1_{2}$ clusters of
$CuAu_{3}$ and $Cu_{3}Au$ is energetically favored. On the other hand, in the
$x_{Cu}=0.5$ case, the formation energy for $L1_{0}$ clusters is close to
$0.11eV$, comparable to that of connected icosahedra for $k\geq 1$. Thus, the
formation of such $L1_{0}$ clusters is much less frequent than of $L1_{2}$
clusters (see Fig. 3).

(a) (b)

Figure 4: $L1_{2}$ ordered clusters of $Cu_{3}Au$ in the supercooled liquid
(a) and in a configuration containing a critical nucleus (b). For clarity,
only the central $Cu$ atom (ochre) and the four closest $Au$ (tan) are shown
as large spheres and liquid-like atoms are translucent. In (a), atoms in
icosahedra are shown as small yellow spheres. In (b), atoms from the critical
nucleus are shown as small blue spheres ($Au$) and small green spheres ($Cu$).

This provides the following picture. For $x_{Cu}=0.25$ and $x_{Cu}=0.75$,
locally favored structures, that are both structurally and chemically ordered,
easily form in the supercooled liquid. An example of a $L1_{2}$ cluster of
$Cu_{3}Au$ is shown on the snapshot of Fig. 4(a). These clusters serve as
precursors for the crystal nucleus, unlike icosahedral clusters in
supercoooled liquid metals. The latter have a structure incompatible with the
crystal, resulting in higher free energy barriers and in the remarkable
stability of supercooled liquid metals. Here, the structural compatibility
between $L1_{2}$ and the stable phase (the structurally ordered and chemically
disordered solid solution) facilitates the nucleation process, resulting in
the decrease in the free energy barrier observed in Fig. 1(a) for both
$CuAu_{3}$ and $Cu_{3}Au$. This synergy continues throughout the nucleation
step: the center of the nucleus converts into the stable phase, while $L1_{2}$
clusters continue to form on its surface (see snapshot of a critical nucleus
of $Cu_{3}Au$ in Fig. 4(b)) and, in turn, assist crystallization. Such a
mechanism is reminiscent of the role played by the body-centered cubic (BCC)
phase during the nucleation of the stable face-centered cubic (FCC) phase from
supercooled Lennard-Jones (LJ) liquids.Oxtoby (1998); ten Wolde et al. (1995)
Indeed, in the LJ case, nucleation starts with the formation of small BCC
clusters that serve as precursors. Then, the center of the LJ nucleus converts
into the FCC structure, while BCC clusters continue to form on its surface to
promote crystallization. On the other hand, the higher formation energy for
$L1_{0}$ clusters implies that the nucleation of the $CuAu$ is not promoted by
ordered precursors and, as a result, the free energy barrier is much higher.

We finally compare the nucleation rates for the three systems. As discussed by
Song et al., the attachment rate ($f^{+}_{n_{c}}$) can, in some cases, be the
dominant factor in the nucleation rate.Song et al. (2018) Thus, following
Frenkel and co-workers,Auer and Frenkel (2001, 2004); Valeriani et al. (2005);
Filion et al. (2010), we determine $f^{+}_{n_{c}}$ for each system (see
supplemental information for details). We find that $f^{+}_{n_{c}}$ are
comparable for the three systems, varying from $2.6\times 10^{-14}~{}s^{-1}$
($x_{Cu}=0.25$) to $2.7\times 10^{-14}~{}s^{-1}$ ($x_{Cu}=0.5$), and finally
to $4.5\times 10^{-14}~{}s^{-1}$ ($x_{Cu}=0.75$), and that the three systems
exhibit kinetic prefactors of the same order. This means that the free energy
barrier is the dominant factor when comparing the nucleation rates of the
three alloys. Thus, given its greater free energy barrier, the nucleation rate
for $CuAu$ is of the order of $10^{32}~{}m^{-3}.s^{-1}$ and is about four
order of magnitude smaller than those obtained both for $CuAu_{3}$and
$Cu_{3}Au$.

The analysis of the nucleation pathway and of its free energy profile allows
us to identify a new nucleation mechanism, in which ordered precursors
facilitate the formation of a critical nucleus of the solid solution. This
result is linked to the presence in the parent phase of small ordered $L1_{2}$
clusters, that are energetically favored over competing low-energy structures
like icosahedra. The $L1_{2}$ structure continues to promote nucleation later
on as clusters keep forming at the solid-liquid interface during the
development of the nucleus, resulting in a lower free energy barrier of
nucleation. On the other hand, the $L1_{0}$ structure does not promote the
nucleation process, leading us to connect the high formation energy of ordered
$L1_{0}$ clusters to a higher free energy barrier of nucleation. As captured
by Ostwald’s step rule,Ostwald (1897) crystallization generally proceeds
through stages that result in a gradual increase in order. Here, however, the
mechanism observed involves an intermediate that is actually more organized
than the end result, since it is both chemically and structurally ordered as
summarized below

$\begin{array}[]{ccccc}\textrm{liquid}&\to&L1_{2}~{}\textrm{precursor}&\to&\textrm{solid
solution}\\\ \textrm{(structural disorder)}&&\textrm{(structural
order)}&&\textrm{(structural order)}\\\ \textrm{(chemical
disorder)}&&\textrm{(chemical order)}&&\textrm{(chemical disorder)}\\\
\end{array}$

This novel pathway suggests new ways to understand and control the formation
of solid solutions, as these have emerged recently as promising candidates for
many applications including water splitting,Maeda et al. (2005) hydrogen
productionYe et al. (2016) and multi-principal elements and high entropy
alloys.Senkov et al. (2015); Wu et al. (2016a)

Acknowledgements Partial funding for this research was provided by NSF through
CAREER award DMR-1052808.

## References

  * Bernstein and Bernstein (2002) J. Bernstein and J. M. Bernstein, _Polymorphism in molecular crystals_ , vol. 14 (Oxford University Press, 2002).
  * Sosso et al. (2016) G. C. Sosso, J. Chen, S. J. Cox, M. Fitzner, P. Pedevilla, A. Zen, and A. Michaelides, Chem. Rev. 116, 7078 (2016).
  * Ostwald (1897) W. Ostwald, Z. Phys. Chem 22, 289 (1897).
  * Oxtoby (1998) D. W. Oxtoby, Acc. Chem. Res. 31, 91 (1998).
  * ten Wolde et al. (1995) P. R. ten Wolde, M. J. Ruiz-Montero, and D. Frenkel, Phys. Rev. Lett. 75, 2714 (1995).
  * Desgranges and Delhommelle (2007a) C. Desgranges and J. Delhommelle, Phys. Rev. Lett. 98, 235502 (2007a).
  * Lechner et al. (2011) W. Lechner, C. Dellago, and P. G. Bolhuis, Phys. Rev. Lett. 106, 085701 (2011).
  * Desgranges and Delhommelle (2011) C. Desgranges and J. Delhommelle, J. Am. Chem. Soc. 133, 2872 (2011).
  * Tanaka and Nishi (1985) H. Tanaka and T. Nishi, Phys. Rev. Lett. 55, 1102 (1985).
  * Desgranges and Delhommelle (2014) C. Desgranges and J. Delhommelle, J. Am. Chem. Soc. 136, 8145 (2014).
  * Zhang et al. (2017) F. Zhang, S. Zhao, K. Jin, H. Xue, G. Velisa, H. Bei, R. Huang, J. Ko, D. Pagan, J. Neuefeind, et al., Phys. Rev. Lett. 118, 205501 (2017).
  * Maeda et al. (2005) K. Maeda, T. Takata, M. Hara, N. Saito, Y. Inoue, H. Kobayashi, and K. Domen, J. Am. Chem. Soc. 127, 8286 (2005).
  * Ye et al. (2016) R. Ye, P. del Angel-Vicente, Y. Liu, M. J. Arellano-Jimenez, Z. Peng, T. Wang, Y. Li, B. I. Yakobson, S.-H. Wei, M. J. Yacaman, et al., Adv. Mater. 28, 1427 (2016).
  * Senkov et al. (2015) O. Senkov, J. Miller, D. Miracle, and C. Woodward, Nat. Commun. 6, 6529 (2015).
  * Wu et al. (2016a) Z. Wu, Y. Gao, and H. Bei, Acta Mater. 120, 108 (2016a).
  * Sun et al. (2000) S. Sun, C. B. Murray, D. Weller, L. Folks, and A. Moser, Science 287, 1989 (2000).
  * Shockley (1938) W. Shockley, J. Chem. Phys. 6, 130 (1938).
  * Ozoliņš et al. (1998) V. Ozoliņš, C. Wolverton, and A. Zunger, Phys. Rev. B 58, R5897 (1998).
  * Luo et al. (2003) S. N. Luo, T. J. Ahrens, T. Cagin, A. Strachan, W. A. Goddard, and D. C. Swift, Phys. Rev. B 68, 134206 (2003).
  * Kart et al. (2005) H. Kart, M. Tomak, and T. Çağin, Modelling Simul. in Mater. Sci. Eng. 13, 657 (2005).
  * (21) See Supplemental Information for a comparison with the experimental data, which includes Refs. [22-24].
  * Kart et al. (2004) S. O. Kart, M. Tomak, M. Uludogan, and T. Cagin, J. Non-Cryst. Solids 337, 101 (2004).
  * Hultgren et al. (1973) R. Hultgren, P. D. Desai, D. T. Hawkins, M. Gleiser, and K. K. Kelley, Tech. Rep., National Standard Reference Data System (1973).
  * Okamoto et al. (1987) H. Okamoto, D. Chakrabarti, D. Laughlin, and T. Massalski, J. Phase Equil. 8, 454 (1987).
  * Desgranges and Delhommelle (2018) C. Desgranges and J. Delhommelle, Phys. Rev. Lett. 120, 115701 (2018).
  * Desgranges and Delhommelle (2016) C. Desgranges and J. Delhommelle, J. Phys. Chem. C 120, 27657 (2016).
  * Torrie and Valleau (1974) G. M. Torrie and J. P. Valleau, Chem. Phys. Lett. 28, 578 (1974).
  * Steinhardt et al. (1983) P. J. Steinhardt, D. R. Nelson, and M. Ronchetti, Phys. Rev. B 28, 784 (1983).
  * Gasser et al. (2001) U. Gasser, E. R. Weeks, A. Schofield, P. N. Pusey, and D. A. Weitz, Science 292, 258 (2001).
  * Auer and Frenkel (2001) S. Auer and D. Frenkel, Nature 409, 1023 (2001).
  * Desgranges and Delhommelle (2007b) C. Desgranges and J. Delhommelle, J. Am. Chem. Soc. 129, 7012 (2007b).
  * Beckham et al. (2008) G. T. Beckham, B. Peters, and B. L. Trout, J. Phys. Chem. B 112, 7460 (2008).
  * Valeriani et al. (2005) C. Valeriani, E. Sanz, and D. Frenkel, J. Chem. Phys. 122, 194501 (2005).
  * Filion et al. (2010) L. Filion, M. Hermes, R. Ni, and M. Dijkstra, J. Chem. Phys. 133, 244115 (2010).
  * Franck (1952) F. C. Franck, Proc. R. Soc. London A 215, 43 (1952).
  * Kawasaki and Tanaka (2010) T. Kawasaki and H. Tanaka, Proc. Natl. Acad. Sci. 107, 14036 (2010).
  * Liang et al. (2014) Y.-C. Liang, R.-S. Liu, Y.-F. Mo, H.-R. Liu, Z.-A. Tian, Q.-Y. Zhou, H.-T. Zhang, L.-L. Zhou, Z.-Y. Hou, and P. Peng, J. Alloy. Comp. 597, 269 (2014).
  * Leocmach and Tanaka (2012) M. Leocmach and H. Tanaka, Nature Commun. 3, 974 (2012).
  * Wu et al. (2016b) Z. Wu, F. Li, C. Huo, M. Li, W. Wang, and K. Liu, Sci. Rep. 6, 35967 (2016b).
  * Desgranges and Delhommelle (2007c) C. Desgranges and J. Delhommelle, J. Chem. Phys. 127, 144509 (2007c).
  * Song et al. (2018) H. Song, Y. Sun, F. Zhang, C.-Z. Wang, K.-M. Ho, and M. I. Mendelev, Phys. Rev. Mater. 2, 023401 (2018).
  * Auer and Frenkel (2004) S. Auer and D. Frenkel, Annu. Rev. Phys. Chem. 55, 333 (2004).


# High-energy Neutrino Source Cross-correlations with Nearest Neighbor
Distributions

Zhuoyang Zhou  Current affiliation: Department of Statistics & Data Science,
Carnegie Mellon University, Pittsburgh, PA, USA  Department of Statistics,
University of Wisconsin-Madison, Madison, WI, USA  Jessi Cisewski-Kehe
Department of Statistics, University of Wisconsin-Madison, Madison, WI, USA
Ke Fang  Department of Physics, Wisconsin IceCube Particle Astrophysics
Center, University of Wisconsin-Madison, Madison, WI, USA  Arka Banerjee
Department of Physics, Indian Institute of Science Education and Research,
Pashan, Pune, India

###### Abstract

The astrophysical origins of the majority of the IceCube neutrinos remain
unknown. Effectively characterizing the spatial distribution of the neutrino
samples and associating the events with astrophysical source catalogs can be
challenging given the large atmospheric neutrino background and underlying
non-Gaussian spatial features in the neutrino and source samples. In this
paper, we investigate a framework for identifying and statistically evaluating
the cross-correlations between IceCube data and an astrophysical source
catalog based on the $k$-Nearest Neighbor Cumulative Distribution Functions
($k$NN-CDFs). We propose a maximum likelihood estimation procedure for
inferring the true proportions of astrophysical neutrinos in the point-source
data. We conduct a statistical power analysis of an associated likelihood
ratio test with estimations of its sensitivity and discovery potential with
synthetic neutrino data samples and a WISE-2MASS galaxy sample. We apply the
method to IceCube’s public ten-year point-source data and find no
statistically significant evidence for spatial cross-correlations with the
selected galaxy sample. We discuss possible extensions to the current method
and explore the method’s potential to identify the cross-correlation signals
in data sets with different sample sizes.

high-energy neutrinos, large-scale structure

††software: The code for the $k$NN-CDFs, joint $k$NN-CDFs, and $k$NN-CDFs
Cross-correlation computation and synthetic event generation is publicly
available at https://github.com/Zhuoyang666/IceCubekNN. Jupyter Notebook files
and code for reproducing some of the key figures are in the paper also
provided. All code is written in Python.

revtex4-1Repair the float revtex4-1Assuming for authors

##  1 Introduction

The existence of extraterrestrial astrophysical neutrino emission has been
confirmed by the IceCube Neutrino Observatory (Aartsen et al., 2013). More
recently, evidence of high-energy neutrino emission from individual sources
has been established (Aartsen et al., 2018, 2020a; Abbasi et al., 2022). The
astrophysical origins of the majority of detected neutrino events, however,
are still largely unknown.

The most commonly used neutrino source search method involves computing the
probability of each neutrino coming from each candidate in a catalog or each
pixel in a source template (Braun et al., 2008). This method has successfully
led to the detection of neutrinos from NGC 1068 (Abbasi et al., 2022) and the
Galactic Plane (Abbasi et al., 2023). It has also been applied to study the
association of neutrinos with nearby galaxies (Aartsen et al., 2020b). Another
method is to search for the cross-correlation between the spatial
distributions of neutrino and source samples (Fang et al., 2020; Urban et al.,
2021; Creque-Sarbinowski et al., 2022; Negro et al., 2023; Guevel & Fang,
2023; Ouellette & Holder, 2024). As the cross-correlation method only relies
on the assumption that the neutrino sources trace the same underlying field as
the chosen catalog, it has the advantage of not requiring a perfect catalog
that contains the actual neutrino sources.

So far the cross-correlation studies are based on the two-point correlation
function (2PCF). Such a correlation analysis could be limiting since 2PCF
cannot give a comprehensive characterization of non-Gaussian spatial
clustering that is present in the distribution of matter at low redshifts and
on small spatial separation scales ($\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower
4.0pt\hbox{$\sim$}}\hss}\hbox{$<$}}}40\,h^{-1}{\rm Mpc}$). One possible
solution is to consider higher $N$-point correlation functions. However, the
associated computational cost will be much higher as the computational
complexity increases exponentially with $N$, which makes it impractical to
measure on large datasets (Banerjee & Abel, 2021a). Given the large
atmospheric background in the full IceCube point source data, a more sensitive
measure of spatial clustering with satisfying computational cost becomes
necessary.

To this end, we adopt the $k$-Nearest Neighbor Cumulative Distribution
Functions (Banerjee & Abel, 2021a), $k$NN-CDFs, in our cross-correlation
study. For a group of data points that fill a volume and a specific $k$, the
$k$NN-CDFs are simply defined as the empirical cumulative distribution
function of the associated $k$-nearest neighbor distances, which can be
estimated through interpolation on the sorted $k$NN distances. In a three-
dimensional space with an Euclidean metric, the value of the $k$NN-CDFs at a
scale $r$ can be interpreted as the estimated probability of finding at least
$k$ data points that are enclosed in a sphere with radius $r$. Based on the
$k$NN-CDFs, a measure of the similarities of spatial clustering of two data
samples called joint $k$NN-CDFs and an associated $k$NN-CDFs cross-correlation
is available (Banerjee & Abel, 2021b).

As mathematically shown in Banerjee & Abel (2021a), the $k$NN-CDFs are related
to the $N$-point correlation functions measured on the same group of data
points but have a reduced computational burden through the use of $k$-d trees
that can be constructed with $N\log N$ operations and a $k$NN searching in
$\log N$ time for each point. The $k$NN-CDFs have demonstrated more
significant detection of spatial clustering in sparse datasets such as the
1000 richest redMaPPer clusters in the SDSS DR8 redMaPPer clusters than the
two-point correlation function (Wang et al., 2022).

In this work, we adapted the $k$NN-CDF method to accommodate the properties of
IceCube point source data, where only the equatorial coordinates are available
for each point source without the depth (redshift) information. We developed
the framework to compute the $k$NN-CDFs cross-correlations and a test of
significance of such correlation between IceCube 10-year point source data
(IceCube Collaboration et al., 2021; IceCube Collaboration, 2021) and a galaxy
catalog based on the Wide-Field Infrared Survey Explorer (WISE; Wright et al.,
2010) and the 2-Micron All-Sky Survey (2MASS; Skrutskie et al., 2006) infrared
databases. By estimating the sensitivity and discovery potential on the
synthetic neutrino datasets generated based on the public ten-year data, we
have shown that the methodology is capable of identifying a highly significant
cross-correlation when the neutrino events and the galaxy sample both trace
the same underlying large-scale density field.

The layout of the paper is as follows: we introduce the backgrounds of $k$NN-
CDFs and related summary statistics in Section 2. The $k$NN-CDFs-related
computations are introduced in Section 3. In Section 4, we summarize the data
engineering for neutrino and galaxy samples that are used in the analysis.
Then we present the main findings based on synthetic samples and the IceCube
public ten-year data in Section 5 and Section 6, respectively. A discussion
and concluding remarks are provided in Section 7.

##  2 Background

The goal of the $k$NN-CDFs is to offer a more sensitive and comprehensive
characterization of the density field from which the data samples (discrete
tracers) are drawn. Thus, a preliminary assumption is that the observed data
points are indeed the discrete tracers of their underlying density field, and
the degrees of spatial clustering of the data points directly reflect their
underlying density fluctuations.

###  2.1 $k$NN-CDFs

We start by defining the $k$NN-CDFs and leave the detailed computation in the
methodology Section 3. Assume there is a group of data points distributed over
a three-dimensional volume $V$, for a specific $k$, the value of $k$NN-CDF at
distance $r$, denoted as $\text{CDF}_{k}(r)$, can be defined as the
probability of finding at least $k$ data points at a distance less than $r$
for a random point in the same volume. With the usual Euclidean metric, we
have the following definition:

$\text{CDF}_{k}(r)=P_{>k-1|V}|_{V=\frac{4}{3}\pi r^{3}},$ (1)

which can be interpreted as the probability of finding spheres of radius $r$
that enclose at least $k$ data points. This gives the $k$NN-CDFs an intuitive
interpretation as measures of the degree of spatial clustering at different
scales with these “probing spheres” on the data sample.

Notice that $k$NN-CDFs can also be defined on two-dimensional surfaces and
other Riemannian manifolds with a corresponding metric, such as the surface of
a unit sphere with a metric defined as the great circle distance. This is
crucial for the application in the context of IceCube point source data since
only the equatorial coordinates are available for each point source without
the depth information or redshift. Thus, a projection onto a spherical surface
is needed.

###  2.2 Joint $k$NN-CDFs

After defining the $k$NN-CDFs on a single data sample, it can be extended to a
joint function called the joint $k$NN-CDFs based on two data samples. Similar
to the definition of $k$NN-CDFs, for specific $k$ (e.g., $k_{1}=k_{2}=k$), the
value of the joint $k$NN-CDF at $r$, denoted as $\text{CDF}_{k_{1},k_{2}}(r)$,
can be defined as the probability of finding a sphere with radius $r$ that
contain at least $k$ data points in both data samples:

$\text{CDF}_{k_{1},k_{2}}(r)=P(>k_{1}-1,>k_{2}-1|V).$ (2)

As shown in Banerjee & Abel (2021b), there is a relationship between the joint
$k$NN-CDFs and the connected $N$-point correlation functions such that the
combined joint $k$NN-CDFs contain the same information as the $N$-point
correlation functions that can be measured on the two density fields from
which the two data samples are drawn. In other words, the joint $k$NN-CDFs
effectively measure the correlation of the two data samples in terms of
spatial clustering.

Before discussing the $k$NN-CDFs cross-correlation, we provide a visual
demonstration in Figure 1 for a more intuitive understanding of the
computation of the $k$NN-CDFs and the joint $k$NN-CDFs within a rectangular
region using the usual Euclidean metric (this can be other surfaces or volumes
with an appropriate metric in practice). The black pluses are reference points
generated from a homogeneous spatial Poisson process. The orange and blue
points, respectively, represent the source catalog (IceCube neutrinos in our
case) and tracers (such as the WISE-2MASS samples). In practice, many more
reference points (black pluses) are generated. We only consider two for
illustration purposes. The reference point at the bottom of Figure 1
illustrates the 3NN-CDF computation. Its radius $r_{3}$ represents one 3NN
distance to the blue points. The 3NN-CDF of the blue points are estimated with
the empirical 3NN-CDF based on many such $r_{3}$ from different reference
points. The reference point at the center illustrates the computation of the
joint $3$NN-CDF. This reference point identifies two 3NN distances, $r_{1}$
and $r_{2}$, to its nearby orange and blue points, respectively. We adopt the
larger of the two, namely $\max(r_{1},r_{2})=r_{1}$, as the joint distance. We
apply such a procedure on all reference points and obtain the empirical CDF
from the joint distances as the joint 3NN-CDF between the blue and orange
points.

Figure 1:  Illustration of the $k$NN-CDFs and joint $k$NN-CDFs computations
for $k=3$ with the Euclidean metric. The orange triangles and blue circles
represent the source catalog (IceCube neutrinos in our case) and tracers (such
as the WISE-2MASS samples), respectively. The two black pluses are reference
points to demonstrate the calculation of $3$NN-CDFs and joint $3$NN-CDFs,
respectively.

###  2.3 $k$NN-CDFs Cross-correlation

To isolate the effect of the clustering from two individual data samples
(density fields), we can subtract the product of the individual $k$NN-CDFs
that are measured from each dataset and obtain a measure of the $k$NN-CDFs
cross-correlation of the two data samples:

$\psi^{k_{1},k_{2}}(r)=\text{CDF}_{k_{1},k_{2}}(r)-\text{CDF}_{k_{1}}^{(1)}(r)\cdot\text{CDF}_{k_{2}}^{(2)}(r).$
(3)

Higher positive values of the cross-correlation defined above imply a higher
level of spatial correlation between the two density fields, while negative
values indicate anti-correlations (Banerjee & Abel, 2021b).

Note that while $k$NN-CDFs and the $k$NN-CDFs cross-correlation can be applied
on samples with different sample sizes by choosing $k_{1}\neq k_{2}$ in
Equation (3). For simplicity, in this paper, we restrict our scope to cases
when $k_{1}=k_{2}=k$ and require that the two data samples have roughly the
same number of data points. To reach this setting, the WISE-2MASS galaxy
sample is down-sampled to match the sample size of the neutrino catalog. A
discussion of how to handle datasets with very different sizes is presented in
Section 7.

##  3 Methodology

In this section, we outline how to compute the $k$NN-CDFs, joint $k$NN-CDFs,
and the $k$NN-CDFs cross-correlation we defined in Section 2 respectively,
specifically on the IceCube point-source data with adaptations of the original
computations outlined in Banerjee & Abel (2021a) and Banerjee & Abel (2021b).

###  3.1 Estimating $k$NN-CDFs

To estimate the $k$NN-CDFs, we first generate a group of reference points from
a homogeneous spatial Poisson process that fills the volume (or surface) for
both the IceCube point-source data and the galaxy catalog. In IceCube’s
context, since all data points lack depth information, we can treat them as
data points distributed on the surface of a unit sphere with right ascensions
(RAs) and declinations (Decs) in radians. The homogeneous spatial Poisson
process is generated on the same unit sphere.

At a specific $k$, instead of performing $k$NN searching within the data
sample, we measure all the $k$-nearest neighbor distances from the reference
points to the IceCube point-source data or the selected tracers. To better
characterize the two tails of the empirical CDF, we use more reference points
than the data points. We can immediately estimate the $k$NN-CDFs by obtaining
the empirical cumulative distribution function by using a 1-D linear
interpolation on the sorted $k$NN distances. As noted by Banerjee & Abel
(2021a), $k$-nearest neighbor searching can be done with high computational
efficiency by constructing a $k$-d tree (Bentley, 1975) on a total of $N$ data
points with $NlogN$ complexity, and the nearest neighbor searching can be
completed in $logN$ time for each reference point.

In our case, instead of using a $k$-d tree and the Euclidean metric, the
nearest neighbor searching is performed using the great-circle distances by
constructing a ball-tree (Omohundro, 1989) on the data points with the Scikit-
learn’s sklearn.neighbors.BallTree implementation (Pedregosa et al., 2011)
using the haversine distance metric (based on great-circle distances). For the
above implementation, the great-circle distance between two points $x,y$ in
the geographic coordinate system is defined as:

$d(x,y)=2r\arcsin{\sqrt{\sin^{2}({\frac{\Delta\phi}{2}})+\cos{\phi_{x}}\cos{\phi_{y}}\sin^{2}({\frac{\Delta\lambda}{2}})}},$

where $r$ is the radius of the sphere ($r=1$ in our case) and
$\Delta\phi=\phi_{x}-\phi_{y}$ and $\Delta\lambda=\lambda_{x}-\lambda_{y}$.
The $\phi$ and $\lambda$ are the latitude and longitude (both in radians),
respectively, of the two data points. We replace the latitudes and longitudes
with an equatorial coordinate system using RA and Dec (J2000) to match the
format of the IceCube public point-source data. The choice of a ball-tree is
due to the fact that while the $k$-d tree partitions the entire space along
the axis ($x$, $y$, and $z$ axis in three-dimensional Cartesian spaces), the
ball-tree only partitions the manifolds on which the data points are located,
and is more suitable when sources are all projected on a spherical surface.

###  3.2 Estimating Joint $k$NN-CDFs and Cross-correlation

When considering two data samples at a specific $k$, each reference point is
separately associated with two $k$NN distances, each from one of the datasets.
To measure the joint term $\text{CDF}_{k_{1},k_{2}}(r)$, for a reference point
at a specific $k=k_{1}=k_{2}$, the associated $k$NN distances used to obtain
the estimated CDF is the maximum of the distances measured from the two data
samples. Similar to the interpretation of $k$NN-CDFs on a single dataset, such
a joint term can be interpreted as the probability of finding spheres of
radius $r$ that enclose at least $k$ data points from both datasets (Banerjee
& Abel, 2021a). This effectively serves as a measure of the similarities in
terms of spatial clustering of the two data samples.

For both the $k$NN-CDFs and joint $k$NN-CDFs, the $k$ corresponding empirical
CDFs are estimated with linear interpolation using Python’s
scipy.interpolate.interp1d (Virtanen et al., 2020) on the middle 90% (between
the 5th and 95th quantiles) of the sorted $k$NN distances to avoid over-
extrapolating into both tails where there may be only a few data points. Then
we generate 50 linearly separated bins (in radians) on the same distance range
for all chosen $k$. Each $k$NN-CDF and joint $k$NN-CDF are represented by
length-50 vectors, where each entry represents the value of the interpolating
function at the corresponding radius. The final $k$NN-CDFs cross-correlations
are calculated directly by plugging in the obtained quantities in Equation
(3). We provide sample code to demonstrate our $k$NN-CDF and joint $k$NN-CDF
calculations in the Software section at the end of Section 7.4.

##  4 Data Samples and Analysis Setup

For a demonstration of our method, we perform a $k$NN-CDFs cross-correlation
analysis using a neutrino sample and a galaxy catalog. In this section, we
describe the neutrino and galaxy samples used in the study and the procedures
for synthetic events generation.

###  4.1 Galaxy and IceCube Neutrino Sample

We use the public all-sky point-source IceCube data from years 2008-2018
(IceCube Collaboration et al., 2021). To avoid the large atmospheric muon
background, we select neutrino events from the northern sky with Dec
$\delta>-5^{\circ}$. We mask the Galactic Plane between Galactic latitudes
$|b|<10^{\circ}$.

The WISE-2MASS galaxy catalog is constructed by combining the photometric data
from the WISE (Wright et al., 2010) and the 2MASS (Skrutskie et al., 2006)
infrared databases following Kovács & Szapudi (2015) as in Fang et al. (2020).
As explained in Section 2.3, we down-sample the galaxy catalog by randomly
selecting $n=451,953$ (about 17% of the original data) galaxies to match the
sample size of the neutrino data after the masking.

###  4.2 Synthetic Neutrino Samples

Synthetic neutrino samples are used to conduct statistical power analyses.
Atmospheric events are generated by randomizing the RAs of the original data.
To generate astrophysical neutrinos, we first randomly sample a fixed number
of sources from the galaxy catalog. The sources are assumed to be equally
bright and each follows a power-law energy spectrum
$dN/dEdAdt=A_{0}\,(E/10^{3}\,{\rm GeV})^{-2}$ between $E_{\rm min}=100\,\rm
GeV$ and $E_{\rm max}=10^{9}\,\rm GeV$. Then for each selected source, we
generate events based on the effective areas and angular resolutions of the
public point-source data. The number of events emitted by a source is a
Poisson random number with the mean being the assumed source flux integrated
over the detector effective area for a given source direction. Finally, we
combine background-only events and astrophysical events from the selected
galaxies to obtain a synthetic sample that contains a total of $N_{\rm ev}$
events and ${f}_{\text{astro}}$ fraction of astrophysical events.

The flux normalization is fixed to $A_{0}=5\times{10}^{-16}\,\rm
GeV^{-1}\,cm^{-2}\,s^{-1}$ in our default cases, which corresponds to a median
of 5 events per source. The source number is determined by the total number of
astrophysical events, $N_{\rm ev}\,f_{\rm astro}$, divided by the event per
source. For example, with ${f}_{\text{astro}}=3.74\%$ and the default $A_{0}$,
the number of sources selected to inject neutrinos is 3500. Our synthetic
event generation code is provided in the Software section at the end of
Section 7.4.

##  5 Likelihood and Sensitivity

Below we define the likelihood function and an associated likelihood ratio
test for the computed $k$NN-CDFs cross-correlation (referred to simply as
“cross-correlation” in what follows). We verify the properties of the
likelihood and estimate the sensitivity and discovery potential of the search.

###  5.1 Likelihood Function and Likelihood Ratio Test

Since the astrophysical and atmospheric neutrino come from unrelated sources
with different density fields, there should not be a contribution to the
correlation with the WISE-2MASS samples. Then the cross-correlation
approximately satisfies:

$\bm{\psi^{(k)}_{\text{g,v}}}(r)=f_{\text{astro}}\cdot\bm{\psi^{(k)}_{\text{g,astro}}}(r)+(1-f_{\text{astro}})\cdot\bm{\psi^{(k)}_{\text{g,atm}}}(r),$
(4)

where $\bm{\psi^{(k)}_{\text{g,v}}}(r)$ denotes the cross-correlation with
$k=k_{1}=k_{2}$ between a neutrino sample of mixed astrophysical and
atmospheric signals with a specific proportion of astrophysical signals,
$f_{\text{astro}}$, and the galaxy sample. The
$\bm{\psi^{(k)}_{\text{g,astro}}}(r)$ denotes the cross-correlation between a
neutrino sample with only astrophysical ($f_{\text{astro}}=1$) events and the
galaxy catalog, and $\bm{\psi^{(k)}_{\text{g,atm}}}(r)$ denotes the cross-
correlation between a neutrino sample with only atmospheric
($f_{\text{astro}}=0$) events and the galaxy catalog. Notice that there should
not be a spatial correlation between the galaxy sample and the atmospheric-
only neutrinos. Hence, $\bm{\psi^{(k)}_{\text{g,atm}}}$ likely captures a
measurement noise due to the spatial distribution of the atmospheric neutrinos
and the declination dependence of the detector effective area, which can be
considered as the benchmark against the detection of astrophysical neutrinos.
The relationship in Equation (4) is verified with synthetic samples in Section
5.2.

As specified in Section 3.2, the cross-correlation is discretized as a vector
with 50 linearly-separated radius bins (in radian). Under the assumption that
the cross-correlation at each radius bin follows a univariate Gaussian
distribution and measurements from different bins are dependent as they are
monotone increasing functions, which makes the cross-correlation function
multivariate Gaussian (verified in Section 5.2). The (log) likelihood function
can be defined as:

$\text{log}\mathcal{L}(f_{\text{astro}}|\bm{\psi}^{(k)}_{\text{g,v}})=-\frac{1}{2}(\bm{\psi}^{(k)}_{\text{g,v}}-\bm{\mu})\bm{\Sigma}^{-1}(\bm{{\psi}^{(k)}_{\text{g,v}}}-\bm{\mu})^{T},$
(5)

where $\bm{\mu}$ is the expected mean of the cross-correlation of a sample of
combined astrophysical and atmospheric events with $f_{\text{astro}}$ and
$\bm{\Sigma}$ is the associated covariance matrix of the cross-correlation
functions. Since we only observed a single collection of neutrino events,
lacking a clear analytical expression, we need to estimate the covariance
matrix with Monte Carlo synthetic samples, which can be practically difficult
as we are working with a high volume of data and a continuous parameter space.
To this end, the expected mean is estimated with Equation (4), and the
covariance matrix is estimated with the covariance matrix of the cross-
correlation between a set of 1000 Monte Carlo simulations of atmospheric-only
events and the galaxy catalog, which only needs to be computed once.

Suppose we generate $n$ Monte Carlo simulations of atmospheric-only events in
total, we denote the cross-correlation at $k$ between the $i$th simulated
atmospheric-only events and the galaxy catalog as
$\bm{\psi^{(k)}_{\text{g,atm,i}}}$. Then an unbiased estimate of the
covariance matrix in Equation (5) is given as:

$\bm{\hat{\Sigma}}=\frac{1}{n-1}\sum_{i=1}^{n}(\bm{\psi^{(k)}_{\text{g,atm,i}}}-\bm{\overline{\psi}})(\bm{\psi^{(k)}_{\text{g,atm,i}}}-\bm{\overline{\psi}})^{T},$
(6)

where $\bm{\overline{\psi}}$ is the sample mean of the $n$ computed cross-
correlations. The inverse is directly taken on the unbiased estimation of the
covariance matrix above. For a specific $f_{\text{astro}}$ and $k$, to
estimate the mean defined with Equation (4), we first generate $n$
atmospheric-only ($f_{\text{astro}}=0$) and astrophysical-only events
($f_{\text{astro}}=1$). The $\bm{\psi^{(k)}_{\text{g,astro}}}$ and
$\bm{\psi^{(k)}_{\text{g,astro}}}$ are estimated with the sample mean of $n$
cross-correlations between the two groups of simulations with the galaxy
catalog, respectively. Then the estimate for $\bm{\psi^{(k)}_{\text{g,v}}}(r)$
is obtained by directly plugging in $f_{\text{astro}}$ and the two sample
means in Equation (4).

The statistical significance of the measured cross-correlation against the
null hypothesis, which is defined as the cross-correlation between the
atmospheric-only events and the source catalog, can be quantified with a
likelihood ratio test with test statistics defined as:

$\text{TS}\equiv
2[\text{log}\mathcal{L}(\hat{f}_{\text{astro}})-\text{log}\mathcal{L}({f}_{\text{astro}}=0)],$
(7)

where $\hat{f}_{\text{astro}}$ is the maximum likelihood estimation for
${f}_{\text{astro}}$. By Wilk Theorem (Wilks, 1938), the above test statistic
asymptotically follows a $\chi^{2}$ distribution with a degree of freedom (df)
equals the number of free parameters in the likelihood function, which is the
${f}_{\text{astro}}$ in our case (df=1).

###  5.2 Verification of Key Assumptions

Here we verify the key assumptions stated in Section 5.1 with $k=8$ at a few
selected ${f}_{\text{astro}}$ values. We choose $k=8$ as it optimizes the
statistical power of the likelihood ratio test and is computationally
efficient. More discussion regarding how to choose $k$ for an analysis is
presented in Section 7.1.

Since we assume a multivariate Gaussian distribution of the estimated cross-
correlation, we test the normality of the cross-correlations measured between
the galaxy sample and 1000 synthetic background-only events
(${f}_{\text{astro}}=0$) with a Henze-Zirkler multivariate normality test
(Henze & Zirkler, 1990). Note that such tests can be performed using synthetic
samples with different ${f}_{\text{astro}}$. The resulting p-value is $0.516$,
so we do not reject the null hypothesis that the background-only cross-
correlations follow a multivariate normal distribution. If the cross-
correlations follow a multivariate normal distribution with mean $\bm{\mu}$
and covariance matrix $\bm{\Sigma}$, $\bm{\psi}^{(k)}_{\text{g,v}}\sim
N_{p}(\bm{\mu},\bm{\Sigma})$ and the covariance matrix $\bm{\Sigma})$ is
positive definite, then we have:

$(\bm{\psi}^{(k)}_{\text{g,v}}-\bm{\mu})^{T}\bm{\Sigma}^{-1}(\bm{\psi}^{(k)}_{\text{g,v}}-\bm{\mu})\sim\chi_{p}^{2},$
(8)

where $p$ is the number of radius bins we selected when estimating the cross-
correlations. In our case, we choose $p=50$. We can visually examine the
conclusion of Equation (8) by plotting the empirical cumulative probabilities
of the calculated chi-squared values against that of a $\chi^{2}$ distribution
with $p=50$ degrees of freedom. We give the visual examination for $k=8$ and
${f}_{\text{astro}}=0$ in Figure 2, which shows a close fit of the empirical
CDF with the expected one.

Then we verify the approximation of Equation (4) with 1000 synthetic neutrino
samples and background-only samples with ${f}_{\text{astro}}=0.0534$ and
${f}_{\text{astro}}=0.107$ for $k=8$ in Figure 3. Notice that the
approximations are fairly good with slight right shifts, which may be
introduced by the uncertainty when estimating
$\bm{\psi^{(k)}_{\text{g,astro}}}(r)$ and $\bm{\psi^{(k)}_{\text{g,atm}}}(r)$.

Figure 2:  The cumulative probabilities of the calculations from Equation (8)
with ${f}_{\text{astro}}=0$ at $k=8$ (the black dashed curve) compared with
that of the $\chi^{2}$ distribution with 50 degrees of freedom (orange solid
curve). Figure 3:  Mean cross-correlations between 1000 synthetic neutrino
samples with ${f}_{\text{astro}}=0.0534$ (blue lines) and $0.107$ (orange
lines) and WISE-2MASS galaxy sample for $k=8$ together with the mean cross-
correlations obtained from the background-only samples (the gray line). The
approximations based on Equation (4) are shown as dashed curves.

###  5.3 Maximum Likelihood Estimation of ${f}_{\text{astro}}$

To evaluate the performance of the likelihood function from Equation (5), we
compute the bias of ${f}_{\text{astro}}$, defined as the differences in its
best-fit and true values, and show the results in Figure 4. The top panel
indicates a likelihood using an approximate covariance matrix based on 1000
Monte Carlo simulations of atmospheric-only events. The bottom panel indicates
a likelihood with an accurate covariance matrix computed using the true
${f}_{\text{astro}}$. When using an accurate covariance matrix, we find no
systematic bias for the estimation across a large range of
${f}_{\text{astro}}$ values. When using the approximate covariance matrix, we
observe a systematic underestimation of ${f}_{\text{astro}}$. The differences
between the estimations and true values are up to the 10% level. Given the
relatively small bias and the unknown ${f}_{\text{astro}}$, we adopt the
likelihood with the approximate covariance matrix. In practice, this could be
improved by maximizing the likelihood using tabulated covariance matrices for
various ${f}_{\text{astro}}$ values and interpolating the tables during the
estimation, rather than assuming a fixed covariance matrix with
${f}_{\text{astro}}=0$.

Figure 4: Differences between the mean of the maximum likelihood estimations
of the ${f}_{\text{astro}}$ of 1000 synthetics neutrino samples and the true
${f}_{\text{astro}}$. The chosen ${f}_{\text{astro}}$ are, from left to right,
0, 0.0186, 0.0534, 0.107, 0.213, 0.320, 0.532. The error bars are based on the
2$\sigma$ of the distributions of the 1000 estimations. The top panel shows
the estimations where the covariance matrices are estimated with the
background-only events. The bottom panel uses a likelihood with covariance
matrices estimated using the true ${f}_{\text{astro}}$. Figure 5:  Cumulative
probabilities of test statistics of the cross-correlation between neutrino and
galaxy samples via the KNN-CDF method. The black curve indicates the null
hypothesis with ${f}_{\text{astro}}=0$, which is consistent with a $\chi^{2}$
distribution with 1 degree of freedom as indicated by the grey dashed curve.
The orange curve indicates neutrino samples with ${f}_{\text{astro}}=1.86\%$
and $A_{0}=5\times 10^{-16}\,\rm GeV^{-1}cm^{-2}s^{-1}$ (about 5 neutrinos per
source). The blue curves are computed using neutrino samples with
${f}_{\text{astro}}=3.74\%$, with dark and light colors indicating $A_{0}=5$
and $1\times 10^{-16}\,\rm GeV^{-1}cm^{-2}s^{-1}$, respectively. The red
dashed line indicates the 50% cumulative probability.

###  5.4 Constraints on Sensitivity and Discovery Potential

Sensitivity and discovery potential are two crucial quantities used to
evaluate the performance of the method. An ideal method should identify
statistically significant cross-correlation signals with a relatively low
${f}_{\text{astro}}$. Following the IceCube convention (e.g., Aartsen et al.,
2020a), we define the sensitivity as the astrophysical neutrino flux, measured
as ${f}_{\text{astro}}$ in a data sample, that yields a statistical power of
0.9 at a statistical significance level of $\alpha=0.5$ (0.5 false-positive
probability). The discovery potential is defined as the astrophysical flux
that results in a statistical power of 0.5 at a $5\sigma$ significance level
($\alpha=2.85\times 10^{-7}$). To evaluate the performance of our framework,
we compute the statistical power of the test at chosen ${f}_{\text{astro}}$
values with the test statistics defined in Equation (7) under the hypothesis:

$H_{0}:{f}_{\text{astro}}=0,\quad H_{1}:{f}_{\text{astro}}>0,$ (9)

where we choose $k=8$ for this calculation. The choice of $k$ depends on
various factors, such as the data samples’ densities and the events’ angular
uncertainties. We will discuss how to choose $k$ and its impact on the
statistical power in Section 7. The ${f}_{\text{astro}}$ values are chosen by
varying the number of galaxies selected to inject synthetic neutrino events.
The numbers we selected are 1500, 1750, 2000, 3500, 4000, 4250, 4500, 4750,
and 5000, which covers the ${f}_{\text{astro}}$ from 1.86% to 5.34%.

Figure 5 presents the cumulative probabilities of the test statistics obtained
by cross-correlating the galaxy sample and synthetic neutrino events. We
verify that the test statistic under the null hypothesis follows the
$\chi^{2}$ distribution with 1 degree of freedom. Among the
${f}_{\text{astro}}$ values we tested, the case with
${f}_{\text{astro}}=1.86\%$ (corresponding to a total of 1750 galaxies
selected to inject synthetic neutrino events) and ${f}_{\text{astro}}=3.74\%$
(corresponding to 3500 galaxies selected) already respectively reached the
$3\sigma$ and $5\sigma$ significance levels with $>0.5$ statistical powers. We
further tested another case with ${f}_{\text{astro}}=2.13\%$ (corresponding to
a total of 2000 galaxies selected to inject synthetic neutrino events) and
found that this failed to reach the $5\sigma$ significance levels with $>0.5$
statistical powers. As we evaluated the statistical power on a grid of $f_{\rm
astro}$, the true discovery potential should be $[2.13\%,3.74\%]$. A more
accurate estimation can be obtained with a binary search on a more refined
grid. We find empirically that the statistical power of the test depends on
the assumed source brightness. Under the same ${f}_{\text{astro}}$, the
statistical power is higher with a lower $A_{0}$, which corresponds to a
scenario with a larger population of dimmer individual sources (see Section
7.1 for more discussion regarding $A_{0}$).

##  6 Results

Figure 6 presents the result obtained by applying our framework to the public
IceCube ten-year point-source data (IceCube Collaboration et al., 2021;
IceCube Collaboration, 2021) and the down-sampled WISE-2MASS galaxy catalog.
The dashed curve shows the $k$NN-CDFs cross-correlation across various angular
sizes for $k=8$. It does not exhibit significant deviation from a cross-
correlation using the background-only events. With the log-likelihood defined
in Equation (5), we obtain a best-fit $\hat{f}_{\text{astro}}=-0.00128$ with a
likelihood ratio test statistic $\text{TS}=0.0589$ (corresponding to a p-value
$0.808$). The negative estimation here is due to our choice of the uniform
prior $f_{\rm astro}\in[-1,1]$. We also take a Bayesian approach and sample
the parameter space of ${f}_{\text{astro}}$ using Markov Chain Monte Carlo
(MCMC) with the uniform prior. With 500 walkers and 2000 steps, the 95%
Bayesian credible interval based on the 5th and 95th quantiles region of the
posterior is $(-0.0115,0.009)$. The results suggest that we find no
statistically significant evidence of astrophysical neutrinos in the IceCube
ten-year point-source data that have spatial cross-correlations with the
selected WISE–2MASS galaxy sample.

Figure 6:  The cross-correlation between the IceCube ten-year point-source
data and our down-sized WISE-2MASS galaxy sample across the scales of interest
at $k=8$ (black dashed line). For comparison, we show the average cross-
correlations from synthetic samples with ${f}_{\text{astro}}=0.0186$ (solid
orange curve) and background-only events (solid gray curve).

##  7 Discussion and Conclusions

In this section, we discuss possible extensions and improvements of our
methods as well as limitations.

###  7.1 The Choice of $k$

The $k$NN-CDFs are ways to measure the degree of spatial clustering at
different scales, which are controlled by the choice of $k$. When computing
the cross-correlations of two data samples, spatial similarities can occur at
some specific scales. This indicates that we can expect to achieve the most
significant cross-correlation signal (with higher power) at some specific $k$.
To demonstrate such a proposition, we compute the statistical power of the
likelihood ratio test formulated in Equation (7) of the cross-correlation
between 1000 synthetic neutrino samples with ${f}_{\text{astro}}=0.0374\text{
and }0.0453$ and the WISE-2MASS galaxy sample under the $5\sigma$ rejection
threshold. We plot the statistical power of the test with $k$ ranging from 1
to 30 in Figure 7.

Figure 7:  Statistical power (1 - Type II error) of the proposed likelihood-
ratio tests, Equation (7), on synthetic neutrino samples with
${f}_{\text{astro}}=0.0374$ (blue) and $0.0453$ (orange) across different $k$.
The rejection threshold is set to be 26.349, which corresponds to a $5\sigma$
rejection threshold with a significance level of $\alpha=2.85\times 10^{-7}$
when the null distribution is a $\chi^{2}$ distribution with 1 degree of
freedom. A good choice of $k$ is the smallest one that maximizes the
statistical power and (ideally) exceeds the 0.5 power threshold.

At low $k$ values, the statistical power is low and reaches a plateau around
$k=8$, and then becomes fairly stable. Lower $k$ corresponds to the angular
uncertainty ranges of the synthetic neutrino data which may lead to reduced
statistical power, or there may simply not be significant correlated spatial
clustering between the two samples at the lower scales. Such structures are
successfully extracted by increasing the scale we are examining by increasing
the $k$. While it is possible to compute at a higher $k$, it suffices to
choose a lower $k$ that gives sufficiently high power to lower the
computational cost. In practice, a reasonable approach for carrying out the
analysis is to consider conducting similar experiments to select the $k$ by
examining the relationship between the statistical power and $k$. An ideal
choice of $k$ will be the smallest one that leads to a cross-correlation with
$r$ that exceeds the angular uncertainties of the neutrino catalog and
maximizes the statistical power and ideally exceeds the 0.5 power threshold.

As noted in Section 5.4, at the same ${f}_{\text{astro}}$, the statistical
power of the test depends on the source brightness (number of events per
galaxy). Figure 8 compares the power of the test at
${f}_{\text{astro}}=0.0374$ under three source brightness assumptions:
$A_{0}=1\times{10}^{-16}$, $5\times{10}^{-16}$, and $1\times{10}^{-15}\,\rm
GeV^{-1}cm^{-2}s^{-1}$, which correspond to a median of 1, 5, and 10 events
per source. The three cases present similar performances as a function of $k$,
though a higher power is obtained for a lower $A_{0}$. This could be because
the sources in the 1 event per source case better match the galaxy catalog.

Figure 8:  Statistical powers (1 - Type II error) of the proposed likelihood-
ratio tests, Equation (7), on synthetic neutrino samples with
${f}_{\text{astro}}=0.0374$ at three different brightness of the sources. The
lines with orange circles, blue diamonds, and green squares correspond to a
median of 1 event, 5 events, and 10 events per source, respectively.

###  7.2 Datasets with Different Sample Sizes

In this paper, we have simplified the analysis by down-sampling the dataset
with a larger sample size, the WISE-2MASS catalog, to match the sparser
IceCube point-source data. Below we discuss how to implement an analysis with
$k_{1}\neq k_{2}$.

One way to handle such a scenario is to adjust $k_{1}$ and $k_{2}$ based on
the ratio of the sample sizes (average densities). This can put the two
samples on equal footing when measuring their spatial similarities via the
joint $k$NN-CDFs. Consider two homogeneous spatial Poisson processes $D_{1}$
and $D_{2}$ with sample sizes $2n$ and $n$, respectively, in the same sample
volume $V$. We expect the average $2$NN distances of the denser sample $D_{1}$
to be roughly equal to the average $1$NN distances of the sparser sample
$D_{2}$. Indeed, in a two-dimensional Euclidean space, the actual ratio of
average 2NN distance to 1NN distances is $\frac{3}{4}\sqrt{2}\approx 1.0606$
(see Appendix A). Thus an intuitive solution would be choosing $k_{1}=2$ for
$D_{1}$ and $k_{2}=1$ for $D_{2}$.

However, such a clear relationship may not exist for all $k$s. For example,
the relationship discussed above for two homogeneous spatial Poisson processes
does not exist for odd values of $k_{1}$. Also, if the sample sizes have a
highly significant difference, it will be unrealistic to adjust the choices of
$k$ simply based on the ratio of the sample sizes (average densities).
Therefore, we propose a more efficient way for framing such problems by
optimizing the areas between $\psi^{(k_{1},k_{2})}_{\text{g,astro}}(r)$ (the
signal) and $\psi^{(k_{1},k_{2})}_{\text{g,atm}}$ (negative control) as
defined in Section 5.1:

$(k_{1},k_{2})=\arg\max_{(k_{1},k_{2})}\int_{0}^{r^{\prime}}\psi^{(k_{1},k_{2})}_{\text{g,astro}}(r)-\psi^{(k_{1},k_{2})}_{\text{g,atm}}(r)\,dr.$
(10)

We suggest choosing $r^{\prime}$ to be slightly larger than the largest $k$NN
distances measured across all selected $k$ for all data samples, which we
adopted in this paper. Another approach can be selecting the range for
integration based on the 5th and 95th quantiles of the measured $k$NN
distances similar to what is outlined in Section 3.2. In practice, this can be
optimized over experimental or synthetic data or with a train-test split on
real data.

Figure 9:  Areas between the two cross-correlations calculated from the
integral defined in Equation (10) across all combinations of $k_{1}$ and
$k_{2}$ up to $k_{1}=k_{2}=30$. $k_{1}$ (x-axis) corresponds to the neutrino
sample and $k_{2}$ (y-axis) corresponds to the galaxy sample. The left panel
demonstrates a case where the sample size of the galaxy sample matches the
sample size of the neutrino sample, while the right panel illustrates a case
where the sample size of the galaxy sample is roughly half of that of the
neutrino sample. Orange stars stand for the $(k_{1},k_{2})$ pairs that
maximize the integrals for a fixed $k_{1}$. The red pluses indicate the
diagonal ($k_{1}=k_{2}$; left) and cells with $k_{2}=k_{1}/2$ for all even
values of $k_{1}$.

We illustrate such an idea with a synthetic neutrino sample with
${f}_{\text{astro}}=0.0186$ and an atmospheric-only sample with
${f}_{\text{astro}}=0$, both with a sample size of $n=451,953$, in Figure 9.
We first match the sample size of the WISE galaxy sample to both the synthetic
neutrino sample and the atmospheric-only sample. Then we compute the integral
defined in Equation (10) across all combinations of $k_{1}$ and $k_{2}$ up to
$k_{1}=k_{2}=30$. The integral results (area between the two cross-
correlations) are displayed on the left panel of Figure 9. Notice that for a
fixed $k_{1}$, the optimal choice of $k_{2}$ that gives the largest integrals
(indicated by the orange stars) consistently occurs around the diagonal
(indicated by the red pluses for even $k_{1}$). This matches our expectation
and intuition that if the two samples have comparable or the same sample
sizes, the maximum cross-correlation should occur (if there is such a
correlation) when $k_{1}$ is roughly equal to $k_{2}$.

Then, we down-sampled the WISE-2MASS galaxy sample so that its sample size is
roughly half of the synthetic neutrino sample and atmospheric-only sample. We
compute the same integral defined in Equation (10) across all combinations of
$k_{1}$ and $k_{2}$ up to $k_{1}=k_{2}=30$, and the results are displayed at
the right panel of Figure 9. Notice that for a fixed $k_{1}$, the optimal
choice of $k_{2}$ that gives the largest integrals (indicated by the orange
stars) consistently occurs around $\frac{1}{2}k_{1}$ (denoted with the red
pluses at entries with $(\frac{1}{2}k_{1},k_{2})$ for even $k_{1}$). This
matches our discussions for the two homogeneous spatial Poisson processes with
sample sizes $2n$ and $n$.

Figure 10:  Same as Figure 9 but with down-sampled neutrino data and a larger
$k$ range. The orange star with $(k_{1},k_{2})=(180,170)$ denotes the entry
with maximum area across all combinations of $k_{1}$ and $k_{2}$.

On the other hand, if there is a significant correlation between the two
samples, we expect to observe higher correlation values at certain scales
rather than uniformly across all scales considered. Recall that $k$NN-CDFs can
be understood as a measure of the degree of spatial clustering at different
scales with “probing spheres” on the data sample. For both very small scale
(smaller than the angular uncertainty of the observatory) and sufficiently
large scale, the clustering of the data sample behaves more like background
noise. For the equal-sample-size case, we do expect that the integral will
increase to a maximum value and then start decreasing as $(k_{1},k_{2})$
increases. The shape of the distribution of the integral defined in Equation
(10) is concave downward for both cases. We illustrate this observation in
Figure 10, where we compute the same areas in the left panel of Figure 9 but
with neutrino data with a sample size of $n=10000$ and
${f}_{\text{astro}}=0.0186$. The decreased sample size enables us to observe
the concave shape at a smaller $k$ with less computation time. Notice that the
area peaks around the orange star (denoting the maximum area) with a concave
downward shape. Especially for very small $(k_{1},k_{2})$, the areas become
smaller and closer to zero. We also expect a decreasing trend of areas for
$(k_{1},k_{2})$ that are larger than the orange star. In practice, we suggest
considering both the optimization results and the power-$k$ curves as in
Figure 5 when selecting the optimal $k_{1}$ and $k_{2}$.

###  7.3 Measurement with Higher-energy Events

We also tested our method’s performance on neutrino samples with smaller sizes
and higher purity of astrophysical events. This can be achieved by using
neutrino events with higher energies. As a demonstration, we compute the
$k$NN-CDFs cross-correlations using a neutrino sample selected from the public
ten-year point source data with energy higher than $10^{4}\,\rm GeV$. The
sample size is reduced to $n=1539$ and we adopt a smaller $k$ for the
analysis, $k=2$. Figure 11 presents the sensitivity of such an analysis.
Similar to the procedures as described in Section 5.4, we tested several
${f}_{\text{astro}}$ values and found that ${f}_{\text{astro}}=10.68\%$ is
needed to reach the $3\sigma$ significance level with $>0.5$ statistical
power.

Figure 11:  Same as Figure 5 but the related measurements are computed with
neutrino events with energy $\geq 10^{4}\,\rm GeV$. The green curve indicates
neutrino samples with ${f}_{\text{astro}}=10.86\%$ and $A_{0}=5.8\times
10^{-16}\,\rm GeV^{-1}cm^{-2}s^{-1}$ (about 1 neutrino per source under the
above energy setting).

###  7.4 Conclusion and Future Work

Compared with previous cross-correlation studies, our framework can better
handle the non-Gaussian features in the dataset with improved computational
efficiency and sensitivity. The $k$NN-CDFs and the associated cross-
correlations also provide an interpretable characterization of the spatial
clustering and spatial similarities between data samples. Such methodologies
can lead to a better understanding of the spatial structure of the detected
neutrinos across different scales.

We find no significant cross-correlation between the IceCube public ten-year
point-source data and the selected WISE-2MASS galaxy sample. This is not
surprising as our analysis does not account for the energy information of the
neutrinos. With the majority of the neutrino events being atmospheric, the
$f_{\rm astro}$ of the sample is well below the percent-level, which is needed
to achieve a detection as shown in Figure 5. An improvement to our study would
be to invoke an energy-dependent analysis, for example by performing the
cross-correlation study in different neutrino energy bins (Guevel & Fang,
2023) or using neutrino samples with higher purity (Abbasi et al., 2021,
2024). The study may be extended by using cosmological survey products with
higher densities, different redshift ranges, and larger sky coverage (Fang et
al., 2020). Different tracers other than galaxies can be analyzed as well.

Our method currently works under the assumption that the neutrino and galaxy
catalogs have roughly the same sample sizes on the spherical surface. As
discussed in section 7.2, an important future direction is to generalize the
framework to catalogs with significantly different sample sizes. An
alternative way to solve such problems can be approximating the denser galaxy
catalog as a continuous field on some HEALPix pixelation, and then apply the
$k$NN version of the cross-correlation between discrete tracers and continuous
fields as formulated in Banerjee & Abel (2023). This potentially retains the
small-scale information in the galaxies while avoiding down-sampling the
catalog. Such a formalism has been applied to study gravitational wave (GW)
transient sources and galaxy samples (Gupta & Banerjee, 2024). Also, we only
select a specific $k$ for the statistical testing and estimation. It is
possible to extend the framework by incorporating multiple $k$NN-CDFs in a
single analysis during the statistical inference.

We thank Michael Larson for his help with using the IceCube public ten-year
point-source data. This research was performed using the computing resources
and assistance of the University of Wisconsin-Madison Center For High
Throughput Computing (CHTC) in the Department of Computer Sciences (Center for
High Throughput Computing, 2006). The CHTC is supported by the University of
Wisconsin-Madison, the Advanced Computing Initiative, the Wisconsin Alumni
Research Foundation, the Wisconsin Institutes for Discovery, and the National
Science Foundation, and is an active member of the Open Science Grid, which is
supported by the National Science Foundation and the U.S. Department of
Energy’s Office of Science. ZZ and JCK acknowledge support from NSF under
Grant Number DMS 2038556. ZZ acknowledges support from the Hilldale
Undergraduate/Faculty Research Fellowship awarded by the University of
Wisconsin-Madison. The work of KF is supported by the Office of the Vice
Chancellor for Research and Graduate Education at the University of Wisconsin-
Madison with funding from the Wisconsin Alumni Research Foundation. KF
acknowledges support from the National Science Foundation (PHY-2110821,
PHY-2238916) and NASA (NMH211ZDA001N-Fermi). This work was supported by a
grant from the Simons Foundation (00001470, KF). KF acknowledges the support
of the Sloan Research Fellowship.

## References

  * Aartsen et al. (2013) Aartsen, M. G., Abbasi, R., Abdou, Y., et al. 2013, Science, 342, doi:10.1126/science.1242856
  * Aartsen et al. (2018) Aartsen, M. G., et al. 2018, Science, 361, 147
  * Aartsen et al. (2020a) —. 2020a, Phys. Rev. Lett., 124, 051103
  * Aartsen et al. (2020b) —. 2020b, JCAP, 07, 042
  * Abbasi et al. (2021) Abbasi, R., et al. 2021, Phys. Rev. D, 104, 022002
  * Abbasi et al. (2022) Abbasi, R., Ackermann, M., Adams, J., et al. 2022, Science, 378, 538.  https://www.science.org/doi/abs/10.1126/science.abg3395
  * Abbasi et al. (2023) Abbasi, R., et al. 2023, Science, 380, adc9818
  * Abbasi et al. (2024) —. 2024, arXiv e-prints, arXiv:2402.18026
  * Banerjee & Abel (2021a) Banerjee, A., & Abel, T. 2021a, Monthly Notices of the Royal Astronomical Society, 500, 5479
  * Banerjee & Abel (2021b) —. 2021b, Monthly Notices of the Royal Astronomical Society, 504, 2911
  * Banerjee & Abel (2023) —. 2023, Monthly Notices of the Royal Astronomical Society, 519, 4856
  * Bentley (1975) Bentley, J. L. 1975, Communications of the ACM, 18, 509
  * Braun et al. (2008) Braun, J., Dumm, J., De Palma, F., et al. 2008, Astropart. Phys., 29, 299
  * Center for High Throughput Computing (2006) Center for High Throughput Computing. 2006, Center for High Throughput Computing, Center for High Throughput Computing, doi:10.21231/GNT1-HW21.  https://chtc.cs.wisc.edu/
  * Creque-Sarbinowski et al. (2022) Creque-Sarbinowski, C., Kamionkowski, M., & Zhou, B. 2022, Phys. Rev. D, 105, 123035
  * Fang et al. (2020) Fang, K., Banerjee, A., Charles, E., & Omori, Y. 2020, The Astrophysical Journal, 894, 112
  * Guevel & Fang (2023) Guevel, D., & Fang, K. 2023, arXiv e-prints, arXiv:2308.03978
  * Gupta & Banerjee (2024) Gupta, K. R., & Banerjee, A. 2024, arXiv preprint arXiv:2404.01428
  * Henze & Zirkler (1990) Henze, N., & Zirkler, B. 1990, Communications in statistics-Theory and Methods, 19, 3595
  * IceCube Collaboration (2021) IceCube Collaboration. 2021, arXiv e-prints, doi:10.21234/sxvs-mt83
  * IceCube Collaboration et al. (2021) IceCube Collaboration, Abbasi, R., Ackermann, M., et al. 2021, arXiv e-prints, arXiv:2101.09836
  * Kovács & Szapudi (2015) Kovács, A., & Szapudi, I. 2015, Monthly Notices of the Royal Astronomical Society, 448, 1305
  * Negro et al. (2023) Negro, M., Crnogorčević, M., Burns, E., et al. 2023, ApJ, 951, 83
  * Omohundro (1989) Omohundro, S. M. 1989, Five balltree construction algorithms (International Computer Science Institute Berkeley)
  * Ouellette & Holder (2024) Ouellette, A., & Holder, G. 2024, arXiv e-prints, arXiv:2405.09633
  * Pedregosa et al. (2011) Pedregosa, F., Varoquaux, G., Gramfort, A., et al. 2011, Journal of Machine Learning Research, 12, 2825.  http://jmlr.org/papers/v12/pedregosa11a.html
  * Skrutskie et al. (2006) Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163
  * Urban et al. (2021) Urban, F. R., Camera, S., & Alonso, D. 2021, A&A, 652, A41
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261
  * Wang et al. (2022) Wang, Y., Banerjee, A., & Abel, T. 2022, Monthly Notices of the Royal Astronomical Society, 514, 3828
  * Wilks (1938) Wilks, S. S. 1938, The Annals of Mathematical Statistics, 9, 60.  http://www.jstor.org/stable/2957648
  * Wright et al. (2010) Wright, E. L., Eisenhardt, P. R., Mainzer, A. K., et al. 2010, Astronomical Journal, 140, 1868

##  Appendix A Average 1NN and 2NN Distances for Homogeneous Spatial Poisson
Processes

The cumulative density function of the nearest neighbor (1NN) distance for
homogeneous spatial Poisson processes is written as Banerjee & Abel (2021a):

$\text{CDF}_{\text{1NN}}(V)=1-e^{-\lambda V}$ (A1)

where $V=\pi r^{2}$ and $\lambda$ is the intensity of the homogeneous spatial
Poisson process. Thus, we have:

$F_{R}(r)_{\text{1NN}}=\text{P}(R\leq r)=1-e^{-\lambda\pi r^{2}}$ (A2)

The probability mass function can be obtained by differentiating the
cumulative density function above:

$f(r)_{\text{1NN}}=\frac{d}{dr}F_{R}(r)_{\text{1NN}}=2\lambda\pi r\cdot
e^{-\lambda\pi r^{2}}$ (A3)

Then the expected 1NN distance can be obtained by integrating the above mass
function on the real line as:

$E[R]_{\text{1NN}}=\int r\cdot 2\lambda\pi r\cdot e^{-\lambda\pi
r^{2}}\,dr=\frac{1}{2\sqrt{\lambda}}$ (A4)

Similarly, with $V=\pi r^{2}$ and intensity as $\lambda$, the cumulative
distribution function of the 2NN distances for homogeneous spatial Poisson
processes can be written as Banerjee & Abel (2021a):

$F_{R}(r)_{\text{2NN}}=\text{P}(R\leq r)=1-e^{-\lambda\pi r^{2}}-(\lambda\pi
r^{2})e^{-\lambda\pi r^{2}}$ (A5)

The probability mass function can be obtained by differentiating the
cumulative density function above:

$f(r)_{\text{2NN}}=\frac{d}{dr}F_{R}(r)_{\text{2NN}}=2\lambda^{2}\pi^{2}r^{3}\cdot
e^{-\lambda\pi r^{2}}$ (A6)

Then the expected 2NN distance can be obtained by integrating the above mass
function:

$E[R]_{\text{2NN}}=\int r\cdot 2\lambda^{2}\pi^{2}r^{3}\cdot e^{-\lambda\pi
r^{2}}\,dr=\frac{3}{4\sqrt{\lambda}}$ (A7)

Consider two homogeneous spatial Poisson processes, $X_{1}$ and $X_{2}$, with
sample sizes $n$ and $2n$, respectively, on a two-dimensional surface. Then we
have intensities $\lambda_{X_{1}}=\frac{1}{2}\lambda_{X_{2}}$, and the average
1NN distance of $X_{1}$ and average 2NN distance of $X_{2}$ can be written as:

$E[R_{X_{1}}]_{\text{1NN}}=\frac{1}{2\sqrt{\lambda_{X_{1}}}}$ (A8)
$E[R_{X_{2}}]_{\text{2NN}}=\frac{3}{4\sqrt{2\lambda_{X_{1}}}}$ (A9)

The ratio of the 2NN distance of $X_{2}$ to the average 1NN distance of
$X_{1}$ is then
$\dfrac{E[R_{X_{2}}]_{\text{2NN}}}{E[R_{X_{1}}]_{\text{1NN}}}=\frac{3}{4}\sqrt{2}\approx
1.0606$.