Exploring the Coronal Magnetic Field with Galactic Cosmic Rays: The Sun Shadow Observed by HAWC

Galactic cosmic rays ( GCRs ) are charged particles that reach the heliosphere almost isotropically in a wide energy range. In the inner heliosphere, the GCR ﬂ ux is modulated by solar activity so that only energetic GCRs reach the lower layers of the solar atmosphere. In this work, we propose that high-energy GCRs can be used to explore the solar magnetic ﬁ elds at low coronal altitudes. We used GCR data collected by the High-Altitude Water Cherenkov observatory to construct maps of GCR ﬂ ux coming from the Sun ’ s sky direction and studied the observed GCR de ﬁ cit, known as Sun shadow ( SS ) , over a 6 yr period ( 2016 – 2021 ) with a time cadence of 27.3 days. We con ﬁ rm that the SS is correlated with sunspot number, but we focus on the relationship between the photospheric solar magnetic ﬁ eld measured at different heliolatitudes and the relative GCR de ﬁ cit at different

12.4 33.4 TeV, which implies that this is the best energy range to study the evolution of magnetic fields in the low solar atmosphere.

Introduction
The solar and heliospheric magnetic field is both dynamic and qualitatively variable as a function of time and heliodistance.At 1 au, there is a multitude of measurements of the field and its fluctuations over many decades.Closer to the Sun, it is less well known, but measurements have been conducted over the last few years with Parker Solar Probe (PSP) and Solar Orbiter (SO).We also know the magnetic field in great detail in the photosphere, through the long-used Zeeman splitting technique, but there is a dearth of data above the photosphere and inside the orbits of the PSP and SO.This is the region where large coronal currents can exist and the solar wind is not the driver of the form and magnitude of the magnetic field as it is the least well known.Our analysis serves to explore the magnetic field in this region by the way it modulates nearby passing Galactic cosmic rays (GCRs) over a wide range of energy on a solar rotation basis.
GCRs consist of predominantly light nuclei that get accelerated to relativistic speeds by violent processes in the galaxy.They are continuously changing their direction during their passage through the ambient interstellar magnetic field.Hence, inside the heliosphere, their flux appears to be isotropic and decreases with the energy following a power-law spectrum with an index of ∼ −2.7 (Antoni et al. 2004).Inside the heliosphere, the low energy particles (E < 100 GeV) are strongly modulated by the solar wind and the magnetic field it carries (Strauss & Potgieter 2014;Martucci et al. 2018), whereas particles of higher energy reach the inner heliosphere without significant deflections.Then, the isotropic and quasiconstant high energy GCR flux may interact with large massive and/or magnetic celestial bodies, and an observer will detect a GCR flux deficit in the direction of such a celestial body.Specifically, in our case, the Moon and Sun cause significant GCR flux deficits, also known as the Moon and Sun shadows (Clark 1957;Amenomori et al. 2013).
The Moon shadow has been used to calibrate the pointing accuracy of several air shower detectors (e.g., Amenomori et al. 1993a;Fiorino et al. 2013, and references therein).On the other hand, different air shower arrays have observed the Sun shadow (SS), e.g., Tibet (Amenomori et al. 2006) (Nan et al. 2022).In general, these studies compare the SS evolution with the solar cycle represented by the sunspot number (Amenomori et al. 1999), the interplanetary magnetic field (Amenomori et al. 1993b;Zhu et al. 2015), and the extrapolated coronal magnetic field (Amenomori et al. 1996;Aartsen et al. 2021).
A common problem of these studies is the large timescale (the integration time used to compute the SS maps) typically of 1 yr, making the comparison with the solar magnetic activity difficult due to the variability of this field at shorter timescales.We overcome this problem with new high temporal resolution data of the HAWC array (Section 2) to construct maps of the Sun and Moon shadows with two integration times: 1 yr and ∼27.3 days (corresponding to one solar rotation known as a Carrington rotation), during the 2016-2021 period, on declining phase until shortly after the past solar minimum.
As the first step, we studied the evolution in terms of the solar cycle, generally characterized by the sunspots number (Section 3).Then, we performed a detailed study of the relationship between the SS and the median photospheric magnetic field computed at different heliolatitudes (Section 4).The discussion of the relationship between the solar magnetic field and the SS is presented in Section 5 and our conclusions in Section 6.

HAWC Array and Data
The HAWC array is a very high-energy detector for cosmic and gamma rays and it is located at an altitude of 4100 m in the state of Puebla, Mexico (N 18°.5904,W 97°.1803).The main array of HAWC consists of 300 water Cherenkov detectors of 7.3 m diameter and 4.5 m height filled with highly purified water and instrumented with four photomultipliers (PMTs).The main array covers an area of about 22,000 m 2 (for a detailed description of the HAWC array see Abeysekara et al. 2023).
With HAWC, it is possible to determine the energy and direction of the primary particle by the characteristics (e.g., number of hit PMTs and time differences between the hits) of the footprint that the corresponding air shower produces on the array.In this way, HAWC efficiently detects gamma and cosmic rays with energies between ∼1 TeV and several hundred TeV, with an angular resolution approaching 0°.2 at high energies (Abeysekara et al. 2017).We have divided the HAWC detection energy into 11 energy bins (Alfaro et al. 2017;Hampel-Arias 2017).In this work, we consider bins 2-10, corresponding to median GCR energies of 2.5-226 TeV (see Table 1 where the median CR energy of each bin is shown); the details on the data acquisition and analysis processes were reported by Abeysekara et al. (2017Abeysekara et al. ( , 2018)).
The HAWC main array has been continuously operational since 2015 (>95% uptime).In this work, we analyze the data observed in the 2016-2021 period.Although, as the detector response depends on the zenith angle of the source, we consider data when the solar zenith angle is between 0°and 25°, and therefore, our annual observation period is reduced to ∼200 days, corresponding to the March-October season of each year.Finally, to better compare our results with the solar magnetic field (as reported in magnetic synoptic maps of one Carrington rotation (CarRot)), the Sun and Moon maps used in this work were integrated over the same time period.The considered CarRots for each year are: 2016: 2175-2182; 2017: 2188-2195; 2018: 2202-2209; 2019: 2215-2222; 2020: 2228-2235; and 2021: 2242-2249, except the rotation 2246 when a data gap took place (see Appendix).
We use standard HAWC software tools (Alfaro et al. 2017) to construct the Moon and Sun shadow maps.The details of the process of computing and quantifying these maps were reported in Abeysekara et al. (2017).We note that we are using equatorial coordinates (α, δ) but, in this case, the maps are centered at the Sun/Moon position (α * = α − α Sun/Moon , δ * = δ − δ Sun/Moon ), meaning that we subtract the coordinates of the Sun/Moon from the coordinates of each event.Also, we have applied a smoothing of 1°in both α and δ to reduce fluctuations in the maps.
The relative deficit in each pixel of the map is calculated as follows: where N on is the number of events in the on-source region and 〈N off 〉 is the average number of events in the off-source or background (see Abeysekara et al. 2017 for details of this process).The significance of the deviation of the data from the isotropic expectation in each bin is calculated using the Li-Ma method (Li & Ma 1983).
As examples of the shadow maps, Figures 1 and 2 show the GCR relative intensity percentage (RI) maps of the Sun and Moon, respectively.Each map was integrated over 1 yr, from 2016 to 2021.To quantify the morphology of the shadows, we fitted a 2D symmetrical Gaussian function for both the RI and significance maps for each CarRot, and the amplitude of the adjusted Gaussian is defined by the Sun/Moon shadow deficit during the corresponding periods.
The contours on each map show the 15%, 30%, 45%, 60%, 75%, and 90% levels of the 2D Gaussian fit to the significance maps in terms of the maximum significance of the whole period; for the Sun maps, this was ∼59σ reached during the CarRot 2232.
Within the time interval analyzed (6 yr), the Moon's shadow remains stable, showing a perfect circular shape and the central position is almost constant (within ∼0°.1).On the other hand, the Sun maps show a large variability, going from circular to oblate shape when the solar activity is low and high, respectively.Also, the SS center position shows variability within ∼0°. 5. Similar changes have been observed and attributed to the solar, interplanetary, and Earth magnetic fields (Amenomori et al. 1993b;Amenomori et al. 2000).Although, until now there is not a clear relationship between these fields and the morphology of the SS.In general, the relative intensity better reflects the SS variations as the solar magnetic cycle evolves, as shown in Table 2 and the bottom panel of Figure 3.The difference between the Sun and the Moon is the existence of the solar magnetic field, which has no counterpart in the Moon.Therefore, we use the SS to explore the magnetic field of the low solar atmosphere.

Time Evolution of the Sun Shadow
As noted in the previous section, the SS shows a large variability over time (see also Aartsen et al. 2021, and references therein).In order to explore such evolution, and as an advantage over previous works where the integration times are of the order of 1 yr, we constructed SS maps integrated over individual Carrington rotations starting with CarRot 2175 (2016 March 16)  The contours mark the 15%, 30%, 45%, 60%, 75%, and 90% levels of the maximum significance (∼ 59σ) reached during the 2020 year.In all frames, the outer contour corresponds to 15%, and the subsequent internal contours correspond to an addition of 15% to each one up to 90%.
and ending with CarRot 2249 (2021 October 21).To minimize the errors caused by large zenith angles of the GCR arrival direction, in this analysis, we only consider the data taken when the Sun is within a zenith angle of 25°at the HAWC site (∼19°o f north latitude), corresponding from the March to October time window of each year, with an integration time of ∼6 hr per day.
The time evolution of the RI during the period of study is shown in Figure 3, where the time is in CarRot number and the corresponding dates are marked in the bottom and top scales, respectively.The bottom panel shows the RI of the Moon (MS RI ; blue line) and Sun (SS RI ; red dots) during the period of study.It is clear that the Moon shadow variation is bounded by a small range of RI, with a mean value of −4.85%.In contrast, the magnitude of SS RI varies in a more extended range following the solar cycle in anticorrelation with the sunspot number (SSN; top panel).The plots show a low absolute value of RI when the SSN was high in mid-2016, corresponding to the decreasing phase of the solar cycle 24.The highest absolute value of the SS RI was reached in 2019, associated with the minimum of the SSN.Then, as the solar activity increases during the beginning of the solar cycle 25, the absolute value of the SS RI decreases until the end of the analyzed period.The photospheric magnetic field shown in the middle panel is more important than the SSN, which is used as a proxy of the solar cycle.See the details in Section 4.
The correlation between the SSN and the SS is important and has been shown by Amenomori et al. (2006) and Aartsen et al. (2021).Although, the SSN represent the strong magnitude and small-scale magnetic fields at low latitudes.Therefore, in order to get a better insight of the interaction of GCRs with the whole Sun, it is important to take into account the global photospheric magnetic field.
In order to understand the differences between the Moon and Sun shadows, it is important to note that the magnetic field behavior is different in the low atmosphere (below the Alfvén surface where strong, closed magnetic fields dominate) and the  1, the contours mark the 15%, 30%, 45%, 60%, 75%, and 90% levels of the maximum significance (∼59σ).The maximum significance was reached during the 2019 year.In this case, the inner contour in the 2016-2019 period corresponds to 90%, whereas the inner contour of the 2021 period reaches only 75% of the maximum significance.Notes: Relative intensity (RI) in %, the Full Width at Half Maximum (FWHM) of the fitted 2D Gaussian in degrees, and the significance for the Sun (columns (2)-( 4)) and Moon (columns ( 5)-( 7)), respectively.
rest of the heliosphere where a weaker magnetic field is transported by the solar wind.For instance, protons with an energy of 5 TeV in a magnetic field of 5 × 10 −5 G (typical value at 1 au, where the Earth and Moon are located) have a gyroradius of ∼22 au.On the other hand, the magnetic field at the solar surface is a few Gauss in quiet regions and reaches 10 3 G in active regions, giving gyroradii of ∼0.3 R e to 10 −4 R e , respectively, for 5 TeV protons, causing the difference between the nature of the SS and Moon shadow.Therefore, the differences between the Moon and Sun shadows must be related to the varying coronal magnetic field.

Solar Magnetic Field and the SS RI
To compare the SS RI with the solar magnetic field, we use synoptic magnetograms, i.e., maps of the photospheric magnetic field observed during a solar rotation acquired by GONG instruments (operated by NISP/National Solar Observatory/AURA/NSF with contribution from NOAA),33 and compute the median magnetic field of each CarRot at different latitudes ranges: low (B LL ; from −30°to 30°) where the closed magnetic field (loops and active regions) dominates in the active region belt; high (B HL ; < − 60°and >60°) in the polar regions where much of the magnetic field is in open coronal holes; and medium (B ML ; from ±30°to ±60°) where coronal holes and active regions are largely absent and the variability is minimum as seen in the photospheric data.
The corresponding scatterplots are shown in Figure 4, where the relationship between the RI and the median magnetic field measured at high latitudes (top panel) shows a good anticorrelation when the magnetic field is lower than ∼4.5 G (the  filled circles mark this anticorrelation).These lower intensity field strengths at high latitude (B HL ) correspond to the decreasing phase of the solar cycle 24, as seen in the middle panel of Figure 3.There is no apparent relationship between the magnetic field at medium latitudes (B ML ) and the RI (middle panel) of Figure 4.The bottom panel of this figure shows the positive correlation between the active region magnetic field (B LL ) and the RI.In this case, the correlation holds for the entire range of the magnetic field (from ∼2.5 to 5 G).

SS RI Energy Dependence
Going further, we make use of the HAWC capability to determine the GCR energy in the ∼2 to ∼200 TeV range and evaluate the degree of correlation between the median photospheric magnetic field computed at low and high latitudes and the RI.In this case, we divided into nine energy bins as shown in the first and second columns of Table 3, which correspond to the number and median energy of each bin, respectively.
As examples of the correlation, Figure 5 shows the scatterplots of the RI as a function of the median magnetic field of low (B LL ; black) and high (B HL ; blue) latitudes for the 4, 6, and 8 energy bins (from top to bottom).We fitted a linear model to both sets, and the resulting slopes (ΔRI/ΔB iL , where i is L, M, or H, for each energy bin) as well as the computed t − s and the reduced χ 2 shown in Table 3, in columns (3), ( 4) and (5) for high latitudes and (6), ( 7) and (8) for low latitudes, respectively.

Discussion and Results
Solar activity produces significant disturbances in the corona, the interplanetary medium, and the Earth's magnetosphere, determining the space weather.Also, and more importantly for the present study, the magnetic fields associated with this activity affect the transport of the GCR in the heliosphere.This effect is known as solar modulation and has been studied for a long time at relatively low energies (<100 GeV).
In this work, we study the effect of solar activity on GCRs in the ∼2 to ∼200 TeV energy range.In particular, we propose using the GCR observed by the HAWC array to explore the low coronal magnetic fields by computing the deficit of the GCR flux (SS RI ) due to the long-term changes of these fields.Note.From left to right columns: the bin number and the corresponding median CR energy, the rate (ΔRI/ΔBiL) with the associated t-statistics and reduced χ2 for high latitudes (HL) and low latitudes (LL), respectively.a In this case, we are using the t-statistics to assess whether the fitted rate is statistically significantly different from zero and basically is the estimated parameter divided by the standard errors.This new technique may be used to evaluate the magnitude of the low coronal magnetic field (although the morphology of the magnetic field is well determined through images at various wavelengths, the magnitude cannot be determined directly).In this way, the relationship between the SS RI and the median magnetic field represents an excellent tool to explore the solar low coronal magnetic field.
The solar cycle 24 started in 2009, reached its maximum in 2014, and ended in 2019.Therefore, part of the time window of our data coincides with the declining phase of this solar cycle.As shown in Figure 3, the SS RI follows the solar cycle (represented by the SSN) with values between −1% and −2% during the decreasing phase of the solar cycle 24 (2016), then both the SSN and the SS RI show a slow decrease until the latter reached a minimum value of −6% during the minimum of the SSN, which took place at the end of 2019.Later, the ascending phase of solar cycle 25 is accompanied by a decrease in the SS RI amplitude reaching a value between −3% and −4% in 2021.
However, the SSN only reflects the phase of the solar cycle; therefore, in order to establish a physical connection between the SS RI and the solar activity, during each CarRot considered in this work we computed the median photospheric magnetic field in three latitudinal zones: the first zone covers the activity belt at low latitudes (from −30°to 30°), where the "closed" magnetic fields associated to active regions is located (B LL ; dark blue curve in Figure 3); the second zone covers midlatitudes (from ±30°to ±60°), where the background magnetic field dominates (B ML ; triangle symbols in Figure 3); and the third zone (from ±60°to ±90°) reflects the bipolar, "open" magnetic field in the polar regions (B HL ; gold curve).The plus symbols correspond to the median magnetic field in all latitudes (B AL ).
The scatterplots of the SS RI as a function of the median magnetic field at low, medium, and high latitudes are shown in the bottom, middle, and top panels of Figure 4, respectively (the filled circles represent the correlations during the descending phase of solar cycle 24).There is an appreciable correlation between SS RI and B LL , as seen in the bottom panel of this figure.This correlation shows that higher absolute values of SS RI are associated with low B LL , reached during the minimum phase of the solar cycle, where there were no active regions implying the absence of strong toroidal magnetic fields.
Both B AL and B ML show similar behavior with small variations during the cycle as shown by the plus and triangle symbols in the central panel of Figure 3. Therefore, there is no correlation between these median magnetic fields and the SS RI , as seen in the middle panel of Figure 4, where we show the scatterplot between SS RI and B ML .
The scatterplot of the SS RI as a function of B HL shows a very strong correlation during the decreasing phase of solar cycle 24 (marked with filled circles in Figure 4).A partial conclusion of the previous analysis is that the SS RI changes according to the evolution of B HL to a high degree and with B LL to a lesser degree.Then, we use the ability of HAWC to determine the energy of the cosmic particles and explore the SS-magnetic field relationship in terms of the GCR energy.
We quantify the relationship between the SS in each energy bin observed at different energies and the median coronal magnetic field by fitting a first-order polynomial function to the SS RI i as a function of the magnetic field, where i = 2, 3, ..., 10 is the energy bin number.
As examples, in Figure 5, we show the SS RI i versus B HL (blue) and B LL (black) fittings for i = 4, 6, and 8 (from top to bottom).In general, the fitting process gives good results for all energy bins, as shown by the t-statistic and reduced χ 2 in Table 3. Although, the ΔRI/ΔB rate is significantly higher for the 5, 6, and 7 energy bins.
In Figure 6, we plotted the rate ΔRI/ΔB slope computed for each energy bin for low-and high-latitude median magnetic fields (black and blue symbols, respectively).Remarkably, the rates of change remain above −2 and below ∼1.6%/G for high and low latitudes, respectively, except for the energy bins 5, 6, and 7.The strongest anticorrelation is reached when we compute the rate of change of the RI as a function of the highlatitude median magnetic field for the energy bin 7.In this case, ΔRI/ΔB = − 4.62%/G, with a Pearson correlation coefficient of 0.95.
The large-scale solar magnetic field changes from a multipolar (disordered) topology during the maximum-activity phase into a bipolar (well-organized) configuration during the minimum-activity phase.The GCRs will face a strong and random magnetic field at low latitudes during the maximum phase and, therefore, suffer large deviations.As the cycle evolves toward the minimum phase, the magnetic field at the active region belt diminishes, allowing the GCRs to pass closer to the solar surface without major deviations.Then, an observer on the Earth will detect a correlation between the (negative) SS RI and the low-latitude magnetic field.
At the same time, during the decreasing phase of the solar cycle, the magnitude of the high-latitude magnetic field increases, and the overall topology approaches a bipolar one.In this way, GCRs will face a weaker and well-ordered magnetic field at all latitudes, except in the polar regions, where the field is stronger but well-ordered.Thus, the observed anticorrelation is because the polar magnetic field magnitude reflects the homogeneity of the overall topology of the solar magnetic field.
A deeper study, taking into account the strength and morphology of the magnetic field above the photosphere, is necessary to understand better the relationship between the solar magnetic activity and the GCRs with energy in the range 10-50 GeV, according to the results obtained with this analysis.

Figure 1 .
Figure1.Relative Intensity maps (in percentage as marked by the color scale) of the GCR flux arriving from the Sun direction and integrated over 1 yr duration periods for the years 2016-2021 as marked in each panel.The contours mark the 15%, 30%, 45%, 60%, 75%, and 90% levels of the maximum significance (∼ 59σ) reached during the 2020 year.In all frames, the outer contour corresponds to 15%, and the subsequent internal contours correspond to an addition of 15% to each one up to 90%.

Figure 2 .
Figure 2. Relative Intensity maps of the GCR flux arriving from the Moon direction.Similar to Figure1, the contours mark the 15%, 30%, 45%, 60%, 75%, and 90% levels of the maximum significance (∼59σ).The maximum significance was reached during the 2019 year.In this case, the inner contour in the 2016-2019 period corresponds to 90%, whereas the inner contour of the 2021 period reaches only 75% of the maximum significance.

Figure 3 .
Figure 3. Evolution of the SS RI (red, bottom panel) as a function of the CarRot during the 2016-2021 period (the corresponding dates are shown at the top horizontal axis).As a reference, the blue line and the gray shadow represent the changes in the MS RI (the Moon shadow).The middle panel shows the evolution of the photospheric median magnetic field during the same period; the plus and triangle symbols represent the magnetic field at all (B AL ) and medium (B ML ) latitudes, respectively, whereas dark blue and yellow lines represent the field at low (B LL ) and high (B HL ) latitudes, respectively.The top panel shows the monthly sunspot number (SSN; https://www.sidc.be/silso/datafiles)during the analysis period.

Figure 4 .
Figure 4.All circles represent the evolution of the SS RI at all energies (from ∼2 to ∼200 TeV) as a function of the median photospheric magnetic field (B) computed at low (B LL ; bottom), medium (B ML ; middle), and high (B HL ; top panel) latitudes for each CarRot considered in this study.The filled circles show the same relation but only during the decreasing phase of the solar cycle 24, from CarRot 2175 to 2205 (declining phase of the solar cycle), when the correlation between the RI and the high-and low-latitude magnetic fields is clear.

Figure 5 .
Figure 5. SS RI as a function of the median photospheric field at low (black circles) and high (blue squares) latitudes for three energy bins (4, 6, and 8 from top to bottom).The lines represent the linear models fitted to both data sets, and the shadow regions represent the 95% confidence intervals of the fitted parameters.The Pearson correlation coefficient (cc) of each linear relationship is shown in the bottom part of the panels; blue for high latitudes and black for low latitudes.

Figure 6 .
Figure6.The rate ΔRI/ΔB given by the linear fits of the RI as a function of the median magnetic field at low (black circles) and high (blue squares) latitudes for each energy bin considered in this study.

Figure 7 .
Figure 7. RI maps integrated during one Carrington rotation, for the years 2016-2021, as marked in each panel.

Table 1
The Bin Number and the Corresponding Median CR Energy

Table 2
Yearly Characteristics of the Sun/Moon Shadow

Table 3
Evaluation of the Degree of Correlation between the Median Photospheric Magnetic Field Computed at Low and High Latitudes and the RI