1 Introduction

It is generally assumed that ultra-high energy cosmic rays (UHECRs) of extragalactic origin are deflected in the Galactic magnetic field (GMF). This assumption for deflections is, on one hand, based on astronomical measurements of Faraday rotation and synchrotron radiation, which indicate magnetic fields of micro-Gauss strengths [1]. On the other hand, measurements of the atmospheric depth of cosmic rays can be explained by a composition of light to medium-heavy nuclei with charge numbers \(Z\ge 1\) [2,3,4]. Together, these measurements predict deflections of the nuclei of several tens of degrees within the galaxy compared to their original extragalactic directions [5,6,7,8,9,10,11]. Most currently used methods for determining the magnetic field strength are limited to measuring the line-of-sight component. Only by studying the cosmic ray deflections, which are induced by the Lorentz force, the component orthogonal to the line-of-sight becomes accessible. Hence, studying the cosmic ray arrival directions on Earth can provide complementary information to the astronomical measurements used for evaluating cosmic magnetic fields.

In previous analyses that aimed to verify such deflections of cosmic rays, local regions of arrival were examined for energy ordering, but no scientific evidence for particle deflections was found for any region [12,13,14]. Recently, we introduced a fit method that determines the most probable extragalactic source directions by inverting the deflections that are caused by a specific GMF model and fitting a particle charge for each cosmic-ray event [15]. The several thousand free parameters are fitted using the backpropagation method developed for neural network training [16].

In this work, we present a novel approach in which all cosmic-ray arrival directions are simultaneously examined for alignment structures without relying on a certain GMF model [17]. The method is independent of energy ordering and analyzes only the arrival directions above a minimum energy threshold. In order to quantify coherent directional deflections, elliptically shaped regions are employed whose orientation is optimized by the frequency of neighboring particles (cf. Fig. 1). Coherence of adjacent ellipses is realized by means of a spherical harmonic expansion which assigns the local orientations of the set of ellipses.

The method is formulated as a likelihood ratio where for each cosmic-ray arrival direction, it is checked whether the cosmic ray is part of a deflection pattern or rather a particle of isotropic arrival directions. As a null hypothesis, circular regions are used instead of elliptical regions to distinguish the effects of coherent deflections from overdensities. The likelihood ratio is employed as the objective function for adjusting the spherical harmonic functions that specify magnetic field deflections. Thus, the test statistics of all measured particles are used to answer the question whether coherent deflections exist in the cosmic-ray arrival directions. If the answer to this question is positive, the orientations of the ellipses indicate the directional deflections caused by the GMF. This novel approach is hereafter referred to as: COherent Magnetic Pattern Alignment in a Structure Search (COMPASS).

Fig. 1
figure 1

Concept of the COMPASS method demonstrated in an astrophysical simulation (cf. Sect. 3). Circular symbols mark UHECR arrival directions, short lines denote the fitted orientation of coherent alignment, and the red stars show the directions of simulated sources. (Left) State of the system at initialization and (right) after the fit. (Top) Shapes of the elliptical signal probability density function in green and the Gaussian background probability density function in red. (Bottom) Local orientation of fitted alignment patterns. The color code in the lower right panel corresponds to the likelihood ratio resulting from the signal and background density functions (see text)

Full size image

The work is structured as follows: first, the analysis strategy is presented, covering the tangent vector field, the definition of the likelihood ratio, and its normalization. Two benchmark simulations are then introduced: one features simplified patterns of point sources to demonstrate the proof of concept and the other one is an advanced astrophysical simulation where UHECR nuclei from uniformly distributed sources are attenuated during propagation in the extragalactic universe. The ability to reconstruct the coherent directional deflections of the GMF and the advantages of using a circular reference model in the likelihood ratio are demonstrated in the following two chapters. Finally, the sensitivity of the method is investigated for both simulations, the simplified patterns, and the astrophysical universe.

2 Analysis strategy

The objective of the COMPASS method is to find alignment patterns in UHECR arrival directions simultaneously across the entire sky. In this approach, an adjustable vector field \(\hat{\varvec{u}}(\vartheta , \varphi )\) tangential to the local celestial sphere determines the orientation of elliptically shaped probability density functions (PDFs) that are centered on each cosmic-ray arrival direction. Here, \(\vartheta \) and \(\varphi \) denote the polar angle and azimuthal angle, respectively, in a spherical coordinate system. The likelihood that the distribution of neighboring arrival directions is better described by an elliptical PDF than by a background hypothesis is then evaluated to optimize the orientation of the ellipses’ major axes for all cosmic-ray events in one single step. In this way, the vector field \(\hat{\varvec{u}}(\vartheta , \varphi )\) locally aligns with elongated structures which are expected to occur from UHECR deflections in the GMF.Footnote 1 Technically, a likelihood ratio (cf. Eq. (14)) serves as an objective function in a minimization based on gradient descent. Additional constraints within the analysis can be accounted for by adding a corresponding penalization term to the objective function.

The basic concept of the COMPASS method is demonstrated in Fig. 1 where the initial state of the system is shown in the left panel and the fitted state in the right panel. Here, the initialization of the vector field \(\hat{\varvec{u}}(\vartheta , \varphi )\) is equal to the local unit vector \(\hat{\varvec{e}}_\vartheta \). The upper panel shows the ellipsoidal PDF (green) and a corresponding Gaussian background PDF (red). One can see that the orientation of the ellipse has changed after the fit where an alignment with a prominent pattern originating from the red marked source is clearly visible. Additionally, the lower panel indicates the orientation of the adjusted vector field \(\hat{\varvec{u}}(\vartheta , \varphi )\) in the vicinity of the pattern. The orientation has changed considerably only for the cosmic-ray events which are part of the pattern, whereas the ellipses of most of the isotropic events have not changed substantially during the fit. This finding is also visualized by the color-coded likelihood ratio where high values are found only for events that are part of the pattern.

The method requires a high number of fit parameters, both for the parameterization of the vector field \(\hat{\varvec{u}}(\vartheta , \varphi )\) and the UHECR model for the likelihood ratio. For the analyzed simulated data set of UHECRs with energies above 40 EeV, the number of free fit parameters is of the order of \({\mathcal {O}}(1000)\). Our method uses the software package TensorFlow [16] to perform a gradient descent-based optimization in this high dimensional parameter space. To enable the computation of gradients within the scope of the backpropagation technique used in the field of machine learning, all operations of this analysis (cf. following subsections) – spherical harmonics expansion, parameterization of density functions, vector algebra operations, likelihood ratio – were written with the TensorFlow API.

2.1 Tangent vector field

A particular challenge is the parameterization of the tangent vector field \(\hat{\varvec{u}}(\vartheta , \varphi )\) which is meant to describe the orientation of deflection patterns caused by the GMF. The large-scale component of the GMF is most likely responsible for patterns of coherent deflection [18, 19]. Thus, the orientation of patterns is expected to vary only slightly within a local domain of the sky.

Here, the adaptable vector field is first realized by a constant vector field \(\hat{\varvec{u}}_0(\vartheta , \varphi )\) which serves as an initialization and is then modified locally by an angle \(\varPsi (\vartheta , \varphi )\). To preserve the local coherence of the resulting vector field \(\hat{\varvec{u}}(\vartheta , \varphi )\), the modification angle is parameterized by a spherical harmonics expansion of order k:

$$\begin{aligned} \varPsi (\vartheta , \varphi ) = \sum _{\ell =0}^{k} \sum _{m=-\ell }^{\ell } a_\ell ^m Y_\ell ^m(\vartheta , \varphi ), \end{aligned}$$
(1)

where \(Y_\ell ^m(\vartheta , \varphi )\) are the spherical harmonics functions and \(a_\ell ^m\) represent a set of free fit parameters to model any continuous differentiable function on the surface of the sphere. Rapid variations on small angular scales can be suppressed by demanding the order k of the spherical harmonics expansion to have an upper limit of \(k=5\). Typical examples for this value k used in the analysis are \(k=4\) and \(k=5\) which yield 25 and 36 free fit parameters, respectively. The resulting modification \(\varPsi \) for a certain direction \(\hat{\varvec{r}}\) is then described by a rotation of \(\hat{\varvec{u}}_0\) around the axis \(\hat{\varvec{r}}\):

$$\begin{aligned} \hat{\varvec{u}}(\vartheta , \varphi ) = R (\hat{\varvec{r}}(\vartheta , \varphi ), \varPsi (\vartheta , \varphi )) \, \hat{\varvec{u}}_0 (\hat{\varvec{r}}(\vartheta , \varphi )), \end{aligned}$$
(2)

where the two arguments of the rotation matrix R are the rotation axis and angle, respectively. Note that here the polar angle \(\vartheta \) is defined as being consistent with the Galactic latitude b; thus, in Cartesian coordinates \(\hat{\varvec{r}}\) is given by:

$$\begin{aligned} \hat{\varvec{r}}(\vartheta , \varphi ) = (\cos \vartheta \, \cos \varphi , \, \cos \vartheta \, \sin \varphi , \, \sin \vartheta )^T. \end{aligned}$$
(3)

For the initialized vector field \(\hat{\varvec{u}}_0 (\hat{\varvec{r}})\), three approaches were investigated in this work. The hairy ball theorem states that there exists no nonvanishing continuous tangent vector field on the surface of the three-dimensional sphere [20]. Thus, \(\hat{\varvec{u}}_0\) always exhibits at least one region on the sphere where the vector field either radially diverges at a certain point or where it circularly curls around it. The following three initializations were used:

  • JF12 GMF: An intuitive approach is to initialize the fit with the best guess of the pattern orientations, e.g. the predictions of the currently most reliable GMF model, namely that developed by Jansson and Farrar [21] (JF12). Here, to obtain the local direction of deflection in the direction of \(\hat{\varvec{r}}\), a magnetically highly rigid particle of \(10^{20}\) eV is backtracked, leaving the Galaxy in direction \(\hat{\varvec{r}}'\). Then, the local tangent vector field is defined as \(\hat{\varvec{u}}_0 (\hat{\varvec{r}}) = \hat{\varvec{r}} \times (\hat{\varvec{r}} \times \hat{\varvec{r}}') / C\) where C is determined by the normalization according to \(|| \hat{\varvec{u}}_0 || = 1\).

  • Galactic meridians: Here, the local tangent vector \(\hat{\varvec{u}}_0 (\hat{\varvec{r}})\) is equal to the local spherical unit vector \(\hat{\varvec{e}}_\vartheta \) in the Galactic coordinate system. The advantage of this initialization is that it is independent of a certain GMF model, while an overall symmetry with respect to the Galactic plane is still maintained. Additionally, certain models favor a general deflection preference towards the Galactic plane [22], which is approximately realized in this case.

  • Equatorial meridians: In analogy to the Galactic meridians, here the initialization is equal to the unit vector \(\hat{\varvec{e}}_\theta \) in the Equatorial coordinate system. This initialization has the advantage that one of the two points of divergence is located in the blind region of a ground-based Observatory. Here, it is used only as a crosscheck for the fit reliability.

Fig. 2
figure 2

Visualization of the utilized tangent vector fields \(\hat{\varvec{u}}_0\) in the Galactic coordinate system. (Top) Local orientation of deflection patterns from simulations in the JF12 field calculated by the displacement of an \(R=10^{20}\) eV particle. (Middle) Meridians in Galactic coordinates along the local \(\hat{\varvec{e}}_\vartheta \) unit vector. (Bottom) Meridians in Equatorial coordinates

Full size image

A visualization of the three tangent vector field initializations \(\hat{\varvec{u}}_0\) is presented in Fig. 2. Regions with curls or divergences of the vector field can be seen in all three initializations. At these locations, the analysis exhibits a decreased sensitivity to find locally aligned structures as the underlying vector field \(\hat{\varvec{u}}_0\) cannot describe them. For the initialization of JF12 GMF, two of these features are visible at Galactic coordinates \((l, b) \approx (90^\circ , 0^\circ )\) and \((l, b) \approx (-60^\circ , -20^\circ )\). The Galactic and equatorial meridians initializations exhibit two divergences at the northern and southern poles of the respective coordinate system. An example of the working principle of the modification function \(\varPsi (\vartheta , \varphi )\) for the Galactic meridians initialization is shown in the lower panel of Fig. 1.

For the JF12 GMF initialization – depending on the reliability of the model predictions – it may be beneficial for the sensitivity to include a penalization term for large model deviations in the objective function. This penalization can be achieved by limiting the integrated squared amplitude of \(\varPsi \) over the entire celestial sphere as:

$$\begin{aligned} F = \frac{1}{4 \pi } \int _0^{2\pi } \int _{-\pi /2}^{\pi /2} \varPsi ^2(\vartheta , \varphi ) \, \cos (\vartheta ) \, d\vartheta \, d\varphi . \end{aligned}$$
(4)

As a potential improvement in scenarios where many directions exhibit alignment patterns, the tangent vector field \(\hat{\varvec{u}}\) may be directly defined in the form of vector spherical harmonics (VSH) [23]. In this way, positions of curls and divergences of the vector field can be shifted dynamically over the sky during the fit. Thus, they will likely stall in sky regions without noteworthy contributions to the likelihood ratio, i.e. in regions without prominent alignment patterns.

2.2 Maximum likelihood ratio

The COMPASS method evaluates the distribution of arrival directions around each cosmic-ray event i in order to search for the existence of an elongated structure. Here, the likelihood is defined in analogy to the approach in [24]: the total UHECR sky model consists of a sum of a signal part \({\mathcal {E}} (\hat{\varvec{r}}) \, {\mathcal {S}}_i (\hat{\varvec{r}})\) (with the contribution \(| f_i |\)) which captures elongated patterns and a purely isotropic part \({\mathcal {E}} (\hat{\varvec{r}})\) which represents the geometrical exposure of the observatory [25]. This likelihood \(\log ({\mathcal {L}}_i^{\mathcal {S}})\) is then compared to a suitable reference model by calculating the likelihood ratio. The fundamental difference with respect to the approach in [24] is that each cosmic-ray event i provides a separate density function \({\mathcal {S}}_i (\hat{\varvec{r}})\) including a fit parameter \(f_i\) which describes the contribution of the respective event to the log-likelihood as:

$$\begin{aligned} \log ({\mathcal {L}}_i^{\mathcal {S}})&= \sum _j^{N_\text {tot}} \log \left[ \, | f_i | \times {\mathcal {E}} (\hat{\varvec{\varTheta }}_j) \, {\mathcal {S}}_i(\hat{\varvec{\varTheta }}_j) \right. \nonumber \\&\left. \quad + (1 - | f_i |) \times {\mathcal {E}} (\hat{\varvec{\varTheta }}_j) \, \right] , \end{aligned}$$
(5)

where \(N_\text {tot}\) is the total number of events in the data set and \(\hat{\varvec{\varTheta }}_j\) the unit vector of the arrival direction of cosmic ray j. Here, the signal and background contributions, \({\mathcal {E}} (\hat{\varvec{r}}) \, {\mathcal {S}}_i (\hat{\varvec{r}})\) and \({\mathcal {E}} (\hat{\varvec{r}})\), are normalized over the surface A of the sphere:

$$\begin{aligned} \iint _A {\mathcal {E}} (\hat{\varvec{r}}) \, \text {d}A = 1 \quad \text {and} \quad \iint _A {\mathcal {E}} (\hat{\varvec{r}}) \, {\mathcal {S}}_i(\hat{\varvec{r}}) \, \text {d}A = 1. \end{aligned}$$
(6)

The set of \(f_i\) represents a total number of \(N_\text {tot}\) free fit parameters which are initialized with a value close to zero. Thus, the total number of free parameters of the COMPASS method is \(n_\text {fit} = (k+1)^2 + N_\text {tot}\), where the \((k+1)^2\) part comes from the spherical harmonics coefficients \(a_\ell ^m\).

The signal part \({\mathcal {S}}_i (\hat{\varvec{r}})\) of Eq. (5) is constructed as an elliptically shaped density function which is centered at cosmic-ray direction \(\hat{\varvec{\varTheta }}_i\). The major axis is aligned with the local direction of the tangent vector field \(\hat{\varvec{u}}_i \equiv \hat{\varvec{u}}(\vartheta _i, \varphi _i)\). Since the GMF is not well known, there is no accurate mathematical description for the expected shape and size of a deflection pattern. Here, the density function \({\mathcal {S}}_i\) is parameterized on the basis of a non-symmetrical Gaussian distribution where the width follows an ellipse equation as

$$\begin{aligned} {\mathcal {S}}_i (\hat{\varvec{r}}) = C \times \exp \left( - \frac{(\hat{\varvec{r}} \cdot \hat{\varvec{u}}_i)^2}{\delta _{\text {max}}^2} -\frac{(\hat{\varvec{r}} \cdot (\hat{\varvec{\varTheta }}_i \times \hat{\varvec{u}}_i))^2}{\delta _{\text {min}}^2} \right) , \end{aligned}$$
(7)

for all directions \(\hat{\varvec{r}}\) located in the same hemisphere as cosmic ray i, i.e. \(\hat{\varvec{\varTheta }}_i \cdot \hat{\varvec{r}} \ge 0\). In the opposite hemisphere of the sky the density function is set to zero. The hyperparameters \(\delta _{\text {max}}\) and \(\delta _{\text {min}}\) denote the angular extent in the direction of the ellipses’ semi-major and semi-minor axes, respectively. C denotes a normalization factor which is investigated in Sect. 2.3.

For every cosmic ray i, the minimal distance \(\delta _{i, j}^{\perp }\) between the direction \(\hat{\varvec{\varTheta }}_j\) of the neighboring cosmic ray j and the orthodrome \(\zeta \), as defined by \(\hat{\varvec{\varTheta }}_i\) and \(\hat{\varvec{u}}_i\), is given by the relation:

$$\begin{aligned} \sin (\delta _{i, j}^{\perp }) = \frac{\hat{\varvec{\varTheta }}_j \cdot (\hat{\varvec{\varTheta }}_i \times \hat{\varvec{u}}_i)}{|| \hat{\varvec{\varTheta }}_i \times \hat{\varvec{u}}_i ||}. \end{aligned}$$
(8)

Since \(\hat{\varvec{\varTheta }}_i\) and \(\hat{\varvec{u}}_i\) are unit vectors, the term \(|| \hat{\varvec{\varTheta }}_i \times \hat{\varvec{u}}_i ||\) is equal to one. Thus, the numerator of the second term in the exponential function of Eq. (7) can be identified as the transverse displacement of cosmic-ray direction \(\hat{\varvec{\varTheta }}_j\) relative to the fitted orientation \(\hat{\varvec{u}}_i\). Likewise, the great-circle distance along the orthodrome \(\zeta \) – and therefore along the pattern orientation – is given by \(\sin (\delta _{i, j}^{\parallel }) = \hat{\varvec{\varTheta }}_j \cdot \hat{\varvec{u}}_i\). Thus, for small angles, i.e. \(\sin \delta _{i, j} \approx \delta _{i, j}\), Eq. (7) can be identified as a two-dimensional Gaussian distribution on the sphere where contour lines of equal function values \({\mathcal {S}}_i\) follow an ellipse equation:

$$\begin{aligned} \left( \frac{\delta _{i, j}^{\parallel }}{\delta _{\text {max}}} \right) ^2 + \left( \frac{\delta _{i, j}^{\perp }}{\delta _{\text {min}}} \right) ^2 = \text {const}. \end{aligned}$$
(9)
Fig. 3
figure 3

Normalized probability density distributions for (top) signal part \({\mathcal {S}}_i (\hat{\varvec{r}})\), (middle) the Gaussian reference part \({\mathcal {G}}_i (\hat{\varvec{r}})\), and (bottom) the geometrical exposure \({\mathcal {E}} (\hat{\varvec{r}})\). For the signal and Gaussian reference parts, the cosmic-ray direction \(\hat{\varvec{\varTheta }}_i\) is centered at Galactic coordinates \(\vartheta = 0^\circ \) and \(\varphi = 0^\circ \) and the ellipse geometry is \(\delta _\text {max}=30^\circ \) and \(\delta _\text {min}=10^\circ \) where the major axis is aligned with the \(\hat{\varvec{e}}_\vartheta \) unit vector. The geometrical exposure is given by the geometry of the Pierre Auger Observatory with a maximum zenith angle of \(80^\circ \)

Full size image

For the reference model of the likelihood ratio, two approaches were explored in this work: a purely isotropic hypothesis following the geometrical exposure \({\mathcal {E}} (\hat{\varvec{r}})\) and a Gaussian reference model with identical signal strength \(f_i\) (compared to the elliptical ones) to cancel out overdensities.

  • Isotropic reference model \({\mathcal {E}} (\hat{\varvec{r}})\): The highest sensitivity to reject an isotropic scenario is obtained by testing explicitly against this hypothesis in the likelihood ratio. Thus, each cosmic ray i provides the same log-likelihood contribution:

    $$\begin{aligned} \log ({\mathcal {L}}_i^{\mathcal {R}})= \sum _j^{N_\text {tot}} \log ( {\mathcal {E}} (\hat{\varvec{\varTheta }}_j) ), \end{aligned}$$
    (10)

  • Gaussian reference model \({\mathcal {G}}_i (\hat{\varvec{r}})\): Here, the reference model is provided by Eq. (7) where both the major axis and minor axis radii are set to \(\sqrt{\delta _\text {max} \times \delta _\text {min}}\) with respect to the ellipse dimensions of the signal distribution. By choosing the geometric average of both dimensions, the effective solid angle of the symmetric reference model is unchanged and Eq. (7) can be written as a symmetric Gaussian-like distribution:

    $$\begin{aligned} {\mathcal {G}}_i (\hat{\varvec{r}}) = C \times \exp \left( - \frac{ \sin ^2( \, \sphericalangle ( \hat{\varvec{r}} , \hat{\varvec{\varTheta }}_i ) \, )}{\delta _{\text {max}} \times \delta _{\text {min}}} \right) . \end{aligned}$$
    (11)

    In the log-likelihood ratio, the same value for the contribution \(f_i\) as in Eq. (5) is chosen to evaluate solely the difference between an elliptically shaped and a symmetrical pattern:

    $$\begin{aligned} \log ({\mathcal {L}}_i^{\mathcal {R}})&= \sum _j^{N_\text {tot}} \log \left[ \, | f_i | \times {\mathcal {E}} (\hat{\varvec{\varTheta }}_j) \, {\mathcal {G}}_i(\hat{\varvec{\varTheta }}_j) \right. \nonumber \\&\left. \quad + (1 - | f_i |) \times {\mathcal {E}} (\hat{\varvec{\varTheta }}_j) \, \right] . \end{aligned}$$
    (12)

    During the TensorFlow fit, the gradient of \(f_i\) is computed only with respect to the signal contribution in Eq. (5) in order to prevent active adaption of the reference model \({\mathcal {R}}\).

Examples of the normalized probability density function of the elliptically shaped signal part \({\mathcal {S}} (\hat{\varvec{r}})\), the symmetric reference part \({\mathcal {G}} (\hat{\varvec{r}})\), and the geometrical exposure \({\mathcal {E}} (\hat{\varvec{r}})\) are visualized in Fig. 3.

For the objective function of the fit, each cosmic ray contributes with a separate log-likelihood ratio \(\text {ts}_i\) as:

$$\begin{aligned} \text {ts}_i= 2 \times [\log ({\mathcal {L}}_i^{\mathcal {S}}) - \log ({\mathcal {L}}_i^{\mathcal {R}}) ]. \end{aligned}$$
(13)

For an isotropic arrival distribution, each individual \(\text {ts}_i\) follows a \(\chi ^2\) distribution with the degrees of freedom corresponding to the free fit parameters according to Wilks’ theorem [26]. Therefore, the corresponding sum of all individual test-statistic contributions is Gaussian distributed according to the central limit theorem. The latter statement holds only if the individual contributions are independent of each other, which is not entirely the case in this application since the density functions of neighboring cosmic rays overlap. Nevertheless, it has been explicitly checked in Monte Carlo simulations of isotropic arrival distributions that the summed test statistics \(\sum _i \text {ts}_i\) approximately follows a Gaussian distribution. Thus, the average test statistic \(\langle \text {ts} \rangle \) provides a well-defined metric to evaluate the pattern alignment over the entire sphere:

$$\begin{aligned} \langle \text {ts} \rangle = \frac{1}{N} \sum _i^{\text {N}} \text {ts}_i. \end{aligned}$$
(14)

To maximize Eq. (14), the negative average test statistic is chosen as the objective function for the gradient descent. Additionally, in the case of the JF12 GMF initialization of the tangent vector field, the objective term F from Eq. (4) is added. Thus, the total objective function J exhibits one hyperparameter \(\lambda _F\) which represents the confidence in the GMF model:

$$\begin{aligned} J = - \langle \text {ts} \rangle + \lambda _F \cdot F. \end{aligned}$$
(15)

For the gradient descent, the fit parameters \(f_i\) and \(a_\ell ^m\) are adapted simultaneously by calculating the derivative of the objective function J with respect to the parameters. As an optimizer for the minimization problem, RMSProp [27] is used, which supports an efficient optimization for adaptive parameters of different magnitudes. This is realized by allowing each adaptive parameter to have a separate step size which is increased (decreased) depending on a consistent (inconsistent) direction of the gradient with respect to the parameter in two consecutive update steps. For stability reasons, the optimizer was complemented by additional conditions for the learning rate adaption.

2.3 Normalization

The normalization factor C in Eq. (7) is generally determined by Eq. (6) and requires a numerical integration for a non-uniform geometrical exposure \({\mathcal {E}} (\hat{\varvec{r}})\). Since the exposure \({\mathcal {E}} (\hat{\varvec{r}}) = {\mathcal {E}} (\delta (\hat{\varvec{r}}))\) depends only on the equatorial declination \(\delta \) [25], the surface integral of the term \({\mathcal {E}} (\delta ) \, {\mathcal {S}}_i (\hat{\varvec{r}})\) depends only on the center direction \(\hat{\varvec{\varTheta }}_i\) of the ellipse and its relative orientation \(| \alpha _i |\) towards the local spherical unit vector \(\hat{\varvec{e}}_\delta \) in equatorial coordinates. As \(\hat{\varvec{u}}_i\) determines the orientation of the ellipse \({\mathcal {S}}_i (\hat{\varvec{r}})\), the angle \(\alpha _i\) is defined by \(\cos \alpha _i = \hat{\varvec{u}}_i \cdot \hat{\varvec{e}}_\delta \).

Fig. 4
figure 4

Integral 1/C of exposure-folded Gaussian probability distribution as a function of the equatorial declination \(\delta \) for the center direction of the ellipse. The semi-major and semi-minor axis are \(\delta _\text {max}=30^\circ \) and \(\delta _\text {min}=10^\circ \), respectively. Colored lines denote different ellipse orientations \(\alpha \). The black dashed line indicates the probability density function value of the normalized exposure \({\mathcal {E}}(\delta )\)

Full size image

The exposure-weighted integral for an ellipse with dimensions \((\delta _\text {max}, \delta _\text {min}) = (30^\circ , 10^\circ )\) is visualized as a function of the equatorial declination \(\delta \) of the center direction and for different values of \(\alpha _i\) in Fig. 4. It can be seen that the inverse normalization factor 1/C approaches the probability density function value of the exposure \({\mathcal {E}}(\delta )\) for intermediate equatorial declinations \(-60^\circ< \delta < 20^\circ \). In the border regions of the geometrical exposure, the integral is larger than the respective exposure values due to the extent of the elliptical distribution. For the gradient descent, the normalization factor was calculated for a grid of equatorial declinations \(\delta \) and orientations \(\alpha \) and then interpolated linearly between the grid points.

3 Benchmark simulations

The sensitivity of the COMPASS fit method is evaluated on two different benchmark simulations. First, a proof of concept is given on the basis of a simple four-source model where each source contributes an equal number \(N_s\) of cosmic rays. This first benchmark simulation is designed to provide a basic understanding of the mechanism for uncovering directions of cosmic ray deflections. Here, we concentrate on light to medium nuclear charges which constitute a subset of actual composition measurements. The second benchmark simulation provides an astrophysical scenario in the vicinity of typical scenarios discussed today, given current data. Here, the simulated charges on Earth are in agreement with the measurements by the Pierre Auger Collaboration [34]. The simulation extends an astrophysical simulation of an extragalactic source population [28] by deflections in the GMF and the observational exposure \({\mathcal {E}}(\hat{\varvec{r}})\). Therefore, given a certain source density, it provides the most reasonable estimate of the sensitivity.

Concretely, the deflections in the GMF in both benchmark simulations are simulated as follows: first, the cosmic rays with magnetic rigidity \(R = E / Z\) are propagated through the regular component of the JF12 model using a magnetic field lens [6, 30]. Next, a rigidity-dependent Gaussian smearing of \(\delta = \delta _0 \times Z/E \,\)[EeV] rad is applied. This procedure fairly approximates complete three-dimensional simulations of the turbulent and striated fields of current magnetic field models as these induce only a beam widening, if the coherence length is much smaller than the propagation distance [30]. This case is expected, as described e.g. in [31] where they adopt coherence lengths of 30–100 pc. They argue that the relevant extent corresponds to the typical size of a supernova remnant, and is thus much smaller than the 20 kpc radius of the Galactic magnetic field model. Since the sources of UHECRs are not yet identified clearly, we cannot obtain the size of the turbulent smearing \(\delta _0\) from UHECR observations themselves and instead rely on current GMF models. The turbulent component of the original JF12 model [21, 32] is likely overestimated as shown by measurements of the Planck collaboration [14, 33]. Therefore we choose \(\delta _0 = 0.5\) which is around \(1/2 - 1/3\) of the median scattering that the JF12 random and striated fields induce [31, 32]. For higher values of \(\delta _0\) we expect a reduced sensitivity with respect to our results presented in Sect. 6.

Both simulations mimic the current data set of the Pierre Auger Observatory for anisotropy studies (e.g. [17]) above energies of 40 EeV with a total event number of \(N_\text {tot} = 1119\) and zenith angles up to \(80^{\circ }\).

3.1 Benchmark 1: distinct source scenario

The first benchmark simulation consists of UHECRs with energies E that follow the parameterized power law from [29] above an energy threshold of 40 EeV. The nuclear charges are assumed to be energy-independent and uniformly distributed between \(Z=1\) and \(Z=8\). While \(N_s\) out of the \(N_\text {tot} = 1119\) UHECR events are assigned to each of four different, randomly placed point-like sources in the sky, the remaining cosmic rays are distributed isotropically, following the geometrical exposure \({\mathcal {E}}(\hat{\varvec{r}})\) of the experiment. A visualization of the arrival directions of this step is shown for \(N_s=20\) source events as black circles in the upper panel of Fig. 5.

Fig. 5
figure 5

Construction of benchmark 1 scenario consisting of four distinct sources, in the Galactic coordinate system. Each of the four sources (black star symbols) contributes \(N_s = 20\) source events with color coded charge number Z. (Top) Construction of implementing GMF uncertainties for the source events. (Bottom) Complete simulations of 80 source events and 1039 isotropically distributed arrival directions

Full size image

If the JF12 GMF initialization is chosen in the fit method, we additionally apply a shift of arrival directions to mimic the uncertainties in the GMF model. Here, the entire pattern of source m is modified by a spherical angle \(\tilde{\varPsi }_m\) performed by a rotation of the individual cosmic-ray arrival directions \(\hat{\varvec{\varTheta }}_i\) from source m around its direction \(\hat{\varvec{r}}_m\) by the angle \(\tilde{\varPsi }_m\). The construction of this rotation is sketched by the dotted line in the top panel of Fig. 5. For the four sources, spherical angles \(\tilde{\varPsi }_m\) are selected in good accordance with uncertainties between existing GMF models [18], as examples \(\tilde{\varPsi }_1 = -30^{\circ }\), \(\tilde{\varPsi }_2 = 0^{\circ }\), \(\tilde{\varPsi }_3 = +45^{\circ }\) and \(\tilde{\varPsi }_4 = +15^{\circ }\) in order of the ascending Galactic longitude l (from right to left in the upper panel of Fig. 5). The colored symbols indicate the resulting arrival directions after this displacement. All arrival directions, consisting of 80 source events and 1039 isotropic events, are shown in the lower panel of Fig. 5.

3.2 Benchmark 2: astrophysical simulation

The second benchmark simulation is based on results obtained in a combined fit of the UHECR observables at the Pierre Auger Observatory [34] and their anisotropy implications for given source densities following [28]. We use the source parameters as given in table 8 and 9 in that work (spectral index \(\gamma =+0.87\)) for the CTG model, as that configuration is also used for our simulations. Here, source candidates are uniformly distributed in the universe following a source density \(\rho _S\) which results in anisotropies due to attenuation effects during the propagation. We test different source densities between \(10^{-1}\)/Mpc\(^{3}\) and \(3\times 10^{-3}\)/Mpc\(^{3}\), a region in which e.g. the typical density of spiral galaxies lies at \(10^{-2}\)/Mpc\(^{3}\) [35]. For more rare source candidates, the density is smaller and the anisotropy measured on Earth hence larger and thus easier to detect.

Fig. 6
figure 6

Example of benchmark 2 scenario consisting of an astrophysical simulation including deflection in the JF12 model for the GMF. Cosmic rays originate from uniformly distributed sources of density \(\rho _S = 10^{-2}\) Mpc\(^{-3}\) and are attenuated by extragalactic photon fields. Gray shaded events denote cosmic rays which originate from a source with at least three contributing events. For this specific source distribution the strongest source contributes 30 cosmic rays, which corresponds to the median value of 1000 simulated universes for this source density. The star symbols denote source directions where the size is proportional to the cosmic-ray contribution. The skymap is shown in the Galactic coordinate system

Full size image

The deflection in the GMF is applied in the same way as for the benchmark 1 scenario. The relative arrival probability for different extragalactic directions and rigidities caused by the GMF (e.g. [22]) are accounted for. The relative observation probability resulting from the geometrical exposure of the observatory is likewise accounted for. Figure 6 shows the resulting arrival directions of the benchmark 2 simulation for a source density of \(\rho _S = 10^{-2}\) Mpc\(^{-3}\). The circular symbols denote cosmic-ray arrival directions, the color scale corresponds to the nuclear charge Z, and the gray shaded events originate from sources which contribute at least three cosmic rays. One can see that patterns may occur in multiple regions of the sky with strongly varying event contributions attributable mainly to the source distance. Some sources situated outside the visible sky of the observatory (e.g. source at Galactic coordinates \(l \approx 110^\circ \) and \(b \approx 50^\circ \)) still contribute a substantial fraction of cosmic rays due to coherent deflection in the GMF.

Additionally, if the tangent vector field is initialized as JF12 GMF, again, an uncertainty angle \(\tilde{\varPsi }\) for the GMF is simulated. To conserve consistent deflection patterns of sources in similar directions of the sky, the uncertainty \(\tilde{\varPsi }(\hat{\varvec{r}}) = \tilde{\varPsi }_a \, \hat{\varvec{r}} \cdot \hat{\varvec{d}}_i\) is modeled as a dipolar function with amplitude \(\tilde{\varPsi }_a = 45^{\circ }\) and random direction of the dipole maximum \(\hat{\varvec{d}}_i\) for each simulated universe i.

This simulation of the UHECR universe exhibits only one single free parameter, the source density \(\rho _S\), which directly determines the degree of anisotropy in the arrival directions. The higher the source density, the more sources are within a horizon where attenuations do not play an important role and, therefore, the more isotropic the sky is.

4 Consistency test for recovering the galactic magnetic field deflections

In this section, proof of concept is provided by showing that the orientation of patterns can be correctly reconstructed based on the benchmark 1 simulation from Sect. 3.

During the minimization process, the modification angle \(\varPsi (\vartheta , \varphi )\) rotates the ellipses of the signal hypothesis of Eq. (1) such that they align with elongated patterns in the cosmic-ray arrival direction distribution. For the JF12 GMF initialized vector field \(\hat{\varvec{u}}_0\), the angle \(\varPsi (\vartheta , \varphi )\) corresponds directly to a correction of the JF12 model in sky regions where a significant pattern is found. Hence, for the benchmark 1 simulation, the final angle \(\varPsi _i \equiv \varPsi (\vartheta _i, \varphi _i)\) of cosmic rays i which originate from one of the simulated sources m is expected to approach the simulated uncertainty \(\tilde{\varPsi }_m\). Here, the necessary correction is modeled by the spherical angles \(\tilde{\varPsi }_m\) for the four sources where \(\tilde{\varPsi }_1 = -30^{\circ }\), \(\tilde{\varPsi }_2 = 0^{\circ }\), \(\tilde{\varPsi }_3 = +45^{\circ }\) and \(\tilde{\varPsi }_4 = +15^{\circ }\) are chosen (cf. Sect. 3). Note that a non-linear deflection behavior in the GMF may disturb the correct values of \(\tilde{\varPsi }_m\). This effect is particularly strong for cosmic rays with a low rigidity \(E_i/Z_i\), i.e. for high absolute deflection angles with respect to their source.

Here, for the first application of the fit, the order of the spherical harmonics expansion of Eq. (1) is defined as \(k = 5\), which corresponds to 36 free fit parameters. In this case, modifications of the GMF model can be performed coherently in sky regions that have angular scales above the order of \(180^\circ / k = 36^\circ \). The degree of modification \(\varPsi \) itself is constrained by the hyperparameter in the objective function (15) where a value of \(\lambda _F = 1\) is chosen for this purpose. For the ellipse geometry, values of (\(\delta _\text {max}, \, \delta _\text {min}) = (10^\circ , \, 5^\circ \)) are chosen in Eq. (7). Furthermore, the Gaussian reference model from Eq. (12) was selected for the likelihood ratio where the effective Gaussian width is \(\sqrt{10^\circ \times 5^\circ } \approx 7.1^\circ \).

The fitted modification function \(\varPsi (\theta , \varphi )\) is visualized in Fig. 7 together with the cosmic rays that originate from the simulated source candidates. As an overall impression the color code in the vicinity of the source candidates m agrees with the simulated uncertainties \(\tilde{\varPsi }_m\). To quantify the method’s reconstruction abilities, for each source m the fitted \(\varPsi _i\) for the 10 closest cosmic rays that originate from the source are averaged. The corresponding averaged values are \(\langle \varPsi _1 \rangle = (-32 \pm 1)^{\circ }\), \(\langle \varPsi _2 \rangle = (6 \pm 2)^{\circ }\), \(\langle \varPsi _3 \rangle = (+47 \pm 5)^{\circ }\) and \(\langle \varPsi _4 \rangle = (+15 \pm 5)^{\circ }\), which are in good agreement with the simulated uncertainties \(\tilde{\varPsi }_m\).

Fig. 7
figure 7

Visualization of the modification angles \(\varPsi (\vartheta , \varphi )\) (color-coded) which is defined relative to the JF12 model (JF12 GMF initialization for \(\hat{\varvec{u}}_0\)) after a fit to benchmark 1 scenario with \(N_s = 20\) source events. Gray circles denote cosmic rays that originate from one of the four sources and the dotted lines mark contours of \(15^{\circ }\)-spacing in \(\varPsi \)

Full size image

The next step is to investigate if orientations of patterns as simulated with the JF12 model can also be captured without including information on the explicit GMF model. For this purpose, we chose the Galactic meridians initialization of \(\hat{\varvec{u}}_0\) where initial ellipse orientations are aligned with the local spherical unit vector \(\hat{\varvec{e}}_\vartheta \) of the latitude in the Galactic coordinate system. Since no information on the simulated GMF is included, the penalization factor of Eq. (15) is not required and is therefore set to \(\lambda _F = 0\). Thus, the tangent vector field can be rotated by the angle \(\varPsi (\theta , \varphi )\) without constraint. For this setup, the degree of the spherical harmonics expansion was decreased to \(k = 4\) to avoid rapid changes of \(\varPsi \) on small angular scales. The order k of the spherical harmonics expansion is the most challenging free parameter since the optimal choice depends on the angular scale of domains with a coherent GMF deflection. While more complex patterns can generally be fitted with an increasing order of k, these structures are more difficult to interpret.

Fig. 8
figure 8

Visualization of the tangent vector field \(\hat{\varvec{u}}(\varTheta _i)\) (black lines) of Eq. (2) using the Galactic meridians initialization of \(\hat{\varvec{u}}_0 (\hat{\varvec{r}})\). Black star symbols show the simulated source directions while the red circular symbols mark the \(N_s = 20\) events per source which have been displaced by the JF12 model without simulated GMF uncertainties

Full size image

A visualization of the fitted tangent vector field \(\hat{\varvec{u}}(\varTheta _i)\) for each individual cosmic-ray arrival direction \(\varTheta _i\) is presented in Fig. 8 where the source events are highlighted in red. All four deflection patterns in this simulation were successfully captured by an alignment of the tangent vector field \(\hat{\varvec{u}}(\hat{\varvec{r}})\) along the local track of source events. There is also structure visible in sky regions without source contribution where fluctuations of isotropically distributed arrival directions were connected along their most prominent patterns. In the sky region at Galactic coordinates \((l, b) \approx (-90^\circ , -75^\circ )\) the isotropic fluctuation was even strong enough to rotate the initialized tangent vector field \(\hat{\varvec{u}}_0\) by up to \(95^\circ \). This suggests that arbitrarily oriented alignment patterns of cosmic-ray arrival directions can be captured even when oriented orthogonally with respect to the chosen initialization \(\hat{\varvec{u}}_0\).

To assess whether a fitted pattern is caused by isotropic fluctuations, the individual cosmic-ray test statistics from Eq. (13) have to be compared to those found in isotropic skies. These sensitivity studies are presented in Sect. 6.

5 Reference model of the likelihood ratio

In this section, two different choices of the reference model (cf. Sect. 2) for the likelihood ratio as defined in Eq. (13) are studied: the isotropic model \({\mathcal {E}}\) which follows the geometrical exposure of an experiment, and the Gaussian model \({\mathcal {G}}_i\) with identical signal contributions \(f_i\) as assigned to the elliptically shaped signal model \({\mathcal {S}}_i\). Again, the ellipse geometry is defined as (\(\delta _\text {max}, \, \delta _\text {min}) = (10^\circ , \, 5^\circ \)), the initialization for \(\hat{\varvec{u}}_0\) is Galactic meridians, the spherical harmonics order is \(k = 4\), and the hyperparameter \(\lambda _F = 0\). To obtain an impression of the performance over the sky, it is useful to investigate the individual test statistics \(\text {ts}_i\) defined in Eq. (13) as well as the anticipated signal contribution \(f_i\) from Eq. (5).

Fig. 9
figure 9

Comparison of fit performance for (left) isotropic \({\mathcal {E}}\) and (right) Gaussian reference model \({\mathcal {G}}_i\), in Galactic coordinates. (Top) Benchmark 1 simulation from Sect. 3 with \(N_s=20\) events per source and deflection in the JF12 model. (Bottom) Gaussian overdensity with a cluster of \(N_s = 20\) events distributed around the star symbols with a Gaussian width of \(10^\circ \) and without GMF deflection. The color code of each event i corresponds to the individual test statistic \(\text {ts}_i\). The size of the circular symbols is proportional to the signal contribution \(f_i\) with identical normalization among all four figures

Full size image

5.1 Isotropic reference model \({\mathcal {E}}\)

The resulting test statistics \(\text {ts}_i\) after the fit using the isotropic reference model (cf. Eq. (10)) is shown in the left panel of Fig. 9. The fit results of the benchmark 1 simulation are displayed in the upper panel with \(N_s = 20\) source events in each of the four patterns. The fitted signal contribution \(f_i\) for the events i is proportional to the size of the circular symbols where a common normalization among all four figures is chosen.

Clearly, for benchmark 1 scenario in the upper panel, the patterns produced by the four simulated sources exhibit cosmic rays with a substantially larger test statistic \(\text {ts}_i\) compared to the isotropically distributed background events. Accordingly, the anticipated signal contribution \(f_i\) of the source events is larger, as shown by the size of the markers. Some local clusters of events with test statistics \(\text {ts}_i> 0\) can also be found in the isotropically distributed background events, however, with considerably smaller values of both the test statistic \(\text {ts}_i\) and the signal contribution \(f_i\). The average test statistic from Eq. (14) is \(\langle \text {ts} \rangle \approx 1.6\). The highest signal fractions reach values of about \(f_i=1.4 \%\), which is in good agreement with the 20 of 1119 injected signal cosmic rays per source. Note that the complete signal contribution of \(20 / 1119 \approx 1.8 \%\) is not necessarily reached even for the innermost cosmic ray of the pattern due to fluctuations in the isotropic background and the ellipse’s limited extent of \(10^{\circ }\) in the semi-major axis, which is mostly less than the extent of the pattern.

To assess the impact of solely overdense but not elongated structures on the test statistic, a new simulation is studied which again consists of four sources each emitting 20 cosmic rays. Instead of simulating deflections in the GMF, the source events are drawn from a Fisher distribution [36] of a width of \(10^{\circ }\) centered on the direction of the source. Here, to enable a better comparison between both scenarios, the source directions were approximately centered within the resulting arrival patterns of the benchmark 1 simulation.

As shown in the lower panel of Fig. 9, the method also responds with high individual test statistics \(\text {ts}_i\) and anticipated signal contributions \(f_i\) due to the event excess of \(N_s = 20\) relative to an isotropic expectation. However, for three of the four Fisher distributions the resulting test statistics are much smaller than in the case of the benchmark 1 simulation.

5.2 Gaussian reference model \({\mathcal {G}}_i\)

The right panel of Fig. 9 again shows the individual test statistic \(\text {ts}_i\) (color coded) and the fitted signal contributions \(f_i\) (size of circles) for the Gaussian reference hypothesis \({\mathcal {G}}_i\) as defined in Eq. (12). For the benchmark 1 simulation in the upper panel, the anticipated signal contribution \(f_i\) is approximately equal to the case where the isotropic reference model \({\mathcal {E}}\) was chosen, as can be estimated from the size of the markers. However, while in the case of the isotropic reference model both the event excess and the elongation of the structure contributed to the test statistic, for the Gaussian reference model only the latter information can be used. Thus, on the one hand, the overall scale of the individual test statistics \(\text {ts}_i\) is much smaller, as reflected by the color scale. Therefore, the average test statistic of Eq. (14) drops to \(\langle \text {ts} \rangle \approx 0.8\). On the other hand, since there is no sensitivity to overdense regions, some of the patterns that were caused by fluctuations of isotropically distributed background events are no longer visible. Thus, the purity of detected patterns is increased compared to when the isotropic reference model was used.

Again, the response to Gaussian overdensities is assessed in the bottom panel of Fig. 9 with four Fisher-distributed event clusters of \(10^\circ \) Gaussian width. Since the Gaussian-shaped event structures are well described by the Gaussian reference hypothesis \({\mathcal {G}}_i\), there is a significant loss in the test statistic for the overdense sky regions compared to when the isotropic reference \({\mathcal {E}}\) is used. In the vicinity of three Gaussian event clusters, there is only barely more fitted signal contribution \(f_i\) compared to the remaining sky. The individual test statistics \(\text {ts}_i\) visibly deviates from natural isotropic fluctuations only for the events of one of the Gaussian clusters, namely at coordinates \((l, b) = (-20^\circ , -20^\circ )\).

As there are already known event excesses in UHECR data, e.g. in data of the Pierre Auger Observatory for this energy threshold (e.g. [24, 37]), there is a risk of detecting these features again rather than new elongated structures when using the isotropic reference model \({\mathcal {E}}\). Therefore, in the following we use the Gaussian-like reference model \({\mathcal {G}}_i\) where the effects of overdensities are mostly canceled out by the likelihood ratio.

6 Sensitivity studies

In this section we investigate the sensitivity of the COMPASS method with respect to its ability to reject isotropic distributions of cosmic-ray arrival directions. According to the findings from Sect. 5, for the following subsections the Gaussian reference model \({\mathcal {G}}_i\) (cf. Eq. (12)) is chosen in the likelihood ratio. In addition, following Sect. 4 and the studies in Sect. 6.3, the tangent vector field \(\hat{\varvec{u}}_0\) is initialized along the Galactic meridians – i.e. \(\hat{\varvec{u}}_0\) is equal to the local spherical unit vector \(\hat{\varvec{e}}_\theta \). Therefore, the penalization term F in Eq. (15) is removed by setting \(\lambda _F = 0\). For a comparison of the sensitivity with a more classical analysis to search for elongated structures refer to [38].

Fig. 10
figure 10

Distribution of individual test statistics \(\text {ts}_i\) as obtained in the benchmark 1 simulation (cf. upper right panel of Fig. 9). The red (gray) histogram shows the contribution of the 80 source (isotropically distributed) events. The dashed vertical black line denotes the average test statistic \(\langle \text {ts} \rangle \) of all events in the sky

Full size image

6.1 Sensitivity for distinct source scenario

The distribution of individual test statistics \(\text {ts}_i\) as obtained in the benchmark 1 simulation from Sect. 3 is presented in Fig. 10. As already suggested by the upper right panel of Fig. 9, most of the events that exhibit high test statistics are attributed to one of the four sources. In total, more than half of the source events show test statistics larger than 3, which, in turn, is only reached for about \(5 \%\) of the isotropic events. Instead, the isotropic distribution peaks close to zero, with about \(50\%\) of events exhibiting test statistics smaller than \(10^{-5}\). One can of course find a statistical measure to reject isotropy based on the evaluation of events with a high test statistic, i.e. based on the tail of the distribution in Fig. 10. However, it was found that the average test statistic \(\langle \text {ts} \rangle \) provides the most stable measure for various simulation setups.

Fig. 11
figure 11

Distribution of average test statistic \(\langle \text {ts} \rangle _\text {iso}\) for \(10^4\) isotropic skies (gray) compared to the values obtained in the benchmark 1 simulation (vertical red lines) with \(N_s = 5\) (solid), \(N_s = 10\) (dashed), \(N_s = 15\) (dash-dotted), and \(N_s = 20\) (dotted) source events

Full size image

In the next step, we evaluated the average test statistic \(\langle \text {ts} \rangle \) for different numbers of source events \(N_s = 5, \, 10, \, 15,\) \(\, 20\) in the benchmark 1 simulation. As shown in Fig. 11, the resulting values for the average test statistics are \(\langle \text {ts} \rangle = 0.23, \, 0.30, \, 0.41, \, 0.72\), respectively. To calculate the chance probability \(p_\text {val}\) of obtaining these average test statistics from an isotropic arrival-direction distribution, the method is additionally applied to \(5 \times 10^4\) isotropic realizations of the sky which follow the geometrical exposure of the observatory. The distribution of the average test statistics \(\langle \text {ts} \rangle _\text {iso}\) for isotropic skies is shown as a gray histogram in Fig. 11. While there is no isotropic sky yielding a higher average test statistic than the scenarios with \(N_s = 20\) source events, the isotropic chance probabilities for the simulations with smaller source events are \(p_\text {val} = 0.40, \, 0.14, \, 0.003\), in order of increasing \(N_s\).

As expected from the central limit theorem (cf. Sect. 2), the gray histogram shows that the average test statistic \(\langle \text {ts} \rangle _\text {iso}\) for an isotropic arrival sky approximately follows a Gaussian distribution. Thus, the sensitivity for the scenario shown in Fig. 5 with \(N_s = 20\) source events can be estimated by fitting a Gaussian distribution to the gray histogram. In this case, the estimated chance probability is about \(3 \times 10^{-11}\) which translates to about \(6.5 \, \sigma \) standard deviations in the normal distribution.

6.2 Sensitivity for the astrophysical universe

While the previous section provided an idea of the sensitivity for comparably clear patterns with a certain signal contribution, this section evaluates the expected implication for an astrophysical universe of uniformly distributed UHECR sources. For the benchmark 2 simulations, 300 simulated universes for each of the source densities \(\rho _S = (10^{-1}, \, 3 \times 10^{-2}, \, 10^{-2}, \, 3 \times 10^{-3})\) Mpc\(^{-3}\) were investigated. The average test-statistic distribution as obtained from the fit exhibits a comparably large spread, which is consistent with the fluctuations in the degree of anisotropy. The median and 68 percentiles of the average test statistics for the four source densities are \(\langle \text {ts} \rangle = 0.27^{+0.16}_{-0.08}, \, 0.38^{+0.39}_{-0.11}, \, 0.54^{+0.31}_{-0.22},\, 0.86^{+1.92}_{-0.30}\), respectively, as visualized in Fig. 12. As expected, the test statistic increases with decreasing source density as the arrival scenarios become increasingly anisotropic.

Fig. 12
figure 12

Distribution of average test statistics \(\langle \text {ts} \rangle _\text {iso}\) for \(5 \times 10^4\) isotropic skies (gray) compared with the values obtained in the benchmark 2 simulation (vertical lines) with source densities of \(\rho _S = 10^{-1} /\text {Mpc}^3\) (blue), \(\rho _S = 3 \times 10^{-2} /\text {Mpc}^3\) (orange), \(\rho _S = 10^{-2} /\text {Mpc}^3\) (green), and \(\rho _S = 3 \times 10^{-3} /\text {Mpc}^3\) (red). The small triangles denote the \(68\%\) quantiles

Full size image

For the isotropic chance probability \(p_\text {val}\), the average test statistic is again compared to the fit results for the \(5 \times 10^4\) isotropic realizations which are shown as a gray histogram in Fig. 12. For the source densities of \(10^{-1}\) Mpc\(^{-3}\) and \(3 \cdot 10^{-2}\) Mpc\(^{-3}\) the isotropic chance probability can be directly determined by the fraction of the gray distribution that is above the corresponding test-statistic values. Here, the chance probabilities yield \(p_\text {val} = 0.229, \, 0.013\) for the two source densities respectively. For the smaller source densities of \(10^{-2}\) Mpc\(^{-3}\) and \(3 \times 10^{-3}\) Mpc\(^{-3}\), the isotropic chance probability can again be estimated by parameterizing the null hypothesis with a Gaussian distribution. In this case, the estimated values are \(8.3 \times 10^{-6}\) and \(1.4 \times 10^{-17}\), respectively, which correspond to a deviation of \(4.3 \, \sigma \) and \(8.5 \, \sigma \) standard deviations in the normal distribution. Thus, in the case of a result on data that is compatible with an isotropic distribution, the density of UHECR sources for this astrophysical model can be estimated.

Fig. 13
figure 13

Example realization of arrival directions from the benchmark 2 simulation with a source density of \(\rho _S = 10^{-2}\) Mpc\(^{-3}\). The strongest of the sources (red star symbols) at Galactic longitude \(l \approx -133^{\circ }\) and Galactic latitude \(b \approx 14^{\circ }\) contributes with 29 cosmic rays. Sizes of the cosmic-ray events (colored circles) correspond to the energy E, the color code to the individual test statistics \(\text {ts}_i\), and the short black lines indicate the tangent vector field \(\hat{\varvec{u}}_0\), i.e. the orientation of the major-axis of the ellipses

Full size image

Figure 13 shows arrival directions in the sky region around the strongest individual test statistic \(\text {ts}_i\) for the scenario that exhibits the median average test statistic \(\langle \text {ts} \rangle \) out of the 300 simulations with a source density of \(10^{-2}\) Mpc\(^{-3}\). Here, the strongest individual test statistic is \(\text {ts}_i= 8.6\) and the corresponding cosmic-ray event (yellow point in the center of the sky patch) is part of the pattern from the source at Galactic coordinates \(l \approx -133^\circ \) and \(b \approx 14^\circ \). Contributing with a total of 29 events, this is the strongest source in this realization. The orientation of the tangent vector field \(\hat{\varvec{u}}(\vartheta , \varphi )\) (short black lines) additionally suggests that the alignment works reasonably well even for patterns that are only separated by about \(20^\circ \). Thus, the vector field \(\hat{\varvec{u}}(\vartheta , \varphi )\) is expected to provide an adequate coherent description of the deflection in the GMF for sufficiently strong signals.

6.3 Optimization of free parameters

In this subsection we evaluate the impact of the free parameters more profoundly based on the astrophysical benchmark 2 simulation with a source density of \(\rho _s = 3 \times 10^{-2}\) Mpc\(^{-3}\) from Sect. 3. Firstly, the initialization method of the tangent vector field \(\hat{\varvec{u}}_0\) and accordingly the free parameter \(\lambda _F\) are addressed. Secondly, the impact of the ellipse geometry, namely the semi-major and semi-minor axes, on the performance is studied.

6.3.1 Confidence in the JF12 GMF initialization

The initialization of the tangent vector field \(\hat{\varvec{u}}_0\) according to the predictions from the JF12 model (JF12 GMF) is visualized in the top panel of Fig. 2. As pointed out in Sect. 2, depending on the reliability of the GMF model it may be beneficial to constrain the allowed deviations \(\varPsi (\vartheta , \varphi )\) with Eq. (4), since this reduces high test statistics from fluctuations in isotropic arrival distributions. Therefore, the isotropic chance probability \(p_\text {val}\) is investigated as a function of the free objective parameter \(\lambda _F\) in Eq. (15) for two reasonable estimates of the uncertainties in GMF models.

The first estimate is obtained by simulating the deflection in the GMF with the model of Pshirkov et al. using an antisymmetric disk field (PT11-ASS) [39] instead of the JF12 model used in Sect. 3. The second estimate is given by a modification of the JF12 model with dipolar distributed modification angles of amplitude \(\varPsi _a=45^\circ \) as described in Sect. 3. The isotropic chance probabilities for both estimates are presented in Fig. 14 as a function of \(\lambda _F\). For high values of \({\mathcal {O}}(\lambda _F) = 10\), the tangent vector field is too stiff and the test statistic is therefore obtained for ellipses which are not aligned with the simulated structures. There is a minimum for both assumptions of GMF uncertainties located consistently at values about \(\lambda _F \approx 0.5\). For an entirely flexible tangent vector field, i.e. for the parameter \(\lambda _F = 0\), the additionally found patterns from isotropic skies reduce the sensitivity slightly; however, the overall isotropic chance probability is still of the same order.

Fig. 14
figure 14

Scan of hyperparameter \(\lambda _F\) for the confidence in the assumed GMF model (here JF12). The initialization of the tangent vector field \(\hat{\varvec{u}}_0\) follows the JF12 GMF procedure and the isotropic chance probability \(p_\text {val}\) is derived for the benchmark 2 simulation with a source density of \(\rho _S = 3 \times 10^{-2}\) Mpc\(^{-3}\) from Sect. 3. The red crosses show deflections in the PT11-ASS model and the blue crosses show deflections in the JF12 model, which is modified by a dipolar modulated uncertainty angle with amplitude \(\varPsi _a\)

Full size image

Since the uncertainties of the GMF models might be even higher than assumed here and particularly uncertain in the Galactic disk region, the advantage of a hyperparameter \(\lambda _F > 0\) may be even smaller. Therefore, in the following ellipse geometry investigation the penalization term F is canceled in Eq. (15) and the Galactic meridians initialization is utilized which features a symmetry with respect to the Galactic disk.

6.3.2 Ellipse geometry

Here, we assess the impact of the ellipse geometry, the semi-major axis width \(\delta _\text {max}\), and the semi-minor axis width \(\delta _\text {min}\), for the astrophysical benchmark 2 simulation. For the two-dimensional scan of the widths, angular bins of (\(3^\circ \), \(5^\circ \), \(7^\circ \), \(10^\circ \), \(15^\circ \), \(20^\circ \)) were chosen where the condition \(\delta _\text {max} > \delta _\text {min}\) is required by design. Thus, there are 15 different scanned ellipse geometries. To calculate the isotropic chance probability \(p_\text {val}\), the analysis is applied to a total of \(10^4\) isotropic skies for each geometry.

Fig. 15
figure 15

Scan of the ellipse geometry (\(\delta _\text {max}\), \(\delta _\text {min}\)) for the Galactic meridians initialization and \(\lambda _F = 0\) evaluated on the astrophysical benchmark 2 simulation with source density of \(\rho _S = 3 \times 10^{-2}\) Mpc\(^{-3}\). The isotropic chance probability is calculated with a total of \(10^4\) isotropic skies for each of the ellipse constellations

Full size image

The resulting median chance probabilities \(p_\text {val}\) of 300 sky realizations are displayed in Fig. 15 where each of the five segments indicate one of the semi-major axes \(\delta _\text {max}\). Generally, larger and less elongated ellipse sizes are beneficial for the sensitivity of the COMPASS method. Consequently, the largest ellipse with values of \(\delta _\text {max} = 20^\circ \) and \(\delta _\text {min} = 15^\circ \) for the semi-major and semi-minor axes, respectively, yields the lowest isotropic chance probability of \(p_\text {val} = 2.8 \times 10^{-3}\). This result is significantly better than the previously considered ellipse geometry of \((\delta _\text {max},\, \delta _\text {min}) = (10^\circ , 5^\circ )\), which exhibits a chance probability of \(p_\text {val} = 1.2 \times 10^{-2}\) in the same benchmark scenario. However, the specific behavior of the sensitivity for the various ellipse geometries may be characteristic for the simulation setup. Since the angular scales of existing structures are unknown for an application to data, it is suggested to scan the ellipse geometry in a reasonable range.

7 Conclusion

In this work we investigated a novel approach to search for structures in the arrival directions of UHECRs induced by cosmic magnetic fields. A dynamic vector field tangential to the local celestial sphere is utilized to fit the orientation of elongated patterns. Thus, elliptically shaped density functions are aligned by the vector field and evaluated in a likelihood ratio with a circular reference model. This work demonstrates that the orientation of the directional deflections of the GMF is detectable by faint signatures of simulated UHECR sources. The sensitivity of the method was investigated by means of an astrophysical simulation of uniformly distributed sources where UHECR nuclei are attenuated during propagation in the extragalactic universe. It was shown that the hypothesis of isotropically distributed arrival directions can be excluded with more than \(4 \, \sigma \) Gaussian significance if a maximum spatial density of UHECR sources of \(\rho _S = 10^{-2}\) Mpc\(^{-3}\) is assumed.