Cosmic-Ray Short Burst Observed with the Global Muon Detector Network (GMDN) on 2015 June 22

K. Munakata; M. Kozai; P. Evenson; T. Kuwabara; C. Kato; M. Tokumaru; M. Rockenbach; A. Dal Lago; R. R. S. de Mendonca; C. R. Braga; N. J. Schuch; H. K. Al Jassar; M. M. Sharma; M. L. Duldig; J. E. Humble; I. Sabbah; J. Kóta

doi:10.3847/1538-4357/aacdfe

1. Introduction

The GRAPES-3 muon telescope observed a "cosmic-ray burst" (CRB) in which the muon count rate increased ∼1% for two hours during 2015 June 22 (Mohanty et al. 2016; hereafter referred to as Paper I, see also Nonaka et al. 2006 for GRAPES-3). This burst was recorded shortly after the arrival of a strong interplanetary shock identified by abrupt increases of the solar wind velocity and the interplanetary magnetic field (IMF) strength. There were two preceding shocks recorded about 15 and 27 hr earlier. These three shocks were formed in front of coronal mass ejections successively ejected from the same active region NOAA 2371. The third shock was followed by a strong enhancement of south directing IMF, which triggered a major geomagnetic storm reaching a K_p index of 8+. Paper I proposed that the CRB was caused by the substantial reduction of the geomagnetic cutoff rigidity and the alteration of asymptotic cosmic-ray orbits in the magnetosphere.

In this paper, we analyze the CRB observed with the Global Muon Detector Network (GMDN) comprising four multidirectional muon detectors located in Japan, Australia, Brazil, and Kuwait, designed for accurate observation of the global cosmic-ray anisotropy. In our previous analyses of the GMDN data, the cosmic-ray anisotropy and density have been derived on an hourly basis and used to analyze the geometry of the cosmic-ray depletion region in individual Forbush decreases (Kuwabara et al. 2009; Kozai et al. 2016). Here, we use 10 minute data to resolve and analyze the rapidly changing anisotropy and density in the CRB to show that the CRB was caused by a rapid enhancement of the cosmic-ray anisotropy and density outside the magnetosphere, and not by the geomagnetic storm.

We describe the GMDN and data analysis in Sections 2.1 and 2.2, respectively, and show the results in Section 2.3. We give a summary and discussion in Section 3.

2. Data Analysis and Results

2.1. Global Muon Detector Network

The GMDN comprises four multidirectional muon detectors, "Nagoya" in Japan, "Hobart" in Australia, "Kuwait City" in Kuwait, and "São Martinho da Serra" in Brazil, recording muon count rates in 60 directional channels viewing almost the entire sky around Earth. Each detector except Kuwait City consists of two horizontal layers separated by 1.73 m of 1 m² plastic scintillators (PSs), each viewed by a 12.7 cm diameter photomultiplier tube. By counting twofold coincidences between pairs of detectors in the upper and lower layers, we record the rate of muons from the corresponding incident direction. Kuwait City consists of four horizontal layers of proportional counter tubes (PCTs), each 5 m long with a 10 cm diameter and with a 50 μm thick tungsten anode along the cylinder axis. The PCT axes are aligned east–west (X) in the top and third layers and north–south (Y) in the second and bottom layers. The top and second layers form an upper pair, while the third and bottom layers form a lower pair. The two pairs are separated vertically by 80 cm. Muon recording is triggered by the fourfold coincidence of pulses from all layers and the incident direction is identified from X–Y locations of the upper and lower PCT pairs.

Table 1 summarizes characteristics of directional channels of the GMDN, while Figure 1 shows the asymptotic viewing directions (corrected for geomagnetic bending of cosmic-ray orbits) of 60 directional channels of the GMDN. The detector configurations in the GMDN are also available at our website.¹² In our calculations of the median primary rigidity (P_m) and the asymptotic viewing direction (ϕ_asymp, λ_asymp) at P_m of each directional channel, we use the response function of the atmospheric muons to the primary cosmic rays given by numerical solutions of the hadronic cascade in the atmosphere (Murakami et al. 1979). We also performed Monte Carlo (MC) simulations of the hadronic cascade by using CORSIKA (HDPM+GHEISHA), calculated the response function and found that both response functions are in good agreement except for minor differences in the high-energy region above 1 TeV. The difference between the calculated median primary rigidities is less than ∼5% even for the most inclined directional channels. We also find that the observed count rates are ∼5% to 15% higher than the calculations. This ∼15% underestimation of our calculations, however, is expected to be constant and does not affect our analyses in this paper based on the fractional intensity (see Section 2.2).

Table 1. Characteristics and Estimated Responses to the Geomagnetic Storm of Four Vertical Channels of GMDN

Detector Name	Detector Type	Count Rate	σ_ci,j	P_m	ϕ_asymp	λ_asymp	P_c	ΔP_c	ΔI_i,j(t)
(No. of Directions)	(Detection Area)	(10⁴/10 minute)	(%)	(GV)	(°)	(°)	(GV)	(GV)	(%)
Nagoya-V	PS	47.8	0.14	58.4	168.9	27.7	12.3	−0.07	+0.08
(17)	(36 m²)	(3.0 ∼ 47.8)	(0.14 ∼ 0.58)	(58.4 ∼ 106.9)	(89.1 ∼ 235.8)	(−24.4 ∼ 64.0)
Hobart-V	PS	23.7	0.21	53.1	171.0	−39.8	1.6	−0.64	+0.00
(13)	(16 m²)	(3.2 ∼ 23.7)	(0.21 ∼ 0.56)	(53.1 ∼ 74.0)	(108.0 ∼ 237.1)	(−76.7 ∼ 5.0)
Kuwait City-V	PCT	24.2	0.20	59.8	78.4	23.8	13.0	−0.10	+0.08
(13)	(18.5 m²)	(4.0 ∼ 24.2)	(0.20 ∼ 0.50)	(59.8 ∼ 94.5)	(24.6 ∼ 127.9)	(−19.2 ∼ 72.0)
São Martinho da Serra-V	PS	42.3	0.15	54.3	331.4	−22.4	9.8	−0.82	+0.29
(17)	(32.0 m²)	(0.7 ∼ 42.3)	(0.15 ∼ 1.22)	(54.3 ∼ 98.4)	(259.4 ∼ 32.7)	(−67.1 ∼ 33.4)

Note. Following the detector name (number of directional channels) and type (detection area), the average 10 minute count rate, count rate error (σ_ci,j), median primary rigidity (P_m), geographic longitude (ϕ_asymp), and latitude (λ_asymp) of asymptotic viewing direction outside the magnetosphere are listed for each of the four vertical channels of the GMDN. Each number in brackets in columns 3–7 indicates a range of corresponding parameter covered by all directional channels available in each detector. The last three columns are the average geomagnetic cutoff rigidity (P_c), its maximum reduction (ΔP_c), and the maximum increase of count rate expected from ΔP_c, which are calculated using the model geomagnetic field during the geomagnetic storm on 2015 June 22 and the response function of each directional channel to primary GCRs (see the text).

Download table as: ASCII Typeset image

A 5 cm thick layer of lead is installed in each detector to absorb the soft component radiation in the air. The muon threshold energy is 300 MeV for the vertical directional channels and 1300 MeV for the most inclined directional channels. In our detectors, coincidence between any pair of layers triggers a muon count. If there are multi-hits on a layer due to delta-rays produced by muons traversing the 5 cm lead layer, the total count in directional channels exceeds the total count of muons. To evaluate this effect, we performed MC simulations for the Nagoya muon detector in which each layer consists of a 6 × 6 horizontal array of 1 m² plastic scintillators. We used GEANT4 (version 10.3 patch-03, 2017 October 20) for our simulations. We randomly generated muons with various energies and incident directions, produced hit-patterns in each layer of 1 m² plastic scintillators and calculated counts in each directional channel. It was found that the ratio of the multi-hit event number to the trigger number increases with increasing muon energy and zenith angle, while the muon flux at Nagoya decreases. By calculating the ratio expected for 17 directional channels available in Nagoya as a function of monitored zenith angle, we found that the contribution from the multi-hit events to the total count rate is about 3% in the vertical channel and about 15% in the most inclined directional channels. Delta-rays affect the absolute count rate, but the contribution is expected to be constant and does not affect our analyses in this paper based on the fractional intensity. The contribution of delta-rays also slightly increases the median primary energy, but the difference from the value in Table 1 is less than a few GV at most and very small when the broad energy response is considered.

2.2. Data Analysis

We analyze each percent deviation of the pressure corrected 10 minute muon count rate, I_i,j(t) in the jth directional channel of the ith detector in the GMDN at universal time t, from the monthly mean in 2015 June. One minute count rates are also available from Nagoya, Hobart, and São Martinho da Serra, but not from Kuwait City, which was being enlarged in 2015. Since the enlargement of Kuwait City, completed in March of 2016, one minute data have become available from all detectors.

We model I_i,j(t) in terms of the cosmic-ray density (or omnidirectional intensity) I₀(t) and three components $({\xi }_{x}^{\mathrm{GEO}}(t),{\xi }_{y}^{\mathrm{GEO}}(t),{\xi }_{z}^{\mathrm{GEO}}(t))$ of the first-order anisotropy vector ( ${{\boldsymbol{\xi }}}^{\mathrm{GEO}}(t)$ ) in the geographic (GEO) coordinate system, as

$\begin{eqnarray}{I}_{i,j}^{\mathrm{fit}}(t)={I}_{0}(t){c}_{0i,j}^{0} & & +{\xi }_{x}^{\mathrm{GEO}}(t)({c}_{1i,j}^{1}\cos \omega {t}_{i}-{s}_{1i,j}^{1}\sin \omega {t}_{i})\\ & & +{\xi }_{y}^{\mathrm{GEO}}(t)({s}_{1i,j}^{1}\cos \omega {t}_{i}+{c}_{1i,j}^{1}\sin \omega {t}_{i})\\ & & +{\xi }_{z}^{\mathrm{GEO}}(t){c}_{1i,j}^{0},\end{eqnarray} \tag{ 1 }$

where t_i is the local time in hours at the ith detector, ${c}_{0i,j}^{0}$ , ${c}_{1i,j}^{1}$ , ${s}_{1i,j}^{1}$ , and ${c}_{1i,j}^{0}$ are coupling coefficients and ω = π/12. In the GEO coordinate system, we set the x-axis to the anti-Sun direction in the equatorial plane, the z-axis to the geographical north perpendicular to the equatorial plane and the y-axis completing the right-handed coordinate system. The coupling coefficients are calculated using the response function of the atmospheric muon intensity to primary cosmic rays mentioned above (Fujimoto et al. 1984). We derive the best-fit set of four parameters $({I}_{0}(t),{\xi }_{x}^{\mathrm{GEO}}(t),{\xi }_{y}^{\mathrm{GEO}}(t),{\xi }_{z}^{\mathrm{GEO}}(t))$ by solving the following linear equations.

$\begin{eqnarray}&&\displaystyle \frac{\partial S}{\partial {I}_{0}(t)}=\displaystyle \frac{\partial S}{\partial {\xi }_{x}^{\mathrm{GEO}}(t)}=\displaystyle \frac{\partial S}{\partial {\xi }_{y}^{\mathrm{GEO}}(t)}=\displaystyle \frac{\partial S}{\partial {\xi }_{z}^{\mathrm{GEO}}(t)}=0,\end{eqnarray} \tag{ 2 }$

where S is the residual defined, as

$\begin{eqnarray}&&S=\displaystyle \sum _{i,j}{({I}_{i,j}(t)-{I}_{i,j}^{\mathrm{fit}}(t))}^{2}{/\sigma }_{{ci},j}^{2}\end{eqnarray} \tag{ 3 }$

and σ_ci,j is the count rate error of I_i,j(t).

The data analysis method based on this equation has been shown to be valid, useful and successful in several previous papers (Kuwabara et al. 2004; Munakata et al. 2005; Okazaki et al. 2008; Fushishita et al. 2010; Rockenbach et al. 2011, 2014; Kozai et al. 2014, 2016). The derived anisotropy vector ${{\boldsymbol{\xi }}}^{\mathrm{GEO}}(t)$ in the GEO coordinate system is then transformed to ${{\boldsymbol{\xi }}}^{\mathrm{GSE}}(t)$ in the geocentric solar ecliptic (GSE) coordinate system for comparisons with the solar wind and IMF data. As seen in Equation (1), the observed I_i,j(t) varies depending on both I₀(t) and ${{\boldsymbol{\xi }}}^{\mathrm{GEO}}(t)$ , which must be determined separately from I_i,j(t). Although the muon count rate in each detector of the GMDN is much smaller than that in GRAPES-3, the global sky-coverage of the GMDN seen in Figure 1 makes it possible to determine the anisotropy and density, separately and accurately. We will discuss this in Section 3.

2.3. Results

Figure 2(b) displays the best-fit density (I₀(t)) during seven days including the CRB on 2015 June 22, together with 10 minute averages of the solar wind velocity and the IMF strength in Figure 2(a). It is seen that I₀(t) decreases in response to each arrival of three shocks indicated by the vertical gray lines, while it starts increasing in the second half of June 22 before the third shock arrival toward the local maximum in the same day and rapidly decreases afterward. The local maximum of I₀(t) on June 22 is lower than the level before the second shock arrival, but higher than the level before the third shock arrival. As will be shown later, this implies that cosmic rays forming the local maximum are transported from the region between the second and third shocks where the cosmic-ray intensity is less depleted by the second shock than at Earth. We note that the MHD simulation of the space weather actually reproduces such a region being formed by the third shock overtaking the second shock, which extends to the north of Earth when the third shock arrived at Earth.¹³

Figures 3(a)–(c) display 1 minute solar wind data from the OMNIWeb data set during the second half of June 22 including the CRB,¹⁴ while solid circles in panels (d)–(g) show I_i,j(t) recorded in four vertical channels of the GMDN for the same period. Following abrupt increases in the solar wind velocity and IMF strength at 18:39 UT in Figure 3(a), indicated by a gray vertical line, the strong southward IMF (B_z < 0 shown by the red curve in Figure 3(c)) discontinuously reduces, while B_x (<0) and B_y (>0) (shown by blue and green curves, respectively) become significant with similar magnitudes at 19:46 indicated by the vertical dotted line. This indicates Earth's crossing of the tangential discontinuity or the heliospheric current sheet (HCS) at this time and the southward IMF changing its orientation to be almost parallel to the nominal Parker field (B_lat ∼ 0°, B_long ∼ 135°) directed away from the Sun along the Archimedean line. By using 10 minute averages of the IMF, we calculate the normal vector ( n^HCS) to this HCS from the vector product $({{\boldsymbol{B}}}^{{\rm{U}}}\,\times \,{{\boldsymbol{B}}}^{{\rm{D}}})$ between the observed southward IMF B^D (B^D_lat = −74 fdg 3, B^D_long = 327 fdg 5) in the down-wind direction and the IMF B^U (B^U_lat = −0 fdg 5, B^U_long = 133 fdg 3) in the up-wind direction (Burlaga & Ness 1969). The GSE-latitude and longitude of the calculated normal vector are −3 fdg 7 and 227 fdg 1, respectively, indicating that the HCS at Earth was nearly perpendicular to the ecliptic plane. Figure 4 illustrates the geometries of the IMF and HCS.

**Figure 3.** Solar wind parameters and muon count rates observed by four vertical channels of the GMDN between 12:00 and 24:00 UT of 2015 June 22. Figures 3(a)–(c) display 1 minute solar wind data from the OMNIWeb data set, while solid circles in panels (d)–(g) show I_i,j(t) recorded in four vertical channels of the GMDN. Each panel displays (a) solar wind velocity (blue curve) on the left vertical axis and IMF strength (red curve) and its dispersion (green curve) on the right vertical axis, (b) proton density (red curve) and temperature (blue curve) on the left and right vertical axes, respectively, (c) three GSE-components of IMF, and (d)–(g) 10 minute muon count rates recorded in four vertical channels of GMDN (solid circles) each with the count rate error. Panels (a)–(c) show 1 minute data, while gray curves in (d), (e), and (g) also display 1 minute count rates (σ_ci,j). Only 10 minute data are available from "Kuwait City" in (f). The dotted curves in panels (d)–(g) display the best-fit count rate ${I}_{i,j}^{\mathrm{fit}}(t)$ , while red and blue curves show contributions from the anisotropy and density to ${I}_{i,j}^{\mathrm{fit}}(t)$ (see the text). Vertical gray solid and dotted lines indicate arrival times of the strong shock and HCS, respectively.
Download figure:
Standard image High-resolution image

**Figure 3.** Solar wind parameters and muon count rates observed by four vertical channels of the GMDN between 12:00 and 24:00 UT of 2015 June 22. Figures 3(a)–(c) display 1 minute solar wind data from the OMNIWeb data set, while solid circles in panels (d)–(g) show I_i,j(t) recorded in four vertical channels of the GMDN. Each panel displays (a) solar wind velocity (blue curve) on the left vertical axis and IMF strength (red curve) and its dispersion (green curve) on the right vertical axis, (b) proton density (red curve) and temperature (blue curve) on the left and right vertical axes, respectively, (c) three GSE-components of IMF, and (d)–(g) 10 minute muon count rates recorded in four vertical channels of GMDN (solid circles) each with the count rate error. Panels (a)–(c) show 1 minute data, while gray curves in (d), (e), and (g) also display 1 minute count rates (σ_ci,j). Only 10 minute data are available from "Kuwait City" in (f). The dotted curves in panels (d)–(g) display the best-fit count rate ${I}_{i,j}^{\mathrm{fit}}(t)$ , while red and blue curves show contributions from the anisotropy and density to ${I}_{i,j}^{\mathrm{fit}}(t)$ (see the text). Vertical gray solid and dotted lines indicate arrival times of the strong shock and HCS, respectively.
Download figure:
Standard image High-resolution image

**Figure 4.** Illustration of the IMF and anisotropy on both sides of the HCS. The left panel displays the up-wind and down-wind IMF vectors $({{\boldsymbol{B}}}^{{\rm{U}}}$ and B^D) by green vectors on both sides of the HCS indicated by a plane with gray frames, while the right panel illustrates the solar wind velocity $({{\boldsymbol{V}}}_{\mathrm{SW}})$ and the velocity of Earth's orbital motion $({{\boldsymbol{v}}}_{{\rm{E}}})$ , which are used for the corrections of the solar wind convection and Compton–Getting effect, by unfilled and gray filled arrows, respectively. Black axes on both panels represent the GSE coordinate system, while blue axes on the right panel represent the HCS coordinate system in which the z-axis directs parallel to the normal vector of the HCS. Three red arrows in the right panel illustrate the observed cosmic-ray streaming orientation (opposite to the anisotropy) rapidly changing across the HCS (see Figure 5(c) and the text).
Download figure:
Standard image High-resolution image

**Figure 4.** Illustration of the IMF and anisotropy on both sides of the HCS. The left panel displays the up-wind and down-wind IMF vectors $({{\boldsymbol{B}}}^{{\rm{U}}}$ and B^D) by green vectors on both sides of the HCS indicated by a plane with gray frames, while the right panel illustrates the solar wind velocity $({{\boldsymbol{V}}}_{\mathrm{SW}})$ and the velocity of Earth's orbital motion $({{\boldsymbol{v}}}_{{\rm{E}}})$ , which are used for the corrections of the solar wind convection and Compton–Getting effect, by unfilled and gray filled arrows, respectively. Black axes on both panels represent the GSE coordinate system, while blue axes on the right panel represent the HCS coordinate system in which the z-axis directs parallel to the normal vector of the HCS. Three red arrows in the right panel illustrate the observed cosmic-ray streaming orientation (opposite to the anisotropy) rapidly changing across the HCS (see Figure 5(c) and the text).
Download figure:
Standard image High-resolution image

The muon count rates in Figures 3(d)–(f) show local maxima around this HCS at similar times to the CRB reported by Paper I, except the rate of Brazilian detector in Figure 3(g), which views in almost the opposite direction to Nagoya-V in Figure 3(d) (see Figure 1). This indicates a significant contribution from the global anisotropy to the CRB. The dotted curves in Figures 3(d)–(g) show ${I}_{i,j}^{\mathrm{fit}}(t)$ reproduced using the best-fit parameters in Equation (1), while the red and blue curves represent contributions to ${I}_{i,j}^{\mathrm{fit}}(t)$ from the cosmic-ray anisotropy ( ${{\boldsymbol{\xi }}}^{\mathrm{GEO}}(t)$ ) and density (I₀(t)), respectively. It is clear that the contribution from the anisotropy (red curves) is out-of-phase for Nagoya in Japan and São Martinho da Serra in Brazil, enhancing (canceling) the common increase due to the density (blue curves) in Nagoya (São Martinho da Serra).

Figure 5 shows the best-fit cosmic-ray density and anisotropy observed in every 10 minutes, over the same period as Figure 3. The shock arrival and the HCS crossing are again indicated by vertical gray and dotted lines, respectively. The error of each parameter is deduced from the muon count rate error (σ_ci,j) (see Figure 6 in Section 3) and the dispersion of 1 minute values of each IMF parameter. It is seen in Figure 5(a) that the cosmic-ray density (I₀(t)) gradually increases toward a local maximum just at the HCS crossing and rapidly decreases afterward. Due to this local maximum, the decrease of I₀(t) appears to start about one hour after the shock arrival. The amplitude of the anisotropy ( ${{\boldsymbol{\xi }}}^{\mathrm{GSE}}(t)$ ) displayed by the green curve in Figure 5(b) also shows a local maximum around HCS and remains large at ∼1% afterward. It is clear that the anisotropy component perpendicular to the IMF (red circles) dominates the total anisotropy (green curve) around the HCS crossing. A striking feature is the anisotropy orientation changing systematically around the HCS in Figure 5(c). This implies that an anisotropy vector with a significant amplitude rapidly passed the field of view (FOV) of GMDN detectors and was observed as CRBs in some directional channels.

**Figure 6.** Errors of the best-fit parameters evaluated from the count rate error (σ_ci,j). From left to right, each panel displays the error of I₀(t), ${\xi }_{x}^{\mathrm{GEO}}(t)$ , ${\xi }_{y}^{\mathrm{GEO}}(t)$ , and ${\xi }_{z}^{\mathrm{GEO}}(t)$ as a function of local time (LT) of Nagoya. Shown in each panel are errors calculated in three cases, (a) best-fitting with only Nagoya data (black curves), (b) best-fitting with only Nagoya data but with count rates virtually enlarged 16 times mimic the GRAPE-3 (blue curves), and (c) best-fitting with the GMDN data (red curves). The vertical gray solid and dotted lines indicate the LT when the asymptotic viewing direction of Nagoya-V directs along the GEO y- and x-axes, respectively (see the text).
Download figure:
Standard image High-resolution image

**Figure 6.** Errors of the best-fit parameters evaluated from the count rate error (σ_ci,j). From left to right, each panel displays the error of I₀(t), ${\xi }_{x}^{\mathrm{GEO}}(t)$ , ${\xi }_{y}^{\mathrm{GEO}}(t)$ , and ${\xi }_{z}^{\mathrm{GEO}}(t)$ as a function of local time (LT) of Nagoya. Shown in each panel are errors calculated in three cases, (a) best-fitting with only Nagoya data (black curves), (b) best-fitting with only Nagoya data but with count rates virtually enlarged 16 times mimic the GRAPE-3 (blue curves), and (c) best-fitting with the GMDN data (red curves). The vertical gray solid and dotted lines indicate the LT when the asymptotic viewing direction of Nagoya-V directs along the GEO y- and x-axes, respectively (see the text).
Download figure:
Standard image High-resolution image

Shown in Figure 5(d) is the anisotropy after subtracting the solar wind convection and the Compton–Getting anisotropy arising from Earth's orbit around the Sun (Amenomori et al. 2004), calculated by using 10 minute average of the solar wind velocity V_SW (blue curve in Figure 3(a)) and the velocity of Earth's orbital motion v_E of 30 km s⁻¹ opposite to the GSE-y axis (see Figure 4). These corrections are made by adding a vector $(2+\gamma )[{{\boldsymbol{V}}}_{\mathrm{SW}}-{{\boldsymbol{v}}}_{{\rm{E}}}]/{\text{}}c$ to ${{\boldsymbol{\xi }}}^{\mathrm{GSE}}(t)$ , where c is the speed of light and γ is the power-law index of the GCR energy spectrum set to 2.7 (e.g., Okazaki et al. 2008). This corrected anisotropy is solely due to the diffusion and the diamagnetic drift, which reflect the spatial distribution of cosmic rays. Seen more clearly in this figure than in Figure 5(b) is a local enhancement of the anisotropy (green curve) around the HCS. This enhanced anisotropy causes the rapid change in the anisotropy orientation seen in Figure 5(c).

The duration of the rapid change of the anisotropy is about two hours and comparable to the timescale (R_L/V_SW ∼ 2.7 hr) for the solar wind to travel across the Larmor radius (R_L ∼ 0.04 au) of 60 GV cosmic rays in IMF (B ∼ 35 nT) with the average velocity (V_SW ∼ 650 km s⁻¹) in the period between 19:00 and 21:00 UT. This implies that the anisotropy and spatial distribution of cosmic rays with large R_L varies only gradually, rather than instantaneously responding to an abrupt change in the IMF orientation across the HCS.

As seen in Figure 5(d), the local enhancement of the anisotropy is dominated by the component perpendicular to the IMF (red circles). By ignoring the contribution from the perpendicular diffusion, i.e., assuming that the perpendicular anisotropy is solely arising from diamagnetic drift, we deduce the spatial density gradient vector (G) perpendicular to the IMF by using the 10 minute average of the observed IMF in Figure 3. Figure 5(e) shows the three GSE components of G. The positive G_z (red circles) increases following the HCS crossing indicate that the "source" from which cosmic rays are transported, probably by the drift along HCS, is located north of Earth. It is noted here that the drift transporting cosmic rays may alter the preexisting spatial distribution of cosmic rays, but it cannot be directly observed as a directional anisotropy along the HCS. Instead, the spatial distribution of cosmic rays produced by the drift is observed by the diamagnetic drift anisotropy, which is perpendicular to the HCS.

Figure 5(f) shows the anisotropy components parallel and perpendicular to the HCS calculated by using the HCS normal vector $({{\boldsymbol{n}}}^{\mathrm{HCS}})$ defined above. The parallel component (blue circles) dominates the anisotropy after the local enhancement of the anisotropy, while the perpendicular component (red circle) becomes dominant around the HCS. We calculate the anisotropy orientation in a coordinate system fixed to the HCS in which the z-axis is parallel to n^HCS, the x-axis is parallel to the up-wind IMF $({{\boldsymbol{B}}}^{{\rm{U}}})$ and the y-axis completes the right-handed coordinate system (see blue axes on the right panel of Figure 4). The calculated longitude and latitude of the anisotropy vector in this coordinate system are shown in Figure 5(g). The anisotropy latitude and longitude are both around zero after the HCS crossing, indicating that the anisotropy is consistent with the streaming along the HCS; nevertheless, the observed local IMF orientation shows rapid and large variations (see Figure 3(c)).

3. Summary and Discussions

We analyzed the CRB observed by the GMDN on 2015 June 22 and found that the CRB was caused by the enhanced cosmic-ray anisotropy and density outside the geomagnetic field around the HCS, where the southward IMF changed its orientation to the nominal Parker field. While the IMF orientation changes almost instantaneously at the HCS, the anisotropy of 60 GV cosmic rays varies gradually over about two hours. This enhanced anisotropy with a significant amplitude rapidly passed the FOV of the GMDN detectors and was observed as "bursts" in some directional channels. This is the first successful result demonstrating that the GMDN is a useful tool for deriving the global cosmic-ray anisotropy for timescales shorter than an hour.

By analyzing the muon count rates recorded in nine directional channels of GRAPES-3, Paper I concluded that the CRB was unlikely to be caused by a cosmic-ray anisotropy outside the magnetosphere. Determining the contribution of the global anisotropy to the observed muon count rate from the observation at a single location on Earth, however, is difficult even for a large detector such as GRAPES-3, particularly when the anisotropy amplitude and orientation rapidly change as indicated in the present paper. The GRAPES-3 is an excellent muon detector with a large detection area (560 m²) and muon count rate, which is about 16 times that of Nagoya (36 m²), while the width of its FOV is similar to that of Nagoya (Nonaka et al. 2006). As seen in Equation (1), the first-order anisotropy vector ${{\boldsymbol{\xi }}}^{\mathrm{GEO}}(t)$ produces broad excess and deficit of the relative intensity, each spreading over 180° of longitude. A single muon detector with a limited FOV, therefore, records similar count rates in all directional channels when it monitors the direction parallel or anti-parallel to ${{\boldsymbol{\xi }}}^{\mathrm{GEO}}(t)$ , being unable to accurately observe the anisotropy. In order to show this quantitatively, we calculate errors of I₀(t) and ${{\boldsymbol{\xi }}}^{\mathrm{GEO}}(t)$ , each as a function of local time (LT) of Nagoya, by propagating the count rate error (σ_ci,j) in Equations (2)–(3). Figure 6 displays errors of I₀(t), ${\xi }_{x}^{\mathrm{GEO}}(t)$ , ${\xi }_{y}^{\mathrm{GEO}}(t)$ , and ${\xi }_{z}^{\mathrm{GEO}}(t)$ calculated in three cases, (a) best-fitting with only Nagoya data (black curves), (b) best-fitting with only Nagoya data but with the count rates virtually enlarged 16 times to mimic the GRAPE-3 (blue curves) and (c) best-fitting with the GMDN data (red curves). It is clear from this figure that the error of each parameter in case (a) of best-fitting with Nagoya data alone is much larger than the error in case (c) with the GMDN data. The longitude of asymptotic viewing direction of Nagoya-V is 168 fdg 9, while the longitude of the detector's location is 137 fdg 0 (see Table 1 and Figure 1). The center of Nagoya's FOV, therefore, directs toward a certain orientation in space ∼2 hr earlier than the corresponding LT at the detector's location and views the direction along the GEO y-axis at 90° (GEO x-axis at 0°) GEO longitude at ∼04:00 LT (22:00 LT), as indicated by the vertical gray solid (dotted) line. It is seen in case (a) that the error of ${\xi }_{y}^{\mathrm{GEO}}(t)$ becomes largest when the FOV of Nagoya directs toward the GEO-y axis, while the error of ${\xi }_{x}^{\mathrm{GEO}}(t)$ becomes smallest at the same local time when Nagoya monitors both the excess and deficit due to ${\xi }_{x}^{\mathrm{GEO}}(t)$ in the FOV. Errors actually become much smaller in case (b) with an enlarged detection area and the reduced σ_ci,j, but still larger than case (c) most of the time. As seen in this figure, errors of ${\xi }_{x}^{\mathrm{GEO}}(t)$ and ${\xi }_{y}^{\mathrm{GEO}}(t)$ in cases (a) and (b) are significantly dependent on the local time, while those in case (c) are much more stable. This is a serious problem for a single detector in observing ${{\boldsymbol{\xi }}}^{\mathrm{GEO}}(t)$ with unknown orientation and amplitude.

Figure 6 demonstrates the difficulty of deriving an accurate anisotropy from muon rates in the FOV of a single detector as a function of time. If the anisotropy can be regarded as constant during a day, on the other hand, even a single detector can observe ${\xi }_{x}^{\mathrm{GEO}}$ and ${\xi }_{y}^{\mathrm{GEO}}$ accurately, by measuring the diurnal variations of muon rates in multidirectional channels (Munakata et al. 2014). In this case, the diurnal variation is generally observed with different phases in different directional channels, because the eastern directional channels observe the anisotropy orientation in space earlier than the western channels according to Earth's spin. If the difference of asymptotic longitudes viewed by the eastern and western channels is Δϕ_asymp, the phase difference of Δϕ_asymp/ω (1 hr/15°) is expected. Paper I concluded that the CRB is not caused by the anisotropy, because no such phase difference is observed in the CRB by GRAPES-3 (see also Mohanty et al. 2018). Their deduction, however, is valid only when the anisotropy is constant and cannot be applied to a CRB in which the anisotropy rapidly changes in a few hours. It is also noted here that, even in the case of constant anisotropy, a single detector cannot observe ${\xi }_{z}^{\mathrm{GEO}}$ because it produces no diurnal variation.

We finally discuss the effect of a major geomagnetic storm recorded in the CRB period. The maximum K_p index of 8+ was recorded during three hours between 18:00 and 21:00 UT of June 22, while the minimum D_st index of −204 nT was recorded at 05:00 UT of June 23.¹⁵ Paper I concluded that the observed CRB is caused by this geomagnetic storm, because the geomagnetic field is weakened during the storm and the geomagnetic cutoff rigidity (P_c) of cosmic rays is reduced, allowing more low energy particles to reach the ground level detectors and possibly increasing the muon count rate. To estimate this effect quantitatively, we performed numerical calculations of cosmic-ray orbits in the latest model of the geomagnetic field (TS05),¹⁶ which is capable of reproducing the geomagnetic field even in a major storm by using the observed solar wind data as inputs (Tsyganenko & Sitnov 2005). In our calculation, we use five-minute averages of these data and reproduce the geomagnetic field every five minutes. This model is also capable of reproducing the variation of D_st index in a good agreement with the observations. Each upper panel of Figure 7 shows the reproduced D_st index (gray curve) and the deviation (ΔP_c) of P_c (black curve) from its nominal value calculated for each vertical channel of GMDN during a period between 18:00 and 22:00 UT on June 22. It is seen that ΔP_c varies roughly in a positive correlation with D_st index and is temporally almost universal and common for the four directional channels, regardless of the location of each detector (IceCube Collaboration et al. 2017). The amplitude of variation of ΔP_c is, on the other hand, significantly different from one location to the other (see the range of the left vertical axis). By comparing one hour average of D_st index reproduced by the model (gray curve) with the observed hourly D_st index (gray diamond), we verified a good correlation with the correlation coefficient of 0.86 and the regression coefficient of 1.05 together with ∼10% off-set of the reproduced minimum D_st index. Based on this, we estimate the uncertainty of the reproduced geomagnetic field to be about 10% which is accurate enough for examining the potential effect of the geomagnetic storm on the CRB.

Listed in the last three columns of Table 1 are the nominal P_c and its maximum reduction ΔP_c during a period in Figure 7. By integrating the response function of the atmospheric muon intensity to primary GCRs with respect to the rigidity above P_c (Murakami et al. 1979; Fujimoto et al. 1984), we calculate the expected increase of muon count rate ΔI_i,j(t) as listed in Table 1. The expected ΔI_i,j(t) is only ∼0.3% at most and much smaller than the observed amplitude of the CRB in Figure 3. The reduction of P_c is −0.64 GV and second largest for Hobart-V, but ΔI_i,j(t) is almost zero. This is because P_c for Hobart-V is already lower than the atmospheric threshold rigidity for producing muons and the reduction of P_c causes no further increase in the muon count rate. Contrary to this, the CRB is also recorded in Hobart-V as seen in Figure 3(e). The largest ΔI_i,j(t) would be expected in São Martinho da Serra-V, but no CRB is observed in this channel, as shown in Figure 3(g). The observed muon count rate in the lower panel of Figure 7 actually shows no clear correlation with ΔP_c in the upper panel. We also confirmed that the asymptotic viewing direction during the storm varies only a few degrees at most, which is insufficient to affect the directional response of muon detectors to primary cosmic-ray intensity and result in the observed CRB. Based on these results, we conclude that the CRB is not caused by the geomagnetic storm.

This work is supported in part by the joint research programs of the Institute for Space-Earth Environmental Research (ISEE), Nagoya University and the Institute for Cosmic Ray Research (ICRR), University of Tokyo. The observations are supported by Nagoya University with the Nagoya muon detector, by CAPES, INPE, and UFSM with the São Martinho da Serra muon detector, and by the Australian Antarctic Division with the Hobart muon detector. The observation with the Kuwait City muon detector is supported by the project SP01/09 of the Research Administration of Kuwait University. A.D.L. thanks CNPq (grant 304209/2014/7). C.R.B. thanks to Sao Paulo Research Foundation (grant 14/24711-6). R.R.S.M. thanks CNPq (grant 152050/2016-7) and the China-Brazil Joint Laboratory for Space Weather. J.K. was partly supported by the Joint Research Program of the Institute for Space-Earth Environmental Research, Nagoya University, Japan. Authors thank Dr. Ryuho Kataoka of the National Institute of Polar Research in Japan for providing us with the information of Tsyaganeno's geomagnetic field models. They also thank Dr. Yoji Hasegawa of Shinshu University for his help in performing MC simulations by GEANT4. The OMNIWeb data set of the solar wind and IMF parameters is provided by the Goddard Space Flight Center, NASA, USA, while the MHD simulation of the space weather during the CRB is provided by the WSA-Enlil solar wind prediction, National Centers of Environmental Information, NOAA, USA. The hourly D_st index is provided by the WDC for Geomagnetism, Kyoto, Japan. We are grateful to the anonymous referee for thoughtful comments that were useful for developing a clearer presentation of this paper.

Cosmic-Ray Short Burst Observed with the Global Muon Detector Network (GMDN) on 2015 June 22

Article metrics

Permissions

Author affiliations

ORCID iDs

Dates

Abstract

1. Introduction