Interplanetary Magnetic Flux Rope Observed at Ground Level by HAWC

The High Altitude Water Cerenkov observatory is located on a plateau between the Sierra Negra and Pico de Orizaba volcanoes in Mexico, at 18°59'41''N, 97°18'306W and at an altitude of 4100 m. HAWC consists of 300 water Cerenkov detectors (WCDs), each 7.3 m in diameter and 4.5 m deep. The WCDs are spread over an area of 22,000 m² and each one is filled with filtered water and instrumented with four photomultiplier tubes (PMTs). A 10 inch PMT is located at the center of the WCD and three 8 inch PMTs are arranged around it, making an equilateral triangle of side 3.2 m (Abeysekara & HAWC Collaboration 2012).

The main data acquisition (DAQ) system measures arrival times and time over threshold of PMT pulses. This information is used to reconstruct the air-shower front and estimate the arrival direction and energy of the primary particles. The electronics are based on time-to-digital converters (TDCs), and the main DAQ also has a scaler system that counts the hits inside a time window of 30 ns for each PMT (R₁ rate from now on) and the coincidences of two, three, and four PMTs in each WCD, called multiplicity M₂, M₃, and M₄, respectively. The TDC–scaler system allows one to measure particles with relatively low rigidity, from the cutoff rigidity at the HAWC site (∼8 GV) to the low limit of the reconstructed showers (∼100 GV). The median energy of observation of the TDC–scaler system and multiplicities are in the range of 40–46 GV. The large area and low cutoff rigidity of HAWC make it a suitable instrument for studying solar modulation in general and space weather in particular. Summarizing, we can say that the TDC–scaler rate used in this work provides information on the primary GCR flux above cutoff rigidity (from 8 GV onward) reaching Earth's atmosphere, which can be measured with high precision.

During long periods of observation, the efficiency of PMTs may vary, so to carry out high-precision studies, it is necessary to correct for these drifts. To perform these corrections, first we have identified relatively stable PMTs during a year by continuously checking their deviations. Using a "singular value decomposition" method, we compute a reference rate that was used to model the slow and small changes in the efficiency of the remaining PMTs and correct for their variation in efficiency. Due to the large collecting area and high altitude of HAWC, the rate of observed particles is high: for instance, during the year 2016 the efficiency-corrected average rate of the HAWC TDC–scaler system was 〈R₁〉 = 23.39 kHz per PMT; similarly for multiplicities the average rates per WCD were 〈M₂〉 = 8.06 kHz, 〈M₃〉 = 5.69 kHz, and 〈M₄〉 = 4.35 kHz. The analysis carried out in this study is using data with "one minute" resolution, and for uniformity the rates were converted to percentage deviation with respect to the mean rate of the year 2016. It is worth mentioning that after the efficiency correction, the measured TDC–scaler rates are in agreement with the expected statistical accuracy of the data. Similar to other air-shower detectors, the TDC–scaler rates show a dependence on barometric pressure. We correct this pressure modulation following the method described by Arunbabu et al. (2019). Finally, due to HAWC's location near the equator and its lower cutoff rigidity of 8 GV, the solar diurnal anisotropy (DA) component is strongly significant in its low-energy observations. The DA was removed from the data using a band rejection filter that removes all the frequencies within the frequency range of ∼1–2 cycles per day.

The efficiency-corrected TDC–scaler rates of all the PMTs were combined to provide the R₁ rate with a statistical accuracy <0.01% per minute. The same efficiency correction process was also applied to the multiplicity rates from all the WCDs, and these were combined to get M₂, M₃, and M₄.

The TDC–scaler data after applying all these corrections are shown in the top panels of Figures 2 and 3, where the mean values of the four channels R₁ (blue), M₂ (black), M₃ (green), and M₄ (red) are shown. The second and third panels of these figures show the solar wind (SW) proton measurements of the speed V (black), the number density N (green), the magnitude of the magnetic field B (black), and its components B_X (blue), B_Y (green), and B_Z (red). The one-minute resolution SW parameters were obtained from the OMNI data service developed and supported by NASA/Goddard's Space Physics Data Facility.²⁹ The flow and magnetic pressures (fourth panel) as well as the D_ST index (bottom panel) were obtained from the World Data Center for Geomagnetism, Kyoto (Nose et al. 2015).

Zoom In Zoom Out Reset image size

In order to show the suitability of the HAWC TDC–scaler system to study the solar modulation on short and medium timescales, Figure 2 shows a time period of approximately three months at the end of 2016, where seven HSSs (marked with yellow shadow areas) and one ICME and its MFR (blue shadow) were observed. From this figure the visible correlation of HAWC rates with the SW velocity resembles the advection effect of the GCRs in the heliosphere. Also the diffusion effects of GCRs inside the heliosphere show correlated decreases in HAWC rate with the magnetic field enhancements. It should be noted that, contrary to expectation, the signature of the MFR structure caused an increase in the rates, which is fairly remarkable (the aim of this analysis is to study the origin of this increase). In a similar way, the changes in the magnetic field, flux, and magnetic pressures and D_ST index caused by this MFR are clearly distinguished from the other SW structures detected over this period.

The GCR modulation associated with the ICME passage during 2016 October observed by HAWC can be seen in Figure 3 (see Section 3.3 for the details of the figure panels). During this time, the TDC–scaler rates increased as a double peak structure, starting around 02:50 UT on 2016 October 14 and lasting 16.8 hr.

Zoom In Zoom Out Reset image size

The HAWC TDC–scaler rates have a background rate of 4.25 × 10⁸ particles per minute. During this event, the integrated counts had an excess of 1.76 × 10⁹ particles above the background level. In order to quantify the noise level of the TDC–scaler system, we computed the mean rates and standard deviation of the 1200 PMT rates for R₁ and 300 WCD rates for multiplicities M₂, M₃, and M₄ every minute. The upper panels of Figure 4 show the distributions of the computed standard deviations during three days of our period of interest. The mean values of these distributions are used as the standard deviation of the observations (σ) for this study and are listed in the second column of Table 1. The magnitude of the observed double peak structure as a percentage is estimated as the difference between the peak and pre-event intensities and is shown in the third and fifth columns of Table 1. In a similar way, the magnitude in terms of standard deviation (σ), i.e., magnitude/σ, which we call the significance of the increase, is shown in the fourth and sixth columns.

Zoom In Zoom Out Reset image size

It should be noted that during the event, the weather at the HAWC site was calm and no electric field disturbances were recorded by the site electric field monitor. The atmospheric pressure at the site was also showing normal behavior, which rules out the possibility of abnormal atmospheric activities. Earthquakes can also cause rate increases in the TDC–scaler system but these last a few minutes and there is no report of an earthquake of magnitude greater than 5.5 during our period of interest in the south-central part of Mexico, where HAWC is sited.

To compare this event with other short-term modulations observed by the TDC–scaler system, we applied a high-pass filter to the rates. This filter removed all the frequencies smaller than 1 day⁻¹ and retained all the fast variations in the data that have timescales less than a day. The bottom panel of Figure 4 shows the high-pass-filtered data during three months, starting on 2016 September 1. The standard deviation of this high-pass-filtered data for the entire year 2016 was estimated as σ_high-pass = 0.068%. From the figure it is clear that the MFR event (marked by the blue shaded area) is significant in comparison to all other short-period modulations, and is much higher than the 5σ_high-pass level. In this figure one can also see atmospheric electric field transients (short-time spikes of tens of minutes). In contrast, the event due to the MFR lasted ∼17 hr. These atmospheric electric field events are listed in Jara Jimenez et al. (2019) except for those on September 21 and 26; these two days were not listed due to short gaps in HAWC TDC–scaler data.

On 2016 October 8 a filament was observed by the Atmospheric Imaging Assembly (AIA, Lemen et al. 2012) on board the Solar Dynamics Observatory (SDO, Pesnell et al. 2012) activated at ∼04:00 UT in the northeast quadrant (see panels (a) and (b) of Figure 5), rising with a speed of ∼20 km s⁻¹ during the next 20 hr. Around 00:00 UT on 2016 October 9 and at an altitude of  ∼3 R_⊙ an acceleration process started as seen by COR1 and COR2, instruments on board the Solar Terrestrial Relations Observatory (STEREO, Kaiser et al. 2008) spacecraft (green and red X symbols, respectively, in panel (e) of Figure 5), and reached a speed of ∼180 km s⁻¹.

Zoom In Zoom Out Reset image size

The filament was part of a halo CME observed on 2016 October 9 at 02:45 UT by the Large Angle and Spectrometric Coronagraph (LASCO) on board the Solar & Heliospheric Observatory (SOHO). This eruption was seen as a limb CME, showing a dark cavity surrounded by bright material (typical structure of the cross section of an MFR, e.g., Dere et al. 1999) by the Sun Earth Connection Coronal and Heliospheric Investigation (SECCHI) on board the STEREO A spacecraft (panels (c) and (f) of Figure 5, respectively).

The ICME propagation from the Sun to the Earth of the 2016 October event has been studied previously by He et al. (2018), who combined remote sensing observations from SDO, STEREO, and SOHO and in situ measurements as extracted from OMNI data. These authors focused their work on understanding the geo-effectiveness associated with the CME magnetic field structure. They reconstructed the CME using the Grad–Shafranov technique to identify the MFR with a strong southward magnetic field, which produced the geomagnetic storm (see Figure 8 of He et al. 2018). They also simulated the SW with the ENLIL-MHD code (Odstrcil et al. 2004) to illustrate the background SW conditions during the ICME propagation. However, they only considered the evolution of the ICME bracketed between two HSSs, without considering that the whole of this system is propagating within a variable SW. We note that in the case of slow CMEs the Sun–Earth travel time may be a few days and therefore they have the opportunity to interact with other SW structures and/or transients that can enhance or ensure their geo-effectiveness (e.g., He et al. 2018; Chen et al. 2019).

To include this ingredient in the scenario (a variable ambient SW), we perform a simulation that includes the observed parcels of SW with different speeds and densities, including the slow ICME and the slow and fast flows (Section 3.4).

When an interplanetary transient is not fast enough, it may still be able to modulate the GCR flux depending on its magnetic field strength and geometry, its dynamics compared with ambient SW, and its interaction with the Earth's magnetic field. This means that, in general, the GCR modulation depends mainly on three factors: the nature of the transient and its intrinsic magnetic field, the resulting conditions after the transient interacts with the ambient SW, and the Earth's magnetosphere.

During 2016 October 12–17, a slow ICME was identified at 1 au using the one-minute resolution in situ measurements (this event has been reported by Nitta & Mulligan 2017; He et al. 2018).

At this time, the TDC–scaler subsystem of the HAWC observed a peculiar GCR flux modulation in the form of a double peak enhancement as shown in the top panel of Figure 3, where HAWC TDC–scaler data—R₁, M₂, M₃, and M₄—are presented. This figure also shows the SW parameters: V, N (second panel), B-field (third panel), the flow and magnetic pressure (fourth panel), as well as the D_ST index (bottom panel). We have marked with a dashed vertical red line the arrival time of the shock at 1 au, whereas the two dotted vertical blue lines mark the starting and ending times of the CME structure at 1 au as discussed by Nitta & Mulligan (2017) and He et al. (2018). The ICME boundaries were selected based on the signatures described on Zurbuchen & Richardson (2006) and Richardson & Cane (2010). There is a shock wave on 2016 October 12 at 22:01 UT and a sheath that lasted for 7.68 hr. The MFR boundaries are very well defined with the magnetic field and plasma parameters, starting on 2016 October 13 at 06:00 UT and ending on 2016 October 14 at 16:00 UT with a bulk speed of 390 km s⁻¹.

Due to the well-defined magnetic structure with clear rotation of B_Y and B_Z (in situ measurements), the slow CME has been identified as an MC or MFR with a strong southward magnetic field (He et al. 2018). As stated above, this MFR was the result of the expulsion of a quiet filament with no obvious eruptive signatures in remote sensing observations (Nitta & Mulligan 2017). At 1 au, the conditions that give rise to the increase in GCRs started with a pre-event HSS (October 1) sweeping the interplanetary medium and imposing very quiet SW conditions (V ∼ 350 km s⁻¹ and N ∼ 7 particles cm⁻³) in front of the ICME. This smooth SW coupled with the low ICME speed (∼450 km s⁻¹) led to a weak shock and small compression region (sheath, marked with red and blue vertical lines around October 13 in Figure 3). The magnetic pressure (${P}_{{\rm{MAG}}}{({\rm{nPa}})=({10}^{-2}/8\pi )B({\rm{nT}})}^{2}$, black) and the SW flow pressure (computed as ${P}_{{\rm{FLOW}}}{({\rm{nPa}})=2\times {10}^{-6}N({{\rm{particles}}{\rm{cm}}}^{-3})V({{\rm{km}}{\rm{s}}}^{-1})}^{2}$, red) are shown in the fourth panel of Figure 3, to make it clear that this slow ICME is related to a strong magnetic field disturbance. It can be seen that within the MFR (between the two dotted vertical blue lines), P_MAG increases while P_FLOW decreases, both considerably in comparison to the pressure conditions before the arrival of the MFR (i.e., before the dashed vertical red line). Moreover, the MFR arrival at the Earth produced a relatively intense geomagnetic storm reaching a minimum D_ST index of ∼−104 nT (He et al. 2018) as shown in the bottom panel of Figure 3.

It is important to note that during 2016 October solar activity was in the declining phase of Cycle 24, therefore the SW conditions that we simulate correspond to an interval of low activity in the inner heliosphere. As the CME propagates 3° east of the Sun–Earth line and 11° north of the ecliptic plane (for further information see Figure 2 and its description from He et al. 2018), we assume that it propagates radially and is directed toward the Earth, and therefore we neglect the deflection and rotation effects (Kay & Opher 2015) as well as its interaction with the SW magnetic field (Schwenn et al. 2005). Under these assumptions, hydrodynamic codes give reliable insights into the dynamics of CMEs in the SW (Niembro et al. 2019).

The simulation was done using the 2D hydrodynamic, adiabatic, and adaptive grid code YGUAZÚ-A, developed by Raga et al. (2000) and used previously to simulate the propagation of ICMEs in the SW (see Niembro et al. 2019, for the details of the code and its application to ICMEs). We assume a standard SW (e.g., Kippenhahn et al. 2012) with a mean molecular weight³⁰ μ = 0.62 to consider a chemical composition with solar mass fractions of 0.7 H, 0.28 He, and 0.02 of the rest of the elements; a specific heat at constant volume of γ = 5/3, and an initial temperature of 10⁵ K (Wilson et al. 2018).

The SW conditions at the injection radius R_inj = 6 R_⊙ (4.2 × 10⁶ km, ∼0.028 au) were estimated by splitting the SW measurements at 1 au into different time periods i = A, B, C, ..., I, in which the SW showed a speed with a relatively stationary behavior. With this, we estimated the median values of the speed $\overline{{V}_{1{\mathrm{au}}_{i}}}$ and the number density $\overline{{N}_{1{\mathrm{au}}_{i}}}$. Each period is characterized with these median values. Then, we computed the mass-loss rate as $\dot{{M}_{i}}=4\mu {m}_{p}\pi {R}_{1\,{\rm{au}}}^{2}\,\bar{{N}_{1{{\rm{au}}}_{i}}}\,\bar{{V}_{1{{\rm{au}}}_{i}}}$ with the proton mass m_p = 1.67 × 10⁻²⁷ kg. With all these parameters obtained from the measurements at 1 au, we obtained the speed ${V}_{{\mathrm{inj}}_{i}}$ at the injection point R_inj (near the Sun) and injection time Δt_i by solving the system of equations described in Cantó et al. (2005). All these values are shown in Table 2, with the MFR corresponding to period C. The computational domain is filled with the conditions described in period A (corresponding to the ambient SW conditions) and the initial injection time is given by the model using the SW speed observed within period B.

Table 2.  Median Values of the SW Speed $\overline{{V}_{1\mathrm{au}}}$ and Number Density $\overline{{N}_{1\mathrm{au}}}$

At 1 au						Near the Sun
(from OMNI data)						(from Cantó et al. 2005)
Period	Initial Datea	End Datea	$\overline{{V}_{1\mathrm{au}}}$	$\overline{{N}_{1\mathrm{au}}}$	$\dot{M}$ × 10−16	Δt	Vinj
			(km s−1)	(cm−3)	(M⊙ yr−1)	(hr)	(km s−1)
A	2016/10/12 00:00	2016/10/12 11:55	372	5.86	2.78	b	363
B	2016/10/12 20:35	2016/10/12 22:15	340	4.82	1.23	29.32	330
C	2016/10/12 23:18	2016/10/13 11:08	431	17.4	5.61	11.83	423
D	2016/10/14 08:48	2016/10/14 16:08	360	7.49	2.02	16.5	351
E	2016/10/14 16:53	2016/10/14 18:33	405	5.53	1.68	1.66	397
F	2016/10/14 21:33	2016/10/14 23:23	358	10.96	2.94	16.79	349
G	2016/10/15 01:03	2016/10/16 06:30	502	4.51	1.69	66.71	495
H	2016/10/16 08:43	2016/10/16 17:03	665	3.87	1.93	44.05	660
I	2016/10/16 18:58	2016/10/18 04:18c	705	2.32	2.25	56.38	700

Notes. The values are extracted from the OMNI data set and the computed mass-loss rate $\dot{M}$ used to reconstruct the SW conditions (V_inj, Δt) near the Sun at the injection radius R_inj = 6 R_⊙ obtained by solving the system of equations described in Cantó et al. (2005). We have tagged with different colors the periods to be identified and followed in Figures 3 and 6.

^a(yyyy/mm/dd hh:mm UT.) ^bPeriod A corresponds to the SW conditions needed to fill in the computational domain in which the different periods will be injected. ^cThe end of the period is not shown in Figures 3 and 6.

Download table as:  ASCIITypeset image

The synthetic speed and number density profiles resulting from the simulation are shown as orange solid lines in the second panel of Figure 3 to be compared with the in situ measurements at 1 au (shown in the figure in black and green, respectively). It can be noticed that our approach accurately resembles the observed speed profile, while the number density profile is consistent with the observations, except for those time ranges in which the 3D configuration of the magnetic field certainly plays a major role in the dynamics (e.g., Cargill & Schmidt 2002), such as in the compression region formed due to the interaction of the MFR with the ambient SW in front of it and delimited by the dashed red vertical line (shock) and the first dotted blue vertical line (initial time of the MFR). The other two enhancements found at later times in the number density profile (2016 October 14 ∼17:00 UT and 2016 October 16 ∼12:00 UT) are related to compression regions formed by the interaction between two flows at different speeds (between periods F and G, and G and H, respectively). It is worth noting that, from the hydrodynamics point of view, the plasma is compressed when two flows at different speeds interact, while a low-density region will be formed when the flows separate from each other,.

An important feature of this particular hydrodynamic code is that it enables us to tag and follow each different plasma parcel. This is shown with the colored symbols used over the synthetic number density profile, which are indicated in Figure 3 following the slightly lighter color code in Table 2. We would like to focus on the evolution of two of the simulated SW parcels colored green (corresponding to period B) and orange (period D). It can be noticed that the SW delimited in these two periods presents low-density cavities that can be easily tracked in Figure 6, where we show the profiles of speed (black solid line) and number density (green solid line with colored symbols) at five different heliospheric distances: 0.1 au, 0.3 au, 0.5 au, 0.7 au, and 1.0 au, from top to bottom.

Zoom In Zoom Out Reset image size

These low-density cavities prevented the interaction of the disturbance with the surrounding ambient SW, preserving the well-defined and strong magnetic structure of the MFR from its source up to 1 au. Moreover, the duration of the CME (48 hr, time delimited by dashed blue lines in Figure 3 and blue shaded area in Figure 6) does not change at different heliospheric distances, thus allowing the guiding of GCRs observed by HAWC.

Besides the spacecraft located around the L1 point, there are several spacecraft orbiting the Earth at different locations that are able to detect changes in the Earth's magnetic field caused by the ICME (shock, sheath, and flux rope) as it arrives in the Earth's vicinity. In this work, we use B-field data from THEMIS-E and THEMIS-C (Angelopoulos 2008), GOES-15 (https://www.goes.noaa.gov/goesfull.html), and MMS-1 (Burch et al. 2016) to characterize the disturbances caused by the 2016 October MFR on the Earth's magnetosphere and its possible effects on the GCR flux.

Figure 7 shows the orbits of these spacecraft in the XY (left) and XZ (right) GSE planes (the geocentric GSE is a right-handed coordinate system with the X-axis pointing toward the Sun and the Z-axis perpendicular to the ecliptic plane and positive pointing north). The triangles on the panels in Figure 7 mark the spacecraft locations at the time when the magnetic field discontinuity associated with the arrival of the ICME-driven shock was observed (column 3 in Table 3 and vertical gray lines in the top panels of Figures 15–19 in Appendix A). The squares in Figure 7 mark the location of the sudden decrease in B_Z associated with the arrival of the leading edge of the MFR at each spacecraft (column 6 in Table 3). The bow shock boundary drawn (continuous line) was calculated from a modified version of Fairfield (1971), while the location of the magnetopause (dashed line) is based on the model of Roelof & Sibeck (1993); both boundaries are shown in Figure 7 only as a reference and are based on the SW average value of P_FLOW = 1.65 nPa measured at 1 au. It is important to note that during the MFR passage these boundaries may be largely modified. According to this and in agreement with the observations, GOES-15, THEMIS-E, and MMS-1 were located inside the magnetosphere during the time interval studied, whereas THEMIS-C was in the magnetosheath region.

Zoom In Zoom Out Reset image size

Table 3.  Parameters Related to the Spacecraft Observations

s/d	Spacecraft	ts/d (UT)	Bup	Bdown	${t}_{{B}_{Z}}$	min(BZ)	Spacecraft
		mm/dd hh:mm:ss	(nT)	(nT)	mm/dd hh:mm:ss	(nT)	Regiona
s	WIND	10/12 21:15:37	3.4	8.2	10/13 05:27:30	−20.0	SW
d	GOES-15	10/12 22:11:51	121	150	10/13 08:00:00	8.0	M-SPH
d	THEMIS-E	10/12 22:13:33	92.5	123	10/13 10:19:21	−60.0	M-SPH
d	MMS-1	10/12 22:13:52	30	47	10/13 10:31:10	−50.0	M-SPH
s	THEMIS-C	10/12 22:24:49	6	10.5	10/13 06:31:37	−22.5	M-SH

Note.

^as/d: shock or discontinuity, SW: solar wind, up: upstream, down: downstream, M-SPH: magnetosphere, M-SH: magnetosheath.

Download table as:  ASCIITypeset image

In general, inside the magnetosphere the magnetic field is dominated by the north–south geomagnetic dipole, i.e., B_Z component. Therefore, we follow and compare the values of this B-field at different locations around the Earth during the MFR passage. In Figure 8, we show the difference ∣B∣ − B_Z normalized to the maximum value of ∣B∣ observed by GOES-15, which, of the spacecraft considered in this work, was the one with the shortest and constant distance to the Earth at the time of the sudden decrease in B_Z. This difference was computed for all the spacecraft (plotted with different colors) during 2016 October 12, 13, and 14 (top, middle, and bottom panels of Figure 8 respectively). The data are in GSE coordinates except for GOES-15, which are presented in EPN (earthward, parallel, normal) coordinates with respect to the spacecraft's orbital plane (see Appendix A). The individual plots of the B-field and its components measured by each spacecraft during the event are shown in Figures 15–19 in Appendix A.

Zoom In Zoom Out Reset image size

In Table 3, we summarize the times at which the spacecraft observed the shock-associated discontinuity of the B-field and their corresponding upstream/downstream B-values (columns 3, 4, and 5, respectively), the times of the sudden decrements of B_Z associated with the MFR and its minimum value (columns 6 and 7) as well as the magnetospheric location (column 8) of each analyzed spacecraft (column 2). The Geo projections of the spacecraft during the MFR passage are plotted in Figure 9.

Zoom In Zoom Out Reset image size

As the MFR enters the magnetospheric region, following the ICME sheath region, the decrease in P_FLOW (see fourth panel in Figure 3) associated with the internal region of the MFR produced an expansion of the Earth's magnetic field. This effect can be seen by comparing the October 13 time intervals 04:00–08:00 UT in Figure 16, 04:00–09:00 UT in Figure 17, and 08:00–10:00 UT in Figure 18 in Appendix A with the same time period during the previous day. This effect is clearer for GOES-15 in Figure 16 due to its geostationary orbit.

This magnetospheric weakening, along with the enhancement of the southward magnetic field component inside the MFR, caused sudden reversals (B_Z < 0) on the Earth's magnetic field observed as step-like signatures by THEMIS-C, THEMIS-E, and MMS-1 as shown in the middle panel of Figure 8, and as changes in the polarity of B_Z in Figures 17–19 in Appendix A.

Because of the unique position of GOES-15, particularly its proximity to the ecliptic plane and its low altitude, the intense local magnetic field makes it difficult to clearly detect the perturbation caused by the MFR, as shown in Figure 8 and in Figure 16 in Appendix A. However, significant disturbances were observed during most of October 13 and continued but with lower amplitude the next day, when the GCR double peak was observed by HAWC. This can be corroborated by observing the middle and bottom panels in Figure 8, showing in this way that the GCR enhancement observed by HAWC was not directly caused by these geomagnetic disturbances.

Although the GCR anisotropy induced by the MFR must be a global event, the response of NMs with cutoff rigidities similar to that of HAWC (6–9 GV), during the passage of the MFR was somewhat poor due to the fact that the sensitivity of NMs is lower than that of HAWC in this energy range. These responses are shown in the inset of Figure 10, where the rates of the selected NMs during a three-day period starting on 2016 October 13 are shown. For the sake of clarity, we have added/subtracted an offset to each time profile, and the M₃ rate of HAWC (red line) is divided by 30 to fit the scale.

Zoom In Zoom Out Reset image size

Even though the behavior of the time profile is different for each monitor, all of them show an increasing trend after day 13 and a more evident decrease after day 14.5. It is interesting to note that the Mexico City NM (magenta line) has a time profile similar to HAWC during and after the second peak, which was caused by low-energy particles (see Section 4.2). The first peak, as mentioned before, was caused by high-rigidity (>15 GV) protons, for which the NM response functions are not well determined (Lockwood 1971).

Due to the fact that neither the initial nor the peak times can be precisely determined, we use the well observed decrease after the second peak to determine the final time of the perturbation on each NM, computed as the mid-time between the maximum and the first change in curvature of the decreasing trend. These times are marked with circles in Figure 10 on the time profiles (inset). The main plot shows the time delay (the fraction of a day after 2016 October 14) of the decrease observed by each NM as a function of its Geo-longitude, and shows a clear relationship between the position of the NM and the final time of the observed interplanetary disturbance. This nonlinear relation is caused by multiple factors such as Earth's rotation, the solar wind velocity at which the MFR is sweeping, the energy of observation, and the detector response function. More detailed study of this effect is required but is out of scope of this paper.

Based on visual inspection, the magnetic configuration of the MFR displays a symmetric magnetic field profile. The structure can be described with a very well-organized single MFR with a south–north (SN) B_Z polarity and positive B_Y. Thereby, the configuration is defined to be a south–east–north configuration and left-handed (SEN-LH) (see Nieves-Chinchilla et al. 2019, and references therein for details). The MFR reconstruction is based on the circular–cylindrical model described by Nieves-Chinchilla et al. (2016). The model assumes an axially symmetric magnetic field cylinder with twisted magnetic field lines of circular cross section. The magnetic field components in the circular–cylindrical coordinate system are described by

$$\begin{eqnarray}&&{B}_{r}=0,\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{B}_{y}={B}_{y}^{0}\left[1-{\left(\displaystyle \frac{r}{R}\right)}^{2}\right],\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{B}_{\phi }=-H\displaystyle \frac{{B}_{y}^{0}}{| {C}_{10}| }\displaystyle \frac{r}{R},\end{eqnarray} \tag{ 3 }$$

where the model estimates B_y⁰, the magnetic field at the center of the MFR, and C₁₀, a measure of the force-freeness along with handedness H (+1 or −1) that indicates whether the MFR is right- or left-handed. The radius of the MFR cross section, R, is a derived parameter (see Equation (5) in Hidalgo et al. 2002) and r is the radial distance from the axis that describes the spacecraft trajectory. The reconstruction technique is based on a multiple regression technique (Levenberg–Marquardt algorithm), which infers the spacecraft trajectory and estimates the orientation of the MFR axis (azimuth and tilt) and the impact parameter (y₀).

The reconstruction of the event corroborates the visual inspection, i.e., the structure corresponds to a symmetric MFR with an axis orientation of longitude ϕ = 99° and latitude θ = −21°. The closest distance of the spacecraft to the axis was y₀ = −0.0054 au and the radius R = 0.146 au, therefore crossing very close to the center of the MFR (y₀/R ≈ 0.04). H is negative, left-handed. The magnetic field strength at the MFR axis is ${B}_{y}^{0}=20.7$ nT and the force-freeness parameter C₁₀ = 1.1 (see Nieves-Chinchilla et al. 2016, for more details).

The twisted reddish curves in Figure 11 represent the fitted MFR at four radial distances from its axis, as quoted in the left side panel. As described in Section 4.2, we used these fitted parameters to model the trajectories (inside the MFR) of the GCRs with different energies and incident angles. The colored curves starting at x = −0.14 au, y = 0 au, and z = 0 au represent these GCR trajectories.

Zoom In Zoom Out Reset image size

As shown in the previous section, the 2016 October ICME had a well-defined magnetic rope topology, ideal for a charged particle to be guided along the MFR axis (governed by the Lorentz force). The B-field inside the MFR reached a maximum of ∼24.64 nT, with a mean value (〈B〉) of ∼18.56 nT. A GCR entering the MFR experiences a Lorentz force that is maximum while it travels perpendicular to the axis and minimum when its direction is parallel to the axis. The ratio of Larmor radius (R_L) to the MFR size (S_MFR) is shown in Figure 12. The low value of this ratio indicates the capacity of the field to redirect the particles and therefore increase the likelihood of an alignment of the GCRs along the B-field (note that this is the field strength along the axis of the MFR).

Zoom In Zoom Out Reset image size

Furthermore, it should be noted that the B-field associated with this MFR had a low turbulence level, which was estimated as σ_turb ≤ 2%, using a running window method (as explained by Arunbabu et al. 2015). The turbulence level enhances the diffusion of particles, but in this event the turbulence level was low. To quantify this effect, we estimated the ratio of diffusion length (L_diff) of the particles to S_MFR using the perpendicular diffusion coefficient described in Snodin et al. (2016) with a turbulence level of σ_turb = 2%. The lower value of this ratio of L_diff/S_MFR as seen in Figure 12 supports our approximation and illustrates the viability of the model.

To model the trajectory of GCRs inside the MFR we define a coordinate system with the origin at the center of the MFR and the Y-axis aligned with the MFR axis. This can be converted from the GSE coordinate system by the rotation operators R_Z and R_Y by the angles θ and ϕ (see Section 3.3), respectively. We assumed that the MFR has a cylindrical cross section, where the boundary of the MFR lies in the XZ-plane with $\sqrt{{x}^{2}+{z}^{2}}\leqslant 0.5{S}_{\mathrm{MFR}}$. The field inside the MFR was modeled using the observations at 1 au and then converted to the coordinate system based on the MFR geometry.

The GCRs entering the MFR are redirected by the Lorentz force ($\tfrac{q}{{m}_{p}}{\boldsymbol{V}}\times {\boldsymbol{B}}$). The net deflection is then transformed (relativistically) into the particle frame. When the particle is inside the MFR, its changes in position and velocity are estimated every 100 m of travel. Figure 11 shows a few examples of the simulated particle trajectories inside the fitted MFR.

To estimate the relevant energies of particles that can be guided along the axis of the MFR we performed a simulation, for which we injected protons with energies ranging from 8 to 100 GeV entering the MFR at an initial position x =−0.5S_MFR, y = 0, and z = 0. These particles enter the MFR with an incident angle in the XY-plane Ψ, ranging from −70° to 70° with an increment of 1°. We tracked the position and direction of each particle inside the MFR by computing the point where it becomes aligned with the axial direction of the MFR (inside the cone with 85° < Ψ < 95°) and the distance it travels parallel to the axis. We found that that a rather narrow range of rigidity behaves in this manner; specifically, particles of energies higher than ∼60 GeV are less likely to be aligned along the axis (see Appendix B for details). Therefore, the enhancement observed by HAWC is due to particles of energies ≲60 GeV in this model.

As an example, in Figure 11 we show the trajectories inside the MFR of protons with energies between 8 and 50 GeV (marked with color shades) entering the MFR with five different incident angles within the range −60° < Ψ < 60° and an increment of 30° (marked by the different colors). The entrance position of this example was x = − 0.5S_MFR, y = 0, and z = 0. The incident angle Ψ = 0 means that the particle is entering the MFR perpendicular to the Y-axis and the initial velocity of the particle was resolved into perpendicular and parallel components with respect to the MFR axis, using this incident angle as ${V}_{\perp }=V\,\cos {\rm{\Psi }}$ and ${V}_{\parallel }=V\,\sin {\rm{\Psi }}$. As expected, the low-energy particles are more affected by the magnetic topology of the MFR. Particles with median energy and larger Ψ are more likely to get aligned parallel to the MFR axis.

For completeness and to cover the maximum possible entries of particles into the MFR, we introduced another angle of incidence Ξ surrounding the MFR in the XZ-plane, varying from 1° to 360° in increments of 1°, and the initial position of entry chosen as $x=-(0.5{S}_{{\rm{MFR}}})\cos {\rm{\Xi }}$, y = 0, and $z\,=-(0.5{S}_{{\rm{MFR}}})\sin {\rm{\Xi }}$. If the particles get aligned and travel in the axial direction for more than 60 R_E (with R_E = 6357 km, the Earth's radius at the equator), the initial position of alignment and the distance traveled parallel to the MFR axis are estimated. The projections of the trajectory of these particles (in the XZ-plane, i.e., in the MFR cross section), with initial energies of 8, 10, 12, 14, 20, and 30 GeV, are shown in Figure 13. The distance traveled parallel to the MFR axis is indicated by the color scale. In this figure, the trajectory of the Earth through the MFR is marked by the brown lines whereas the red and green circles mark the times when the double peak structure was observed by HAWC. Therefore, the number of points inside these circles shows the trajectory of the particles that are likely to be heading toward the Earth at these times. As seen in this figure, our model suggests that the first enhancement detected by HAWC was mainly caused by protons in the energy range 12–30 GeV, whereas the second enhancement was caused by low-energy protons in the range 8–12 GeV.

Zoom In Zoom Out Reset image size

To reach the lower atmosphere and be detected by HAWC, the GCR anisotropy caused by the MFR (parallel to its axis) must be matched by the HAWC asymptotic directions. Hence, we have estimated the asymptotic direction of these particles using a back-trace method (Smart & Shea 2005) based on IGRF12. The computed directions are shown with colored triangles in Figure 9, where the color code corresponds to the energy. As expected, the low-energy particles are subjected to large deviations inside the geomagnetic field and therefore their asymptotic directions are almost opposite to that of the particles with median energy.

The alignment angle (Λ) between the normal vector of the asymptotic direction ($\hat{N}$) and the interplanetary B-field can be represented as

$$\begin{eqnarray}&&\cos ({\rm{\Lambda }})\,=\,\displaystyle \frac{\hat{N}\cdot {\boldsymbol{B}}}{| B| }.\end{eqnarray} \tag{ 4 }$$

The cosine of Λ reaches the value 1 when B is parallel to the HAWC asymptotic direction, allowing the detection of particles of the specific energy (marked with colors in Figure 14). On the other hand, when Λ < 1 the particles are not able to reach the detector. It is important to note that during the local minimum at the middle of the double peak structure, $\cos ({\rm{\Lambda }})\leqslant 0$ and the parallel distance traveled was  ≤ 2000R_E (Figure 14) for all particle energies. At this point, the MFR axis was perpendicular to the observing directions of HAWC. Considering this alignment, the first peak of the double peak structure was dominated by particles with energy in the range ∼14–30 GeV, whereas the second peak was dominated by lower-energy particles. From this analysis (and as seen in Figure 14) we conclude that the double peak structure observed by HAWC is a geometric effect due to the alignment between the MFR axis and the HAWC asymptotic direction.

Zoom In Zoom Out Reset image size