Effects of the Galactic Magnetic Field on the UHECR Correlation Studies with Starburst Galaxies

Ryo Higuchi; Takashi Sako; Toshihiro Fujii; Kazumasa Kawata; Eiji Kido

doi:10.3847/1538-4357/acc739

1. Introduction

Cosmic rays (CRs) are high-energy nuclei that exist throughout the universe. Specifically, ultrahigh-energy cosmic rays (UHECRs) with energies of about 100 EeV (10²⁰ eV) are observed, but their origin is not yet known. There are two leading experiments that observe UHECRs: the Telescope Array (TA) experiment (Kawai et al. 2008; Sagawa 2020; it is located in the U.S.A., 39.3 deg N, 112.9 deg W) and the Auger experiment (Aab et al. 2015a; it is located in Argentina, 35.2 deg S, 69.5 deg W), which cover the sky in the northern and southern hemispheres, respectively. Thanks to a large number of UHECR events observed through these experiments in the past decade, the arrival directions of UHECRs are found to be anisotropic in the intermediate angular (i.e., ∼10° to ∼20°) scale (Abreu et al. 2012; Abbasi et al. 2014; Aab et al. 2015b), which is believed to give us keys to knowing the UHECR origins. Due to the energy loss through the photopion production or photonuclear dissociation, most UHECRs cannot propagate more than 30–100 Mpc (Greisen 1966; Zatsepin & Kuz'min 1966; the GZK limit). Consequently, as astronomical candidates of UHECR sources, it is natural to consider nearby extragalactic high-energy objects. Previous studies have investigated the correlation between the arrival directions of UHECRs and their source candidates (Abreu et al. 2007, 2010; He et al. 2016; Aab et al. 2018; Abbasi et al. 2018). Especially, recent studies (Abreu et al. 2010; Aab et al. 2018; Abbasi et al. 2018) try to explain the arrival direction of UHECRs by the weighted sum of events originating from sources and isotropic backgrounds (see also the review in Batista et al. 2019). In these studies, the flux of UHECRs (CR flux model) is composed of the source-originated flux (source-flux model) and the isotropic backgrounds (isotropic flux model). Based on a catalog of source candidates, the source-flux model is constructed as a superposition of the Gaussian-smeared angular distributions of an individual point source. Previous studies introduced the anisotropic fraction f_ani as the fraction of the source contribution to the total CR flux and the separation angular scale θ as the scale of the Gaussian smearing. The separation angular scale θ is considered to reflect the deflections and scattering by the Galactic and extragalactic magnetic fields (EGMF). The parameters (f_ani, θ) are searched to best fit the observed CR angular distribution (see Section 3 for details).

One of the possible source-flux models suggested by previous studies is that UHECRs originate from nearby starburst galaxies (SBGs). Aab et al. (2018) investigated the correlations between UHECRs observed by the Auger experiment and the CR flux models constructed with nearby extragalactic high-energy objects such as SBGs, active galactic nuclei, and gamma-ray bursts. They reported that the best correlation was found above 39 EeV with the SBG source-flux model (SBG model), and the best-fit anisotropic fraction ${f}_{\mathrm{ani}}^{\mathrm{Auger}}$ and the separation angular scale θ^Auger were estimated to be 9.7% and 12 fdg 9, respectively (we call this best-fit model the Auger best-fit model).

The TA experiment studied the correlation between the UHECR arrival directions observed in the northern sky and the source-flux model with SBGs using the best-fit parameters reported by Aab et al. (2018). It concluded that the UHECR arrival directions are compatible with both the isotropic distribution (f_ani = 0%) and the Auger best-fit model (Aab et al. 2018) with current statistics (Abbasi et al. 2018). While the nearby SBGs belong to the group of attractive candidates for the UHECR sources, it is surprising that even this best-fit source model can explain only ∼10% of the observed UHECRs. This question motivated us to study possible biases that cause an underestimate of f_ani and to estimate a realistic constraint of f_ani deduced from the observations.

The isotropic scattering approximation in previous studies does not reflect the actual structure of the GMF, which deflects UHECR trajectories in a certain direction (coherent deflection). Current analyses treating the coherent deflection as a part of isotropic scattering may result in a smaller f_ani and a larger θ than the true values. In this study, we call these systematic effects on the parameter estimation caused by the GMF "the GMF bias." To consider the GMF bias (in addition to introducing the GMF model), we take the following two components into account: (1) The dependence on the arrival direction. The coherent deflection around the Galactic center (GC) and the Galactic plane (GP) is much larger than in other regions of the sky. This dependence also affects analyses whose samples are divided into the northern sky (TA experiment) and the southern sky (Auger experiment). Generally, UHECRs observed in the southern sky should be more affected by the GMF than those in the northern sky. These effects, caused by the limitation of the sky coverage of each experiment, are not evaluated in the previous studies. (2) The dependence on the rigidity R = pc/Ze ∼ E/Ze. Here pc, E, and Ze represent the particle momentum, energy, and electric charge, respectively. The approximation is valid in the ultrarelativistic regime considered here. Because the deflection angle is proportional to 1/R of each CR, the magnetic field effect is strongly coupled with the energy spectrum and the mass composition. The Auger experiment suggests that the mass composition of UHECRs is higher at higher energy (Aab et al. 2017; Batista et al. 2019; Heinze & Fedynitch 2019).

In this study, we investigate the GMF bias in previous studies by applying a commonly used GMF model (Jansson & Farrar 2012a, 2012b) by means of Monte Carlo (MC) simulations. We generate mock event data sets assuming the true parameters $({f}_{\mathrm{ani}}^{\mathrm{true}},{\theta }^{\mathrm{true}})$ , taking coherent deflections by the GMF and a mass model into account (a mixed-mass spectrum model proposed in Heinze & Fedynitch 2019). For some samples of mock events, we demonstrate that the event arrival directions are apparently displaced from the real source directions. The size of the displacements strongly depends on the direction of the sky and the rigidity. Then, we applied a maximum likelihood analysis to the mock events in the same manner as in previous studies. In order to focus on the bias in the previous analysis (Aab et al. 2018; Abbasi et al. 2018), we fix the source-flux model to be the SBG model. Based on the analyses of the mock data sets, we discuss the biases in the parameter estimation separately in the northern and southern hemispheres and in the all-sky regions. We also develop an analysis technique to reduce the GMF bias.

2. Mock Data-set Production

The mock event data sets are generated under a set of assumptions (the SBG model, the GMF model, and the mass-dependent energy spectrum) with the flow shown in Figure 1. The construction of the CR flux models before and after considering the GMF deflection is described in Sections 2.1 and 2.2, respectively. To generate mock events taking the coherent deflections caused by the GMF into account, one needs to assign a rigidity R, i.e., energy and mass, to each event. We use a mixed-mass assumption to reflect a realistic situation as described in Section 2.3. Results with single-mass assumptions are summarized in Appendix A.1. In all cases, we generate 1000 data sets, each of which contains 4000 mock events across the whole sky (all-sky data set). To compare the data sets with the observed UHECRs by the TA and Auger experiments, we select in Section 2.4 the north-sky and south-sky data sets from the all-sky data sets, taking the sky coverage of each experiment into account.

**Figure 1.** A schematic for a mock data-set generation.
Download figure:
Standard image High-resolution image

2.1. CR Flux Model that Originated from the SBG Model

With the source model fixed to the SBG model, the CR flux model is constructed with the assumption of a two-parameter set (f_ani, θ). In Aab et al. (2018) and Abbasi et al. (2018), the CR flux model F_org( n , θ) is determined as the superposition of the von Mises-Fisher function (Fisher 1953; the Gaussian distribution on the sphere) of each source,

$\begin{eqnarray}&&{F}_{\mathrm{org}}({\boldsymbol{n}},\theta )=\displaystyle \frac{{\sum }_{i}{f}_{i}\exp ({{\boldsymbol{n}}}_{i}\cdot {\boldsymbol{n}}/{\theta }^{2})}{{\int }_{4\pi }{\sum }_{i}{f}_{i}\exp ({{\boldsymbol{n}}}_{i}\cdot {\boldsymbol{n}}/{\theta }^{2})d{\rm{\Omega }}}.\end{eqnarray} \tag{ 1 }$

Note that i indicates each SBG, and n _i and f_i mean its direction and relative flux (contribution from each source), respectively. Table 1 is the list of SBGs in our SBG model defined in Aab et al. (2018), which contains the values of f_i and n _i. The relative flux of SBGs f_i in Table 1 is determined by their continuum radio flux. The source directions and relative contributions in Table 1 are visualized in Figure 2. We can see that most SBGs are located along the supergalactic plane (SGP) and the top-4 contributions of SBGs dominate with ∼60% of the total flux.

**Figure 2.** Directions and contributions of SBGs in Table 1 from Aab et al. (2018; in equatorial coordinates). Circles show the direction of SBGs. The color and area of each marker scale their relative flux contribution to the CR flux models. The gray dots (circles) indicate the GP (the SGP).
Download figure:
Standard image High-resolution image

Table 1. Catalog of SBGs in Aab et al. (2018)

ID ^a	l (°) ^b	b (°) ^b	D (Mpc) ^c	f (%) ^d
NGC 253	97.4	−88	2.7	13.6
M82	141.4	40.6	3.6	18.6
NGC 4945	305.3	13.3	4	16
M83	314.6	32	4	6.3
IC 342	138.2	10.6	4	5.5
NGC 6946	95.7	11.7	5.9	3.4
NGC 2903	208.7	44.5	6.6	1.1
NGC 5055	106	74.3	7.8	0.9
NGC 3628	240.9	64.8	8.1	1.3
NGC 3627	242	64.4	8.1	1.1
NGC 4631	142.8	84.2	8.7	2.9
M51	104.9	68.6	10.3	3.6
NGC 891	140.4	−17.4	11	1.7
NGC 3556	148.3	56.3	11.4	0.7
NGC 660	141.6	−47.4	15	0.9
NGC 2146	135.7	24.9	16.3	2.6
NGC 3079	157.8	48.4	17.4	2.1
NGC 1068	172.1	−51.9	17.9	12.1
NGC 1365	238	−54.6	22.3	1.3
Arp 299	141.9	55.4	46	1.6
Arp 220	36.6	53	80	0.8
NGC 6240	20.7	27.3	105	1
Mkn 231	121.6	60.2	183	0.8

Notes.

^aSBG Names.^bSBG directions (galactic coordinates).^cDistances from the Earth.^dRelative flux contributions normalized by the radio flux at 1.4 GHz.

Download table as: ASCII Typeset image

An example of the SBG model F_org( n , θ = 10°) is presented in Figure 3. As expected from Table 1 and Figure 2, a small number (<10) of sources dominates the distribution. The normalized CR flux model F_norm is defined as the weighted sum of the SBG model F_org and the isotropic flux model F_iso,

$\begin{eqnarray}\begin{array}{rcl}{F}_{\mathrm{norm}}({\boldsymbol{n}},{f}_{\mathrm{ani}},\theta ) & = & {f}_{\mathrm{ani}}{F}_{\mathrm{org}}^{{\prime} }({\boldsymbol{n}},\theta )+(1-{f}_{\mathrm{ani}}){F}_{\mathrm{iso}}\\ {F}_{\mathrm{org}}^{{\prime} }({\boldsymbol{n}},\theta ) & = & \displaystyle \frac{{F}_{\mathrm{org}}({\boldsymbol{n}},\theta )}{{\displaystyle \int }_{4\pi }{F}_{\mathrm{org}}d{\rm{\Omega }}},{F}_{\mathrm{iso}}=1/4\pi .\end{array}\end{eqnarray} \tag{ 2 }$

**Figure 3.** An example of the SBG model with θ = 10° (in equatorial coordinates). The dotted white lines represent the GP and the SGP (SGP). The top six contributing SBGs are noted as gray stars.
Download figure:
Standard image High-resolution image

2.2. Flux Mapping

At this stage, the flux model does not include the GMF deflections. To include them, we use a back-propagation technique using the CR propagation code CRPropa3 (Batista et al. 2016). We adopt the JF12 model (Jansson & Farrar 2012a, 2012b) as the GMF model. We calculate the trajectories of antiprotons of energy E emitted from the Earth to a sphere of 20 kpc radius from the GC (the galaxy sphere). The position of the Earth is defined to be 8.5 kpc away from the GC, following the JF12 model. These trajectories represent the trajectories of particles with the same rigidity that can arrive on the Earth through the Galactic sphere. The trajectory of a higher-mass particle with the charge Ze is replaced with the trajectory of a proton with the rigidity R = E/Ze.

Based on the CR trajectories obtained through the back-propagation, we convert the CR flux model on the galaxy sphere into that on Earth through the GMF model (flux mapping). We define the original CR flux model as F_org( n _org, θ), where n _org indicates the direction on the galaxy sphere. The conversion from the directions on the Earth n _earth to those on the galaxy sphere n _org of particles with rigidity R is expressed as

$\begin{eqnarray}&&{{\boldsymbol{n}}}_{\mathrm{org}}={A}_{\mathrm{BT}}({{\boldsymbol{n}}}_{\mathrm{earth}},R),\end{eqnarray} \tag{ 3 }$

where A_BT indicates the conversion function. As Liouville's theorem tells that the flux value along each CR trajectory remains constant (Bradt & Olbert 2008), we can determine the CR flux on Earth F_earth as

$\begin{eqnarray}\begin{array}{rcl}{F}_{\mathrm{earth}}({{\boldsymbol{n}}}_{\mathrm{earth}},\theta ,R) & = & {F}_{\mathrm{org}}({{\boldsymbol{n}}}_{\mathrm{org}},\theta )\\ & = & {F}_{\mathrm{org}}({A}_{\mathrm{BT}}({{\boldsymbol{n}}}_{\mathrm{earth}},R),\theta ).\end{array}\end{eqnarray} \tag{ 4 }$

Because of the nearby source contributions, the photonuclear interaction during intergalactic propagation is not taken into account in this study. Some examples of F_earth for different R based on F_org in Figure 3 are shown in Figure 4. As shown in Figure 4 (top left), at the highest rigidity (R = 100 EV), the GMF does not affect the CR flux. As the rigidity R decreases, the peaks of CR flux around NGC1068 and NGC253 become displaced from the true source directions (top right panel in Figure 4) for $\mathrm{log}(R/\mathrm{EV})=1.5$ ). At lower rigidity (bottom left panel in Figure 4 for R = 10 EV), the displacements become larger and the peak around NGC4945 splits along the GP. Through visual inspections, it is clear that the GMF bias is stronger in the southern hemisphere.

**Figure 4.** Examples of the SBG model as seen from Earth when R = 10^2.0, 10^1.5, and 10^1.0 EV and θ = 10° (the JF12 model). The color scale is the same as in Figure 3. The rigidity R is shown in the top right corner in each panel in the log scale.
Download figure:
Standard image High-resolution image

2.3. Generation of Mock Data Sets

To quantitatively discuss the GMF biases in Section 3, we generate mock events as follows. We adopt a best-fit function and the parameters given in Heinze & Fedynitch (2019) based on the observed UHECRs from the Auger experiment (Aab et al. 2017). In Heinze & Fedynitch (2019), the energy spectrum of each mass (A) at the source is assumed by the following function J_A, and fitting was performed for the mass composition observed on Earth,

$\begin{eqnarray}&&{J}_{A}(E)={{ \mathcal J }}_{A}{f}_{\mathrm{cut}}(E,{Z}_{A},{R}_{\max }){n}_{\mathrm{evol}}(z){\left(\displaystyle \frac{E}{{10}^{9}\ \mathrm{GeV}}\right)}^{-\gamma }.\end{eqnarray} \tag{ 5 }$

The cutoff function f_cut is given as

$\begin{eqnarray}{f}_{\mathrm{cut}}=\left\{\begin{array}{ll}1 & (E\lt {Z}_{A}{R}_{\max })\\ \exp \left(1-\tfrac{E}{{Z}_{A}{R}_{\max }}\right) & (E\gt {Z}_{A}{R}_{\max }).\end{array}\right.\end{eqnarray} \tag{ 6 }$

Because we only focus on the nearby sources, the redshift-evolution term n_evol(z) is approximated to be 1. We also assume that the mass composition observed on Earth and the composition at the source are the same. The fractions of elements are defined as f_A = ${{ \mathcal J }}_{A}/{{\rm{\Sigma }}}_{A}{{ \mathcal J }}_{A}$ at 10 EeV in Heinze & Fedynitch (2019). We adapt the best-fit parameters from Heinze & Fedynitch (2019) as γ = −0.80 and ${R}_{\max }=1.6\,\mathrm{EV}$ . We also adapt the values of f_A as (¹H,⁴He,¹⁴N, ²⁸Si,⁵⁶Fe) = (0.0, 82.0, 17.3, 0.6, 2.0 · 10⁻²)[%]. According to this mass fraction and to the spectra, we determine the mass and energy, i.e., R, of each mock event. The energy E of a mock event is randomly sampled from the spectrum, and the rigidity R of the event is calculated through the formula R = E/Ze. We adapt the minimum energy ${E}_{\min }=40\,\mathrm{EeV}$ according to previous studies (Aab et al. 2018, ${E}_{\min }=39\,\mathrm{EeV}$ ). Using the selected rigidity R, we determine the arrival direction of the event based on the CR flux model defined in Section 2.2.

An example of the distribution of the mock event arrival directions is provided in Figure 5. The distribution in Figure 5 is similar to the distribution of pure-carbon case (Figure 10 in Appendix A.1). Although we can see clustering around M82 and NGC4945, the centers of the distributions are displaced from the source directions. The events that originated from NGC1068 and NGC253 are mostly deflected. This suggests that the real UHECR distribution cannot be reproduced with a single isotropic smearing, and the deflections by the GMF depend on the arrival directions.

**Figure 5.** Example of the distribution of the mock events ( $({f}_{\mathrm{ani}}^{\mathrm{true}},{\theta }^{\mathrm{true}})=(100 \% ,10^\circ )$ ). The gray dots show the arrival directions of 4000 mock events. The directions of the SBGs whose contribution exceeds 5% are shown by black stars and their name. The area of each star indicates its relative contribution in Tabel 1. Black dots (circles) present the GP (SGP).
Download figure:
Standard image High-resolution image

2.4. The Sky Coverage of the Experiments

To make the comparison with the analysis of the observed UHECRs (Aab et al. 2018; Abbasi et al. 2018), the sky coverage of the TA and Auger experiments is considered based on the equations given by Sommers (2001). In Sommers (2001), the sky coverage ω( n _CR) depends on the decl. δ,

$\begin{eqnarray}&&\omega (\delta )\propto \cos ({a}_{0})\cos (\delta )\sin ({\alpha }_{{\rm{m}}})+{\alpha }_{{\rm{m}}}\sin ({a}_{0})\sin (\delta ).\end{eqnarray} \tag{ 7 }$

$\begin{eqnarray}{\alpha }_{{\rm{m}}}=\left\{\begin{array}{ll}0\ & (\xi \gt 1)\\ \pi \ & (\xi \lt -1)\\ {\cos }^{-1}(\xi )\ & (-1\lt \xi \lt 1)\end{array}\right.\end{eqnarray} \tag{ 8 }$

$\begin{eqnarray}&&\xi =\displaystyle \frac{\cos ({\theta }_{{\rm{m}}})-\sin ({a}_{0})\sin (\delta )}{\cos ({a}_{0})\cos (\delta )},\end{eqnarray} \tag{ 9 }$

where θ_m is the maximum zenith angle and a₀ is the latitude of the experimental site. We adopt the latitude a₀ = 39 fdg 3 (−35 fdg 2) and the maximum zenith angle θ_m = 55° (60°) for the TA (Auger) experiment. From the all-sky data sets, we randomly select mock events with the probability of each sky coverage. Out of the 4000 mock events in each data set, approximately 1000 mock events each are selected by the TA and Auger coverage. We define the data set selected by the sky coverage of TA (Auger) as the north-sky (south-sky) data set.

3. Analysis

In order to investigate how much GMF deflections affect the estimated parameters, we conduct the same maximum likelihood analysis as in Aab et al. (2018) on the mock data sets. We test two hypotheses for the CR flux models. One is a flux with nonzero f_ani with the SBG model (F_norm), and the other is the isotropic flux (F_iso), i.e., f_ani = 0. The test statistics TS are calculated as the log-likelihood ratio,

$\begin{eqnarray}&&{TS}=2\mathrm{ln}(L({F}_{\mathrm{norm}})/L({F}_{\mathrm{iso}})).\end{eqnarray} \tag{ 10 }$

The likelihood of each model L(F) is given as

$\begin{eqnarray}&&L(F)=\displaystyle \prod _{\mathrm{CR}}\displaystyle \frac{F({{\boldsymbol{n}}}_{\mathrm{CR}})\omega ({{\boldsymbol{n}}}_{\mathrm{CR}})}{{\int }_{4\pi }F({\boldsymbol{n}})\omega ({\boldsymbol{n}})d{\rm{\Omega }}},\end{eqnarray} \tag{ 11 }$

where F, ω( n _CR), and n _CR are the normalized CR flux model, the sky coverage of each experiment (Equations (7)), and the arrival directions of the observed UHECRs, respectively.

By scanning the set of parameters (f_ani, θ), the best-fit parameters that maximize the TS in Equation (10) are determined.

4. Results

We show the distribution of the best-fit parameters for 1000 mock event data sets in Figures 6. The different panels show the results for the different true parameters $({f}_{\mathrm{ani}}^{\mathrm{true}},{\theta }^{\mathrm{true}})$ at the top of each panel and with the gray star in each panel. From the left to right column, the true parameter ${f}_{\mathrm{ani}}^{\mathrm{true}}$ is given as 20%, 40%, and 60%. From the top to the bottom row, the true parameter θ is given as 10°, 20°, and 30°. The black, blue, and red contours indicate the 68% and 95% percentile containment for the distributions for the all-sky, north-sky, and south-sky data sets, respectively. The black cross, blue circle, and red triangle show the most frequent values $({\tilde{f}}_{\mathrm{ani}},\tilde{\theta })$ of the best-fit parameters of the all-sky, north-sky, and south-sky data sets, respectively.

**Figure 6.** Distributions of the best-fit parameters for the 1000 mock event data sets. From the left to right column, the true parameter ${f}_{\mathrm{ani}}^{\mathrm{true}}$ is given as 20%, 40%, and 60%. From the top to the bottom row, the true parameter θ is 10°, 20°, and 30°. The true parameters are marked by the gray stars. The black, blue, and red contours indicate the 68% and 95% percentile containment for the all-sky, north-, and south-sky data sets, respectively. The black cross, blue circle, and red triangle show the most frequent values $({\tilde{f}}_{\mathrm{ani}},\tilde{\theta })$ for the all-syk, north-, and south-sky data sets, respectively. The distributions of the best-fit parameters are smoothed with a kernel-Gaussian distribution. The best-fit parameter $({f}_{\mathrm{ani}}^{\mathrm{Auger}},{\theta }^{\mathrm{Auger}})=(9.7 \% ,12.9\,\deg )$ in Aab et al. (2018) is shown as the black triangle.
Download figure:
Standard image High-resolution image

**Figure 6.** Distributions of the best-fit parameters for the 1000 mock event data sets. From the left to right column, the true parameter ${f}_{\mathrm{ani}}^{\mathrm{true}}$ is given as 20%, 40%, and 60%. From the top to the bottom row, the true parameter θ is 10°, 20°, and 30°. The true parameters are marked by the gray stars. The black, blue, and red contours indicate the 68% and 95% percentile containment for the all-sky, north-, and south-sky data sets, respectively. The black cross, blue circle, and red triangle show the most frequent values $({\tilde{f}}_{\mathrm{ani}},\tilde{\theta })$ for the all-syk, north-, and south-sky data sets, respectively. The distributions of the best-fit parameters are smoothed with a kernel-Gaussian distribution. The best-fit parameter $({f}_{\mathrm{ani}}^{\mathrm{Auger}},{\theta }^{\mathrm{Auger}})=(9.7 \% ,12.9\,\deg )$ in Aab et al. (2018) is shown as the black triangle.
Download figure:
Standard image High-resolution image

The distributions of the all-sky, north-sky, and south-sky data sets do not agree with each other, especially for a higher anisotropic fraction ${f}_{\mathrm{ani}}^{\mathrm{true}}$ (the right panels). When the anisotropic fraction ${f}_{\mathrm{ani}}^{\mathrm{true}}$ is larger and separation angular scale θ^true is smaller (the four upper right panels), the estimated separation angular scale θ becomes larger due to the deflection of the GMF. When the anisotropic fraction ${f}_{\mathrm{ani}}^{\mathrm{true}}$ is smaller (left panels), the distribution of the all-sky, north-, and south-sky data sets are similar due to the low contrast between the SBG model and the isotropic backgrounds.

Focusing on the north–south difference, the GMF affects the results of south-sky data sets more than north-sky data sets. This can be explained with two reasons: first, the GMF deflections are larger around the GC. In the rest of the sky, the GMF deflection is larger in the Galactic South (see also Figure 11 in Jansson & Farrar 2012a). NGC253 and NGC1068, which contribute to the south-sky data sets, are located in the region of the sky in which the GMF deflections are large. In any parameter set $({f}_{\mathrm{ani}}^{\mathrm{true}},{\theta }^{\mathrm{true}})$ , it is found that the most frequent values of the anisotropic fraction ${\tilde{f}}_{\mathrm{ani}}$ for the south-sky data sets (red triangles) are largely underestimated, namely below 50% of the true ones ${f}_{\mathrm{ani}}^{\mathrm{true}}$ (gray stars). Regardless of the true parameters $({f}_{\mathrm{ani}}^{\mathrm{true}},{\theta }^{\mathrm{true}})$ , the distributions for the south-sky data sets include the best-fit parameters $({f}_{\mathrm{ani}}^{\mathrm{Auger}},{\theta }^{\mathrm{Auger}})=(9.7 \% ,12.9\,\deg )$ , which are indicated by the black triangles in the 2σ contour. We discuss this tendency in Section 5.2.

5. Discussions

5.1. The Uncertainty of the GMF Models

In this section, we discuss the uncertainty effect on the GMF models. To test the effects caused by the uncertainty of the halo components of the JF12 model, we conduct the same analysis, but change the halo components in the model within 1σ uncertainty, generate the mock event data sets, and repeat the analysis. It is found that the uncertainty of the halo components does not have a large effect on the most frequent parameters $({\tilde{f}}_{\mathrm{ani}},\tilde{\theta })$ .

For an independent comparison with the JF12 model, we also refer to the Pshirkov & Tinyakov (2011) model (PT11; Pshirkov et al. 2011). We generate the mock event data sets based on the PT11 model. Except for the GMF model, the other assumptions (the SBG model and the mixed-mass composition) are the same. Although the separation on the angular scale θ in the south-sky data sets becomes smaller than with the JF12 model, the most frequent values of the anisotropic fractions ${\tilde{f}}_{\mathrm{ani}}$ are also reduced by more than 50% compared to the true value ${f}_{\mathrm{ani}}^{\mathrm{true}}$ (see also Figure 16 in Appendix B).

5.2. Comparison with the Best-fit Parameters in Aab et al. (2018)

In this section, we search for a set of the true parameters $({f}_{\mathrm{ani}}^{\mathrm{true}},{\theta }^{\mathrm{true}})$ that is compatible with the best-fit parameters $({f}_{\mathrm{ani}}^{\mathrm{Auger}},{\theta }^{\mathrm{Auger}})$ . We generate 4000 south-sky data sets in the same manner as in Section 2.3, but each has the same number of events (894 events) as was used in the analysis of Aab et al. (2018). These mock event data sets are analyzed in the same manner, and the best-fit parameters are obtained. Examples of the distributions of 4000 best-fit parameters are shown in Figure 7. Because of the GMF bias and the statistical fluctuation regardless of the true parameters $({f}_{\mathrm{ani}}^{\mathrm{true}},{\theta }^{\mathrm{true}})$ , marked by the gray star, the estimated parameters tend to be distributed around the Auger best-fit parameters $({f}_{\mathrm{ani}}^{\mathrm{Auger}},{\theta }^{\mathrm{Auger}})$ , marked by the black triangle within the 68% or 95% containment levels. We classify the true parameters $({f}_{\mathrm{ani}}^{\mathrm{true}},{\theta }^{\mathrm{true}})$ according to whether $({f}_{\mathrm{ani}}^{\mathrm{Auger}},{\theta }^{\mathrm{Auger}})$ is contained in the 68% or 95% contours. Figure 8 shows the result of this classification. From Figure 8, except in the bottom right corner, a wide range of parameters is still compatible with $({f}_{\mathrm{ani}}^{\mathrm{Auger}},{\theta }^{\mathrm{Auger}})$ . Considering the GMF effect and the mass-dependent energy spectrum, a large contribution of SBGs to the UHECR flux is still possible.

**Figure 7.** Examples of the distributions of the best-fit parameters estimated from the mock event data sets. Red contours show the 68 and 95% percentile containments of the best-fit parameters. The red triangle shows the most frequent value of the best-fit parameters $({\tilde{f}}_{\mathrm{ani}},\tilde{\theta })$ of the mock event data sets. The best-fit parameter $({f}_{\mathrm{ani}}^{\mathrm{Auger}},{\theta }^{\mathrm{Auger}})=(9.7 \% ,12.9\,\deg )$ in Aab et al. (2018) is shown as the black triangle. Note that the mock event data sets are generated with a set of assumptions for the source model (SBG), the GMF model (the JF12), and the mass-composition model (Heinze & Fedynitch 2019).
Download figure:
Standard image High-resolution image

**Figure 7.** Examples of the distributions of the best-fit parameters estimated from the mock event data sets. Red contours show the 68 and 95% percentile containments of the best-fit parameters. The red triangle shows the most frequent value of the best-fit parameters $({\tilde{f}}_{\mathrm{ani}},\tilde{\theta })$ of the mock event data sets. The best-fit parameter $({f}_{\mathrm{ani}}^{\mathrm{Auger}},{\theta }^{\mathrm{Auger}})=(9.7 \% ,12.9\,\deg )$ in Aab et al. (2018) is shown as the black triangle. Note that the mock event data sets are generated with a set of assumptions for the source model (SBG), the GMF model (the JF12), and the mass-composition model (Heinze & Fedynitch 2019).
Download figure:
Standard image High-resolution image

**Figure 8.** Excluded region for the best-fit parameter in Aab et al. (2018). Circles show the searched-for true parameters $({f}_{\mathrm{ani}}^{\mathrm{true}},{\theta }^{\mathrm{true}})$ . White (gray) indicates the true parameters $({f}_{\mathrm{ani}}^{\mathrm{true}},{\theta }^{\mathrm{true}})$ that reproduce the best-fit parameter $({f}_{\mathrm{ani}}^{\mathrm{Auger}},{\theta }^{\mathrm{Auger}})=(9.7 \% ,12\buildrel{\circ}\over{.} 9)$ in Aab et al. (2018) within the 68% (95%) percentile. The best-fit parameter $({f}_{\mathrm{ani}}^{\mathrm{Auger}},{\theta }^{\mathrm{Auger}})=(9.7 \% ,12\buildrel{\circ}\over{.} 9)$ in Aab et al. (2018) is shown as the black triangle. The parameters are scanned with resolutions of ${\rm{\Delta }}{f}_{\mathrm{ani}}^{\mathrm{true}}=10 \%$ and Δθ^true = 5°.
Download figure:
Standard image High-resolution image

**Figure 8.** Excluded region for the best-fit parameter in Aab et al. (2018). Circles show the searched-for true parameters $({f}_{\mathrm{ani}}^{\mathrm{true}},{\theta }^{\mathrm{true}})$ . White (gray) indicates the true parameters $({f}_{\mathrm{ani}}^{\mathrm{true}},{\theta }^{\mathrm{true}})$ that reproduce the best-fit parameter $({f}_{\mathrm{ani}}^{\mathrm{Auger}},{\theta }^{\mathrm{Auger}})=(9.7 \% ,12\buildrel{\circ}\over{.} 9)$ in Aab et al. (2018) within the 68% (95%) percentile. The best-fit parameter $({f}_{\mathrm{ani}}^{\mathrm{Auger}},{\theta }^{\mathrm{Auger}})=(9.7 \% ,12\buildrel{\circ}\over{.} 9)$ in Aab et al. (2018) is shown as the black triangle. The parameters are scanned with resolutions of ${\rm{\Delta }}{f}_{\mathrm{ani}}^{\mathrm{true}}=10 \%$ and Δθ^true = 5°.
Download figure:
Standard image High-resolution image

5.3. Maximum Likelihood Analysis Method with the CR Flux Models on the Earth

To calculate the likelihood (Equation (11)) and TS (Equation (10)), we use the original CR flux model F_org instead of the CR flux model on the Earth F_earth (Equation (2)). This causes the GMF bias in the parameter estimations. To reduce the GMF bias in the previous parameter estimation, it is necessary to replace F_org with F_earth in Equation (2). Note that this analysis is valid only when the GMF and mass-dependent spectrum models are correct. In other words, we need to test a set of assumptions together. We rewrite Equation (2) as follows:

$\begin{eqnarray}\begin{array}{rcl}{F}_{\mathrm{norm}}({\boldsymbol{n}},{f}_{\mathrm{ani}},\theta ,R) & = & {f}_{\mathrm{ani}}{F}_{\mathrm{earth}}^{{\prime} }({\boldsymbol{n}},\theta ,R)\\ & & +\,(1-{f}_{\mathrm{ani}}){F}_{\mathrm{iso}}\\ {F}_{\mathrm{earth}}^{{\prime} } & = & \displaystyle \frac{{F}_{\mathrm{earth}}({\boldsymbol{n}},\theta ,R)}{{\displaystyle \int }_{4\pi }{F}_{\mathrm{earth}}d{\rm{\Omega }}},{F}_{\mathrm{iso}}=1/4\pi .\end{array}\end{eqnarray} \tag{ 12 }$

Here, F_earth( n _earth, θ, R) is obtained using Equation (4). Thus, we can rewrite the CR flux models from the sources ${F}_{\mathrm{earth}}^{{\prime} }({\boldsymbol{n}},{f}_{\mathrm{ani}},\theta ,R)$ as

$\begin{eqnarray}&&{F}_{\mathrm{earth}}^{{\prime} }({{\boldsymbol{n}}}_{\mathrm{CR}},{f}_{\mathrm{ani}},\theta ,{R}_{\mathrm{CR}})=\displaystyle \frac{{F}_{\mathrm{org}}({A}_{\mathrm{BT}}({{\boldsymbol{n}}}_{\mathrm{CR}},{R}_{\mathrm{CR}}),\theta )}{{\int }_{4\pi }{F}_{\mathrm{org}}({A}_{\mathrm{BT}}({\boldsymbol{n}},R),\theta )d{\rm{\Omega }}}.\end{eqnarray} \tag{ 13 }$

The denominator ∫_4π F_org(A_BT( n , R), θ)dΩ in Equation (13) is derived by integrating F_earth.

We conduct the improved maximum likelihood analysis following Equation (12) to the same data sets as in Section 4. Note that the new analysis is carried out assuming that we know the event-by-event mass of each mock event. Figure 9 illustrates the results in the same manner as Figure 6, but for the estimation following Equation (12). In all cases, the analysis improves the estimates of the parameters (f_ani, θ). Specifically, the GMF bias is reduced when the separation angular scales θ are small. For a larger separation angular scale, there is still a significant difference between the estimated parameters and the true parameters $({f}_{\mathrm{ani}}^{\mathrm{true}},{\theta }^{\mathrm{true}})$ . Because the GMF bias caused by the regular component of the GMF is effectively reduced within the 1σ contour, the origins of the dispersion should be the statistical fluctuation. Future observations with large statistics are expected to reduce the dispersion. In this analysis, we assumed a perfect event-by-event energy and (rigidity) resolution, but the effect of realistic resolutions will be discussed in a future publication.

**Figure 9.** Same as Figure 6, but with the improved analysis method in Section 5.3. The analysis is applied for the same mock event data sets with a mixed-mass assumption and the JF12 model in Figure 6.
Download figure:
Standard image High-resolution image

5.4. Discussions of the Random GMF and EGMF

In this section, we refer to the components of the magnetic fields that are not fully considered in this study. To reveal the effects of coherent deflection by the GMF, we only focus on the regular component of the GMF. We assume that the separation angular scale θ includes both the EGMF and a random component of the GMF (random GMF). Aab et al. (2018) and Abbasi et al. (2018) also included regular components of the GMF. In general, the deflections by the EGMF and the random GMF also should have a rigidity dependence. Bray & Scaife (2018) suggested that the upper limit to the EGMF lies at the nano-Gauss scale. Although the upper limit to the EGMF is smaller than that of the GMF, the distance between each source and the Earth is much larger than the radius of our galaxy. We need to consider this distance dependence (see also Anchordoqui 2019). A random component of the GMF has an arrival-direction dependence that is the same as for the regular component of the GMF. Pshirkov et al. (2013) investigated the random deflections in the GMF. They suggested that the random component of the GMF deflects 40 EeV protons by less than 1°–2° in most of the sky and ∼5° along the GP.

Although the physically correct description of random components of the GMF and EGMF is important, it is beyond the scope of this study. The effect of random components in the GMF and EGMF should also have an arrival-direction dependence and rigidity dependence. We will take them into account in a future realistic model.

6. Summary

We estimate the biases caused by the coherent deflections due to the Galactic magnetic field in searches for UHECR sources in the literature. We generated mock event data sets with a set of assumptions for the source model (Aab et al. 2018), coherent deflection by the GMF model (Jansson & Farrar 2012a, 2012b), and a mass-composition model (Heinze & Fedynitch 2019), and conducted a maximum likelihood analysis on the data sets neglecting the GMF in the same manner as in previous studies. Our main results are listed below.

1.
The distributions of the estimated parameters (f_ani, θ) are displaced from the true parameters $({f}_{\mathrm{ani}}^{\mathrm{true}},{\theta }^{\mathrm{true}})$ . This confirms the existence of the GMF bias.
2.
The distributions of the estimated parameters (f_ani, θ) in the all-sky, north-sky, and south-sky data sets do not agree with each other. The directional or sky dependence of the GMF bias is also confirmed.
3.
We find that the estimated f_ani is systematically reduced by more than 50% in the south-sky data sets.
4.
We search for the true parameters $({f}_{\mathrm{ani}}^{\mathrm{true}},{\theta }^{\mathrm{true}})$ that are compatible with the best-fit parameters reported in Aab et al. (2018), taking the number of events used in their study into account. Except for the narrow region with the large anisotropic fraction ${f}_{\mathrm{ani}}^{\mathrm{true}}$ and small separation angular scale θ^true, a wide parameter space is still compatible with the experimental result within the 95% confidence limit.
5.
We develop a maximum likelihood analysis that takes the GMF deflections into account and confirm that the parameters would be correctly estimated within the 1σ contour under the ideal condition that we know the event-by-event energy and mass of each UHECR event and the GMF structure.

Note that again this study is conducted under a specific set of assumptions: the source model (Aab et al. 2018), the magnetic field model (Jansson & Farrar 2012a, 2012b), the energy spectrum, and the mass composition (Tsunesada et al. 2017; Heinze & Fedynitch 2019).

Although these models and the assumptions are to be tested and updated regularly, the technique in Section 5.3 can be applied to future models and updated observational data sets. The extension of the TA and Auger experiments (Castellina & Pierre Auger Collaboration 2019; Abbasi et al. 2021, TA × 4 and AugerPrime) and next-generation UHECR observation (Hörandel 2021) will play an important role. The improvement of the GMF and CR propagation models also leads us to more realistic source searches (Boulanger et al. 2018).

We thank the members of the Telescope Array collaboration for fruitful discussions. We are grateful to Peter Tinyakov, Anatoli Fedynitch, and Federico Urban for fruitful discussions and revisions. This work was supported by JSPS KAKENHI grant Nos. JP19J11429 and JP19KK0074 and the joint research program of the Institute for Cosmic Ray Research (ICRR) at the University of Tokyo. E.K. is grateful for support from the "Pioneering Program of RIKEN for Evolution of Matter in the Universe (r-EMU)."

Appendix A: Analysis with a Single-mass Assumption

A.1. Mock Data Sets with a Single-mass Assumption

As single-mass assumptions, we test the pure-proton, He, C, Si, and Fe cases. We fix the energy spectrum as a broken power law with spectral indexes γ = −2.69 (E < 10^1.81 EeV) and γ = −4.63 (E > 10^1.81 EeV), which is as reported by the TA experiment (Tsunesada et al. 2017). In the same manner as in Section 2.3, we choose the arrival direction of the anisotropic event based on the generated CR flux models F_earth. Examples of the distribution of the mock event arrival directions with the single-mass assumption and with $({f}_{\mathrm{ani}}^{\mathrm{true}},{\theta }^{\mathrm{true}})=(100 \% ,10^\circ )$ are shown in Figure 10. For the light masses of proton and helium, for instance, the distributions are similar to that of Figure 3, which means that the GMF bias is small. On the other hand, for the higher masses, the distortion due to the GMF is significant. In the pure-Fe assumption, the clustering of events around the top four contributing SBGs (M82, NGC4945, NGC1068, and NGC253) are not seen.

**Figure 10.** Examples of the distribution of the mock event arrival directions with the single-mass assumptions ( $({f}_{\mathrm{ani}}^{\mathrm{true}},{\theta }^{\mathrm{true}})=(100 \% ,10^\circ )$ and 4000 events over the all-sky region).
Download figure:
Standard image High-resolution image

A.2. Estimated Parameters in Previous Studies with Single-mass Assumption

The results of the likelihood analysis for single-mass models are shown in Figures 11–15. In the pure-proton (Figure 11) and pure-He (Figure 12) cases, the most frequent best-fit values $({\tilde{f}}_{\mathrm{ani}},\tilde{\theta })$ fall near the true values $({f}_{\mathrm{ani}}^{\mathrm{true}},{\theta }^{\mathrm{true}})$ , which means that the GMF bias is small. However, when the separation angular scale θ^true is larger, the dispersion of the estimated parameters, especially for f_ani, becomes larger. This dispersion gives an intrinsic statistical uncertainty in the f_ani estimation. In the pure-C case (Figure 13), a discrepancy in the distributions between the north-sky and south-sky data sets can be seen. This tendency becomes larger for the heavier single-mass cases. In the pure-Si and pure-Fe cases, the most frequent value of ${\tilde{f}}_{\mathrm{ani}}$ in the south-sky data sets becomes 0% with any values of the true parameters (Figures 14 and 15). Although both are separated from the true parameters, the distributions of north-sky data sets are closer to those of the all-sky data sets in any case. The single and dominant source contribution of M82 and smaller deflection by the GMF in the northern sky can explain this tendency.

**Figure 13.** Same as Figure 11, but for the pure-C case.
Download figure:
Standard image High-resolution image

**Figure 14.** Same as Figure 11, but for the pure-Si case.
Download figure:
Standard image High-resolution image

**Figure 15.** Same as Figure 11, but for the pure-Fe case.
Download figure:
Standard image High-resolution image

Appendix B: Analysis with the PT11 Model

For an independent comparison with the JF12 model, we also refer to the Pshirkov & Tinyakov 2011 model (PT11) in this study (Pshirkov et al. 2011). The mock event data sets are generated in the same manner as in Section 2.3, except for the GMF model. Figure 16 shows the distributions of the best-fit parameters in the same manner as Figure 6. Although there is a quantitative difference, both results show the same tendency due to the GMF bias.

**Figure 16.** Same as Figure 6, but for the mock event data sets generated with the PT11 model and the mixed-mass assumption.
Download figure:
Standard image High-resolution image

Effects of the Galactic Magnetic Field on the UHECR Correlation Studies with Starburst Galaxies

Article metrics

Author e-mails

Author affiliations

ORCID iDs

Dates

Abstract

1. Introduction

2. Mock Data-set Production