.tifpng.pngconvert #1 .tiffpng.pngconvert #1 .tif, .tiff

# A Pre-training Oracle for Predicting Distances in Social Networks

Gunjan Mahindre  Electrical and Computer Engineering, Colorado State
UniversityFort CollinsUSA ,  Randy Paffenroth  Mathematical Sciences,
Worcester Polytechnic InstituteWorcesterUSA ,  Anura Jayasumana  Electrical
and Computer Engineering, Colorado State UniversityFort CollinsUSA and  Rasika
Karkare  Data Science, Worcester Polytechnic InstituteWorcesterUSA

###### Abstract.

In this paper, we propose a novel method to make distance predictions in real-
world social networks. As predicting missing distances is a difficult problem,
we take a two-stage approach. Structural parameters for families of synthetic
networks are first estimated from a small set of measurements of a real-world
network and these synthetic networks are then used to pre-train the predictive
neural networks. Since our model first searches for the most suitable
synthetic graph parameters which can be used as an “oracle” to create
arbitrarily large training data sets, we call our approach “Oracle Search Pre-
training” (OSP).

For example, many real-world networks exhibit a Power law structure in their
node degree distribution, so a Power law model can provide a foundation for
the desired oracle to generate synthetic pre-training networks, if the
appropriate Power law graph parameters can be estimated. Accordingly, we
conduct experiments on real-world Facebook, Email, and Train Bombing networks
and show that OSP outperforms models without pre-training, models pre-trained
with inaccurate parameters, and other distance prediction schemes such as Low-
rank Matrix Completion. _In particular, we achieve a prediction error of less
than one hop with only 1% of sampled distances from the social network._ OSP
can be easily extended to other domains such as random networks by choosing an
appropriate model to generate synthetic training data, and therefore promises
to impact many different network learning problems.

social network, autoencoder, pre-training, prediction, network sampling,
missing data

##  1\. Introduction

Shortest node-pair distances in graphs pack a lot of insight regarding the
network as a whole and thus have been used in several graph analysis
algorithms. In social network analysis, for example, hop-distances are used
for predicting friendships (Zhang and Pang, 2015; Wahid-Ul-Ashraf et al.,
2019), extracting sensor network topology (Jayasumana et al., 2019; Dhanapala
and Jayasumana, 2013), and modeling social networks (Sarkar and Moore, 2005;
Boguná et al., 2004) to name a few. However, missing measurements in sampled
datasets adversely affect the structural properties of social networks
(Kossinets, 2006; Sharma et al., 2015) and the original network
characteristics are misrepresented. Predicting these missing measurements
would solve this problem.

Neural networks have proved to be useful in solving many prediction problems
but they often require a large amount of training data to perform well. Pre-
training has alleviated this issue with transfer learning, especially when the
second task is much more difficult for training (Radford et al., 2018; Erhan
et al., 2010). However, relatively little is known about the pre-training
behavior related to neural networks used for social network analysis as social
networks are difficult to measure or sample completely due to their cost of
measurement, privacy, or storage issues. The ever increasing size of social
networks today only adds to this problem.

We demonstrate a prediction model based on efficient pre-training for distance
inference in social networks when only sparse samples are known. We use
artificially simulated training data to compensate for the lack of sufficient
real-world measurements. For accurate prediction, this artificial pre-training
data has to be faithful with the characteristics of the target social network.
However, we do not know the characteristics of the target network. All we have
is a small fraction of randomly measured node-pair distances. A way to
estimate the right pre-training parameters is required. The central idea of
this paper is to build an “Oracle” that points us towards these suitable set
of pre-training parameters.

Fig. 1 gives a preview of our results for a 1133 node real-world Email network
when only 0.5% of the total node-pair distances are known. We can see that
accurate pre-training parameters improve the prediction performance by 33%
over pre-training with inaccurate parameters and by 53% over models without
pre-training.

Figure 1. Comparison of proposed method with Low-rank Matrix Completion,
neural network without pre-training, and pre-training with inaccurate
parameters on a real world Virgili Emails network with only 0.5% of distances
sampled. OSP improves the performance by 53% from a neural network without
pre-training and by 33% from when inaccurate pre-training parameters are used.

We list the main contributions of our work in Section 2. Section 3 discusses
relevant literature. Required theoretical background is provided in Section 4
with the complete methodology in Section 5. Sections 6 and 7 provide
performance analysis and sensitivity study of the model. We summarize our work
in Section 8.

##  2\. Contribution

The main objective of this paper is to present a sophisticated method of
efficiently pre-training an autoencoder (AE) for prediction when sufficient
target measurements are not available for training. A multi-stage setup for
estimating network characteristics from partial distance measurements is
proposed. We integrate the idea of pre-training with that of best parameter
selection for optimal prediction.

These optimally pre-trained autoencoders are demonstrated to predict missing
shortest node-pair distances in social networks. We show the effectiveness of
OSP on real-world social networks such as Virgili Emails, Facebook, and Train
Bombing network which vary in their size and areas of function. The prediction
performance is also compared against an existing distance inference model
based on Low-rank Matrix Completion (Candès and Recht, 2009; Jayasumana et
al., 2019; Mahindre et al., 2019) (Section 6).

We improve with respect to the existing work in the following ways:

  1. (1)

Artificially generated training data is used to pre-train an autoencoder when
sufficient real-world data is not available (Section 5.1).

  2. (2)

A novel two stage oracle is presented to predict suitable parameters for
optimal pre-training. The oracle points us towards the suitable parameters
solely using the sampled distances and without knowing any characteristics of
the target network (Section 5.5).

  3. (3)

OSP does not require large amounts of measurements from the target network
that needs to be predicted. Our experiments show results for networks sampled
as low as 0.1%, 0.5%, and 1% while also covering up to 80% sampling (Section
6).

  4. (4)

Sensitivity of the predictive model is studied towards different pre-training
parameters. This helps us understand pre-training behaviour of social networks
and identify best set of pre-training parameters for the given test network
parameters (Section 7).

##  3\. Related Work

In this section, we discuss the most relevant literature to this paper:
prediction techniques for node-pair distances, completing missing network
measurements from partial measurements, and using pre-training for predicting
distances.

Low-rank Matrix Completion (LMC) (Candès and Recht, 2009) has been used
effectively in reconstructing missing distances in graphs such as IoT or
social networks when partial distance measurements are known, either randomly
or from a selected set of landmarks. It has demonstrated accurate prediction
up to 20-40% of sampled measurements (Jayasumana et al., 2019; Khuller et al.,
1996; Mahindre et al., 2019). This deterministic technique works well due to
the low-rankness of the distance matrices. Thus we use Low-rank Matrix
Completion as one of the reference points for comparing our results.

Deep learning techniques have found its way into social network analysis for
distance prediction. Graph convolutional neural networks have been used to
predict missing links in recommender systems (Berg et al., 2017). Preserved
distances of each node from a few landmark nodes are successfully used (Qi et
al., 2020) to train a multi-layer perceptron (MLP) and predict distance
between two vertices without a high space cost. WiDE, (Liu et al., 2018), uses
unsupervised stacked autoencoders to pre-train a deep neural network and
predicts person-to-person distances using surrounding Wifi signals. However,
these techniques rely on a significant amount of measurements within the real-
world networks.

Our work is inspired by the methods proposed in (Mahindre et al., 2020;
Karkare et al., 2021) where Hadamard Autoencoders are trained on synthetic
data when only sparse measurements are available from the target network.
Though this helps us train on artificial data, all social networks are not the
same and vary greatly in their characteristics. This leaves room for improving
the pre-training process to be customized for each social network being
reconstructed.

Thus, we hypothesize that building a pre-training “Oracle” that guides us into
a narrower range of the pre-training data being used will help significantly
improve the prediction performance by generating customized training data that
is more faithful to the given social network.

##  4\. Theory

In this section, we provide necessary background by giving a brief
introduction of the autoencoders and the data used. We also briefly talk about
the concept of pre-training which forms the foundation of this work.

###  4.1. Supervised Autoencoders

An autoencoder is a special type of artificial neural network that tries to
reproduce the input at the output while learning its representation in a low
dimensional space. It consists of two parts, namely, encoder and decoder. The
encoder takes the input, x, and maps it to a hidden representation h, such
that:

(1) $\textbf{h}=\sigma(\textbf{W}\textbf{x}+\textbf{b}),$

where W is the weight matrix, b is the bias vector and $\sigma$ is the
activation function for the encoder. We use Leaky ReLU as the activation
function throughout the autoencoder. It is defined as:

(2) $\sigma(z)=\begin{cases}z,&\textit{if }z>0\\\ \alpha z,&z<=0,\end{cases}$

where $z$ is the input and $\alpha$ is a small non-zero, constant gradient.
Leaky ReLU helps to solve the “vanishing gradient” problem caused by ReLU
where the function outputs 0 if $z<0$, seizing the update of the
hyperparameters and thus the learning.

The decoder reconstructs this hidden representation h to the output,
$\hat{\textbf{x}}$, given by:

(3)
$\hat{\textbf{x}}=\sigma(\textbf{W}^{\prime}\textbf{h}+\textbf{b}^{\prime}),$

where $\textbf{W}^{\prime}$ is the weight matrix, $\textbf{b}^{\prime}$ is the
bias vector and $\sigma$ is the activation function for the decoder.

The input and output layers have the same number of neurons in order to
reconstruct the given input. Each neuron operates on the incoming signal and
generates an output based on its activation function, weight and bias values.
The weight and bias values are found by training so that the difference
between $\hat{\textbf{x}}$ and x is minimized.

###  4.2. Training Protocols

We use supervised training method to train the hyperparameters, i.e., while
training, the expected output is known. Back-propagation is used to train the
weights and biases for each layer during training. Back-propagation consists
of two phases. During the first phase, input is fed to the autoencoder with
weights and biases initialized to random values and the output layer
reconstructs the input. In the second phase, the error between expected and
predicted output is calculated using a mean squared error loss function given
by Eq. 4. The weights and biases are updated iteratively using this error
value.

(4)
$\mathcal{L}(\textbf{x},\hat{\textbf{x}})=||\textbf{x}-\hat{\textbf{x}}||^{2},$

where x is the known expected output, $\hat{\textbf{x}}$ is the predicted
output as reconstructed by the output layer. Thus from Eq. (1), (3), and (4),
we can say that

(5)
$\mathcal{L}(\textbf{x},\hat{\textbf{x}})=||\textbf{x}-\sigma(\textbf{W}^{\prime}(\sigma(\textbf{W}\textbf{x}+\textbf{b}))+\textbf{b}^{\prime})||^{2}.$

In other words, the loss function is the mean
$(\frac{1}{n}\sum\limits_{i=1}^{n})$ of the squared errors
$({x}_{i}-\hat{x}_{i})^{2}$, which is an easily computable quantity for the
given data sample, and can be expressed as
$MSE=\frac{1}{n}\sum\limits_{i=1}^{n}({x}_{i}-\hat{{x}}_{i})^{2}$, where $n$
is the total number of samples.

The size of the hidden layer can be adjusted in proportion to the size of the
target data. We used learning rate 0.001, batch size 1, and hidden layer size
50 for the Facebook network and 20 for the Virgili Emails and Train Bombing
network, as determined from the observed measurements.

###  4.3. Pre-training

The concept of pre-training is inspired from the way humans learn. While
learning new things, we re-use our knowledge of the past. Pre-training a
neural network refers to training the model during one task and using the
trained parameters as a head start for another task. This way, the model does
not have to learn from scratch for the second task. It is especially useful if
the two tasks are similar or the final task is much more difficult and can be
broken down in simpler pre-training tasks (Radford et al., 2018; Erhan et al.,
2010) . For instance, training an autoencoder to predict missing measurements
of sparsely available data can be difficult. Instead, we pre-train the
autoencoder on readily available simulated data and try to predict the
parameters for this pre-training data so that it is similar in characteristics
to the target data.

Analysis of both natural and human-created real-world networks have shown that
such networks closely follow specific rules, such as Power law behaviour in
their degree distributions $P(k)\sim k^{-\gamma}$, with the exponent varying
between 2 and 3. The characteristic that only a few nodes have an extremely
large number of connections while most of the nodes have very few links is due
to ‘preferential attachment’ (Albert and Barabási, 2002; Udny Yule, 1925;
Barabási et al., 2016). Furthermore, social networks like Facebook and WeChat
are observed to have Power law growth dynamics (Zang et al., 2018). Thus, for
a given real-world network, it is possible to identify a closely related
network type that can be used to generate synthetic data for training. Social
networks that we are interested in for this paper, e.g., Emails, Facebook, and
Train Bombing networks, are well documented to follow Power law.

Training the autoencoder on Power law networks with similar parameters helps
the model treat the target real-world network as just a variation of the
training network. Another advantage is that we can generate sufficient data as
required for proper training and effective performance of the neural network.

###  4.4. Data

Given a graph $G=\left\\{V,E\right\\}$, where $V$ is the set of nodes in the
network and $E$ is the set of undirected edges, we use shortest hop-distances
between node-pairs as our data. For a given graph $G$, its _shortest hop-
distance matrix_ H can be given as

(6) $\textbf{H}=\left[h_{ji}=h_{ij}=\textit{shortest hop-distance from
}i\textit{ to }j\right],$

The shortest _Hop-Distance_ Matrix (HDM), $\textbf{H}\in\mathbb{N}^{N\times
N}_{0}$ can represent the entire network, where $\mathbb{N}_{0}$ denotes the
non-negative integers $\mathbb{N}\cup\\{0\\}$ and $N$ denotes number of nodes
in $G$.

Note that nodes and links may have different interpretation in different
networks, e.g., in the Facebook network, nodes are users and links represent
friendships, where as in a software system network, nodes represent software
modules and links may represent dependency between the two modules. While our
work focuses only on _connected, unweighted, and undirected graphs_ , there
are several commercially important social networks in this domain such as
Friendship, Email, crime, and interaction networks to name a few.

###  4.5. Sampling Scheme

We aim to predict missing distances in networks when the target network is
only partially sampled.

Figure 2. Visual illustrations of P (sampled H) for random sampling scheme
with $N=8$. Known entries are denoted by dark blue and unknown values by light
blue.

In our experiments, we measure a small set of distances between random pairs
of nodes from the distance matrix H to achieve distributed, unbiased sampling.
The randomly sampled H results in an $N\times N$ matrix P, as shown in Fig. 2,
where the sampled entries are depicted by solid squares and ‘?’s represent
missing entries. Other sampling schemes such as random walk or snow-ball
sampling tend to be biased towards higher degree nodes, producing a
misrepresentation of typical large scale, real-world, unstructured networks
with non-uniform degree distributions, e.g., Power law degree distribution
(Awan et al., 2006).

The conventional approach for network reconstruction from a sparse set of
random samples are based on Low-rank Matrix Completion which requires a
sampling scheme to measure at least one non-zero distance value in each row
for the test network (Jayasumana et al., 2019). However, our approach
overcomes this restriction.

###  4.6. Low-rankness of Data

Low-rankness in a dataset implies redundancy which allows reconstruction of
missing measurements from a small set of entries, i.e., partial measurements.
We leverage this while reconstructing distance matrices.

Here, we demonstrate low-rankness of the real-world social networks used. Fig.
3 shows the logarithm of all the singular values of H. Only a small set of
singular values are dominant with the later magnitudes falling down rapidly
towards ‘0’, indicating low-rankness of the H matrices for shown networks.

Figure 3. Low-rankness of hop-distance matrices of social networks used in
this paper.

Note that low-rank data lies on low-dimensional linear subspace whereas low-
dimensional data lies on low-dimensional non-linear subspace.

##  5\. Methodology

An overview of the proposed simple approach for predicting distances from a
small set of measurements is shown in Fig. 5. In this section, we discuss how
pre-training data is artificially generated to train a neural network, which
parameters are tuned to make the synthetic data similar to real-world data we
are predicting on, and how the ideal values of this parameter are estimated
using our proposed pre-training “Oracle”.

Figure 4. Illustration of selecting windows for parameter values for
generating artificial Power law networks. Neural network is trained and tested
to evaluate effectiveness of each window on prediction. Figure 5. Oracle
Search Pre-training (OSP): Model architecture. STAGE 0 helps narrow down
suitable parameter values from a very broad range and STAGE I helps to select
the three most suitable parameters for the given sampled network. The neural
network is pre-trained on artificial data generated using these selected
parameters and fine tuned using sampled distances form the target network in
STAGE II to predict missing distances in the real-world social network.

###  5.1. Generation of Artificial Training Data

We learnt in Section 4.3 that our target networks, i.e., social networks, are
Power law in nature. Thus, we use synthetic Power law networks as a base for
generating our pre-training data.

Power law cluster model, in Python’s Networkx library, creates graphs for a
given network size and edges using preferential attachment with the function
powerlaw_cluster_graph(N, m, p), where a graph of $N$ nodes is grown by
attaching new nodes, each with $m$ edges that are preferentially attached to
existing nodes with high degree (Jin et al., 2001; Hagberg et al., 2008), and
$p$ is the probability of forming a triangle after adding a random edge. The
model generates a network whose degree distribution follows a Power law curve
and thus, represents social networks such as protein-protein interactions,
World-Wide-Web, citations network, Facebook, and others. However, depending on
the given parameter values, the generated graphs could be vastly different
from the specific network under consideration. Even for two networks of the
same type, e.g., social networks, the parameters could be significantly
different. Next we address the method for determining the appropriate
parameters that need to be tuned, corresponding to the specific network of
interest.

###  5.2. Tuning Parameter for Synthetic Data Generation

Node degree is the simplest way of describing how nodes are connected to each
other. Yet, it has a prodigious effect on the network characteristics such as
clustering coefficient, betweenness centrality, diameter of the network, etc.,
proving to be a simple but significant network feature.

As $m$ represents number of edges connected to each new node, in the function
powerlaw_cluster_graph(N, m, p), it is analogous to the degree of the node.
Hence, we focus on finding a suitable range of values for parameter $m$ for
the given target network. We generate various training networks with fixed $N$
\- to generate networks of the size same as that of the target network,
consciously selected values of $m$ \- to study the effect of different $m$
values on prediction performance, while using widely spread values of
probability $p$, i.e., from 0.1 to 0.9 - to generate networks with different
extent of clustering. This allows us to train the autoencoder on different
clustering densities, as we do not know the clustering coefficient of our
target network.

Next, we use window method, as described below, to systematically vary $m$ and
study its effect on prediction performance.

###  5.3. Window Method for Varying the Parameter Value

Our goal is to find a narrow range of values for $m$ to generate training data
so that it comes closer to the characteristics of our target network.

If we assume that the average node degree of the target network is less than
some value $n_{d}$ and that the network has a Power law node degree
distribution, then the distribution is highly right-skewed and we can safely
say that the set of suitable $m$ values, for training data generation, are
more likely to fall on the left side of average node degree as most of the
nodes have node degree less than the average node degree of the network. We
also demonstrate this in Section 7.

As $m>0$ for connected networks, we partition the range 1 to $n_{d}$ in small
windows as shown in Fig. 4, where $n_{d}=19$ for illustration purposes. An
autoencoder is trained on networks generated for the $m$ values in each window
and is tested. This allows us to evaluate effectiveness of each parameter
window with respect to its prediction performance. The goal is to find an
ideal range of parameters that will give us the best prediction for the given
target network.

Note that the assumed value of $n_{d}$ may vary as per the type of the target
network. For example, we know that in a road network, a node (an intersection)
is less likely to have more than 10 edges (roads) so $n_{d}$ can be set to a
lower value. Setting a higher $n_{d}$ has no major drawbacks in the model as
it only produces a few more windows for analysis but underestimating $n_{d}$
might make us miss the ideal set of values of $m$ for the target data.

###  5.4. Evaluation Parameters

We now describe the performance metric used to evaluate prediction performance
of our model. We aim to calculate the cumulative error in prediction. Thus, we
measure and compare mean error and absolute hop-distance error to evaluate
prediction performance and to demonstrate effectiveness of our model. Note
that we measure these error values for the set of node-pairs under test. This
can be a set of observed node-pairs or unobserved node pairs depending on the
stage of the model and will be mentioned accordingly.

####  5.4.1. Mean error

We define mean error $M_{e}$ as follows:

(7) $M_{e}={\left[\sum\limits_{\forall\
i,j\in\check{\Omega}}{|\check{h_{ij}}(f)-h_{ij}|}\right]}\bigg{/}{\left[\sum\limits_{\forall\
i,j\in\check{\Omega}}{h_{ij}}\right]},$

where, $\check{\Omega}$ is a set of all node pairs under test,
$\check{h_{ij}}(f)$ refers to the predicted shortest hop-distance from node
$i$ to node $j$ when $f$ percentage of random measurements are sampled and
$h_{ij}$ refers to its original expected value. Mean error gives a percentage
value of the prediction error with respect to sum of the original hop-
distances in the network.

####  5.4.2. Absolute hop-distance error

Absolute hop-distance error (AHDE) is defined as

(8) $H_{e}=\bigg{\\{}\sum\limits_{\forall\
i,j\in\check{\Omega}}|\check{h_{ij}}(f)-h_{ij}|\bigg{\\}}\bigg{/}\\{card(\check{h})\\},$

where $card(\check{h})$ gives the total number of node-pair entries in H under
test. AHDE captures average deviation in prediction in the magnitude of hop-
distances. For instance, the absolute hop-distance error of 1 implies that on
average a predicted path length will be off by 1 hop, i.e., for the original
hop-distance of 10, predicted value will be somewhere from 9 to 11.

###  5.5. Model

A visual illustration of the model architecture can be seen in Fig. 5. A
detailed explanation of various stages of the proposed model follows where we
explain the ideal parameter selection (STAGE 0 and STAGE I), systematic pre-
training and final prediction of missing node-pair distances (STAGE II).

####  5.5.1. STAGE 0-Broad range parameter selection

In initial stages of building the Oracle, the autoencoder is evaluated for a
range of $m$ values as we are not aware of the ideal values for the tunable
parameter. These values are chosen using a sliding window to include a wide
range (Section 5.3 and Fig. 4). As we saw in Section 5.3, a parameter range
can be assumed for $m$ to start. If the range is broad, it can be divided into
smaller windows. This can be viewed as STAGE 0.

For example, being a communication network, Virgili Emails network might have
a high average node degree. So we first test the autoencoder performance on a
set of windows $\mathcal{W}_{0}$ of $m$ values of [1,5,10], [10,15,20], and so
on up to [90,95,100] (as shown in Figs. 6[a]). The best performing window of
[1,5,10], say $\mathcal{W}_{0\\_k}$, is selected from this set for the next
stage, i.e., STAGE I.

Note that STAGE 0 is optional and can be used when a quick narrow down of a
possible wide range of $m$ values is required.

####  5.5.2. STAGE I-Narrow range parameter selection

The best performing window from STAGE 0 is further split into a narrower set
of windows, $\mathcal{W}_{I}$ and $\mathcal{W}_{I}\subset\mathcal{W}_{0}$ (as
shown in Fig. 5), to tune the parameter values further. The autoencoder is
trained for each of these windows, $\mathcal{W}_{I\\_k}$, independently and is
tested after each session on observed distances from the target network.

For example, we selected the range [1,5,10] from STAGE 0 for the Virgili
Emails network, thus, in STAGE I, we have trained the autoencoder on narrower
windows of $m$ values as shown in Fig. 6[b] from [1,2,3], [3,4,5], and so on
up to [9,10,11].

(a) STAGE 0 evaluation

(b) STAGE I evaluation

Figure 6. Parameter selection for pre-training in the Oracle for the Virgili
Emails network: (a) STAGE 0: Each broad range window is evaluated for its
prediction effectiveness on sampled distances. The best parameter range is
selected for the next stage. Here, parameter range [1,5,10] is seen to perform
the best. (b) STAGE I: The braod range of (1,5,10) selected from STAGE 0 is
split into narrower windows and each window is evaluated.

As we are evaluating the system over partial and presumably a very small
quantity of sampled data, and because no single window performed significantly
better than the rest in STAGE I, we select the top three parameter windows
that performed well for the sampled measurements. You can see from Figs. 7 and
8 that windows [5,6,7], [7,8,9], and [9,10,11] perform the best for the
Virgili Emails and Facebook network. The Train Bombing network being small in
size, we can safely pick a narrower range and train the final model on the
best performing window from STAGE 1.

####  5.5.3. STAGE II-Final Prediction

The best parameter values selected from STAGE I are eventually used to
generate pre-training data. Once trained, this model can be fine tuned on
observed distances from the target network. We examine the results in both
cases (with and without fine tuning) in Section 6. The trained autoencoder
then reconstructs the complete distance matrix from only a fraction of sampled
distances from the test network.

Note that STAGE 0 and STAGE I act as an “Oracle” to search the ideal pre-
training parameters corresponding to the target network while STAGE II pre-
trains AE on artificial data generated with these parameters.

###  5.6. Implementation Details

We used TensorFlow, Google’s Machine Learning Library, to build the
autoencoder and wrote the script in Python 3.6. We ran the experiments on
Google Colab.

Data used and the code for our model are available on Github:
https://github.com/anonymous/anonymous

##  6\. Performance evaluation of the proposed method

We measure mean error (Section 5.4.1) and absolute hop-distance error (Section
5.4.2) to evaluate the prediction performance of our model and compare it with
the Low-rank Matrix Completion (Jayasumana et al., 2019) based approach. Here,
the errors are measured only over the unobserved node-pair distances. We know
that all diagonal values in a distance matrix are ‘0’ and no off-diagonal
entries are ‘0’, so we set all diagonal values of $\hat{\textbf{P}}$ to 0 and
round-up all off-diagonal values between 0 and 1 to 1.

(a) Mean error for prediction in Virgili Emails Network

(b) AHDE for prediction in Virgili Emails Network

Figure 7. Virgili Emails network: Prediction performance. The plots show that
pre-trained neural networks (yellow, pink and blue) outperform prediction
especially when lower percentages of entries are sampled from the real-world
network. Graceful degradation is also observed for pre-training models. Fine
tuning along with pre-training predicts with higher accuracy than a neural
network trained only on sampled entries.

(a) Mean error for prediction in Facebook Network

(b) AHDE for prediction in Facebook Network

Figure 8. Facebook network: Prediction performance. The plots show that pre-
trained neural networks (yellow, pink and blue) outperform prediction
especially when lower percentages of entries are sampled from the real-world
network. Graceful degradation is also observed for pre-training models. Fine
tuning along with pre-training predicts with higher accuracy than a neural
network trained only on sampled entries.

(a) Mean error for prediction in Train Bombing Network

(b) AHDE for prediction in Train Bombing Network

Figure 9. Train Bombing network: Prediction performance. The plots show that
pre-trained neural networks (yellow, pink and blue) outperform prediction
especially when lower percentages of entries are sampled from the real-world
network. Graceful degradation is also observed for pre-training models. Fine
tuning along with pre-training predicts with higher accuracy than a neural
network trained only on sampled entries.

The model has been evaluated for different variations and are labeled as shown
in Figs. 7, 8, and 9. These labels are explained below:

  * •

Trivial 0: all missing values are replaced by 0

  * •

Trivial 1: all missing values are replaced by 1
Trivial 0 and Trivial 1 cases provide a reference to the worse case model
performance.

  * •

MC: prediction using Low-rank Matrix Completion

  * •

Trained on observed: the AE has been trained only on the sampled node-pair
distances from the real-world network and no pre-training has been done at all

  * •

[x,y,z] - Powerlaw q: the STAGE II AE, i.e., trained on artificial Power law
networks with parameter $m$ values x,y, and z, was tested on a same size Power
law network artificially generated with $m=q$

  * •

[x,y,z] - Network: the STAGE II AE, i.e., trained on artificial Power law
networks with parameter $m$ values x,y, and z, was tested on the real-world
social network

  * •

[x,y,z] - Network, fine tuned: the STAGE II AE, i.e., trained on artificial
Power law networks with parameter $m$ values x,y, and z, was fine tuned with
the observed entries of the real-world network and then tested on the
unobserved entries of the same real-world network

Now let us take a look at the performance of the model for three real-world
social networks.

Virgili Emails network: This is a 1133 nodes network with average node degree
of 9.6.

Facebook network: This is a 4039 node network with average node degree of
43.6.

Train Bombing network: This is a 64 node network with average node degree of
7.5.

From STAGE 0 and STAGE I of the Virgili Emails network, as shown in Fig. 6,
three best windows were selected, [5,6,7], [7,8,9], and [9,10,11]. In STAGE
II, the AE is thus trained on [5,7,9] to cover these three windows. This
behavior is also seen in the Facebook network. The STAGE 0 and STAGE I results
for the Facebook and Train Bombing networks are not shown here due to space
restrictions. However, can be found at Github:
https://github.com/anonymous/anonymous

Note that the parameter values selected by the “Oracle” are not the same as
the average node degree of the test network but tend to be on the lower side
of the average node degree value. We will explore this relation and the
sensitivity of the model towards parameter value $m$ in Section 7.

The STAGE II (Virgili Emails: Fig. 7, Facebook: Fig. 8, Train Bombing: Fig. 9)
show that AE performs far better than MC for lower percentages of sampled
measurements. We also see that AE trained on artificial data performs better
than AE trained only on the observed entries as the neural network has much
more data to train on from synthetic Power law networks. Note that the network
performs best for a variation of Power law network (see label
“[x,y,z]-Powerlaw”) which is obvious as the network was trained on very
similar data. The error, when tested on the real-world network is slightly
higher than for “[x,y,z]-Powerlaw” as the real-world network is different than
the synthetic Power law networks. However, we notice that fine tuning the AE
almost always bridges this gap. Especially when we have a higher percentage of
sampled node-pair measurements from the network, fine tuning gives lower error
than mere Power law trained AE for all three networks.

It is worthwhile to say that OSP performs well not only for large networks but
also for smaller networks such as the Train Bombing network here, showing
independence of the model performance with respect to the network size.

##  7\. Discussion on optimum parameter selection for pre-training

This section furthers the intuition behind parameter selection for optimal
pre-training when the average node degree of the test network is (tentatively)
known.

(a) Prediction performance with various window sizes.

(b) Prediction performance with various window sizes: Visual comparison.

Figure 10. Sensitivity analysis: The plots show that our model predicts more
accurately when pre-trained on a more suitable and narrow range of pre-
training parameters. More so, the accuracy is higher for parameters on the
lower side of the average node degree as Power law networks have a right
skewed node degree distribution with a high number of nodes with low node
degree.

We saw that the parameter $m$ is analogous to node degree as it represents
number of edges attached to the new node while the artificial network is being
constructed. We conduct a sensitivity analysis on the prediction performance
towards parameter range. This helps us understand the model behavior and
select best training parameters for the given average node degree of the test
network.

Mean error over the unobserved values of the Virgili Emails network, when
autoencoder is pre-trained over various parameter windows, is shown in Fig.
10[a]. A more comparative visual display of the results can be seen in Fig.
10[b]. We can see that for an average node degree of 9.6, the autoencoder
performs better when trained on narrower windows as compared to wider. Thus,
the performance improves from [2,10,18] to [9,10,11]. We also observed that a
narrow window placed on the lower side of the target node degree gives a
slightly better performance than window placed on the higher side of the
spectrum, i.e., [7,8,9] performs better than [10,11,12], even if both the
windows are equally narrow. This is because our target network is Power law in
nature and a lot of nodes have degrees on the lower side of the average node
degree.

##  8\. Conclusion

In this paper, we have shown how to optimally pre-train a neural network when
only sparse samples are available from the target network. We call this model
“Oracle Search Pre-training” (OSP) as it helps us search the optimal pre-
training parameter values specifically for the given social network. We use
artificially generated data to compensate for the scarcity of real-world
measurements to effectively train a neural network. The low-rankness of the
distance matrices is leveraged.

Prediction performance is evaluated on three real-world networks, namely,
Virgili Emails, Facebook, and Train Bombing network, under different
variations of the model, and is also compared with a state-of-the-art Low-rank
Matrix Completion approach. The model infers distances within one hop of the
original value even when only 1% of measurements are sampled. We also study
sensitivity of the model towards the tuning parameter values and its relation
with average node degree of the social network. It is worth mentioning that
though we illustrate results on social networks, the method can be generalized
to other network domains by adapting the synthetic data.

In the future, we would like to test our model for directed networks and
networks from other domains such as crime networks. The effect of other
network parameters, such as clustering coefficient, on reconstruction can also
be studied.

## References

  * (1)
  * Albert and Barabási (2002) Réka Albert and Albert-László Barabási. 2002\.  Statistical mechanics of complex networks.  _Reviews of modern physics_ 74, 1 (2002), 47\.
  * Awan et al. (2006) Asad Awan, Ronaldo A Ferreira, Suresh Jagannathan, and Ananth Grama. 2006.  Distributed uniform sampling in unstructured peer-to-peer networks. In _Proceedings of the 39th Annual Hawaii International Conference on System Sciences (HICSS’06)_ , Vol. 9. IEEE, 223c–223c.  https://doi.org/10.1109/HICSS.2006.126
  * Barabási et al. (2016) Albert-László Barabási et al. 2016\.  _Network science_.  Cambridge University Press.
  * Berg et al. (2017) Rianne van den Berg, Thomas N Kipf, and Max Welling. 2017\.  Graph convolutional matrix completion.  _arXiv preprint arXiv:1706.02263_ (2017).
  * Boguná et al. (2004) Marián Boguná, Romualdo Pastor-Satorras, Albert Díaz-Guilera, and Alex Arenas. 2004\.  Models of social networks based on social distance attachment.  _Physical review E_ 70, 5 (2004), 056122.  https://doi.org/10.1103/PhysRevE.70.056122
  * Candès and Recht (2009) Emmanuel J Candès and Benjamin Recht. 2009.  Exact matrix completion via convex optimization.  _Foundations of Computational mathematics_ 9, 6 (2009), 717\.  https://doi.org/10.1007/s10208-009-9045-5
  * Dhanapala and Jayasumana (2013) Dulanjalie C Dhanapala and Anura P Jayasumana. 2013.  Topology preserving maps—Extracting layout maps of wireless sensor networks from virtual coordinates.  _IEEE/ACM Transactions on Networking_ 22, 3 (2013), 784–797.  https://doi.org/10.1109/TNET.2013.2263254
  * Erhan et al. (2010) Dumitru Erhan, Aaron Courville, Yoshua Bengio, and Pascal Vincent. 2010.  Why does unsupervised pre-training help deep learning?. In _Proceedings of the thirteenth international conference on artificial intelligence and statistics_. JMLR Workshop and Conference Proceedings, 201–208.
  * Hagberg et al. (2008) Aric A. Hagberg, Daniel A. Schult, and Pieter J. Swart. 2008\.  Exploring Network Structure, Dynamics, and Function using NetworkX. In _Proceedings of the 7th Python in Science Conference_ , Gaël Varoquaux, Travis Vaught, and Jarrod Millman (Eds.). Pasadena, CA USA, 11 – 15\.
  * Jayasumana et al. (2019) Anura P Jayasumana, Randy Paffenroth, Gunjan Mahindre, Sridhar Ramasamy, and Kelum Gajamannage. 2019\.  Network topology mapping from partial virtual coordinates and graph geodesics.  _IEEE/ACM Transactions on Networking_ 27, 6 (2019), 2405–2417.  https://doi.org/10.1109/TNET.2019.2953921
  * Jin et al. (2001) Emily M Jin, Michelle Girvan, and Mark EJ Newman. 2001\.  Structure of growing social networks.  _Physical review E_ 64, 4 (2001), 046132.  https://doi.org/10.1103/PhysRevE.64.046132
  * Karkare et al. (2021) Rasika Karkare, Randy Paffenroth, and Gunjan Mahindre. 2021\.  Blind Image Denoising and Inpainting Using Robust Hadamard Autoencoders.  arXiv:2101.10876 [eess.IV]
  * Khuller et al. (1996) Samir Khuller, Balaji Raghavachari, and Azriel Rosenfeld. 1996\.  Landmarks in graphs.  _Discrete applied mathematics_ 70, 3 (1996), 217–229.  https://doi.org/10.1016/0166-218X(95)00106-2
  * Kossinets (2006) Gueorgi Kossinets. 2006\.  Effects of missing data in social networks.  _Social networks_ 28, 3 (2006), 247–268.
  * Liu et al. (2018) Wenping Liu, Yufu Jia, Guoyin Jiang, Hongbo Jiang, Fan Wu, and Zhicheng Lv. 2018\.  WiFi-Sensing Based Person-to-Person Distance Estimation Using Deep Learning. In _2018 IEEE 24th International Conference on Parallel and Distributed Systems (ICPADS)_. IEEE, 236–243.  https://doi.org/10.1109/PADSW.2018.8644549
  * Mahindre et al. (2019) Gunjan Mahindre, Anura P Jayasumana, Kelum Gajamannage, and Randy Paffenroth. 2019.  On sampling and recovery of topology of directed social networks–a low-rank matrix completion based approach. In _2019 IEEE 44th Conference on Local Computer Networks (LCN)_. IEEE, 324–331.  https://doi.org/10.1109/LCN44214.2019.8990707
  * Mahindre et al. (2020) Gunjan Mahindre, Rasika Karkare, Randy Paffenroth, and Anura Jayasumana. 2020.  Inference in Social Networks from Ultra-Sparse Distance Measurements via Pretrained Hadamard Autoencoders. In _2020 IEEE 45th Conference on Local Computer Networks (LCN)_. IEEE, 256–266.  https://doi.org/10.1109/LCN48667.2020.9314769
  * Qi et al. (2020) Jianzhong Qi, Wei Wang, Rui Zhang, and Zhuowei Zhao. 2020\.  A learning based approach to predict shortest-path distances.  (2020).  https://doi.org/10.5441/002/edbt.2020.34
  * Radford et al. (2018) Alec Radford, Karthik Narasimhan, Tim Salimans, and Ilya Sutskever. 2018.  Improving language understanding by generative pre-training.
  * Sarkar and Moore (2005) Purnamrita Sarkar and Andrew W Moore. 2005.  Dynamic social network analysis using latent space models.  _ACM SIGKDD Explorations Newsletter_ 7, 2 (2005), 31–40.  https://doi.org/10.1145/1117454.1117459
  * Sharma et al. (2015) Rajesh Sharma, Matteo Magnani, and Danilo Montesi. 2015\.  Investigating the types and effects of missing data in multilayer networks. In _2015 IEEE/ACM International Conference on Advances in Social Networks Analysis and Mining (ASONAM)_. IEEE, 392–399.  https://doi.org/10.1007/s13278-016-0384-3
  * Udny Yule (1925) G Udny Yule. 1925\.  A mathematical theory of evolution, based on the conclusions of Dr. JC Willis, FRS.  _RSPTB_ 213 (1925), 21–87.  https://doi.org/10.1098/rstb.1925.0002
  * Wahid-Ul-Ashraf et al. (2019) Akanda Wahid-Ul-Ashraf, Marcin Budka, and Katarzyna Musial. 2019\.  How to predict social relationships—Physics-inspired approach to link prediction.  _Physica A: Statistical Mechanics and its Applications_ 523 (2019), 1110–1129.  https://doi.org/10.1016/j.physa.2019.04.246
  * Zang et al. (2018) Chengxi Zang, Peng Cui, Christos Faloutsos, and Wenwu Zhu. 2018\.  On Power Law Growth of Social Networks.  _IEEE Transactions on Knowledge & Data Engineering_ 30, 09 (sep 2018), 1727–1740.  https://doi.org/10.1109/TKDE.2018.2801844
  * Zhang and Pang (2015) Yang Zhang and Jun Pang. 2015.  Distance and friendship: A distance-based model for link prediction in social networks. In _Asia-Pacific Web Conference_. Springer, 55–66.


11institutetext:  Colorado School of Mines, Department of Physics, Golden, CO,
USA 22institutetext: Bergische Universität Wuppertal, Department of Physics,
Wuppertal, Germany 33institutetext: Karlsruhe Institute of Technology (KIT),
Institute for Astroparticle Physics, Karlsruhe, Germany 44institutetext:
Observatorio Pierre Auger, Av. San Martín Norte 304, 5613, Malargüe,
Argentina.
Full author list: https://www.auger.org/archive/authors_2022_10.html

# Update on the indication of a mass-dependent anisotropy above 1018.7 eV in
the hybrid data of the Pierre Auger Observatory

Eric Mayotte 1122<EMAIL_ADDRESS>Thomas Fitoussi for the Pierre Auger
Collaboration 3344<EMAIL_ADDRESS>
###### Abstract

We test for an anisotropy in the mass of arriving cosmic-ray primaries
associated with the galactic plane. The sensitivity to primary mass is
obtained through the depth of shower maximum, $X_{\text{max}}$, extracted from
hybrid events measured over a 14-year period at the Pierre Auger Observatory.
The sky is split into distinct on- and off-plane regions using the galactic
latitude of each arriving cosmic ray to form two distributions of
$X_{\text{max}}$, which are compared using an Anderson-Darling 2-samples test.
A scan over roughly half of the data is used to select a lower threshold
energy of $10^{18.7}$ eV and a galactic latitude splitting at
$|b|=30^{\circ}$, which are set as a prescription for the remaining data. With
these thresholds, the distribution of $X_{\text{max}}$ from the on-plane
region is found to have a $9.1\pm 1.6^{+2.1}_{-2.2}$ g/cm2 shallower mean and
a $5.9\pm 2.1^{+3.5}_{-2.5}$ g/cm2 narrower width than that of the off-plane
region and is observed in all telescope sites independently. These differences
indicate that the mean mass of primary particles arriving from the on-plane
region is greater than that of those from the off-plane region. Monte Carlo
studies yield a $5.9\times 10^{-6}$ random chance probability for the result
in the independent data, lowering to a $6.0\times 10^{-7}$ post-penalization
random chance probability when the scanned data is included. Accounting for
systematic uncertainties leads to an indication for anisotropy in mass
composition above $10^{18.7}$ eV with a $3.3\,\sigma$ significance.
Furthermore, the result has been newly tested using additional FD data
recovered from the selection process. This test independently disfavors the
on- and off-plane regions being uniform in composition at the $2.2\,\sigma$
level, which is in good agreement with the expected sensitivity of the dataset
used for this test.

##  1 Introduction

The Pierre Auger Observatory is currently the largest observatory dedicated to
studying cosmic rays with energies in the EeV range, so-called ultra-high-
energy cosmic rays, UHECR PierreAuger:2015eyc . To do this, it uses both an
array of particle detectors on the Earth’s surface, the Surface Detector, SD
PierreAuger:2007kus , and an array of fluorescence telescopes monitoring the
atmosphere above the SD, the Fluorescence Detector, FD PierreAuger:2009esk .
The highest quality data set of the Observatory is made up of UHECR events
which have simultaneously been measured by both the FD and SD, so-called
hybrid events. In a hybrid event reconstruction, the geometry of the shower
axis is highly constrained by combining the triggered pixel geometry/timing
from the FDs and the high confidence core location/timing provided by the SDs.
This results in an angular resolution for the pointing direction of the shower
axis of better than $0.5^{\circ}$ Bonifazi:2009ma , and a resolution on the
location of the shower core of $50$ m Mostafa:2006id .

The low uncertainty geometric reconstruction provided by the hybrid method
allows the evolution of the intensity of UV fluorescence light measured by the
FD to be inverted to model the shower profile, which is the number of charged
particles in the air shower as a function of the amount of matter it has
traversed, the slant depth, $X$. From the shower profile, the slant depth at
which the maximum development of the shower occurs, $X_{\text{max}}$, can be
extracted. $X_{\text{max}}$ is closely related to the mass of the primary
cosmic ray which induced the air-shower, but is subject to large fluctuations
meaning that it can not be used on a shower-by-shower basis to determine
primary mass. However, if collected with sufficient statistics, the first and
second moments of distributions of $X_{\text{max}}$, $\langle
X_{\text{max}}\rangle$ and $\sigma(X_{\text{max}})$ respectively, can be used
to make high certainty estimations of the mean mass of the UHECR events used
to form that particular $X_{\text{max}}$ distribution Aab:2014kda .

Up until recently PierreAuger:2021jlg , $X_{\text{max}}$ derived from hybrid
measurements has been used to study the average composition of the cosmic-ray
sky as whole, rather than being used to compare the mean compositions of
different subsets of the sky. This choice was likely driven by the relatively
sparse statistics available in hybrid studies due to the upper-limit of
exposure available to them being set by the relatively low 14 % up-time of
FDs. The possibility of splitting the data set into subsamples was therefore
limited by the need to maintain sufficient statistics to say something useful
about primary composition. However, the Pierre Auger Observatory has now
collected more than 14-years of FD data, and tens of thousands of high-quality
hybrid measurements. With this quantity of data, the sky can be split into
different regions and the mean mass of UHECRs arriving from them can be
studied. This new reality then prompts two questions:

  1. 1.

Are the systematic uncertainties which trend with event arrival direction for
hybrid reconstruction low enough to allow the mean mass arriving from
different regions of the sky to be meaningfully compared?

  2. 2.

Is there a reason to expect that different regions of the sky may display
differing compositions due to astrophysical causes?

Question 1) will be explored in section 4. For question 2) it is clear that
the opportunity exists. The flux above the ankle at ${\sim}5$ EeV
PierreAuger:2020kuy  is mixed in composition and has long been thought (now
confirmed) to be extragalactic in origin Linsley:1963bk . Furthermore, it
definitively displays anisotropy above 8 EeV Aab:2017tyv . Additionally, as
was nicely put by Alan Watson in 1990:

  * “… the Larmor radius of a proton of $10^{18}$ eV in a 3 $\mu$G field is about 400 pc, comparable to the thickness of the galactic disk. It follows, therefore, that, if the bulk of cosmic rays are protons, anisotropies associated with the magnetic field structure of the galactic disk might appear as the energy increases.” \- Watson:1990fj

Indeed, there were hints of such a spectral feature starting somewhere around
$10^{18.5}$ eV for directions within $30^{\circ}$ of the galactic plane
Szabelski:1986rx ; Watson:1990fj . Unfortunately, such an excess so far does
not appear to be significant in the data of current experiments and therefore
has not been given much attention since those initial publications.

Now it is known that the flux above 1 EeV is best described as an evolving mix
of light-, intermediate-, and high-mass primaries pierre2014aab ;
PierreAuger:2021mmt . Due to the galactic magnetic field, GMF, the different
mass components present at any given energy will be deflected to different
degrees as they travel from their extragalactic sources to Earth. This mass
dependent deflection suggests that an anisotropy associated with the structure
of the GMF would kick in for increasingly heavier components as energies
climb. It is therefore distinctly possible that an anisotropy associated with
the galactic plane could arise in the higher mass components of the flux at
some energy in the EeV range.

What follows below are specific tests for such a mass-dependent anisotropy
associated with the galactic plane using hybrid data of the Pierre Auger
Observatory collected between 2004 and the end of 2018. To avoid repeating the
contents of the ICRC 2021 proceedings on this result PierreAuger:2021jlg  in
its entirety, the contents of this proceeding will aim to include components
of the analysis that could not be fit in the eight pages allotted in that
publication. Therefore, while this proceeding will cover the details of the
analysis, increased space will be given to the cross checks and studies of the
systematics uncertainties of the analysis. Additionally, a new test on an
independent hybrid data set, recovered from the quality cuts, will also be
discussed.

##  2 Reconstruction and selection

The reconstruction methods, selection cuts, and core analysis of
$X_{\text{max}}$ distributions used here are largely the same as those
described in Yushkov:2020nhr . A rigorous description of these methods can be
found in Aab:2014kda . Other than the fiducial field-of-view cut, FidFoV,
which is treated below, a detailed description of the methods will not be
provided here.

The important differences from Yushkov:2020nhr  are that the minimum energy
for inclusion in the data set has been raised to $E>10^{18.4}$ eV, and that,
as described below in subsection 2.1, the event acceptance, reconstruction
bias, and $X_{\text{max}}$ resolution are now treated based on the arrival
direction of each event. The lower limit of $10^{18.4}$ eV has been chosen as
above this energy, the composition is well mixed and expected to be primarily
of extragalactic origin Aab:2016zth . The period over which data has been
collected has also been slightly expanded to span events observed between
2004-12-01 and 2018-12-31, yielding 7572 high-quality events. A further
subdivision of the data is necessary to test for the hypothesized anisotropy.
Following the results of the scan shown in Figure 5, the data is split into
the events with a galactic latitude,
$\mathinner{\\!\left\lvert\,b\,\right\rvert}\leq 30^{\circ}$, the on-plane
sample, and, $\mathinner{\\!\left\lvert\,b\,\right\rvert}>30^{\circ}$, the
off-plane sample.

#### Fiducial field-of-view selection

To ensure a good reconstruction of $X_{\text{max}}$ with the FD,
$X_{\text{max}}$ itself should be directly observed. The telescopes of the FD
have a field-of-view, FoV, which is vertically constrained. Additionally, some
showers in more vertical geometries will not reach $X_{\text{max}}$ before
impacting the ground. These factors together lead to a geometric and
$X_{\text{max}}$ dependence for what events end up in the analyzed data set.
If unaccounted for, this $X_{\text{max}}$ acceptance, will inevitably bias a
composition study based on FD $X_{\text{max}}$ data. At the Observatory, this
$X_{\text{max}}$ acceptance is addressed primarily through mitigation of the
effect with the fiducial field-of-view, FidFoV, cut. The FidFoV cut constrains
the FD detector volume to only event geometries where the expected range of
$X_{\text{max}}$ values would be visible in the FD FoV. As can be seen in
Figure 1, this changes the natural FD $X_{\text{max}}$ acceptance (gray) to
one which is unbiased from ${\sim}600$ to ${\sim}900$ g/cm2, which spans the
typical range of observed event $X_{\text{max}}$ values. Remaining effects on
rare events outside this range are corrected using parameterizations as a
function of $X_{\text{max}}$ and primary energy.

Figure 1: FD $X_{\text{max}}$ acceptance before (gray) and after (black)
FidFoV cuts. A 4-variable parameterization of the post-FidFoV acceptance is
shown in red. The range of $X_{\text{max}}$ values with unbiased sampling is
shown in blue.

###  2.1 Distributions of $\mathbf{X_{\text{max}}}$ and arrival direction

After measurement, reconstruction, and selection, the observed
$X_{\text{max}}$ distribution does not quite represent the true
$X_{\text{max}}$ distribution of all cosmic rays landing within the
Observatory. This is due to energy dependent biases on the reconstruction of
$X_{\text{max}}$ ($B$), the resolution on $X_{\text{max}}$ of the hybrid
reconstruction method ($R$), and the residual effects of the $X_{\text{max}}$
acceptance ($A$). Since the location of an event and its inclination with
respect to the observing fluorescence telescope plays a role in the magnitude
of $A$, $R$, and $B$, they have an inherent geometric dependence. If
unaccounted for, this could bias this study. Distributions of the key
geometric relationships between events and the FD can be seen in Figure 2 for
on- and off-plane regions. It is clear that there is little difference between
the two regions, so $A$, $R$, and $B$ are expected to also be similar.

Figure 2: Geometries of the on- and off-plane samples.

To explicitly verify $A$, $R$, and $B$ similarity on and off the galactic
plane, CONEX Bergmann:2006yz  is used to generate showers with Sibyll-2.3c
Riehn:2015oba . The showers are then isotropically thrown into detector
simulations which include the time-dependent state of the FD and SD from 2004
through 2018. These simulations, therefore, mimic the measurement conditions,
trigger efficiency, and up-time of the real data, accurately modelling the
exposure and geometries of events arriving from all parts of the sky
Abreu:2010aa . Two sets of these simulations are produced, one formed from
showers generated with a flat sampling of $X_{\text{max}}$ between 300 and
1500 g/cm2, the flat-MC, and one formed from an equal number of proton,
helium, nitrogen, and iron simulations which are then weighted to their
abundances observed in data as reported in Bellido:2017cgf , the mixed-MC.
These simulated event sets are then subjected to the same reconstruction and
selection techniques used on the real data, so that $A$, $R$, and $B$ are
accurately included in them.

The flat-MC is split into on- and off-plane subsamples which are then used to
extract the functional form of $A$ in 0.1 $\log_{10}(E/\text{eV})$ energy bins
using the method illustrated in Figure 1. The form of $A$ is extracted by
leveraging the flatly sampled $X_{\text{max}}$ generation of the Monte Carlo,
as, once the plateau of the distribution is normalized to one, the height of
each bin represents the acceptance of events with $X_{\text{max}}$ values in
that range. The acceptance is then fit with the 4-component parameterization
illustrated in Figure 1. The energy evolution of $x_{1},\lambda_{1},x_{2}$,
and $\lambda_{2}$ is then parameterized separately with 2D polynomials for the
on- and off-plane subsamples, resulting in Figure 3. From Figure 3, it is
clear that above $10^{18.4}$ eV there is no statistically significant
difference in $X_{\text{max}}$ acceptance between the on- and off-plane
regions. Even so, the region-specific parameterizations of $A$ are used to
correct the 1.4 % of events with partial $X_{\text{max}}$ acceptance using the
up-weighting technique outlined in Aab:2014kda . Uncertainties in these
parameterizations result in the systematic uncertainties on the first and
second moments specified in Table 1.

Figure 3: The $X_{\text{max}}$ acceptance parameterizations for the on- and
off-plane sky regions from Monte Carlo. Figure 4: The $X_{\text{max}}$
reconstruction bias and resolution parameterizations for the on- and off-plane
sky regions from Monte Carlo. Note: only the detectors and reconstruction
$X_{\text{max}}$ resolution is shown. Other effects lowering the resolution
are included as specified in Aab:2014kda .

The mixed-MC is likewise split into on- and off-plane subsamples which are
then used to extract the energy evolution of $B$ and $R$ for each region. $B$
and $R$ are extracted by forming a distribution of the difference between the
FD reconstructed value of $X_{\text{max}}$ and the Monte Carlo truth value in
0.1 $\log_{10}(E/\text{eV})$ energy bins. From these, the mean reconstruction
bias, $B=\langle X_{\rm max}^{\rm FD}-X_{\rm max}^{\rm MC}\rangle$, and the
$X_{\text{max}}$ resolution, $R=\sigma\left(X_{\rm max}^{\rm FD}-X_{\rm
max}^{\rm MC}\right)$, are extracted in each energy bin. Again, the evolution
of each is parameterized with a 2D polynomial. Figure 4 shows that $B$ and $R$
for the two regions are found to agree within errors. These however are also
corrected for separately. Fit uncertainties again result in the systematic
uncertainties outlined in Table 1. At this point, the $X_{\text{max}}$
distributions and moments from the on- and off-plane regions in each energy
bin can be compared without bias from selection and reconstruction.

##  3 Testing for anisotropy

The specific hypothesis to be tested is whether, above some energy threshold,
$E_{\rm th}$, the mean composition of UHECRs coming from directions near to
the galactic plane is significantly higher in mass than those arriving further
from it. This is to be tested using $X_{\text{max}}$ as a mass sensitive
parameter. Typically, $X_{\text{max}}$ based composition analyses leverage the
first two moments of $X_{\text{max}}$ distributions binned in energy, to
comment on primary mass. This approach, however, does not lend itself well to
quantifying the significance of a result testing the above statement. Instead,
a test statistic, $TS$, which quantifies the degree of dissimilarity between
the $X_{\text{max}}$ distributions in the two regions in a single value is
preferred. For this, the returned value from the Anderson-Darling two-sample
homogeneity test andersondarling , AD-test, has been selected as it scales
with the dissimilarity of the tested distributions. The AD-test has good
sensitivity to the full width of a distribution scholz1987k , and has more
power than the Kolmorogov-Smirnov test while remaining robust against false
positives engmann2011comparing .

To use the AD-test and $X_{\text{max}}$ for this purpose, two modifications
are required. First, a single $TS$ comparing all events in each region above
$E_{\rm th}$ is desired. So, all events with $E\geq E_{\rm th}$ in the on- and
off-plane samples separately need to be collected into a common on-plane
distribution and a common off-plane distribution. To do this, the natural
evolution of $X_{\text{max}}$ with energy needs to be removed so that spectral
features in the flux do not influence the result. Therefore, we define an
energy-normalized $X_{\text{max}}$ value

$X_{\text{max}}^{{}^{\prime}}=X_{\text{max}}-\underbrace{\left(649+63.1\,Z_{18}+1.97\,Z_{18}^{2}\right)}_{\text{EPOS-
LHC elongation rate for iron}},$ (1)

where $Z_{18}=\log_{10}\left(E_{\text{rec}}/\,\text{EeV}\right)$. The last
term in Equation 1 is the natural energy evolution of mean $X_{\text{max}}$
for iron primaries as predicted by EPOS-LHC Pierog:2013ria ***Choice of
hadronic interaction model varies result by $\sim 0.02$ g/cm2.. Second, the
$X_{\text{max}}^{\prime}$ distribution of an on-plane sample populated with
primaries which are on average heavier than those in the off-plane sample will
display a lower mean and a narrower width than that of the off-plane
$X_{\text{max}}^{\prime}$ distribution. Since the null hypothesis is that
there is either no composition difference or a heavier off-plane sample, a
$TS$ sensitive to the ordering of the $X_{\text{max}}^{\prime}$ distributions
is required†††Modifying the test to also require
$\sigma(X_{\text{max}}^{\prime})^{\rm on}<\sigma(X_{\text{max}}^{\prime})^{\rm
off}$ would be more restrictive, but conservatively has not been applied.. The
AD-test is insensitive to ordering, so it is modified to

$TS=\begin{cases}AD\mathrel{\mathop{\mathchar 58\relax}}\langle
X_{\text{max}}^{\prime}\rangle^{\rm on}<\langle
X_{\text{max}}^{\prime}\rangle^{\rm off}\\\ -3\hskip
2.84526pt\mathrel{\mathop{\mathchar 58\relax}}\text{else}\end{cases},$ (2)

where $AD$ is the result of the AD-test comparing the on- and off-plane
distributions, and $-3$ is selected as it is well below the minimum of the AD-
test.

#### Scan for energy and galactic latitude thresholds

A scan has been used to select the optimal on/off splitting latitude, $b_{\rm
split}$, and minimum energy, $E_{\rm th}$, as uncertainties in GMF models and
source distributions make other approaches impractical. In this scan, each
trial [$E_{\rm th}$, $b_{\rm split}$] pair is used to form on- and off-plane
subsets and the $TS$ is extracted. To preserve the statistical strength of the
sparse FD data set, a coarse scan of $5^{\circ}$ steps in
$\mathinner{\\!\left\lvert\,b\,\right\rvert}$ from $20^{\circ}$ to
$35^{\circ}$ and 0.1 $\log_{10}(E/\text{eV})$ steps in energy from $18.4$ to
$19.4$ $\log_{10}(E/\text{eV})$ is used. The scan is performed on the data set
from Aab:2014kda , which includes events through Dec 31st 2012\. At the time
of writing, this scan data set represents $54\,\%$ of the analyzed events. The
remaining $46\,\%$ of events, the post-scan data set, is reserved as blind.

Figure 5: Parameter scan over 54% of the data.

Interestingly, as shown in Figure 5, all tested pairs result in $\langle
X_{\text{max}}^{\prime}\rangle^{\rm on}<\langle
X_{\text{max}}^{\prime}\rangle^{\rm off}$. An optimal [$E_{\rm th}$, $b_{\rm
split}$] of [$10^{18.7}$ eV,$30^{\circ}$] was found with a $TS=8.4$. The
selected [$E_{\rm th}$, $b_{\rm split}$] is applied as a prescription to the
post-scan data set, which independently confirms the result with a $TS=12.6$,
for a total $TS=21.0$ for the full data set.

#### Statistical significance

The chance probability of the observed TS occurring with in an isotropic sky
is tested using Monte Carlo methods on randomized skies derived from the real
data. To form each randomized sky, the arrival direction is first decoupled
from the energy and $X_{\text{max}}^{\prime}$ values of each event. These are
then randomly re-paired to create a new sky which maintains the real
$X_{\text{max}}$, energy, and sky exposure distributions, but has a scrambled
arrival direction/composition pairing. The above analysis is then used to
extract a $TS$ from each sky which is compared to the result in data. Skies
which display more extreme on-/off-plane differences than those observed in
data are tallied and used to calculate the probability of an isotropic sky
generating the observed $TS$. The results of this procedure are shown in
Figure 6.

Figure 6: The Monte Carlo determination of the post-scan (red) and all-data
(blue) significance with 1 and 10 billion randomized skies, respectively.

For the blind, post-scan data set, the prescribed [$E_{\rm th},b_{\rm split}$]
pair is used to split each randomized sky into on- and off-plane samples and a
$TS$ is extracted. In one billion random skies, only 5865 resulted in a more
extreme $TS$ than the 12.6 observed in data. This indicates a chance
probability of $5.87\times 10^{-6}$ which corresponds to 4.4 $\sigma$.

To calculate the significance of the result when the scan- and post-scan data
sets are combined, the entire analysis chain, including the scan, is
duplicated. In each random sky, 54 % of the data is used to scan for the
[$E_{\rm th},b_{\rm split}$] pair which results in the most extreme result,
fully penalizing for the scan. These values are then used to split all data in
the random sky into on- and off-plane subsamples and the $TS$ for the sky is
extracted. From 10 billion random skies, only 5964 resulted in a more extreme
$TS$ than the 21.0 observed in data. This indicates to a chance probability of
$5.96\times 10^{-7}$ which corresponds to 4.9 $\sigma$. The strong
penalization of the scanned data is evident as the additional 54 % of the data
(with $\Delta\langle X_{\text{max}}^{\prime}\rangle$ $=8.5$ g/cm2) only
resulted in an 11 % increase of the significance of the observation.

####  $X_{\text{max}}$ moments and trends

To illustrate the difference in composition on and off the plane, the first
two moments of the $X_{\text{max}}$ distribution in each 0.1
$\log_{10}(E/\text{eV})$ energy bin has been plotted in section 3 for both
regions. Above $10^{18.7}$ eV there is a clear separation in $\langle
X_{\text{max}}\rangle$ for all energy bins. Most energy bins also display a
separation in $\sigma(X_{\text{max}})$. Heavier primaries are expected to, on
average, have a shallower $X_{\text{max}}$ and lower shower-to-shower
fluctuations. Therefore the correlated difference seen here indicates that,
for this data sample, primaries from the on-plane region have a higher mean
mass than that of the off-plane region above $10^{18.7}$ eV.

To evaluate the degree to which fluctuation plays a role in the observed
result, the growth of the $TS$ over time has been plotted in Figure 8. The
time evolution of the signal is consistent with linear growth at a rate of 1.3
$TS$ yr-1. This behavior is in line with expectations for a real difference in
mean mass between the subsamples. The shaded region of Figure 8 shows
preliminary data from 2019. These reconstructions were not subject to a
validated reconstruction chain and may change. Still, when added, a 3.7/4.4
$\sigma$ (post-scan/all data) statistical significance is expected. The best
fit rate of growth of 1.3 $TS$ yr-1 remains unchanged.

Figure 7: The first (left) and second (right) moments of the $X_{\text{max}}$
distributions from on- and off-plane regions.

Figure 8: The time evolution of the TS with significance indicated on the
right. The shaded region is preliminary data.

##  4 Systematic uncertainties

Figure 2 shows that there is very little difference between the on- and off-
plane regions in the local reference frame of the detector. The data sets from
the two regions also consist of events measured with the same instrumentation
in the same location, and have been reconstructed with the same methods.
Because systematic uncertainties are derived from measurement effects in the
local frame, these similarities between the on- and off-plane samples mean
that the majority of uncertainty sources outlined in Aab:2014kda  will apply
equally to both regions and thereby cancel in a comparison. Furthermore, from
the acceptance, resolution, and bias studies in subsection 2.1, the two
regions are also free from selection and reconstruction biases.

Table 1: Systematic uncertainties on the difference in the means, $\Delta\langle X_{\text{max}}\rangle$, and widths, $\Delta\,\sigma(X_{\text{max}})$, of the on- and off-plane $X_{\text{max}}$ distributions. Source | Uncertainty [g/cm2] of
---|---
$\Delta\langle X_{\text{max}}\rangle$ | $\Delta\,\sigma(X_{\text{max}})$
$A$ correction | ${}^{+1.14}_{-0.71}$ | ${}^{+2.37}_{-1.61}$
$B$ correction | $\pm 0.36$ | $\pm 0.01$
$R$ correction | 0 | ${}^{+1.78}_{-0.24}$
Seasonal | ${}^{+1.00}_{-1.53}$ | ${}^{+1.19}_{-1.23}$
Instrumentation | $\pm 1.41$ | $\pm 1.41$
Sum in Quadrature | ${}^{+2.10}_{\mathbf{-2.23}}$ | ${}^{+3.49}_{\mathbf{-2.48}}$

To test for potential systematic effects derived from uncertainties in the
$A$, $B$, and $R$ corrections, possible seasonal effects, and small
differences between the instrumentation at different fluorescence telescope
sites, FD-sites, several studies were performed. All permutations of the
uncertainties in the $A$, $B$, and $R$ corrections were evaluated and the
maximum changes in $\Delta\langle X_{\text{max}}^{\prime}\rangle$ and
$\Delta\sigma(X_{\text{max}}^{\prime})$ were recorded. To look for on-/off-
plane biases associated with instrumentation differences, events seen by two
or more FD-sites were used to compare the $X_{\text{max}}$ reconstructions of
each site. No significant biases were found, and the instrumentation-derived
systematic uncertainties from Aab:2014kda  were adopted. To check for
systematics derived from a combination of the seasonal dependence of the FD
exposure and normal yearly variation of the atmospheric quality, $\langle
X_{\text{max}}^{\prime}\rangle$ and $\sigma(X_{\text{max}}^{\prime})$ for the
on- and off-plane regions were tracked over the course of the year, and is
shown in Figure 9. The maximum exposure-weighted difference between the two
regions at any point of the year has been adopted as the systematic
uncertainty due to seasonal effects. These uncertainties and the total
systematic uncertainty on $\Delta\langle X_{\text{max}}^{\prime}\rangle$ and
$\Delta\sigma(X_{\text{max}}^{\prime})$ are listed in Table 1.

Figure 9: Seasonal fluctuation of the first (left) and second (right)
$X_{\text{max}}$ moments for the on-plane (blue) and off-plane (red) samples.
The double black line shows the difference between these curves. The maximum
and minimum values of these differences are taken as the systematic
uncertainty on $\Delta\langle X_{\text{max}}^{\prime}\rangle$ and
$\Delta\sigma(X_{\text{max}}^{\prime})$.

#### Confidence level considering systematic uncertainties

The observed $\Delta\langle X_{\text{max}}^{\prime}\rangle$ of $9.1\pm 1.6$
g/cm2 is 4.1 times larger than the 2.2 g/cm2 systematic uncertainty listed in
Table 1. The observed $\Delta\sigma(X_{\text{max}}^{\prime})$ of $5.9\pm 2.9$
g/cm2 is 2.4 times larger than its 2.5 g/cm2 systematic uncertainty. This
means it is unlikely that the result could be entirely due to systematic
effects. However, the systematic uncertainties in $\Delta\langle
X_{\text{max}}^{\prime}\rangle$ and $\Delta\sigma(X_{\text{max}}^{\prime})$
may increase the likelihood of an extreme result occurring in data. To
quantify the result significance taking possible systematic effects into
account, a two step approach is taken. First, the on-/off-plane difference is
reduced by adding a value sampled from a Gaussian distribution with $\mu=2.2$
g/cm2 and $\sigma=2.5$ g/cm2 to the on-plane sample. Then, the AD-test is
applied to the resulting on- and off-plane distributions. Repeating this
process 1 million times results in a mean TS of $11.3\pm 0.5$ ‡‡‡Treating the
other side of the systematic errors in the same way results in a $TS$ of
$31.8\pm 1.1$ ($6.3\,\sigma$).. If these values are converted to significances
using the data from Figure 6, this corresponds to at least 3.3 $\sigma{}$.
Conservatively, to include systematic effects, this lower bound of 3.3
$\sigma{}$ is adopted as the confidence level of the result.

#### Cross-check: Results by zenith angle and FD-site

If astrophysical in nature, the difference in composition of UHECR arriving
from the on- and off-plane regions should be independently observed by each
FD-site and in all zenith angle ($\theta$) ranges. Figure 10 shows that the
difference in $X_{\text{max}}$ on and off the plane is indeed present in all
zenith angles bins and is also observed by all FD-sites independently. Even
more stringently, when the response of each FD-site is split in
$\cos^{2}\theta$ bins, it appears in 22 out of 28 bins. This independent
observation at all sites and zeniths is a strong confirmation of the stereo
study described above, showing that detector systematics can not play a large
role in the result. Furthermore, because the FD-sites have FoVs differing by
$90^{\circ}$ on average, each sees the galactic plane at a different local
geometry and time during the year, making it unlikely that some unidentified
detector, reconstruction, or atmospheric effect is causing the observed
anisotropy.

Figure 10: $\Delta\langle X_{\text{max}}^{\prime}\rangle$ by FD-site and
zenith bin.

##  5 Independent Test of the Result

To cross-check the validity of the previous results, an independent data set
is needed. The selection presented in section 2 aims to select only high-
quality events. In particular, the FidFoV cut is needed to ensure that the
distribution of selected events covers the expected range of $X_{\text{max}}$
values with an unbiased acceptance. The events rejected by the FidFoV cut are
still high-quality and perfectly well reconstructed. In this section, a second
data set is built from the events removed in this FidFoV cut to form an out-
FidFoV data set which is 82.9 % as large as the data set already used in this
analysis, the in-FidFoV data set.

Figure 11: Top panel: acceptance for out-FidFoV data set split in 0.1
$\log_{10}(E/\text{eV})$ energy bins. The energy in the legend is the lower
energy edge of the bin, excepting 19.7 which corresponds to
$\log(E/\text{eV})\geq 19.7$. Lower panel: the $X_{\text{max}}$ distribution
of events for the full energy range.

Figure 11 represents the $X_{\text{max}}$ acceptance of this out-FidFoV data
set in different energy ranges. For all energies, the efficiency maximized for
low $X_{\text{max}}$ values and decreases as $X_{\text{max}}$ increases.
Indeed, without the FidFoV selection, deep showers tend to be under-
represented (Aab:2014kda, ). Thus, the new data set is made out of shallower
$X_{\text{max}}$ events. As consequence, the $\langle X_{\text{max}}\rangle$
distribution of the new data set is on average $10$ g/cm2 shallower. Despite
this bias, if the difference between on- and off-plane is astrophysical, it
should also appear in the out-FidFoV data set.

To test this hypothesis, this new data set is divided into on- and off-plane
samples using the $E_{\rm th}=10^{18.7}$ eV, $b_{\rm split}=30^{\circ}$
splitting determined by the scan. The resulting on/off distributions show a
somewhat smaller $\Delta\langle X_{\text{max}}\rangle$ of ${\sim}5$g/cm2,
which is only 55 % of what was obtained with the in-FidFoV data set. The AD-
test returns a $TS=1.8$. The post-trial significance of this test is re-
evaluated with the method in section 3, using two million randomized MC trials
generated from the out-FidFoV data set. The corresponding $TS$ distribution is
represented by the orange histogram in Figure 12. The red dashed line depicts
the AD-test obtained for the data set. The corresponding significance is
${\sim}2.2\sigma$, which gives a probability of 0.03 that this would result
from an isotropic sky. The significance seen in this sample is, therefore,
much lower than with the in-FidFoV data set.

Figure 12: Anderson-Darling test from randomized skies (orange) and from
forward folding (blue) for the out-FidFoV data set. The red dashed line
indicates the value obtained for data in the on/off test.

The lower significance of the out-FidFoV sample begs the question, why is
there such a large difference between the in-FidFoV and out-FidFoV data sets?
The bottom panel of Figure 11 shows the distribution of $X_{\text{max}}$
values from events. Clearly the bulk of the distribution lies in a region
where the acceptance is decreasing quickly. It is possible that this could
reduce the sensitivity to a difference in composition. To test this
hypothesis, two million mock data sets have been generated from the on- and
off-plane $X_{\text{max}}$ distributions of the in-FidFoV data set. The non-
flat acceptance and lower $X_{\text{max}}$ resolution of the out-FidFoV data
set has then been forward folded onto these samples to create out-FidFoV mock
data sets which assume the difference in-FidFoV is real. The AD-test is then
computed for each mock data set. The blue histogram represents the
corresponding $TS$ distribution. It peaks exactly at the value observed with
the out-FidFoV data set, showing that the out-FidFoV is indeed less sensitive
to the tested anisotropy. Overall, the independent test with the out-FidFoV
shows that the on/off separation is present in both data sets, and that the
difference seen in the out-FidFoV sample is consistent with the magnitude of
the difference seen in the in-FidFoV sample.

##  6 Mapping the UHECR sky in composition

To aid interpretation of the latitude-dependent difference in composition, a
test statistic quantifying the relative difference in $X_{\text{max}}$ between
different parts of the sky is mapped in Figure 13, for UHECR primaries with
$E\geq 10^{18.7}$ eV. To produce this map from the in-FidFoV data set, first
the requirement $E\geq 10^{18.7}$ eV is imposed. Then, because small portions
of the sky are analysed, in contrast to the on/off study, each event has its
$B$, $R$, and $A$ corrected for based on its arriving declination instead of
its arriving galactic latitude. This is because local geometry has a time-
independent relationship with arrival declination. This means $B$, $R$, and
$A$ can be corrected equally well for each direction in the sky111Using
declination-dependent corrections changes the on/off comparison only by $+0.1$
g/cm2 and increases systematic uncertainties..

At this point, a top-hat sampling is used to collect all events with arrival
directions within $30^{\circ}$ of a point $(\ell,b)$ into an in-hat sample.
All other events are placed in an out-hat sample. The distributions of
$X_{\text{max}}^{\prime}$ for the in-hat and out-hat samples are then compared
using Welch’s t-test welch1938significance :

$TS=\frac{\langle X_{\text{max}}^{\prime\,\text{in}}\rangle-\langle
X_{\text{max}}^{\prime\,\text{out}}\rangle}{\sqrt{\left(\sigma\left(X_{\text{max}}^{\prime\,\text{in}}\right)/\sqrt{N^{in}}\right)^{2}+\left(\sigma\left(X_{\text{max}}^{\prime\,\text{out}}\right)/\sqrt{N^{out}}\right)^{2}}},$
(3)

where $N^{in}$ and $N^{out}$ are the event counts for the in- and out-hat
samples respectively§§§Because Welch’s t-test considers event statistics, the
FD arrival direction-dependent exposure is naturally treated through its use..
This procedure is repeated for top-hats centered on each point in a
$5^{\circ}$ by $5^{\circ}$ galactic latitude and longitude grid. The result is
shown in Figure 13, which illustrates the relative composition of UHECRs with
$E\geq 10^{18.7}$ eV arriving from each point in the sky.

In Figure 13, positive $TS$ values (red) indicate that events within
$30^{\circ}$ of that point have a lighter mean mass than the rest of the sky.
Negative values (blue) indicate that events within $30^{\circ}$ of that point
have a heavier mean mass than the rest of the sky. An excess of heavy
particles within $30^{\circ}$ of the galactic plane is visible. This can not
be due to detector systematics as they would be declination dependent and
appear as radial patterns centered on $\ell=-57^{\circ},b=-27^{\circ}$ due to
the geographic location of the Observatory.

Figure 13: Sky map of comic ray composition for $E\geq 10^{18.7}$ eV

##  7 Conclusions and Outlook

There is an apparent difference in the mean mass of primaries with energies
greater than $10^{18.7}$ eV that arrive from within $30^{\circ}$ of the
galactic plane. This has been observed at least at the $3.3\,\sigma$ level in
the standard hybrid data set used for $X_{\text{max}}$-based composition
analyses. It has now been independently confirmed with an additional
$2.1\,\sigma$ significance in a second hybrid data set formed from high
quality events cut by a selection aimed at reducing the bias caused by the
$X_{\text{max}}$-dependent event acceptance. The combined significance of
these two results has not yet been evaluated. Further tests of the on-/off-
plane difference are being planned using analyses of data from the SD and will
be reported elsewhere.

Currently, this result should be considered to primarily provide a new
verification of a mixed composition above the ankle as it is clear no such
difference could be observed in a flux with a single mass component. Though
the analysis provides a possible indication that the galactic magnetic field
may have an observable impact on mass-dependent anisotropies, the result found
in this analysis does not necessarily support a causal relationship with
galactic structures. The differing horizons of different nuclear species at a
given energy could also result in composition-dependent anisotropic patterns
Ding:2021emg . It is important, however, to note there is significant tension
with models Allard:2021ioh . Alternative scenarios are being explored along
these lines of thought.

## References

  * (1) A. Aab et al. (Pierre Auger), Nucl. Instrum. Meth. A 798, 172 (2015), 1502.01323
  * (2) I. Allekotte et al. (Pierre Auger), Nucl. Instrum. Meth. A 586, 409 (2008), 0712.2832
  * (3) J. Abraham et al. (Pierre Auger), Nucl. Instrum. Meth. A 620, 227 (2010), 0907.4282
  * (4) C. Bonifazi (Pierre Auger), Nucl. Phys. B Proc. Suppl. 190, 20 (2009), 0901.3138
  * (5) M. Mostafa, Nucl. Phys. B Proc. Suppl. 165, 50 (2007), astro-ph/0608670
  * (6) A. Aab et al. (Pierre Auger), PRD 90, 122005 (2014)
  * (7) E. Mayotte et al. (Pierre Auger), PoS ICRC2021, 321 (2021)
  * (8) A. Aab et al. (Pierre Auger), Phys. Rev. Lett. 125, 121106 (2020), 2008.06488
  * (9) J. Linsley (1963), ICRC 8 77
  * (10) A. Aab et al. (Pierre Auger), Science 357, 1266 (2017)
  * (11) A.A. Watson, Nucl. Phys. B Proc. Suppl. 22, 116 (1991)
  * (12) J. Szabelski, J. Wdowczyk, A.W. Wolfendale, J. Phys. G 12, 1433 (1986)
  * (13) Pierre Auger Collaboration, PRD 90, 122006 (2014)
  * (14) P. Abreu et al. (Pierre Auger), PoS ICRC2021, 311 (2021)
  * (15) A. Yushkov (Pierre Auger), PoS ICRC2019, 482 (2020)
  * (16) A. Aab et al. (Pierre Auger), JCAP 04, 038 (2017)
  * (17) T. Bergmann et al., Astropart. Phys. 26, 420 (2007)
  * (18) F. Riehn et al., PoS ICRC2015, 558 (2016)
  * (19) P. Abreu et al. (Pierre Auger), Astropart. Phys. 34, 368 (2011)
  * (20) J. Bellido (Pierre Auger), PoS ICRC2017, 506 (2018)
  * (21) T.W. Anderson, D.A. Darling, Ann. Math. Statist. 23, 193 (1952)
  * (22) F.W. Scholz, M.A. Stephens, JASA 82, 918 (1987)
  * (23) S. Engmann, D. Cousineau, JAQM 6 (2011)
  * (24) T. Pierog et al., PRC 92, 034906 (2015)
  * (25) B.L. Welch, Biometrika 29, 350 (1938)
  * (26) C. Ding, N. Globus, G.R. Farrar, Astrophys. J. Lett. 913, L13 (2021), 2101.04564
  * (27) D. Allard, J. Aublin, B. Baret, E. Parizot, Astron. Astrophys. 664, A120 (2022), 2110.10761