11institutetext:  Department of Physics & ITCP, University of Crete, GR-70013,
Heraklion, Greece 22institutetext: Institute of Astrophysics, Foundation for
Research and Technology-Hellas, Vasilika Vouton, GR-70013 Heraklion, Greece
33institutetext: Ludwig Maximilian University of Munich, Geschwister-Scholl-
Platz 1, 80539 Munich, Germany 44institutetext: Max Planck Institute for
Astrophysics, Karl-Schwarzschild-Straße 1, 85748 Garching, Germany

# Reconstructing Galactic magnetic fields from local measurements for
backtracking ultra-high-energy cosmic rays

Alexandros Tsouros<EMAIL_ADDRESS>Gordian Edenhofer , 3344
Torsten Enßlin , 3344 Michalis Mastorakis , 1122 Vasiliki Pavlidou , 1122

(Received ; accepted )

###### Abstract

Context. Ultra-high energy cosmic rays (UHECRs) are highly energetic charged
particles with energies exceeding $10^{18}$ eV. These energies are far greater
than those achieved in Earth-bound accelerators, and identifying their sources
and production mechanism can shed light into many open questions in both
astrophysics as well as high energy physics. However, due to the presence of
the Galactic magnetic field (GMF) they are deflected and hence the location of
their true source on the plane of the sky (PoS) is concealed. The
identification of UHECR sources is an open question, excacerbated by the large
uncertainties in our current understanding of $3$-dimensional structure of the
GMF. This difficulty arises from the fact that currently all GMF observations
are integrated along the line-of-sight (LoS). However, thanks to upcoming
stellar optopolarimetric surveys as well as Gaia data on stellar parallaxes,
we expect that local measurements of the GMF in the near future will become
available.

Aims. The question then arises: given such a set of (sparse) local GMF
measurements, what is the optimal way to use them in backtracking UHECRs
through the Galaxy? In this paper, we evaluate the reconstruction of the GMF,
in a limited region of the Galaxy, through Bayesian inference using principles
of Information Field Theory.

Methods. We employ methods of Bayesian statistical inference in order to
estimate the posterior distribution of the GMF configuration within a certain
region of the Galaxy from a set of sparse simulated local measurements. Given
the energy, charge, and arrival direction of a UHECR, we can backtrack it
through GMF configurations drawn from the posterior, and hence calculate the
probability distribution of the true arrival directions on the PoS, by solving
the equations of motion for each case.

Results. We show that, for a weakly turbulent GMF, it is possible to correct
for its effect on the observed arrival direction of UHECRs to within $\sim
3^{\circ}$. For completely turbulent fields, we show that our procedure can
still be used to significantly improve our knowledge on the true arrival
direction of UHECRs.

###### Key Words.:

Galactic magnetic field – Ultra high energy cosmic ray sources – Interstellar
turbulence

##  1 Introduction

The identification of ultra high energy cosmic ray (UHECR) sources remains one
of the central questions in modern high-energy astrophysics. A possible
resolution to this problem would shed light on the astrophysical mechanism
that produces UHECRs, as well as their composition (either protons or heavier
nuclei), which at high energies still remains a subject of debate. Knowledge
of their composition may, in turn, provide insights into high-energy particle
phenomenology at energies not yet probed by particle accelerators (Pavlidou &
Tomaras 2019; Romanopoulos et al. 2022a; Romanopoulos et al. 2022b).

Despite the fact that several theoretical candidates for UHECR sources have
been proposed (Bhattacharjee & Sigl 2000; Torres & Anchordoqui 2004), the
identity of their sources remains elusive. The primary reason for this is that
UHECRs are charged particles and are hence deflected by the Galactic magnetic
field (GMF) as well as the intergalactic magnetic field. Even if several of
the detected UHECR originate from a single nearby cosmic-ray-bright source di
Matteo et al. 2023, their arrival directions would be very spread out on the
sky, and any residual clustering would be centered off-source, due to magnetic
deflections. This is unlike the case of photons or neutrinos, where positional
identifications events with their likely sources can be attempted even with
very low-number statistics.

The main difficulty in resolving the GMF stems from the fact that $3D$
tomographic realizations of the intervening magnetic fields are notoriously
difficult to acquire. Specifically for the GMF, most observables that are
currently available are integrated along the line of sight (LoS). Due to this
limitation, and specifically for the GMF, the principal approach (Takami &
Sato 2010) is relying on GMF models that are acquired by parameter fitting
three separate components - a toroidal, a poloidal, and a Gaussian random
field (Jansson & Farrar 2012a; Jansson & Farrar 2012b; Sun et al. 2008; Sun &
Reich 2010).

Nonetheless, it is possible to acquire direct information regarding the $3D$
structure of the GMF. Gaia data on stellar distances have localised more than
a billion stars in the Galaxy by accurately determining stellar parallaxes
(Gaia Collaboration et al. 2016; Gaia Collaboration et al. 2021; Bailer-Jones
et al. 2021). This data - in addition to other available spectroscopic data -
has been used in order to construct $3D$ tomographic maps of the dust density
distribution of certain Galactic regions (Green et al. 2019; Lallement et al.
2019; Leike & Enßlin 2019a; Leike et al. 2020). These reconstructions,
however, do not constrain the magnetic field which is of primary interest in
cosmic ray physics.

Nonetheless, there do exist probes that provide information regarding the
structure of the GMF in $3D$. For instance, the linear polarization of
starlight is of particular interest; while starlight starts off unpolarized
from its source, it will usually acquire a linear polarization until it is
observed on Earth, due to dichroic absorption by dust particles aligned with
the ambient magnetic field (Andersson et al. 2015).

Upcoming optopolarimetric surveys, such as PASIPHAE (Tassis et al. 2018;
Maharana et al. 2021; Maharana et al. 2022; Magalhães 2012), are expected to
provide a large number of high quality stellar polarization measurements on
more than a million of stars. With stellar distances also known from the Gaia
survey, this data can be used to provide direct tomographic measurements of
the plane of the sky (PoS) component of the GMF at the location of the dust
clouds (Panopoulou et al. 2017; Pelgrims et al. 2022; Davis 1951;
Chandrasekhar & Fermi 1953; Skalidis et al. 2021). Used jointly with available
LoS information, we can expect to have local and sparse data of the GMF in the
near future, that could be used in order to provide a $3D$ tomographic map of
a particular region of interest. Given such a map, one can then backtrack
UHECRs through that region, thus improving the localisation of their source on
the sky, modulo the effects of the intergalactic magnetic field. To be more
specific, one might be interested in reconstructing the region of the GMF
through which the UHECRs ‘hotspots’ (Kawata et al. 2019;Abbasi et al.
2014;Pierre Auger Collaboration et al. 2017) have to travel through.

In this paper, we will address the problem of reconstructing the posterior
density function (PDF) for the true arrival directions of UHECRs, using local,
sparse, and (statistically) uniformly distributed GMF measurements. In
essence, we are interested in solving the inverse problem in a Bayesian
setting, wherein one is given local data and is tasked with calculating the
posterior distribution for the GMF configurations in the region of interest.
For Bayesian field inference, information field theory (IFT) was developed
Enßlin et al. 2009, and demonstrated in a number of contexts Enßlin 2022.
Here, we will use an adapted version of the algorithm used in Leike & Enßlin
2019b, after we generalise it in order to make it applicable to divergence-
free vector fields. Using the reconstructed posterior distribution over GMF
realizations, we correct for the effect of the GMF on the observed arrival
direction of the UHECR.

The work is structured as follows: in section 2 we describe the Bayesian
setting that will be employed, and motivate the principal components of the
algorithm that forms the core of our reconstruction scheme, as well as
describe our assumptions. In section 3 we present the main results of this
paper, by applying our reconstruction scheme on magnetic fields that are
locally and sparsely sampled, for different relative strengths of the
turbulent and the uniform component, as well as different sampling intervals.
In section 4 we summarize and discuss our conclusions.

##  2 Methodology

In Bayesian inference for continuous signals, we are in general interested in
reconstructing the posterior probability distribution $\varphi(\mathbf{x})$
defined over a domain $\mathcal{V}$ from a given data set $d$:

$P(\varphi|d)\propto P(d|\varphi)P(\varphi)$ (1)

The distribution $P(d|\varphi)$ is the _likelihood_ , and quantifies how
likely acquiring $d$ as our data is, given the configuration is
$\varphi(\mathbf{x})$. The function $P(\varphi)$ is the _prior probability_ ,
containing our information on the $\varphi(\mathbf{x})$ before $d$ is taken
into account. The proportionality factor in (1) is determined by normalising
the posterior to unity.

Mathematically, $\varphi(\mathbf{x})$ is a continuous function, and the
distributions $P(\varphi)$, $P(\varphi|d)$, and $P(d|\varphi)$ are
functionals. For example, an important case is the Gaussian distribution,

$\mathcal{G}(\varphi-m,\varPhi)\equiv\frac{\exp\left[-\frac{1}{2}(\varphi-m)\varPhi^{-1}(\varphi-m)\right]}{|2\pi\varPhi|^{\frac{1}{2}}},$
(2)

where $m(\mathbf{x})=\langle\varphi(\mathbf{x})\rangle$ is the mean. The
quantity
$\varPhi(x_{1},x_{2})=\langle\delta\varphi(x_{1})\delta\varphi(x_{2})\rangle$
is the covariance of the distribution, where we have defined
$\delta\varphi(x)\equiv\varphi(x)-m(x)$. $|\varPhi|$ denotes the determinant
of $\varPhi$. In equation (2), note that there is also an implicit
integration, that is,

$\displaystyle(\varphi-m)\varPhi^{-1}$ $\displaystyle(\varphi-m)\equiv$
$\displaystyle\int
d^{3}xd^{3}x^{\prime}[\varphi(\mathbf{x})-m(\mathbf{x})]\varPhi^{-1}(\mathbf{x},\mathbf{x}^{\prime})[\varphi(\mathbf{x}^{\prime})-m(\mathbf{x}^{\prime})].$

In order to avoid cluttered notation, this integration will always be written
implicitly.

###  2.1 Prior

The signal that we are interested in reconstructing is the GMF,
$\mathbf{B}(\mathbf{x})$, which is a vector quantity. Following equation (2),
the prior is written as

$\mathcal{G}(\delta\mathbf{B}(\mathbf{x}),M)\equiv\frac{\exp\left[-\frac{1}{2}\delta
B_{i}M_{ij}^{-1}\delta B_{j}\right]}{|2\pi M|^{\frac{1}{2}}},$ (3)

where $\delta B_{i}(\mathbf{x})=B_{i}(\mathbf{x})-\langle
B_{i}(\mathbf{x})\rangle$, with $\langle B_{i}(\mathbf{x})\rangle$ denoting
the mean of the distribution.

Here the indices correspond to individual components of the vector
$\mathbf{B}(\mathbf{x})$, and the Einstein summation convention is assumed,
and $\mathbf{B_{0}}$ is the mean field. Here, the covariance matrix inherits
two indices since by definition

$M_{ij}(\mathbf{x},\mathbf{x}^{\prime})\equiv\langle\delta
B_{i}(\mathbf{x})\delta B^{*}_{j}(\mathbf{x}^{\prime})\rangle,$ (4)

In lack of information regarding the geometry and statistics of the GMF in the
region of interest, we model the prior as being statistically isotropic and
homogeneous. In Fourier space, this is equivalent to assuming that the two-
point correlation function (4) takes the form

$\langle\delta\hat{B_{i}}(\mathbf{k})\delta\hat{B}^{*}_{j}(\mathbf{k}^{\prime})\rangle=\frac{(2\pi)^{3}}{2}\mathcal{P}_{ij}(\mathbf{k})\delta^{(3)}(\mathbf{k}-\mathbf{k}^{\prime})P(k),$
(5)

where

$\mathcal{P}_{ij}(\mathbf{k})\equiv\delta_{ij}-\hat{k}_{i}\hat{k}_{j},$ (6)

is the transverse projection operator, $\hat{k}_{i}$ denotes the unit
$\mathbf{k}$-vector in the $i$-th direction, $\delta_{ij}$ is the kronecker
delta, $\delta^{(3)}(\mathbf{k}-\mathbf{k}^{\prime})$ is the three dimensional
Dirac delta function, and $\delta\hat{B}_{j}^{*}(\mathbf{k})$ denotes the
complex conjugate of $\delta\hat{B}_{j}(\mathbf{k})$. Further, the norm of the
$\mathbf{k}$-vector is henceforth denoted as $k$. The function $P(k)$ is the
magnetic power spectrum, and needs to be inferred as well. Equation (5)
assumes that no prior expectation on the presence of helicity of any sign
exists111Helicity would add a complex, antisymmetric term of the form
$i\varepsilon_{ijl}k_{l}H(k)$, with $H(k)$ the real helicity spectrum, which
satisfies $-P(k)<H(k)<P(k)$.. We note that magnetic helicity in the
reconstructed magnetic field is not excluded thereby, it just needs to be
requested by the data and will not be enforced by the prior.

Instead of inferring $\mathbf{B}(\mathbf{x})$ directly, it is more practical
to work with a latent vector field $\bm{\varphi}(\mathbf{x})$ that does _not_
satisfy the divergence free condition, and thus contains more degrees of
freedom than needed. The Fourier space correlation structure for
$\bm{\varphi}$ will be assumed to be isotropic and homogeneous, and this will
take the form

$\langle\delta\hat{\varphi}_{i}(\mathbf{k})\delta\hat{\varphi}^{*}_{j}(\mathbf{k}^{\prime})\rangle=(2\pi)^{3}\delta_{ij}\delta^{(3)}(\mathbf{k}-\mathbf{k}^{\prime})P(k),$
(7)

with the understanding that $\mathbf{B}(\mathbf{x})$ and
$\bm{\varphi}(\mathbf{x})$ related to each other by an application of
$\mathcal{P}_{ij}(\mathbf{k})$ by

$\hat{B}_{i}(\mathbf{k})\equiv\frac{3}{2}\mathcal{P}_{ij}(\mathbf{k})\hat{\varphi}_{j}(\mathbf{k})=\frac{3}{2}\left(\hat{\varphi}_{i}(\mathbf{k})-[\hat{\varphi}_{j}(\mathbf{k})k_{j}]\frac{k_{i}}{k^{2}}\right).$
(8)

in Fourier space. This definition is such that the resulting field
$\mathbf{B}$ is guaranteed to be divergence free, since then
$k_{i}\hat{B}_{i}(\mathbf{k})=0$ in harmonic space, implying
$\nabla\cdot\mathbf{B}=0$ as required. The $3/2$ factor compensates for the
loss of power that the subtraction of degrees of freedom causes. This is
justified by the original a priori assumption of statistical isotropy for
$\bm{\varphi}$, which leads to equipartition along the three directions in
Fourier space. Equation (8) is in accordance with the correlation structure
assumed in equation (5) (Jaffe et al. 2012).

The power spectrum $P(k)$ is unknown and needs to be inferred as well. It is
modeled as the sum of a power law component and an integrated Wiener process
component (for details regarding the generative model for power spectra, the
reader should refer to Arras et al. 2020). The parameters that define it, as
well as their respective prior PDFs are

  1. 1.

The total field offset for all field components. This controls the mean value
around which the random field fluctuates, $\mathbf{B_{0}}$. We set a value of
$0$, which means that in the absence of data, the mean overall field is
assumed to be zero-centered.

  2. 2.

The standard deviation of the total offset. This is a random variable with a
prior log-normal distribution. We set its mean at $3$ $\mu$G and its own
standard deviation at $1$ $\mu$G. This reflects our current understanding of
typical GMF values in the interstellar medium that a (randomly oriented) mean
magnetic field of this typical strength within the reconstructed volume is
conceivable, but not enforced.

  3. 3.

The total spectral energy222Here, the term ‘energy’ is used in the context of
signal processing, and it does not refer to physical energy, although the two
are related.. This controls the amplitude of the fluctuations in configuration
space. Its prior probability distribution is log-normal with a mean and
standard deviation set to unity. This is essentially the standard deviation of
the field values calculated over the complete set of voxels.

  4. 4.

The spectral index; the exponent of the pure power law component. Its prior
distribution is normal, with a mean set as the Kolmogorov index, $-11/3$, and
we assume unit standard deviation.

  5. 5.

Amplitude of the integrated Wiener process component, controlling the
deviations from a pure power-law.

Its prior distribution is log-normal, with mean and standard deviation equal
to $1.5$ and $1$ respectively. These parameters are chosen so as to allow for
deviations from a pure power law without destroying the approximately power-
law overall behaviour of the power spectrum.

The parameters of the amplitude model are assumed to be relative to a unit-
less power spectrum, i.e. the parameters are assumed to be agnostic to changes
in the volume of the target subdomain $\mathcal{V}^{3}$. In table 1, we
summarise our choice for the probability distributions for the parameters of
the power spectrum, as well as their respective means and standard deviations.

Table 1: Parameters that define the power spectrum prior Parameter | Distribution | Mean | Standard deviation
---|---|---|---
Total offset ($\mathbf{B_{0}}$) | Not-applicable | $0$ | Not-applicable
Total offset st. dev. | Log-normal |  $3$ $\mu$G |  $1$ $\mu$G
Total spectral energy | Log-normal |  $1$ $\mu$G |  $1$ $\mu$G
Spectral index | Normal | $-\frac{11}{3}$ | $1$
Int. Wiener process amplitude | Log-normal | $1.5$ | $1$

In closing the discussion on the prior, we consider the implementation of
anisotropy. In general, the Galactic magnetic field in the region within which
we wish to reconstruct it, will be assumed to admit the two-component
structure

$\mathbf{B}(\mathbf{x})=\mathbf{B}_{0}+\mathbf{b}(\mathbf{x}),$ (9)

where $\mathbf{B}_{0}$ is a uniform field parallel to the PoS. Since the
observed arrival velocity of the UHECR will be chosen to be parallel to the
LoS (see section 3), the PoS component of $\mathbf{B}_{0}$ will dominate the
UHECR’s deflection and the LoS component will be hardly constrained.
Additionally, $\mathbf{b}(\mathbf{x})$ is a fluctuating field with zero mean,
which can be physically interpreted as the turbulent field. _Our prior is
agnostic to the direction and magnitude of $\mathbf{B}_{0}$, and it is a
central task to infer it from the data_. The relative strength of the two
components will be quantified by the turbulent-to-uniform ratio $\lambda$,
defined as

$\lambda\equiv\frac{b_{\text{rms}}}{B_{0}},$ (10)

where $B_{0}$ is the norm of $\mathbf{B}_{0}$, and $b_{\text{rms}}$ is the
root mean square value of the fluctuating component’s magnitude, in the domain
$\mathcal{V}$. For $\lambda\gg 1$ we have _strong turbulence_ , for
$\lambda\ll 1$ we have _weak turbulence_ , and for $\lambda\simeq 1$ we have
_intermediate turbulence_. As we will demonstrate, $\lambda$ is the main
parameter that controls how well the sources of UHECRs can be localized by the
use of local GMF data along, of course, with the sampling rate.

###  2.2 Likelihood

In order to construct the likelihood, we consider how the data is acquired
from the true signal $\bm{B}$. We will assume that the $i$-th datapoint is
given by a measurement process of the form

$\mathbf{d}^{(i)}=\int
R(\mathbf{x},\mathbf{x}_{i})\mathbf{B}(\mathbf{x})d^{3}x+\mathbf{n}^{(i)}.$
(11)

The operator $R$ is a map from the signal space to the data space. For this
proof of concept work, we will make the simplifying assumption that local
measurements of field components are possible, i.e.
$R(\mathbf{x},\mathbf{x}_{i})=\delta^{(3)}(\mathbf{x}-\mathbf{x}_{i})$. In
practice, most astrophysical measurements of magnetic fields (Faraday
rotation, UHECR deflection, etc.) provide at least a LoS, particle trajectory,
or even sub-volume averaged field strength. This only complicates the
numerical reconstruction of fields, but does not imply conceptual changes in
the formalism.

Thus, since we are assuming local measurements, $R$ is here assumed to be a
mask operator. In practice, since the mock signal is discretised into voxels,
$R$ is simply a very sparse matrix. The full vector $\bm{d}$ is therefore a
concantenation of a number of $3D$ vectors known on a finite number of points
inside the domain wherever we have measurements, and undefined anywhere else.
Simply put, then, the operator $R$ merely acts a selection operator that picks
out the values of $\mathbf{B}(\mathbf{x})$ where we have measurements.
Physically, the assumption is that at the positions where GMF measurements can
be obtained, all three components of $\mathbf{B}(\mathbf{x})$ can be measured.
This can be accomplished, e.g., by Zeeman observations for the LOS component
of the magnetic field. The POS component of the magnetic field can be
recovered through a combination of stellar optopolarimetry and stellar
distance measures. Optopolarimetric measurements of stars of known distances
can be decomposed to the contributions of individual clouds along the line of
sight Pelgrims et al. 2022. Then, the average direction of polarization angles
induced on starlight by dichroic absorption due to a single cloud reveals the
direction of the magnetic field in that cloud. The dispersion of these
directions can yield a measurement (within a factor of 2) of the magnetic
field strength in the same cloud (Chandrasekhar & Fermi 1953; Davis 1951;
Skalidis et al. 2021; Skalidis & Tassis 2021).

Additionally, the vector $\mathbf{n}^{(i)}$ is random noise added to the
$i$-th _noiseless_ data point. It is assumed to be drawn from a Gaussian
distribution with standard deviation equal to half the root-mean-square
magnitude of the mock signal field. Physically, $\mathbf{n}^{(i)}$ represents
the measurement error, assuming an average signal-to-noise ratio (SNR) of $2$.

The likelihood is then calculated by marginalizing over the noise:

$\displaystyle P(\bm{d}|\bm{B})$
$\displaystyle=\int\mathcal{D}\bm{n}P(\bm{d},\bm{n}|\bm{B})=\int\mathcal{D}\bm{n}P(\bm{d}|\bm{n},\bm{B})P(\bm{n}|\bm{B})$
$\displaystyle=\int\mathcal{D}\bm{n}\delta(\bm{d}-(R\bm{B}+\bm{n}))P(\bm{n}|\bm{B})$
$\displaystyle=\mathcal{G}(\bm{d}-R\bm{B},N).$ (12)

where $N$ is the noise covariance, and we used the shorthand notation
$R\bm{B}\equiv\int R(\mathbf{x},\mathbf{x}_{i})\mathbf{B}(\mathbf{x})d^{3}x$.
We can write the likelihood in terms of $\bm{\varphi}$ by absorbing
$\mathcal{P}_{ij}$ into $R$, that is, $R^{\prime}\equiv R\mathcal{P}_{ij}$.
Then, the likelihood becomes

$P(\bm{d}|\bm{\varphi})=\mathcal{G}(\bm{d}-R^{\prime}\bm{\varphi},N).$ (13)

###  2.3 Approaching the posterior: Geometric variational inference

In the previous two sections, we described how we model the prior and the
likelihood for our inference setting. From equation (1), the posterior is
proportional to the product of the two.

Due to the fact that the magnetic power spectrum $P(k)$ needs to be inferred
along with the configuration of the GMF, this inference problem is nonlinear.
A way to see this, is that the magnetic field $\mathbf{B}(\mathbf{x})$ couples
to $P(k)$ through equations (3) and (5). Moreover, the prior PDFs for the
parameters that determine $P(k)$ are not Gaussian (with the exception of the
spectral index), which is an additional source of nonlinearity. Further, there
is no small parameter that could be used in a perturbative analysis about a
linear inference case. For this reason, a non-perturbative scheme, called
geometrical variational inference (geoVI) developed by Frank et al. 2021 will
be utilised. In this section we motivate the basic premises of geoVI.

The idea is to approximate the true posterior, $P$, with an approximate one,
$Q$. The approximate posterior $Q$ is chosen such that the Kullback-Leibler
divergence (Kullback & Leibler 1951)

$D_{KL}(Q,P)\equiv\int dQ\log\left(\frac{Q}{P}\right)$ (14)

between the actual posterior $P$ and an approximate posterior $Q$ is
minimized. The main idea of geoVI is to achieve this minimization in a new
coordinate system, chosen such that $P$ \- in the new coordinate system -
locally closely resembles a normalized standard distribution. Once this is
done, the approximating posterior $Q$ is chosen to be of the form (3). Then,
the mean and covariance are chosen as the parameters with respect to which the
KL divergence is minimized.

A few more details on geoVI can be found in Appendix A and in Frank et al.
2021. Essentially, the algorithm provides approximate posterior samples that
can follow the non-Gaussian structure of the posterior to a certain degree.
geoVI can be invoked by the Numerical Information Field Theory (NIFTy) package
in Python (Selig et al. 2013; Steininger et al. 2017; Arras et al. 2019). The
input that is required is the likelihood and the prior of the original
physical model, as described in sections 2.2 and 2.1 respectively.

##  3 Results

###  3.1 General Procedure

We are now ready to construct a number of representative examples, in which we
will apply the procedure described in the previous section, in order to
reconstruct various assumed GMF geometries from a set of local, sparse, noisy,
$3D$ GMF mock observations. All of the examples that we will study are created
according to the following scheme:

  1. 1.

Define the domain: We choose a cube $\mathcal{V}$ of side length $L=3$ kpc
with periodic boundary conditions (topologically a $3$-torus), and evenly
divide it into $N^{3}$ voxels. The number $N$ defines our resolution.

  2. 2.

Produce a power spectrum $P(k)$ by randomly sampling each of its defining
parameters from its respective distribution (see section 2.1).

  3. 3.

Produce the latent field $\bm{\varphi}$ with a correlation structure dictated
by equation (7), with it taking a constant value over each voxel.

  4. 4.

Create the fluctuating part of the synthetic signal in Fourier space, by
acting on $\bm{\varphi}$ with the transverse projection operator;

$b_{i}(\mathbf{k})=\mathcal{P}_{ij}(\mathbf{k})\varphi_{j}(\mathbf{k}).$

After transforming into configuration space and computing the RMS value of the
field, $b_{\text{rms}}$, we add the uniform component that lies along the
$x-y$ plane (serving as our PoS) of a magnitude $\lambda^{-1}b_{\text{rms}}$,
for a desired value of $\lambda$. Rescale the total field such that the total
RMS magnitude is $5$ $\mu$G.

  5. 5.

Act on the mock signal with a mask operator, which uniformly chooses a
specified small fraction of voxels from the subdivided domain. This is done as
follows: if $N_{d}$ is the (approximate) number of data points we wish to
have, then we act on the original array with an operator which masks each
voxel with a probability of $1-N_{d}/N$. Roughly, then, the number of voxels
that survive is $N_{d}$, by construction.

However, it is more useful to refer to the _mean_ sampling ‘rate’; since by
construction the distribution of data points is spatially homogeneous, the
mean distance between the data points, as a function of $N_{d}$, is

$\ell=\frac{L}{N_{d}^{1/3}}.$

The mean sampling rate, is thus defined as

$k_{\text{sample}}\equiv\ell^{-1}=\frac{N_{d}^{1/3}}{L}.$ (15)

In physical reality, the data points will be located wherever HI clouds exist
in the region under study, which are not positioned uniformly within
$\mathcal{V}$, and so the sampling rate will vary with vertical distance from
the Galactic plane. Furthermore, there are also non-local measurements which,
for example, average over a LoS. In future work, we shall also consider such
effects.

  6. 6.

Add noise: to each data vector created in the previous step, we add a random
vector sampled from a multivariate Gaussian distribution meant to represent
observational error. The variance of the distribution is

$\sigma^{2}=\frac{B^{2}_{\text{rms}}}{4}=\frac{25}{4}(\mu\text{G})^{2}$ (16)

where $B_{\text{rms}}$ is the RMS value of the total magnetic field’s norm
within $\mathcal{V}$, which is set to $5$ $\mu$G (see step 4). The denominator
on the right hand side is chosen such that the average signal to noise ratio
is $2$.

  7. 7.

Apply the geoVI algorithm on the data set (see section 2.3). The output of the
algorithm are samples of the posterior distribution.

An illustrative example is shown in Fig. 1, where we choose
$k_{\text{sample}}=(600\text{ }pc)^{-1}$ as the mean sampling rate, and a
turbulent-to-uniform ratio $\lambda=0.2$ (see equation (10)). In this case,
the posterior mean essentially identifies the uniform component (Fig. 1(c))
provided the sparse and local data shown in Fig. 1(b). The dominance of the
ordered component results in a large scatter in the power spectrum (Fig. 2).
To see why this is the case, note that the zero mode of the field is

$\hat{\mathbf{B}}(0)=\int\mathbf{B}(\mathbf{x})d^{3}\mathbf{x}\propto\mathbf{B}_{0}.$
(17)

In the scenario of a strong uniform component, the modes
$\hat{\mathbf{B}}(\mathbf{k})$ with $k>0$ are poorly constrained by the data,
which mainly inform with respect to the zero mode.

It is important to note, however, that the UHECR’s deflection will primarily
be influenced by the uniform component (or the zero mode), and will be
affected much less by the small scale, turbulent fluctuations. Therefore, for
the purposes of our inquiry, the recovery of the zero mode is enough to
provide the leading order correction on the original arrival direction of
UHECRs. In the next section, we will apply the geoVI algorithm as exemplified
in Fig. 1 specifically to the case of backtracking UHECRs back to their
sources by utilising local data.

(a)

(b)

(c)

Figure 1:  Reconstruction of a 3D magnetic field within a cube of side $L=3=$
kpc; Top: The ground truth; a uniform field parallel to the $x-y$ plane with
magnitude $5$ $\mu$G, plus a turbulent, random, field with RMS magnitude of
$1$ $\mu$G, corresponding to $\lambda=0.2$ in equation (10). Middle: Local
data sampled randomly, with a constant mean sampling rate $(600\text{
pc})^{-1}$. The colormap is saturated at the maximum magnitude that appears in
Fig. 1(a). Bottom: The mean of the approximating posterior distribution
attained via the geoVI algorithmn based on the data provided in Fig. 1(b).
Figure 2: The power spectrum $P(k)$ corresponding to the example showcased in
Fig. 1. The black line is the power spectrum of the mock signal, while the
orange lines correspond to the posterior samples, providing the variance for
the posterior mean power spectrum estimate. The wavevector $k$ is given in
units of $2\pi/N$. The large scatter is due to the fact that we have few and
noisy datapoints, in a signal that is dominated by its uniform component, and
thus the modes with $k>0$ are poorly constrained.

(a)

(b)

Figure 3: The path of a test charged particle backtracked through the true
magnetic field shown in Fig. 1(a) (red path), and samples of the posterior GMF
given the sparse and local data, shown in Fig. 1(b) (magenta paths). For
visibility purposes, the paths through the posterior samples are cut short in
Fig. 1(a). The mean sampling rate of the local data is
$k_{\text{sample}}=(600\text{ }\text{pc})^{-1}$. Top: the two paths over the
whole cube $\mathcal{V}$ within which the field was reconstructed, and bottom:
the same as on the top, but with a truncated $x-y$ plane so that the
deviations between the paths are visually amplified.

###  3.2 Correcting the arrival directions of UHECRs

In the previous section we demonstrated the results of the geoVI algorithm
applied to sparse and local measurements of a magnetic field for the case of a
strong uniform component. However, the quality of the reconstruction should be
judged in context, as different degrees of resolution are required in
different applications. In fact, in this context, the GMF reconstruction is a
means to an end, the latter being correcting for the effect of the GMF on the
arrival direction of the UHECRs, given their arrival direction and energy as
observed on Earth.

Stated explicitly, _given an observed arrival direction
$\mathbf{\hat{v}}^{\text{obs}}$ of a charged particle of known lab frame
energy and charge, what is the posterior distribution for its original,
extragalactic direction assuming we are provided with sparse and local
measurements of the GMF within a region $\mathcal{V}$ within which it
traveled?_ Since the geoVI algorithm returns as output sample configurations
of the GMF drawn from the posterior, we can reconstruct the UHECR’s path
through each sample by solving the equations of motion backwards (for details
on how this is carried out in practice, see Appendix B). In the end, we are
left with a distribution for possible paths, all converging to the velocity
observed on Earth.

An example can be seen in Fig. 3. For the situation presented in Fig. 1, a
particle of charge $e$ with final position $\mathbf{r}_{i}=(L/2,L/2,0)^{T}$
(Earth’s location inside $\mathcal{V}$), velocity direction
$\mathbf{\hat{v}}^{\text{obs}}=(0,0,-1)^{T}$, and energy in the lab frame
$E=5\times 10^{19}\text{ }eV$ is backtracked through the true GMF (red line),
as well as samples of the posterior distribution of the GMF given the local
and sparse data shown in Fig. 1(b) (magenta lines). A kernel density
estimation (KDE) of the posterior distribution for the true arrival directions
of the UHECR, based on $100$ posterior samples, is also shown in Fig. 4. It
can be seen, at least in this case, that we are able to substantially correct
for the effect of the GMF. In order to quantify this improvement, we use the
Mahalanobis distance (Mahalanobis 1936) between a given arrival direction
$\hat{\mathbf{v}}$, and the posterior distribution $P(\hat{\mathbf{v}}|d)$
acquired from the geoVI algorithm,

$d_{M}[\hat{\mathbf{v}};P]=\sqrt{(\hat{\mathbf{v}}-\hat{\mathbf{\mu}})^{T}S^{-1}(\hat{\mathbf{v}}-\hat{\mathbf{\mu}})},$
(18)

where $S^{-1}$ and $\hat{\mathbf{\mu}}$ are the inverse covariance matrix and
mean of $P(\hat{\mathbf{v}}|d)$, respectively. Essentially, $d_{M}$ measures
how many standard deviations $\hat{\mathbf{v}}$ is located away from
$\hat{\mathbf{\mu}}$. If $\mathbf{\hat{v}}^{\text{true}}$ and
$\mathbf{\hat{v}}^{\text{obs}}$ are the true and observed arrival directions
respectively, then for the specific example shown in Figs 1, 3, 4, we
calculate $d_{M}[\hat{\mathbf{v}}^{\text{true}};P]=0.2\sigma$ and
$d_{M}[\hat{\mathbf{v}}^{\text{obs}};P]=35.0\sigma$.

Figure 4: Kernel density estimation of the posterior PDF for the UHECR arrival
directions on the PoS given local and sparse GMF data based on $100$ posterior
samples (magenta colormap). The GMF reconstruction problem is the one shown in
1. The red star denotes the true arrival direction, while the black star
denotes the observed arrival direction, and they are found $0.2\sigma$ and
$35.0\sigma$ from the posterior’s mean, respectively. This is the same data
shown in Fig. 3, projected on the PoS. The black ellipses signify the
$1\sigma$, $2\sigma$, and $3\sigma$ confidence levels. The inset plot zooms in
at the location of the posterior’s peak. The colormap is in arbitrary units.

We also quantify our results using a different metric. Let
$\mathbf{\hat{v}}^{\text{rec}}$ denote the arrival direction as obtained by
backtracking through a reconstructed magnetic field. If we define the angle

$\theta^{\prime}=\cos^{-1}(\mathbf{\hat{v}}^{\text{true}}\cdot\mathbf{\hat{v}}^{\text{rec}}),$
(19)

we are interested in the posterior mean of $\theta^{\prime}$, henceforth
referred to as $\theta\equiv\langle\theta^{\prime}\rangle$.

###  3.3 Effect of the sampling rate

In order to see how local sampling of the GMF might help us resolve the UHECR
initial directions in general, the backtracking procedure that led to Fig. 3,
is carried out for an ensemble of $100$ GMF configurations, for three
different values of the mean sampling rate; $k_{\text{sample}}=(300\text{
pc})^{-1}$, $(425\text{ pc})^{-1}$, and $(600\text{ pc})^{-1}$, and each time
a different field serves as ground truth and is reconstructed, via the process
described in section 3.1. The angle $\theta$ defined in the end of the last
section is calculated for each case, and the histograms for each sampling rate
considered are shown in blue, in Fig. 5. For reference, we also calculate the
angle between $\mathbf{\hat{v}}^{\text{true}}$ and the final velocity
direction, $\mathbf{\hat{v}}^{\text{obs}}$,

$\theta_{0}=\cos^{-1}(\mathbf{\hat{v}}^{\text{true}}\cdot\mathbf{\hat{v}}^{\text{obs}}).$
(20)

The histogram for $\theta_{0}$ is shown in red for each sampling rate
considered. At zeroth order, without attempting to reconstruct the GMF at all,
the UHECR’s path would not be altered, and so in this case
$\mathbf{\hat{v}}^{\text{true}}\simeq\mathbf{\hat{v}}^{\text{obs}}$ or
$\theta\simeq\theta_{0}$, trivially. Therefore, $\theta_{0}$ is the zeroth
order approximation, without any data points considered. As can be clearly
seen, the usage of the local and sparse data $\mathbf{d}$ significantly
corrects for the effect of the GMF on our knowledge of the UHECR arrival
direction. Further, as expected, increasing sampling rates tend to move the
PDF for $\theta$, $P(\theta)$, more towards $\theta=0$. In the extreme case of
a perfectly known GMF, $P(\theta)$ would be a delta function centered at zero.

These results are for the set value of $\lambda=0.2$. As $\lambda$ becomes
larger, or as the fluctuating/turbulent component $\mathbf{b}(\mathbf{x})$ in
eq. (9) becomes larger, our results should become increasingly worse. The
reason for this is that if $\lambda\ll 1$ then the uniform component
dominates, and we are required to only infer that, as the fluctuating part
will only play a subdominant role in the deflection of the UHECR. As $\lambda$
becomes larger, then the process will need to discern more irregular
structures that play an increasingly important role in deflecting the UHECR.
This task is more difficult to achieve with sparse data. Additionally, the
larger $\lambda$ becomes the less the UHECR path is expected to be deflected,
as for dominating $\mathbf{b}(\mathbf{x})$ (which has zero mean) the UHECR
will travel through regions that will partly cancel each other’s effect on the
path resulting on a random walk, rather than a systematic deflection, about
the true position of the source on the sky. Therefore, in Fig. 5, as $\lambda$
increases, we expect $P(\theta)$ (blue histogram) to diffuse towards larger
$\theta$ (owing to the increasing resolution that is required in that case).
Conversely, the PDF for the angular distance between the final velocity and
the true initial velocity (red histograms) is expected to shift towards
smaller $\theta$.

In Fig. 6 we test this expectation, by performing the same calculation as in
Fig. 5, but this time with increasing $\lambda$ and keeping
$k_{\text{sample}}=(600\text{ }\text{pc})^{-1}$ fixed. The results confirm our
intuition. We notice however, that even in the limit where the turbulence is
completely isotropic, $\lambda\rightarrow\infty$, the PDF resulting from
considering the reconstruction is narrower than the the one where only
observed directions are used.

To see this, consider raising $k_{\text{sample}}$ to $(300\text{
}\text{pc})^{-1}$ which is not an unrealistic expectation for the density of
HI clouds relatively close to the Galactic disk. Then, for a strongly
turbulent case $\lambda\rightarrow\infty$, our reconstruction can be seen in
Fig. 7. The mean of the posterior as found by the geoVI algorithm, given the
data of Fig. 7(b), can be seen in Fig. 7(c). In this example, the mean is
essentially a ‘combed’ version of the true field, in that it does contain all
the main structures and larger-scale features, but completely misses out the
fluctuations below a certain length scale. It should be made clear, however,
that the result of the reconstruction - as is the case with any inference
problem - is not just the mean but rather the whole posterior distribution.
The mean is drawn as a representative example of the posterior distribution,
but in general samples from the whole posterior are utilised in backtracking
the cosmic rays and thus obtaining the posterior for the initial arrival
directions.

In Fig. 8, the power spectrum defined via equation (5). We notice that in the
case of isotropic turbulence, the modes with nonzero wavevector are
constrained much better from the data compared to the anisotropic case of a
dominating zero mode.

(a) $k_{\text{sample}}=(600\text{ }\text{pc})^{-1}$

(b) $k_{\text{sample}}=(425\text{ }\text{pc})^{-1}$

(c) $k_{\text{sample}}=(300\text{ }\text{pc})^{-1}$

Figure 5: Histograms of the deviation angle between the true arrival direction
of UHECRs, and the one obtained by backtracking through the magnetic field as
reconstructed using local data sampled at three different mean rates
$k_{\text{sample}}=(600\text{ pc})^{-1}$, $(425\text{ pc})^{-1}$, and
$(300\text{ pc})^{-1}$ (blue). For each mean sampling rate, we consider $100$
independent GMFs. The observed arrival direction, the charge, and lab frame
energy are assumed to be known and are $(0,0,-1)^{T}$, $+e$, and $5\times
10^{19}\text{ }eV$, respectively. The turbulent-to-uniform magnitude ratio
(eq. (10)) is set at $\lambda=0.2$. For reference, we also plot the same angle
by completely neglecting the effect of the GMF on the particle’s trajectory
(red). It can be seen that local and homogeneously distributed GMF data
significantly improve our estimate of the true arrival directions of UHECRs,
and that there is successive improvement as $k_{\text{sample}}$ becomes
larger, as expected.

As was done with the weakly turbulent example, in Fig. 10 we show the actual
path of a UHECR with charge $e$, observed arrival direction
$\mathbf{\hat{v}}^{\text{obs}}=(0,0,-1)^{T}$, and energy in the lab frame
$E=5\times 10^{19}\text{ }eV$ through this strongly turbulent field (red
path). The magenta paths correspond to those that the same UHECR would have
taken if it travelled through samples of the posterior distribution of the
GMF. The black path is the zeroth order approximation - unaltered by the
presence of the GMF. Especially in Fig. 10(b), it can clearly be seen that in
the extreme case of completely isotropic turbulence, there is still a
substantial deviation due to the presence of the GMF, and the zeroth order
approximation (that is, assuming
$\hat{\mathbf{v}}^{\text{true}}\simeq\hat{\mathbf{v}}^{\text{obs}}$) is much
further away from the posterior’s mean than the true arrival direction. In
Fig. 9 we plot a KDE of the initial velocity’s posterior on the PoS, based on
$100$ posterior samples As with Fig. 4, the red star denotes the true arrival
direction, while the black star denotes the observed arrival direction. The
Mahalanobis distance between these directions and the posterior (magenta
contours) is computed at $2.0\sigma$ and $7.6\sigma$ respectively.

(a) $\lambda=0.5$

(b) $\lambda=1$

(c) $\lambda=\infty$

Figure 6: Histograms of the deviation angle between the true arrival direction
of UHECRs, and the one obtained by backtracking through the magnetic field as
reconstructed using local data sampled at three turbulent-to-uniform ratios
$\lambda=0.5$, $\lambda=1$, and $\lambda=\infty$, for $\lambda$ defined in eq.
(10). For value of $\lambda$, we consider $100$ independent GMFs. The observed
arrival direction, the charge, and lab frame energy are assumed to be known
and are $(0,0,-1)^{T}$, $+e$, and $5\times 10^{19}\text{ }eV$, respectively.
The sampling rate (eq. (15)) is set at $k_{\text{sample}}=600\text{ pc}$ . For
reference, we also plot the same angle by completely neglecting the effect of
the GMF on the particle’s trajectory (red). It can be seen that as $\lambda$
becomes larger, the deviations of the original direction as obtained by
backtracking through the reconstructed GMF and the true arrival direction
become more dispersed, while the observed discrepancy diffuses towards smaller
angles.

(a)

(b)

(c)

Figure 7: Reconstruction of a isotropically turbulent 3D magnetic field within
a cube of side $L=3\text{ }kpc$. Top: the ground truth; an isotropically
turbulent GMF with RMS magnitude of $5\text{ }\mu$G. Colormap denotes the
vector’s magnitude. Middle: Local data sampled uniformly. Mean distance
between each datapoint is $300\text{ }\text{pc}$. The colormap is saturated at
the maximum magnitude that appears in Fig. 7(a). Bottom: The mean of the
approximating Gaussian posterior distribution attained via the geoVI
algorithmn based on the data provided in Fig. 7(b). Figure 8: The power
spectrum corresponding to the example showcased in Fig. 7. The black line is
the power spectrum of the mock signal, while the orange lines correspond to
the posterior samples, providing the variance for the posterior mean power
spectrum estimate. The wavevector $k$ is given in units of $2\pi/N$. In this
case, all modes are relevant and well constrained, as opposed to the case of
Fig. 2. Further, note that as $k$ becomes larger the posterior of the power
spectrum starts deviating from that of the signal, implying a loss of
information at smaller lengthscales, as observed in Fig. 7. Figure 9: Kernel
density estimation of the posterior PDF for the UHECR arrival directions on
the PoS given local and sparse GMF data based on $100$ posterior samples
(magenta colormap). The GMF reconstruction problem is the one shown in Fig. 7.
The red star denotes the true arrival direction, while the black star denotes
the observed arrival direction, and they are found $2.0\sigma$ and $7.6\sigma$
from the posterior’s mean, respectively. This is the same data shown in Fig.
10, projected on the PoS. The black ellipses signify the $1\sigma$, $2\sigma$,
and $3\sigma$ confidence levels.

(a)

(b)

Figure 10: The path of a charged particle backtracked through the true
magnetic field shown in Fig. 7(a) (red path), and samples of the posterior GMF
given the sparse and local data, shown in Fig. 7(b) (magenta paths). For
visibility purposes, the paths through the posterior samples are cut short in
Fig. 10(a). The mean sampling rate of the local data is
$k_{\text{sample}}=(300\text{ }\text{pc})^{-1}$. Top: the two paths over the
whole cube $\mathcal{V}$ within which the field was reconstructed, and bottom:
the same as on the top, but with a truncated $x-y$ plane so that the
visibility of the deviations between the paths is amplified.

##  4 Summary, Conclusions, and Outlook

###  4.1 Summary

In this work, we used $3D$ local and sparse mock observations of the GMF
scattered across a cubic domain of $3$ kpc side within the Galaxy that are
statistically uniformly distributed, in order to obtain the posterior
distribution for the UHECR arrival directions before they entered the domain
of influence of the GMF. We considered a two component field, which was
comprised of a uniform and random (turbulent) part, for various relative
strengths between the two components. We used techniques from information
field theory and geoVI in order to construct a systematic way for calculating
sample configurations from the posterior PDF. We have adapted the correlated
field model of IFT to include the case divergence-free vector field such as
the GMF, so that the divergence-free condition is also taken into account
during the computation of the posterior samples.

Once the above framework was established, we used it to test the quality of
our reconstruction by comparing it to the ground truth, for various cases of
sampling rates and turbulent-to-uniform ratios. By backtracking a UHECR of
known arrival direction and rigidity (energy per unit charge) equal to
$5\times 10^{19}\text{ }V$ through sample configurations drawn from the
posterior PDF for the GMF given the sparse and local mock data, we calculated
their reconstructed direction, thus producing a posterior distribution for the
true arrival directions. For an ensemble of GMFs sharing the same turbulence-
to-uniform ratios and the same mean sampling rate, we were able to estimate
the PDF of the mean angle between the true arrival direction and the ones
obtained by backtracking through posterior samples, $\theta$, _and we were
able to show that for the case of weak turbulence very modest sampling rates
are needed in order for the predictions to not deviate from the true arrival
direction by more than a few degrees in relative angle_. In particular, for a
mean distance of $600$ pc between each measurement, the vast majority of the
relative angles between the inferred arrival directions and the true ones are
below $\sim 3^{\circ}$, for the cases considered. For increasing sampling
rates our results become increasingly better, as should be expected. For
reference, completely neglecting the information provided by the local and
sparse data, thus also neglecting the effect of the GMF, results in a relative
angle between the observed arrival direction and the true one ($\theta_{0}$)
with mean value $\theta_{0}\sim 15^{\circ}$.

As soon as the turbulent component starts dominating over the uniform one, our
results deteriorate, and in the extreme case of fully isotropic turbulence
there is a significant overlap between the PDF for $\theta$ and $\theta_{0}$.
However, even in this extreme and unrealistic case, the local and sparse GMF
data still provide a significant improvement in determining the true arrival
directions, and by extension the location UHECR sources.

###  4.2 Conclusions

In conclusion, our results are briefly summarised as follows:

  1. 1.

We have developed a systematic framework, based on information field theory,
that infers the posterior distribution of the true arrival directions of
UHECRs, given sparse and local data for the GMF scattered within a region of
the Galaxy.

  2. 2.

For weakly turbulent fields and by using uniformly sampled local data of the
GMF with an average signal to noise ratio of $2$, we are able to correct for
the effect of the GMF to within $\theta\sim 3^{\circ}$, where $\theta$ is the
angle between the inferred arrival direction and the true arrival direction.
The required mean distance between each measurement can be as large as $600$
pc, in this case. In applications to real data, the locations of the data
points will coincide with the locations of molecular clouds, and a mean
distance of $\sim 600$ pc is a rather conservative estimate.

  3. 3.

For the extreme, but unrealistic, case of an isotropically turbulent GMF our
results are worse compared to the case of a dominating structured field, but
local data can still provide substantial improvement of our knowledge on the
true arrival directions of UHECRs, and hence on the identity of their sources.

###  4.3 Final Comments & Outlook

In this work, the data points were chosen to have a constant mean sampling
rate throughout the domain of study, with the mean distance between the
measurements serving as an adjustable parameter. However, in real applications
this parameter is not set by us, but rather by Nature, as our measurements
will be localised wherever HI clouds exist inside the domain within which we
wish to reconstruct the GMF. In particular, their spacing will depend on the
distance from the Galactic plane. It might be the case that the assumption of
statistical homogeneity of data introduces biases into our results. In future
work, we will consider the case of inferring UHECRs true arrival directions,
by using data from MHD simulations of Galactic evolution in Milky Way-like
galaxies, where our mock observations will not be chosen at random by a
predetermined distribution, but rather from the distribution of clouds as
given by the simulation data. It should also be noted that $3D$ information on
each synthetic data point is a rather idealised situation; real data will
mainly consist of PoS information alone, often even integrated (with some
structured weighting) along the LoS. However, in the small deflection case
that is relevant for our problem, the PoS component is expected to dominate
the UHECR path’s deflection, since the observed velocity is parallel to the
LoS. This issue will be adressed in detail in future work. Finally, we like to
add that for any UHECR source that could be identified, the difference between
the observed and initial travel directions of the corresponding UHECRs will
then provide information on the galactic magnetic field.

###### Acknowledgements.

A.T. and V.P. acknowledge support from the Foundation of Research and
Technology - Hellas Synergy Grants Program through project MagMASim, jointly
implemented by the Institute of Astrophysics and the Institute of Applied and
Computational Mathematics. A.T. acknowledges support by the Hellenic
Foundation for Research and Innovation (H.F.R.I.) under the “Third Call for
H.F.R.I. Scholarships for PhD Candidates” (Project 5332). V.P. acknowledges
support by the Hellenic Foundation for Research and Innovation (H.F.R.I.)
under the “First Call for H.F.R.I. Research Projects to support Faculty
members and Researchers and the procurement of high-cost research equipment
grant” (Project 1552 CIRCE). A.T. would like to thank Vincent Pelgrims,
Raphael Skalidis, Georgia V. Panopoulou, and Konstantinos Tassis for helpful
tips and stimulating discussions. G.E. acknowledges the support of the German
Academic Scholarship Foundation in the form of a PhD scholarship
(”Promotionsstipendium der Studienstiftung des Deutschen Volkes”).

## References

  * Abbasi et al. (2014) Abbasi, R. U., Abe, M., Abu-Zayyad, T., et al. 2014, ApJ, 790, L21
  * Amari (2016) Amari, S.-i. 2016, Invariant Geometry of Manifold of Probability Distributions (Tokyo: Springer Japan)
  * Andersson et al. (2015) Andersson, B. G., Lazarian, A., & Vaillancourt, J. E. 2015, ARA&A, 53, 501
  * Arras et al. (2019) Arras, P., Baltac, M., Ensslin, T. A., et al. 2019, Astrophysics Source Code Library
  * Arras et al. (2020) Arras, P., Frank, P., Haim, P., et al. 2020 [arXiv:2002.05218]
  * Bailer-Jones et al. (2021) Bailer-Jones, C. A. L., Rybizki, J., Fouesneau, M., Demleitner, M., & Andrae, R. 2021, AJ, 161, 147
  * Bhattacharjee & Sigl (2000) Bhattacharjee, P. & Sigl, G. 2000, Phys. Rep, 327, 109
  * Chandrasekhar & Fermi (1953) Chandrasekhar, S. & Fermi, E. 1953, ApJ, 118, 113
  * Davis (1951) Davis, L. 1951, Physical Review, 81, 890
  * di Matteo et al. (2023) di Matteo, A., Anchordoqui, L., Bister, T., et al. 2023, arXiv e-prints, arXiv:2302.04502
  * Enßlin (2022) Enßlin, T. 2022, Entropy, 24, 374
  * Enßlin et al. (2009) Enßlin, T. A., Frommert, M., & Kitaura, F. S. 2009, Phys. Rev. D, 80, 105005
  * Frank et al. (2021) Frank, P., Leike, R., & Enßlin, T. A. 2021, Entropy, 23
  * Gaia Collaboration et al. (2021) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2021, A&A, 649, A1
  * Gaia Collaboration et al. (2016) Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1
  * Green et al. (2019) Green, G. M., Schlafly, E., Zucker, C., Speagle, J. S., & Finkbeiner, D. 2019, ApJ, 887, 93
  * Jaffe et al. (2012) Jaffe, T., Waelkens, A., Reinecke, M., Kitaura, F. S., & Ensslin, T. A. 2012, Hammurabi: Simulating polarized Galactic synchrotron emission, Astrophysics Source Code Library, record ascl:1201.014
  * Jansson & Farrar (2012a) Jansson, R. & Farrar, G. R. 2012a, ApJ, 757, 14
  * Jansson & Farrar (2012b) Jansson, R. & Farrar, G. R. 2012b, ApJ, 761, L11
  * Kawata et al. (2019) Kawata, K., di Matteo, A., Fujii, T., et al. 2019, in International Cosmic Ray Conference, Vol. 36, 36th International Cosmic Ray Conference (ICRC2019), 310
  * Kullback & Leibler (1951) Kullback, S. & Leibler, R. A. 1951, Ann. Math. Statist., 22, 79
  * Lallement et al. (2019) Lallement, R., Babusiaux, C., Vergely, J. L., et al. 2019, A&A, 625, A135
  * Leike & Enßlin (2019a) Leike, R. H. & Enßlin, T. A. 2019a, A&A, 631, A32
  * Leike & Enßlin (2019b) Leike, R. H. & Enßlin, T. A. 2019b, A&A, 631, A32
  * Leike et al. (2020) Leike, R. H., Glatzle, M., & Enßlin, T. A. 2020, A&A, 639, A138
  * Magalhães (2012) Magalhães, A. M. 2012, in Science from the Next Generation Imaging and Spectroscopic Surveys, 7
  * Mahalanobis (1936) Mahalanobis, P. C. 1936, Proceedings of National Institute of Sciences (India), 2, 6
  * Maharana et al. (2022) Maharana, S., Anche, R. M., Ramaprakash, A. N., et al. 2022, Journal of Astronomical Telescopes, Instruments, and Systems, 8, 038004
  * Maharana et al. (2021) Maharana, S., Kypriotakis, J. A., Ramaprakash, A. N., et al. 2021, Journal of Astronomical Telescopes, Instruments, and Systems, 7, 014004
  * Nocedal & Wright (2006) Nocedal, Jorge, P. & Wright, S. J. 2006, Large-Scale Unconstrained Optimization (New York, NY: Springer New York), 164–192
  * Panopoulou et al. (2017) Panopoulou, G. V., Psaradaki, I., Skalidis, R., Tassis, K., & Andrews, J. J. 2017, MNRAS, 466, 2529
  * Pavlidou & Tomaras (2019) Pavlidou, V. & Tomaras, T. 2019, Phys. Rev. D, 99, 123016
  * Pelgrims et al. (2022) Pelgrims, V., Panopoulou, G. V., Tassis, K., et al. 2022, arXiv e-prints, arXiv:2208.02278
  * Pierre Auger Collaboration et al. (2017) Pierre Auger Collaboration, Aab, A., Abreu, P., et al. 2017, Science, 357, 1266
  * Romanopoulos et al. (2022a) Romanopoulos, S., Pavlidou, V., & Tomaras, T. 2022a, arXiv e-prints, arXiv:2206.14837
  * Romanopoulos et al. (2022b) Romanopoulos, S., Pavlidou, V., & Tomaras, T. 2022b, in 37th International Cosmic Ray Conference. 12-23 July 2021. Berlin, 475
  * Selig et al. (2013) Selig, M., Bell, M. R., Junklewitz, H., et al. 2013, aap, 554, A26
  * Skalidis et al. (2021) Skalidis, R., Sternberg, J., Beattie, J. R., Pavlidou, V., & Tassis, K. 2021, A&A, 656, A118
  * Skalidis & Tassis (2021) Skalidis, R. & Tassis, K. 2021, A&A, 647, A186
  * Steininger et al. (2017) Steininger, T., Dixit, J., Frank, P., et al. 2017, ArXiv e-prints [arXiv:1708.01073]
  * Sun & Reich (2010) Sun, X.-H. & Reich, W. 2010, Research in Astronomy and Astrophysics, 10, 1287
  * Sun et al. (2008) Sun, X. H., Reich, W., Waelkens, A., & Enßlin, T. A. 2008, A&A, 477, 573
  * Takami & Sato (2010) Takami, H. & Sato, K. 2010, ApJ, 724, 1456
  * Tassis et al. (2018) Tassis, K., Ramaprakash, A. N., Readhead, A. C. S., et al. 2018, arXiv e-prints, arXiv:1810.05652
  * Torres & Anchordoqui (2004) Torres, D. F. & Anchordoqui, L. A. 2004, Reports on Progress in Physics, 67, 1663

##  Appendix A Geometric Variational Inference (geoVI)

In this appendix we provide a brief step-by-step overview of the geoVI
algorithm, which is the main algorithm used to approximate the posterior
distribution of the magnetic field, given sparse and local data. The idea is
to approximate the true posterior, $P$, with an approximate one, $Q$. The
approximate posterior $Q$ is chosen such that the Kullback-Leibler divergence
(Kullback & Leibler 1951)

$D_{KL}(Q,P)\equiv\int dQ\log\left(\frac{Q}{P}\right)$ (21)

between the actual posterior $P$ and an approximate posterior $Q$ is
minimized. The main idea of geoVI is to achieve this minimization in a new
coordinate system, chosen such that $P$ \- in the new coordinate system -
locally closely resembles a normalized standard distribution. Once this is
done, the approximating posterior $Q$ is chosen to be of the form (3). Then,
the mean and covariance are chosen as the parameters with respect to which the
KL divergence is minimized.

  1. 1.

First, a coordinate transformation $\bm{\varphi}=f(\bm{\xi})$ is performed,
such that the new _prior_ is Gaussian with unit covariance and zero mean, the
_standardized coordinate system_. Henceforth, the measure
$\mathcal{D}\bm{\xi}$ signifies integration over all possible configurations
of the vector field $\bm{\xi}$.

  2. 2.

In this new coordinate system, we calculate the Fisher information metric
(Amari 2016) $\mathcal{M}(\bm{\xi})$ for the likelihood, marginalising the
data and joining it with the unit prior metric, $\mathbb{1}$. Intuitively,
this may be regarded as a metric over the statistical manifold associated with
the likelihood $P(d|\bm{\xi})$.

  3. 3.

We seek another transformation $g$ on top of the original $f$, that turns
$\mathcal{M}+\mathbb{1}$ into the Euclidean metric, locally. The motivation is
that in this coordinate system, since the geometry of the statistical manifold
is as simple as possible, the original posterior is more likely to be
accurately described by a Gaussian. The accuracy of the step depends on the
choice of the expansion point for the local transformation. It can be computed
locally around the mean of the approximating Gaussian in this new coordinate
system.

  4. 4.

The KL-divergence (21) is minimised in the coordinate system to which $g$ maps
to with respect to the mean referred to in the previous step, using a second-
order quasi-Newton method, called Newton Conjugate Gradient (NewtonCG)
(Nocedal & Wright 2006). This is achieved by drawing sample configurations and
using them to compute the KL divergence, minimising it with respect to the
mean.

##  Appendix B Back-propagating the UHECRs through the GMF

The equations of motion for a relativistic charged particle of charge $q$ in a
static magnetic field in the lab frame $\mathbf{B}=\mathbf{B}(\mathbf{x})$ are

$\frac{d(\gamma m\mathbf{v})}{dt}=q\mathbf{v}\times\mathbf{B},$ (22)

and

$\frac{d(\gamma mc^{2})}{dt}=0,$ (23)

where $\gamma$ is the particle’s Lorentz factor, and $\mathbf{v}$ its
velocity. Equation (23) follows from the absence of an electric field.
Substituting equation (23) into (24), we get

$\frac{d\mathbf{v}}{dt}=\frac{qc^{2}}{E}\mathbf{v}\times\mathbf{B},$ (24)

where $E$ is the particle’s lab-frame energy.

If $\delta t$ is a small time interval, then the change in the velocity during
the interval $\delta t$ is

$\delta\mathbf{v}=\mathbf{v}(t)-\mathbf{v}(t-\delta t).$ (25)

Using equation (24), we may write

$\frac{\delta\mathbf{\hat{v}}}{\delta
t}=\frac{qc^{2}}{E}\mathbf{\hat{v}}\times\mathbf{B},$ (26)

where we divided both sides by $|\mathbf{v}|\simeq c$, and
$\mathbf{\hat{v}}\simeq\mathbf{v}c^{-1}$ is the velocity’s direction at any
given time. Substituting equation (25) into (26) and solving for
$\mathbf{v}(t-\delta t)$, we obtain

$\mathbf{\hat{v}}(t-\delta
t)=\mathbf{\hat{v}}(t)-\frac{Zec^{2}}{E}(\mathbf{\hat{v}}(t)\times\mathbf{B})\delta
t,$ (27)

where $q=Ze$, with $e$ the electron charge and $Z$ the atomic number of the
UHECR.

If we are also given the position of the UHECR at time $t$, and we wish to
calculate it at time $t-\delta t$, then we may write

$\mathbf{r}(t-\delta t)=\mathbf{r}(t)-\mathbf{\hat{v}}(t)c\delta t,$ (28)

where we once again made the assumption $|\mathbf{v}|\simeq c$ throughout the
particle’s path.

Therefore, if we are given the position, charge, lab-frame energy, and
observed arrival direction of a UHECR, we can use equations (27) and (28)
iteratively in order to solve the equations of motion numerically. We choose
$\delta t$ in the iterative process such that the length $c\delta t$ is equal
to the our resolution; the total domain is subdivided into voxels of side
length $c\delta t$, within which the GMF is assumed constant. For this work,
this amounts to setting $c\delta t=60$ pc, while the side of the total cubic
domain is $3$ kpc.

Finally, once the initial arrival direction is obtained - this happens when
the coordinates of the UHECR’s location exceeds the boundaries of the domain -
it is translated into galactic coordinates via the equations

$\displaystyle b$ $\displaystyle=\sin^{-1}(\hat{v}_{z}),$ (29)
$\displaystyle\ell$
$\displaystyle=\text{sign}(\hat{v}_{y})\cos^{-1}\left(\frac{\hat{v}_{x}}{\cos(b)}\right),$

where if $\ell<0$ add $2\pi$ to avoid negative angles.