Influence of Earth-directed Coronal Mass Ejections on the Sun's Shadow Observed by the Tibet-III Air Shower Array

Coronal mass ejections (CMEs) are large magnetized clouds/blobs of plasma ejected by solar eruptions spreading into the interplanetary space. A fast CME forms/drives a strong shock and a magnetic sheath containing the strong and turbulent magnetic field behind the shock. Some CMEs also have rope-like magnetic structures behind the magnetic sheath (Burlaga et al. 1981; Klein & Burlaga 1982). A CME propagating away from the Sun affects (the otherwise steady) galactic cosmic rays (GCRs) in a number of ways. The best-known consequence is the so-called "Forbush decrease," which is observed when the detector enters a region of suppressed GCR density behind the shock driven by a CME (Cane 2000). It has been known that the CMEs sometimes may trigger major geomagnetic storms when they reach the Earth (Gosling et al. 1990; Gosling 1993).

While in-situ measurements on board a satellite give the 1D distributions of the magnetic field in the CME along the satellite's path, it is still difficult to precisely derive the 3D large-scale magnetic structure of a CME. TeV and sub-TeV GCRs with large Larmor radii can sense the remote CMEs and give useful additional information on their structure. Observations of these high-energy GCRs are also relevant to space weather, because these particles travel nearly at the speed of light, sense the CME, and escape into the upstream region of the shock. They can easily overtake the much slower CME, bringing an advance warning of the Earth-directed CME (ECME) approaching the Earth (Munakata et al. 2000; Rockenbach et al. 2014). The solar modulation of sub-TeV GCRs has been already reported (Munakata et al. 2010), but the influence of ECMEs has not yet been addressed above 1 TeV.

The Sun blocks high-energy GCRs arriving from the direction behind the Sun and casts a shadow (Sun's shadow) in the GCR intensity, which is influenced by the magnetic field of the Sun (Clark 1957). The small air shower (AS) array (Tibet-I) consisting of 45 scintillation detectors with a 15 m spacing started operation at Yangbajing (4300 m above sea level) in Tibet, China in 1990 and successfully observed the first evidence of the magnetic field's influence on the Sun's shadow in the 10 TeV energy range (Amenomori et al. 1993, 1996). Afterwards, the large AS array (Tibet-II) consisting of 221 scintillation detectors with the same 15 m spacing started operation in 1996 and successfully observed a clear solar-cycle variation of the Sun's shadow during a full activity-cycle from 1996 to 2009 (Amenomori et al. 2013). To clarify the physical implications of the observed solar-cycle variation, detailed numerical simulations of the Sun's shadow have been developed based on the potential field models of the coronal magnetic field (Wiegelmann et al. 2015). It was found that the simulated intensity deficit in the Sun's shadow is very sensitive to the coronal magnetic field models, and the observed solar-cycle variation is best reproduced by the current sheet source surface (CSSS) model. This was the first successful attempt of evaluating the coronal magnetic field models by using the Sun's shadow observed in the TeV GCR flux (Amenomori et al. 2013). A higher density array with a 7.5 m spacing then started operation in late 1999 (Amenomori et al. 2003). This array (Tibet-III) can observe lower energy GCRs by lowering the energy threshold below 3 TeV (Amenomori et al. 2018). In this paper, we analyze the Sun's shadow observed by Tibet-III, and report the influence of ECMEs on the Sun's shadow observed in the 3 TeV GCR intensity.

We analyze AS events observed from the direction of the Sun by Tibet-III in 2000–2009. The overall angular resolution and the modal energy of primary GCRs recorded with the array are estimated to be 09 and 3 TeV, respectively. We define the angular resolution as an angular radius containing 50% of AS events produced by primary GCRs arriving from an identical direction. We calculate the number of on-source events (N_on) as the number of events arriving from a given direction within a circular window centered at a certain point on the celestial sphere. The number of background or off-source events ($\langle {N}_{\mathrm{off}}\rangle$), on the other hand, is calculated as the average number of events within each of the eight off-source windows which are located at the same zenith angle as the on-source window, but apart by ±96, ±128, ±160 and ±192 in the azimuthal direction (Amenomori et al. 2009). The window radius of 09 (angular resolution for 3 TeV GCRs) is adopted for calculating both N_on and $\langle {N}_{\mathrm{off}}\rangle$. We calculate N_on and $\langle {N}_{\mathrm{off}}\rangle$ on each gridpoint on 01 × 01 mesh of the geocentric solar ecliptic (GSE) longitude and latitude surrounding the optical center of the Sun. We then estimate the deficit relative to the number of background events as ${D}_{\mathrm{obs}}=({N}_{\mathrm{on}}-\langle {N}_{\mathrm{off}}\rangle )/\langle {N}_{\mathrm{off}}\rangle$ at every grid. We also apply the same method to calculate D_obs in the Moon's shadow in the equatorial coordinate system.

Figures 1(a) and (b) show yearly maps of D_obs in the Sun's and the Moon's shadows, respectively, observed at 3 TeV in 2000–2009, except for 2006, which is excluded from Figure 1(a) because of the insufficient statistics. The Moon's shadow provides us with a good reference of the detector stabilities, because the Moon has no magnetic field and its optical size viewed from the Earth is almost the same as the Sun (Amenomori et al. 2009). It is seen in Figure 1(a) that the Sun's shadow is darker (with larger negative D_obs) around 2008 when the solar activity was close to the minimum, while it becomes faint (with smaller negative D_obs) around 2001 when the activity was close to the maximum. In contrast, the observed Moon's shadow shown in Figure 1(b) is quite stable, ensuring the instrumental stability during the same period.

Zoom In Zoom Out Reset image size

For the primary GCRs in our Monte Carlo (MC) simulations of the Sun's shadow, we assume the energy spectra and the elemental compositions modeled by Shibata et al. (2010), which compile various measurements between 0.3 TeV and 1000 TeV. We set the minimum energy of primary particles at 0.3 TeV, well below the threshold energy for triggering our AS detection. We perform the simulations as follows: (1) we first generate AS events randomly at the top of the atmosphere along the Sun's orbit using the CORSIKA code (Heck et al. 1998) with the EPOS-LHC interaction model (Pierog et al. 2015) so that AS cores uniformly distribute within a 300 m radius from the center of array. This area within 300 m from the center sufficiently covers detectors actually hit by AS particles; (2) we then distribute these simulated events among detectors configuring Tibet-III by the GEANT4 code (Agostinelli et al. 2003) and calculate the output of each detector, which can be analyzed for AS reconstruction and event selections in the same way as the experimental data; (3) we assign an opposite charge to each primary particle at the top of the atmosphere and shoot it in random directions within 4° from the Sun, each called as the initial shooting direction; (4) we trace the orbital motion of each particle back to the Sun in the model magnetic fields (see the next section) by the fourth-order Runge–Kutta method; (5) if a particle hits the photosphere, its initial shooting direction is regarded as the direction contributing to the intensity deficit in the Sun's shadow. After smearing initial shooting directions mimic the angular resolution event by event, we finally get the expected Sun's shadow equivalent to the observation (Amenomori et al. 2013).

For the solar magnetic field model in the MC simulation, we adopt the CSSS model which is the potential field model most successfully reproducing the temporal variation of the Sun's shadow observed with Tibet-II at 10 TeV (Amenomori et al. 2013). The potential field models describe the coronal magnetic field based on the optical measurements of the photospheric magnetic field. In our MC simulations, we use the photospheric field observed with the spectromagnetograph of the National Solar Observatory at Kitt Peak (Jones et al. 1992) in each Carrington rotation (CR) period (∼27.3 days). The CSSS model (Zhao & Hoeksema 1995) involves four free parameters, the radius R_ss of the spherical source surface (SS) where the supersonic solar wind starts blowing radially, the order n of the spherical harmonic series describing the observed photospheric field, the radius R_cp (=1.7 R_⊙) of the sphere where the magnetic cusp structure in the helmet streamers appears, and the length scale l_a of horizontal coronal electric currents. In the present paper, we set l_a to be one solar radius (l_a = R_⊙) and examine two different cases with R_ss = 2.5 R_⊙ and R_ss = 10 R_⊙. The former R_ss is a standard value used in the original paper (Zhao & Hoeksema 1995), while the latter gained support from some recent evidences (Balogh et al. 1995; Zhao et al. 2002; Schüssler & Baumannk 2006). We set n = 10 which is sufficient to describe fine structures relevant to the orbital motion of high-energy particles with large Larmor radii. The radial component of the coronal magnetic field at R_ss is then stretched out forming the Parker's spiral interplanetary magnetic field (IMF; Parker 1958). For the radial solar wind speed needed in the Parker's model, we use the solar wind speed synoptic chart estimated from the interplanetary scintillation measurement in each CR and averaged over the Carrington longitude (Tokumaru et al. 2010).²⁷ We adopt a dipole model for the geomagnetic field.

In the present paper, we analyze the temporal variation of D_obs, measured at the center of the yearly mean 2D map in Figure 1(a). For reference, the solid curve in Figure 2(a) displays the monthly mean sunspot number representing the solar activity on the right vertical axis, while the solid squares in Figure 2(b) show the temporal variation of D_obs observed by Tibet-II at 10 TeV. It is seen in Figure 2(b) that the magnitude, $| {D}_{\mathrm{obs}}|$, shows a clear solar-cycle variation decreasing with increasing solar activity. It is indicated by our MC simulations that, during the solar maximum, GCRs passing near the solar limb are scattered by the strong and complicated coronal magnetic field and may appear from the direction of the optical solar disc reducing $| {D}_{\mathrm{obs}}|$. Also shown in Figure 2(b) by the red open circles and green open triangles are predictions by two different MC simulations using the CSSS models with R_ss = 2.5 R_⊙ and R_ss = 10.0 R_⊙, respectively. It is seen that these models successfully reproduce the temporal variation over an entire period (Amenomori et al. 2013).

Zoom In Zoom Out Reset image size

Solid squares in Figure 2(c) show the temporal variation of D_obs observed at 3 TeV by Tibet-III. The solar-cycle variation is again clearly seen, but with a larger amplitude when compared with Figure 2(b) at 10 TeV, due to the larger magnetic deflection expected for lower energy GCRs in the solar corona. This energy dependent feature of the solar-cycle variation is overall reproduced by the MC simulations using the CSSS models. We also note, however, that $| {D}_{\mathrm{obs}}|$ in Figure 2(c) is significantly less than those predicted by MC simulations particularly in 2000–2002 around the solar maximum with the statistical significance 4.6σ and 5.4σ for the CSSS models with R_ss = 2.5 R_⊙ and R_ss = 10.0 R_⊙, respectively. Such deviations from MC simulations are not seen in Figure 2(b) for 10 TeV GCRs.

The χ² test of the agreement between the observed and simulated D_obss in Figure 2 are summarized in Table 1. The solid squares in Figure 2(f) show the temporal variation of D_obs in the Moon's shadow observed by Tibet-III array, while the dashed curve indicates the variation of D_obs expected from the small variation of the distance between the Earth and the Moon. We estimate the systematic error of D_obs in Table 1 to be 0.178% from the deviation of D_obs from the dashed curve in this figure. It is clear in Table 1 that the deviation of D_obs from the MC simulation during an entire period in Figure 2(c) is exceeding the acceptable error range.

In order to examine the influence of the ECMEs that occur most frequently during the solar maximum period, we use the catalog of CMEs compiled by Richardson & Cane (2010) that includes 228 ECMEs recorded in 2000–2009, each with a given arrival time at the Earth. While an ECME arrives at the Earth 2–4 days after the solar eruption, 3 TeV GCRs take only ∼8 minutes to reach the Earth after passing near the Sun. GCRs arriving at the Earth from the direction of the Sun, therefore, can be affected by the ECME during a transit period between the solar eruption time observed by the Solar and Heliospheric Observatory (SOHO)/Large Aperture Solar Coronagraph (LASCO) and the arrival time at the Earth. In this study, we examine the influence of ECMEs by analyzing D_obs observed during a period with/without these transit periods of ECMEs. For ECMEs lacking relevant eruption time available from the SOHO/LASCO observation, we assume 4 days for the transit period, which is an average of all ECMEs covered by the SOHO/LASCO observation.

The histograms in Figure 2(a) show the number of ECMEs recorded in each year. The gray histogram displays all 228 ECMEs listed in Richardson & Cane's catalog, while the blue histogram shows only 118 ECMEs each with a transit period covered by the observation of D_obs by Tibet-III. It is seen that the number of ECMEs varies roughly in a positive correlation with the sunspot number. Due to this correlation, almost half of the analysis period is occupied by 20–35 ECMEs per year in 2000–2002 around the solar maximum, while more than 90% of the analysis period contains no ECME in 2007–2009 around the solar minimum. The solid squares in Figure 2(d) display D_obs observed during periods without ECMEs transit periods, while the solid squares in Figure 2(e) show D_obs observed during the ECME transit periods that are excluded in Figure 2(d). It is seen in these figures that the deviation of D_obs from the MC simulations in 2000–2002 is significantly reduced in Figure 2(d) than in Figure 2(c), while it is increased in Figure 2(e), as confirmed by the χ² tests in Table 1. Interestingly, the D_obs observed at 10 TeV by the Tibet-II array during the ECME transit period also seem to be marginally inconsistent with the MC simulation at the chance probability of 0.044 as seen in this table.

Due to insufficient statistics and limited angular resolution, the current observation of the Sun's shadow with Tibet-III does not allow us to analyze each individual ECME event by event. Hence, we superposed the temporal variations of D_obs on a daily basis with the CME eruption time set at t = 0 in Figure 3. In order to reduce the influence of the overlapping ECMEs, we excluded an ECME from this plot when it has another ECME eruption(s) within ±5 days from its eruption time. The blue solid circles show D_obs observed in the entire period including 2000–2002 when the transit periods were recorded with significant overlapping in the table by Richardson & Cane (2010), while the red triangles show D_obs observed in 2003–2009 when the transit periods were less overlapping. Although the variations in this plot may not be statistically compelling, one can still see the tendency that D_obs has a maximum at t = 0, and gradually decreases during the following couple of days. This may imply that the ECME continues to affect the Sun's shadow even after getting out to interplanetary space. This feature is better seen in the 2003–2009 period (red triangles) that is less contaminated by overlapping CMEs. We also note that D_obs starts to gradually increase before the ECME eruption (t < 0). This again may suggest that the Sun's shadow can be a useful tool for space weather forecasting. It is also noted in this figure that the temporal variation of D_obs due to the ECME is separate from the solar-activity-cycle variation that appears in the difference between the background (average) levels in 2000–2009 and 2003–2009.

Zoom In Zoom Out Reset image size

These results give a supporting evidence of the influence of the ECMEs on D_obs observed at 3 TeV. The magnetic deflection of GCR orbit in the CME is proportional to L/R_L, where L is the scale length of the magnetic field and R_L is the Larmor radius of GCRs. Due to the self-similar expansion of the CME, L increases with the radial distance r from the Sun as L ∝ r, while R_L also increases as R_L ∝ r² because of the transverse field component, B_t, responsible to the orbital deflection, decreasing as B_t ∝ 1/r² in the expanding magnetic flux rope. So, the deflection power of the CME varies as 1/r. The same deflection is, however, more effective in reducing D_obs when it occurs farther from the Sun, resulting in another factor proportional to r. The net influence on D_obs, therefore, may remain nearly constant during the transit period of the ECME. If we assume B_t = 10 nT and L = 0.1 au at 1 au, we obtain ∼2° for the deflection angle which is larger than the optical diameter (∼05) of the Sun viewed from the Earth and enough to decrease D_obs.

As R_L of 3 TeV cosmic rays is as large as 0.1 R_⊙ even in a 1 G magnetic field, the Sun's shadow is expected to be more sensitive to the magnetic field with L comparable to (or larger than) R_L and less sensitive to the smaller scale field near the photosphere, even if the field is stronger. $| {D}_{\mathrm{obs}}|$ is, therefore, expected to decrease significantly when such a large-scale magnetic field appears over the active regions on the solar limb. In order to check how a small-scale field in the solar corona can affect D_obs, we performed the MC simulations with the CSSS model by increasing the maximum order of the harmonic series from n = 10–30. The smaller scale in coronal magnetic field is expressed by the higher order harmonic component. As a result, we found that D_obs is insensitive to the smaller scale field with n > 10. Since the nth order harmonic component of the field at r is proportional to ( R_⊙/r)⁽ⁿ⁺²⁾, a smaller and stronger field on the photosphere represented with a larger n diminishes faster with r and cannot occupy a large volume (∼L³) sufficient to deflect TeV cosmic-ray orbits. We also performed our simulations by changing the photospheric magnetic field observations at Kitt Peak (Synoptic Optical Long-term Investigations of the Sun (SOLIS)/Kitt Peak Vacuum Telescope (KPVT)) to those by the Michelson Doppler Imager (MDI)²⁸ and Global Oscillation Network Group (GONG),²⁹ which are known to have the different spatial resolution and dynamic range of the observed field strength (Riley et al. 2014). The average strength of the photospheric magnetic field by the MDI (GONG) observation is 1.5–2 times higher (lower) than that of SOLIS/KPVT, but we found no significant difference in the calculated solar-cycle variations of D_obs. This indicates that the strong magnetic fields around the active regions on the photosphere, which largely contribute to the average field strength, are mostly closed in the lower corona and incapable of occupying a large volume sufficient to affect D_obs. These results indicate that D_obs is insensitive to small-scale magnetic structures in the potential field model.

The discussion above is only qualitative and we need the detailed magnetic field models of the ECME for quantitative studies. It is noted in Figures 2(c)–(e) that the prediction of our MC simulation is almost insensitive to the ECME transit periods. This is because the potential field models used in our simulation are capable of reproducing only the average magnetic field in each CR period and incapable of properly representing the dynamic variations of the field in association with the ECMEs. Shiota et al. (2010) and Den et al. (2015) have recently developed the magnetohydrodynamic (MHD) simulations to reproduce the propagation and evolution of the ECMEs in the Sun–Earth space aiming at the space weather prediction. In the future, we will incorporate the magnetic fields by these MHD simulations into our MC simulations of the Sun's shadow. Comparisons of such simulations with the observed Sun's shadow will provide us with an alternative tool for evaluating the MHD models and for improving the space weather prediction. Due to the insufficient statistics and the limited angular resolution, the current observation of the Sun's shadow with Tibet-III does not allow us to analyze each individual ECME event by event. The future huge AS arrays will enable us to analyze the Sun's shadow during an individual ECME transit period on a daily or even on hourly basis.

We confirmed a clear solar-cycle variation of the Sun's shadow observed at 3 TeV by Tibet-III in 2000–2009. By comparing this observed variation with the MC simulations, assuming the CSSS models for the coronal magnetic field, we find that the magnitude ($| {D}_{\mathrm{obs}}|$) of GCR intensity deficit in the Sun's shadow is significantly less than those of the MC simulations, particularly during years around the solar maximum when a significant number of ECMEs are recorded at the Earth. The χ² test of the agreement between the observations and the MC simulations shows that $| {D}_{\mathrm{obs}}|$ is considerably increased by excluding the ECME transit periods from the analysis, while it is reduced in the analyses of the ECME transit periods. This is the first experimental result indicating the influence of the ECMEs on the Sun's shadow observed at 3 TeV.

The collaborative experiment of the Tibet AS arrays has been performed under the auspices of the Ministry of Science and Technology of China (No. 2016YFE0125500) and the Ministry of Foreign Affairs of Japan. This work was supported in part by a Grant-in-Aid for Scientific Research on Priority Areas from the Ministry of Education, Culture, Sports, Science, and Technology, by Grants-in-Aid for Science Research from the Japan Society for the Promotion of Science in Japan, and by the National Natural Science Foundation of China (Nos. 11533007 and 11673041) and the Chinese Academy of Sciences and the Key Laboratory of Particle Astrophysics, Institute of High Energy Physics, CAS. This work is also supported by the joint research program of the Institute for Cosmic Ray Research (ICRR), the University of Tokyo. K.K. is supported by the Toray Science Foundation. The authors thank Dr. Xuepu Zhao of Stanford University for providing the usage of the CSSS model. The authors thank Dr. József Kóta of University of Arizona for his useful comments and discussions.