Cosmic ray north-south anisotropy: rigidity spectrum and solar cycle variations observed by ground-based muon detectors
M. Kozai Y. Hayashi K. Fujii K. Munakata C. Kato N. Miyashita A. Kadokura R. Kataoka S. Miyake M. L. Duldig J. E. Humble K. Iwai



Masayoshi Kozai
kozai.masayoshi@nipr.ac.jp

0000-0002-3948-3666]M. Kozai
Polar Environment Data Science Center (PEDSC), Joint Support-Center for Data Science Research (ROIS-DS), Research Organization of Information and Systems, Tachikawa, Tokyo, Japan

0000-0002-0890-0607]Y. Hayashi
Physics Department, Shinshu University, Matsumoto, Nagano, Japan

0000-0002-4315-6369]K. Fujii
National Institute of Informatics (NII), Research Organization of Information and Systems, Chiyoda-ku, Tokyo, Japan

0000-0002-2131-4100]K. Munakata
Physics Department, Shinshu University, Matsumoto, Nagano, Japan

0000-0002-4913-8225]C. Kato
Physics Department, Shinshu University, Matsumoto, Nagano, Japan

0009-0006-3569-7380]N. Miyashita
Physics Department, Shinshu University, Matsumoto, Nagano, Japan

0000-0002-6105-9562]A. Kadokura
Polar Environment Data Science Center (PEDSC), Joint Support-Center for Data Science Research (ROIS-DS), Research Organization of Information and Systems, Tachikawa, Tokyo, Japan

0000-0001-9400-1765]R. Kataoka
National Institute of Polar Research (NIPR), Research Organization of Information and Systems, Tachikawa, Tokyo, Japan

0000-0002-3067-655X]S. Miyake
National Institute of Technology (KOSEN), Gifu College, Motosu, Gifu, Japan

0000-0001-7463-8267]M.L. Duldig
School of Natural Sciences, University of Tasmania, Hobart, Tasmania, Australia

0000-0002-4698-1671]J.E. Humble
School of Natural Sciences, University of Tasmania, Hobart, Tasmania, Australia

0000-0002-2464-5212]K. Iwai
Institute for Space-Earth Environmental Research (ISEE), Nagoya University, Nagoya, Aichi, Japan

§ ABSTRACT 

The north-south (NS) anisotropy of galactic cosmic rays (GCRs) is dominated by a diamagnetic drift flow of GCRs in the interplanetary magnetic field (IMF), allowing us to derive key parameters of cosmic-ray propagation, such as the density gradient and diffusion coefficient.
We propose a new method to analyze the rigidity spectrum of GCR anisotropy and reveal a solar cycle variation of the NS anisotropy's spectrum using ground-based muon detectors in Nagoya, Japan, and Hobart, Australia.
The physics-based correction method for the atmospheric temperature effect on muons is used to combine the different-site detectors free from local atmospheric effects.
NS channel pairs in the multi-directional muon detectors are formed to enhance sensitivity to the NS anisotropy, and in this process, general graph matching in graph theory is introduced to survey optimized pairs.
Moreover, Bayesian estimation with the Gaussian process allows us to unfold the rigidity spectrum without supposing any analytical function for the spectral shape.
Thanks to these novel approaches, it has been discovered that the rigidity spectrum of the NS anisotropy is dynamically varying with solar activity every year.
It is attributed to a rigidity-dependent variation of the radial density gradient of GCRs based on the nature of the diamagnetic drift in the IMF.
The diffusion coefficient and mean-free-path length of GCRs as functions of the rigidity are also derived from the diffusion-convection flow balance.
This analysis expands the estimation limit of the mean-free-path length into ≤200 GV rigidity region from <10 GV region achieved by solar energetic particle observations.







§ INTRODUCTION
 
The cosmic-ray intensity from ∼0.1 to ∼100 GeV energy region shows a clear anti-correlation with the solar activity typically represented by the sunspot number, explicitly indicating that galactic cosmic rays (GCRs) are influenced by the helispheric state, or solar modulation, through which they propagate before approaching Earth.
Measurements of GCR properties, including the intensity, energy spectrum, composition, anisotropy, and their temporal variations, on the ground and in the space have been driving the elucidation of the GCR propagation and space environment.
The anisotropy, or momentum-space distribution of GCRs, is expected to be the next key to reveal solar modulation;
although it is theoretically proved to have a close relation with the interplanetary magnetic field (IMF) structure and GCR diffusion coefficients <cit.>, observed solar cycle variations of the anisotropy are nearly unexplored in quantitative reproductions by heliospheric simulations, unlike the GCR density variations.
Even in the broader areas of cosmic-ray physics, including galactic propagation, a consistent interpretation of the anisotropy with other properties has been an open question <cit.>.

Ground-based observations of cosmic rays measure secondary cosmic rays produced by GCR interactions with atmospheric nuclei and the subsequent reactions.
They feature long-term stability, improved statistics from large detection areas, and large angular acceptances by a worldwide network, all essential to study solar modulation phenomena consisting of diverse temporal variations from hourly (e.g., interplanetary shocks) to >10-year (solar cycles) scales.
Among secondary cosmic-ray species on the ground, muons have the highest flux (statistics) and relatively accurate response of the incident angle to primary cosmic rays because they directly originate from pion decays immediately occurring after the atmospheric nuclear spallations by primary cosmic rays typically around the upper edge of troposphere.
These advantages make the ground-based muon detector a unique means for measuring the anisotropy whose magnitude is only ∼0.1% of GCR intensity.

However, anisotropy observation by muon detectors has suffered from an atmospheric temperature variation that perturbs muon counting rates.
Conventional methods <cit.> have corrected the temperature effect by subtracting muon counting rates of directional channels from each other to cancel out the effect.
The subtracted directional channels must be from the identical location where the temperature effect is expected to be similar.
Additionally, in some analyses, the statistics are degraded by this subtraction.
On the other hand, <cit.> proved the validity of a more physics-based approach in which the temperature effect is directly corrected in each directional channel by using meteorological reanalysis data.
Contrary to the conventional method, this correction method allows for a direct combination of observations from different sites free from the local temperature effect.
<cit.> applied this correction to derive the rigidity spectrum of the anisotropy in a solar eruption event, demonstrating the new temperature correction method for anisotropy studies.
However, there is still a problem in that power-law indices of the rigidity spectra of the derived anisotropy have unreasonably unstable fluctuations, especially in the case of a small anisotropy amplitude.
The power-law function supposed to reconstruct the rigidity spectrum in their study is likely too definitive for the dynamically varying spectrum, causing this problem.

In this study, we propose a new method applying Bayesian estimation for deriving the rigidity spectrum.
It does not pre-suppose an analytical function such as the power-law function for the spectrum, while it uses the Gaussian process to confine the smoothness of the spectrum instead.
This high tolerance of varying spectral shapes enables us to trace the dynamic variation of the spectrum in solar modulation phenomena and provides a reliable analysis result.
Additionally, we introduce general graph matching <cit.> in graph theory to survey optimal channel combinations in the inter-station network, leveraging the advantage of the temperature correction method that allows us to combine multiple station data directly.
We demonstrate these analysis ideas by analyzing the rigidity spectrum of the north-south (NS) anisotropy observed by ground-based muon detectors in Nagoya, Japan, and Hobart, Australia.
The NS anisotropy, or NS flow of GCRs, shows a polarity reversal according to the local IMF sector assigned as Toward or Away each time <cit.>.
This phenomenon proves that the NS anisotropy is dominated by a diamagnetic drift flow of GCRs induced by a combination of the gyro-motion and radial (anti-sunward) density gradient of GCRs in the IMF.
Based on this mechanism, we derive rigidity-dependent variations of the density gradient and, subsequently, the radial diffusion coefficient of GCRs, which are essential parameters to elucidate solar modulation.
We also report a plan for the NS conjugate observation in the polar region, which will vastly improve the sensitivity of the muon detector network to the NS anisotropy.



§ METHOD
 



 §.§ Preprocessing of muon counting rates
 
Table <ref> lists characteristics of the Nagoya and Hobart muon detectors used in this paper, compiled by <cit.>.
It summarizes the geographical latitude, longitude, and altitude of the observation sites (λ_D, ϕ_D, and alt.), the number of directional channels (ch-no.) each counting incident muons, the geomagnetic cutoff rigidity P_ cut, the hourly muon counting rate (cph), its statistical error σ (0.01%), the median rigidity P_ median of primary cosmic rays producing muons measured in each directional channel, and the asymptotic trajectory direction (λ_ asymp, ϕ_ asymp) of the primary cosmic rays with the median rigidity outside the magnetosphere.
Nagoya and Hobart stations started their observations in 1970 and 1992, respectively, and currently form a part of the Global Muon Detector Network (GMDN).
We analyze their muon counting rates from 1995 when both the muon detectors and the solar wind observation data used in this paper have sufficient duty cycles.
Archive data of the muon counting rates in GMDN are published in Shinshu University's institutional repository <cit.>.


Our processing of muon counting rates described below partially refers to the conventional analysis method of the NS anisotropy, known as Nagoya-GG <cit.>.
It uses muon counting rates I_i uncorrected for the temperature effect in Nagoya station and takes differences

    GG =  (I_ N2 - I_ S2) + (I_ N2 - I_ E2) /2

where factor 2 is introduced in this paper for a quantitative comparison with our method.
The channel identifiers N2, S2, and E2 indicate secondary-inclined channels whose central viewing directions are north-, south-, and eastward, respectively, with a zenith angle of 49^∘.
The first NS differential term, I_ N2 - I_ S2 ensures a high response of Nagoya-GG to the NS anisotropy.
The second term I_ N2 - I_ E2 is introduced to cancel out the temperature effect remaining in the first term.
Nagoya-GG is proven to have sufficient sensitivity to the NS anisotropy and has contributed to revealing it <cit.>.

In our method, on the other hand, the atmospheric temperature effect on muons is corrected in each directional channel of Nagoya and Hobart detectors by a physics-based method using the meteorological reanalysis data, as described in Appendix <ref>.
The corrected muon counting rate I_i(t) for each directional channel i is coupled to the GCR anisotropy in space as <cit.>

    I_i(t) ∼∫_P=0^∞[    ξ_0(t,P) c_0,i^0(P) + ξ_z(t,P) c_1,i^0(P) +
       ξ_x(t,P) c_1,i^1(P) cosω t_ st - s_1,i^1(P) sinω t_ st
                + ξ_y(t,P) s_1,i^1(P) cosω t_ st + c_1,i^1(P) sinω t_ st] .

Here, the anisotropy is expanded into the zeroth and first-order spherical harmonics in momentum space with rigidity P.
The expansion coefficients ξ_0, ξ_z, ξ_x and ξ_y are defined in a geocentric (GEO) coordinate system whose z-axis points geographic north pole and x-axis points to midnight in the equatorial plane.
Indices n and m attached to coefficients c_n,i^m and s_n,i^m are degree and order of the spherical harmonics, respectively, representing components of the zeroth and first-order harmonics.
These coefficients are known as the differential coupling coefficients and are derived from numerical calculations of GCR propagation in the magnetosphere, atmosphere, and detector <cit.>.
The local time t_ st of each station to which the directional channel i belongs is related to universal time t as ω t_ st = ω t + ϕ_ st where ϕ_ st is the station's geographic longitude and ω=π/12.
The higher order (≥2nd) harmonics are omitted because they have negligible amplitudes compared to the zeroth and first-order anisotropy in the usual state.
This equation describes the coupling between the observed muon counting rate I_i and the space anisotropy (ξ_0, ξ_z, ξ_x and ξ_y) via the differential coupling coefficients at each rigidity P.
The zeroth order anisotropy ξ_0, or isotropic component, represents a variation of the GCR density from its Bartels rotation average because the counting rate I_i is converted into a deviation from its Bartels rotation average as shown in Appendix <ref>.
Equatorial components of the first-order anisotropy (ξ_x,ξ_y) represent a diurnal anisotropy or the GCR flow in the equatorial plane.
The remaining component ξ_z is the NS anisotropy representing the NS component of the GCR flow, which is the main target of this study.

The hourly counting rate I_i(t) is averaged for each day in which all counting rates and hourly IMF data mentioned below are available for at least 20 hours.
This daily average is expected to minimize a contribution from the diurnal anisotropy (ξ_x, ξ_y) in the counting rates.
Then the daily counting rates are sorted into days assigned to each IMF sector, Toward or Away, and averaged over each sector's period in each year, as I_i^T for the Toward sector or I_i^A for the Away sector.
In this process, only Bartels rotations each with ≥5 Toward, ≥5 Away, and total ≥15 available days are used to secure the sector reversal in each Baterls rotation.
Errors of I_i^T and I_i^A, σ_i^T and σ_i^A respectively, are derived as standard deviations of the daily counting rates.

The IMF data is obtained from OMNIWeb service <cit.> in the geocentric solar ecliptic (GSE) coordinate system whose x-axis points toward Sun and y-axis opposes Earth's orbital motion.
In the picture of the diamagnetic drift dominating the NS anisotropy <cit.>, the IMF parallel to the ecliptic plane is assumed, and its GSE-y component (B_y) is essential for determining the anisotropy.
Therefore, we omit days with |B_z|>|B_y| for the daily average IMF, and the IMF sector is identified daily as a Toward (Away) sector when B_y < 0 (B_y ≥ 0).
It is also noted that a definition of each year is slightly modified as described in Table <ref>, so that its start and end dates correspond to boundaries of the Bartels rotation periods.


Expanding the concept of Nagoya-GG into our temperature-corrected dataset, we take a difference between counting rates of directional channels i and j averaged for Toward (Away) sector in each year as

    η_ij^T(A) = I_i^T(A) - I_j^T(A).

Its standard error is derived as

    (σ_ij^T(A))^2 = (σ_i^T(A))^2 + (σ_j^T(A))^2 .

The channel pair ij is defined to have sufficient sensitivity to the NS anisotropy in Section <ref>.
The sector reversal of the NS anisotropy is extracted by taking a difference between η_ij^T and η_ij^A, and from equation (<ref>), it is expressed as

    η_ij^TA = η_ij^T - η_ij^A/2∼∫_P=0^∞( ξ_z^TAc_ij^z + ϵ_0^TAc_ij^0 + ϵ_c^TAc_ij^d + ϵ_s^TAs_ij^d )

where

    c_ij^z = c_1,i^0 - c_1,j^0,  c_ij^0 = c_0,i^0 - c_0,j^0,  c_ij^d = c_1,i^1 - c_1,j^1,ands_ij^d = s_1,i^1 - s_1,j^1.

Rigidity P as a variable of ξ and c (s) is omitted in this equation for simplicity.
A standard error of η_ij^TA is derived as

    (σ_ij^TA)^2 = (σ_ij^T)^2 + (σ_ij^A)^2/4 .

The parameter ξ_z^TA represents a sector reversal of the NS anisotropy ξ_z, defined as ξ_z^TA = (ξ_z^T - ξ_z^A)/2 where ξ_z^T(A) denotes the Toward (Away) sector average of ξ_z in each year.
Hereafter, we call ξ_z^TA(P) the rigidity spectrum of the NS anisotropy.
The parameters ϵ_0^TA and (ϵ_c^TA,ϵ_s^TA) also correspond to Toward - Away sector differences of GCR density ξ_0 and diurnal anisotropy (ξ_x,ξ_y) respectively in each year.
Contrary to the NS anisotropy, these parameters have no physical reason to systematically depend on the IMF polarity in sufficiently long period averages.
Therefore these are expected to play only as perturbations on η_ij^TA, while ξ_z^TA has a substantial contribution on η_ij^TA, leading to an approximation of equation (<ref>) as

    η_ij^TA∼∫_P=0^∞ξ_z^TAc_ij^z.

This equation indicates that the observed value η_ij^TA is a convolution of the NS anisotropy spectrum ξ_z^TA(P) through its differential response c_ij^z(P)/dP.

Unlike Nagoya-GG described above, the channel pair ij can be formed from any two of all channels in Nagoya and Hobart muon detectors, thanks to the physics-based temperature correction in each channel.
Nagoya and Hobart muon detectors have 17 and 13 directional channels, respectively;
the total 30 channels contain 30×29/2=435 kinds of channel pairs, and we need to select some pairs from them to analyze the NS anisotropy.
Simultaneously forming multiple pairs makes this problem more complex.
In the next section, we solve this problem using graph theory algorithm.



 §.§ Optimization of channel pairing
 
Our optimization problem is defined as searching for a pattern of the channel pairing (ij's) which maximizes the sum of sensitivity

    A = ∑_ij a_ij

where a_ij is a sensitivity of each channel pair ij to the NS anisotropy, derived in Appendix <ref>.
The set of the pairs has to meet a condition that each channel (i or j) does not duplicate in the set to ensure that derived η_ij^TA's are independent of each other.
This problem is equivalent to the maximum matching problem in general graph, or general graph matching in graph theory, and solved by the  function of the Python  library <cit.>.

Figure <ref> visualizes this problem and its solution.
Each node in this network corresponds to a directional channel (i or j) of each pair passing the threshold described in equation (<ref>), and each edge connecting the nodes represents each channel pair.
The node positions are determined by the Fruchterman-Reingold algorith <cit.>.
The line width of each edge is proportional to the sensitivity a_ij of each pair, and the edges with red colors represent the optimized pairing pattern which maximizes the total sensitivity A.

Each node's code “NAG” or “HOB” represents Nagoya or Hobart station.
Each character string below the station code indicates each directional channel, where V, N, S, E, and W indicate the vertical channel and north-, south-, east-, and westward inclined channels, respectively.
The 1st inclined channels (N, S, E, W) have central viewing direction with a zenith angle of 30^∘, while the 2nd and 3rd inclined channels denoted by 2 and 3 in the channel name have the zenith angles of 49^∘ and 60^∘ respectively.
Table <ref> lists the derived channel pairs, ij's, along with their coupling coefficient c_ij^z for a flat spectrum and sensitivity a_ij each derived as described in Appendix <ref>.
All left side channels i's are those in Nagoya station in the northern hemisphere, while the other side channel j subtracted from the channel i in equation (<ref>) is all in Hobart station in the southern hemisphere.
This is reasonable as a result of the optimization maximizing the NS anisotropy sensitivity.



A solid line in Figure <ref> displays a differential response c_ij^z/P for each channel pair ij in table <ref>.
This differential response function has a positive value when the asymptotic direction of channel i at rigidity P points to a higher latitude (more northward) direction than channel j.
In other words, its negative value in the lower rigidity region in Figure <ref> indicates that the NS relation between channels i and j is swapped for low rigidity GCRs.
The low-rigidity GCRs passing through the mid- or low-latitude magnetosphere, such as those observed by Nagoya or Hobart detectors, are largely deflected from the original viewing angle of each directional channel by the geomagnetic field, causing this negative response.
This problem is solved by placing detectors in Arctic and Antarctic regions where the geomagnetic effect is minimized, and such an observation plan is mentioned in Section <ref>.
The difference of the differential response between channel pairs displayed in Figure <ref> allows us to derive the rigidity spectrum ξ_z^TA(P) of the NS anisotropy, as described in the next section.

Finally, the dashed line in Figure <ref> displays the differential response of Nagoya-GG for comparison.
The channel pairs ij's used in this paper (solid lines) generally have higher response peaks than Nagoya-GG, thanks to the optimization of the channel pairing.



 §.§ Bayesian estimation of the rigidity spectrum
 
In conventional approaches, certain analytical functions are proposed for the rigidity spectrum of the anisotropy.
Most typically use a power-law function,

    ξ_z^TA(P) = β P^γ,

and searches for parameters (β, γ) fitting to the observed counting rates <cit.>.
Although the power-law function is inapplicable for a spectrum crossing zero or with an insignificant magnitude, such a case can happen in the anisotropy component because the anisotropy is a vector quantity, unlike the absolute density spectrum.
<cit.> suggests that this is a cause of the unstable fluctuations of the power-law index γ in their results.
A single analytical function such as the power-law spectrum is probably too restrictive a premise to express all spectral shapes appearing in solar modulation phenomena. 
On the other hand, this paper derives the NS anisotropy spectrum and its error range by the Bayesian approach as described below.
The Gaussian process is introduced as a prior probability distribution in Bayesian estimation, confining acceptable ranges of the smoothness and magnitude of the spectrum instead of the analytical function.
This approach makes the spectrum estimation tolerant of varying spectral shapes, allowing us to trace dynamic variations of the spectrum.

Hereafter, the GCR rigidity P [GV] is expressed by its logarithm q = log_10(P) in equations.
We approximate the rigidity spectrum by a step function with a finite number (N) of rigidity bins as

    ξ_z^TA(q) = θ_k for q_k-1≤ q < q_k

where k=1,…,N indicates each rigidity bin and q_k = log_10 P_k represents an upper limit of the k-th rigidity bin.
The NS anisotropy θ_k in each rigidity bin k is approximated to be constant for the interval q_k-1≤ q < q_k.
The rigidity spectrum ξ_z^TA(P) is now expressed by an N-dimension vector parameter θ = ( θ_1,θ_2,…,θ_N)^⊤, and its mean and error can be expressed by a probability distribution in the N-dimension space in which each dimension corresponds to each rigidity bin k.
Based on Figure <ref>, our responsive rigidity range is from ∼10 to ∼400 GV, and we split this range into N=8 bins.
Table <ref> lists each rigidity bin with its index k and upper boundary q_k.
Median rigidity P_k^m [GV] of each rigidity bin satisfies log_10 P_k^m = (q_k-1 + q_k)/2.


The differential response c_ij^z/P is also discretized and expressed by a matrix

    C_lk = ∫_P=P_k-1^P_kc_ij^z(P)

where each channel pair ij is replaced with l = 1,2,…,8 in the index of Table <ref>.
From equation (<ref>), expected value of the observable η_l^TA = η_ij^TA from the parameter θ is derived as

    η̃_l^TA = ∑_k=1^N C_lkθ_k = (Cθ)_l.

Therefore, the conditional probability of the observed values η_l^TA's for the parameter θ is modeled as

    𝒫(η | θ)
               = ∏_l 𝒩(η_l^TA | η̃_l^TA, (σ_l^TA)^2)
       ∝∏_l
            exp[- η_l^TA - (Cθ)_l^2/2 (σ_l^TA)^2]
            = exp[-1/2∑_l η_l^TA - (Cθ)_l^2/(σ_l^TA)^2]

where η is a vector consisting of η_l^TA's for all channel pairs and σ_l^TA=σ_ij^TA in equation (<ref>).
The symbol 𝒩(x | μ,σ^2) expresses a normal distribution for a variable x with its mean μ and variance σ^2.

On the other hand, the Gaussian process is introduced by considering the probability distribution of the spectrum ξ_z^TA(q) as a multivariate normal distribution in which each dimension corresponds to each value of q; the distribution is defined in an infinite-dimensional space in principle.
In this practical case, which represents the spectrum by a finite number (N) of rigidity bins, the Gaussian process is expressed by a multivariate normal distribution for the parameter θ in the N-dimension space, as

    𝒫(θ)
        = 𝒩(θ | Θ_G, Σ_G)

where Θ_G and Σ_G are the mean vector and covariance matrix of the Gaussian process, respectively.
We adopt the radial basis function kernel for the covariance matrix as

    Σ_G,kk' = σ_G^2 exp( -|q_k - q_k'|^2/b^2 ).

The diagonal component of Σ_G is provided as Σ_G,kk = σ_G^2, representing a variance of the parameter θ_k.
The non-diagonal component divided by the variance σ_G^2, i.e., exp(-|q_k - q_k'|^2 / b^2), represents a correlation coefficient between θ_k and θ_k', confining an acceptable range of the smoothness of the spectrum.
Detailed descriptions of the conditional probability 𝒫(η | θ) and Gaussian process 𝒫(θ) are provided in Appendix <ref> and <ref>.

Implementing the Gaussian process 𝒫(θ) as a prior probability distribution in Bayesian estimation, the posterior distribution of the parameter θ is derived as

    𝒫( θ | η ) ∝𝒫( η | θ ) 𝒫(θ).

The Gaussian process 𝒫(θ) forms a conjugate prior distribution of the normal distribution 𝒫( η | θ ) in this Bayesian estimation.
Therefore, the posterior distribution is also a multivariate normal distribution written as

    𝒫( θ | η ) =
        𝒩( θ | Θ, Σ).

This equation can be analytically solved, providing another advantage of our method compared to the conventional approach, which generally requires a numerical fitting of parameters, such as the index γ for the power-law spectrum.
From Appendix <ref>, the mean vector and covariance matrix in equation (<ref>) are derived as

    Θ   =
            Σ_G ( Σ_L + Σ_G )^-1Θ_L
            + Σ_L ( Σ_L + Σ_G )^-1Θ_G and
    Σ   =
            Σ_L ( Σ_L + Σ_G )^-1Σ_G.

The vector Θ_L and matrix Σ_L are parameters of the conditional probability 𝒫(η | θ) and derived as

    Θ_L = Σ_L pandΣ_L = Q^-1wherep = C^⊤WηandQ = C^⊤WC.

The weight matrix W is a diagonal matrix with its diagonal elements W_ll = 1/(σ_l^TA)^2.
In the hyperparameters of the Gaussian process, the mean vector Θ_G is determined every year, while parameters for the covariance matrix Σ_G are set as common values σ_G = 0.1% and b = 0.4 [log_10(GV)] for all years, as described in Appendix <ref>.
The marginal distribution of 𝒫(θ | η) in each rigidity bin k is derived as

    𝒫(θ_k | η)
        = ∫𝒫(θ | η) θ_(k)'
        = ∫𝒩(θ | Θ, Σ) θ_(k)'
        = 𝒩(θ_k  | Θ_k, Σ_kk)

where θ_(k)' is a subset of the vector θ from which an element θ_k is removed.
Therefore the 1-σ error σ_k attached to the mean spectrum Θ_k is derived from Σ as

    σ_k^2 = Σ_kk.

The next section presents the results of the mean spectrum Θ and its error σ_k derived for each year. 



§ RESULTS


Figure <ref> represents derived rigidity spectra of the NS anisotropy in (a) 2009 on left panels and (b) 2015 on right panels, as sample years around solar activity minimum and maximum, respectively.
A solid black line in each upper panel is the mean spectrum Θ derived by equation (<ref>), and dashed lines above and below the mean spectrum show its ±1σ error range derived by equation (<ref>) from observed values η_l^TA's in each year.
A horizontal blue line is the mean vector Θ_G of the prior probability distribution.
The mean spectrum represented by the solid black curve seems biased to the prior distribution in the outer rigidity bins (∼10 GV and ∼300 GV), resulting in some kinks of the mean spectrum near these rigidities.
This is due to the less responses in these rigidities as shown in Figure <ref>, while it negligibly affects scientific discussions, covered by large errors in these rigidity bins.

Solid circles in each lower panel in Figure <ref> show a scatter plot between the observable η_l^TA's of all channel pairs (vertical axis) and their expected values η̃_l^TA's from the mean spectrum (horizontal axis) in each year.
The expected value η̃_l^TA is calculated by replacing θ in equation (<ref>) with the mean spectrum Θ.
They are consistent with observed values within error bars in each year, ensuring that our new method successfully reconstructs the observed rigidity spectrum of the NS anisotropy with muon counting rates.
An open circle in each lower panel represents the observed and expected values of Nagoya-GG.
The Nagoya-GG is not used to reconstruct the NS anisotropy spectrum in this paper, but its expected value from the derived spectrum shows agreement with its observed value in a comparable range with solid circles.
This result demonstrates that our new analysis method is consistent with the conventional method while providing advanced insights, such as the yearly variation of the anisotropy spectrum.

It is also found that the standard error of Nagoya-GG is smaller than solid circles by a factor less than ∼1/2.
One of the causes is that Nagoya-GG double-counts channel pairs in equation (<ref>), reducing the statistical error by a factor of ∼1/√(2).
The remaining factor is likely a local effect different in different stations and remained uncorrected <cit.>, representing the limit of our method of directly combining multiple stations.
On the other hand, the solid circles show larger magnitudes than the open circle in each panel, enlarging their significance, thanks to a higher response of each channel pair than Nagoya-GG as shown in Figure <ref>.


The harder slope of the mean spectrum in the solar maximum (2015) than the minimum (2009) in Figure <ref> represents a common feature of the solar cycle variation of the NS anisotropy, as demonstrated by Figure <ref>.
It displays the mean spectrum every year from 1997 to 2020, split into the solar-activity ascending and descending phases in the solar cycle 23 or 24 in each panel.
The edge rigidity bins around ∼10 GV and ∼300 GV are truncated in this figure, considering their relative unreliability mentioned above.
Each year is denoted by a legend below each panel, and the gradation of line colors indicates the solar activity, where the red and blue colors correspond to years around the activity maximum and minimum, respectively.
Overall, the blue lines have a softer slope than the red lines in each panel, indicating that the NS anisotropy exhibits a softening of its rigidity spectrum in the solar activity minimum.

<cit.> also suggested a variation of the NS anisotropy's rigidity spectrum according to the solar activity, and since then, quantitative estimation of the spectrum has been rarely reported.
<cit.> provided only two average spectra around the maximum and declining solar activity phases, each in 1969-1970 and 1971-1973, and concluded an insignificant difference of the spectrum between the two periods.
In that analysis, multiple types of observations, including underground muon detectors, were required, and the power-law rigidity spectrum of the anisotropy was presumed.
These limitations probably prevented the study from revealing the dynamic variation of the anisotropy with a better temporal resolution.
On the other hand, our approach successfully reconstructs the spectrum every year, as shown in Figure <ref>, thanks to the better tolerance of varying spectral shapes.
A physical origin of the revealed variation of the spectrum is discussed in Section <ref>.



§ DISCUSSIONS



 §.§ Modulation parameters

The sector reversal of the NS anisotropy is expected as a consequence of the diamagnetic drift of GCRs with the radial density gradient g(P) in the IMF, expressed as <cit.>

    ξ_z^TA(P)/cosδ_E∼ R_L  g(P) sinψ = P/cB g(P) sinψ.

The radial density gradient g(P) is related to the GCR density as

    g(P) = 1/U(r_E,P). ∂ U(r,P)/∂ r|_r=r_E

where U(r) is a GCR density at the distance r from Sun and r_E is the position of Earth.
In equation (<ref>), it is approximated that the density gradient of GCRs in the ecliptic plane is dominated by its radial component, especially in the Bartels rotation mean, which averages the azimuthal (GSE-y) component out to zero.
The gyro-radius R_L of GCRs is derived as R_L = P/(cB) from the GCR rigidity P, IMF magnitude B and light speed c.
Derivation of the NS anisotropy in Section <ref> is based on the GEO coordinate system, and ξ_z^TA is expected to be a projection of the sector reversal of the GSE-z component of the anisotropy onto the GEO-z axis.
Therefore, the NS anisotropy is divided by the factor cosδ_E, where δ_E is the inclination angle of Earth's rotation axis from GSE-z axis, in equation (<ref>) to convert the NS anisotropy ξ_z^TA from the GEO to GSE coordinate system.

Referring to the derivation procedure of η_ij^TA in Section <ref>, GSE-x and y components of the IMF vector are calculated from OMNIWeb data as

    B_x(y)^TA = B_x(y)^T - B_x(y)^A/2

where B_x(y)^T and B_x(y)^A are averages of the GSE-x(y) component of the IMF in the Toward and Away sectors respectively in each year.
The IMF magnitude and spiral angle ψ in equation (<ref>) are consequently derived as

    B = √((B_x^TA)^2 + (B_y^TA)^2)andsinψ = - B_y^TA/B.


Replacing ξ_z^TA(P) in equation (<ref>) with the mean spectrum Θ_k derived in Section <ref>, the most probable value of the radial density gradient in each rigidity bin k is derived as

    g_k ∼ - c/P_k^m cosδ_E·B^2/B_y^TAΘ_k

where the rigidity P is represented by its median P_k^m in each rigidity bin.
A square of its error is estimated as

    σ^2(g_k)    = ∂ g_k/∂ B_x^TAσ(B_x^TA) ^2 +
                   ∂ g_k/∂ B_y^TAσ(B_y^TA) ^2 +
                   ( ∂ g_k/∂Θ_kσ_k )^2
       = ( c/P_k^m cosδ_E B_y^TA )^2
        [  2 B_x^TAΘ_k σ(B_x^TA) ^2
        + ( B_y^TA - (B_x^TA)^2/B_y^TA ) Θ_k σ(B_y^TA) ^2
        + ( B^2 σ_k )^2 ]

where σ(B_x^TA) and σ(B_y^TA) are errors of B_x^TA and B_y^TA respectively.
Similarly to the error of η_ij^TA in equation (<ref>), these errors are derived as

    σ^2(B_x(y)^TA) = σ^2(B_x(y)^T) + σ^2(B_x(y)^A)/4

where σ(B_x(y)^T) and σ(B_x(y)^A) are derived from standard deviations of the GSE-x(y) component of the daily IMF in the Toward and Away sectors respectively.


To reveal a solar cycle variation of the radial density gradient along with its rigidity dependence, we extract g_k's in k=3 and k=6 bins corresponding to the rigidity of P∼30 GV and P∼130 GV, respectively.
Figure <ref>c displays temporal variations of them in a unit of %/AU and a logarithmic vertical axis as a function of year.
The solar activity and interplanetary states are represented by Figures <ref>a and <ref>b;
the panel (a) displays yearly averages of the IMF magnitude B (blue line) and spiral angle ψ (red line) derived in equation (<ref>), red line in the panel (b) shows yearly averages of the solar wind speed V, and black points in the panel (b) are occurrence rate of the X-class solar flares every week.
In the same manner as the IMF, the yearly average of the solar wind speed is derived as V = (V^T + V^A)/2 where V^T and V^A are averages of the solar wind speed in the OMNIWeb for the Toward and Away sectors, respectively, every year.
The IMF spiral angle ψ and solar wind speed V show an anti-correlation with each other, following the Parker spiral picture.
The X-class flare occurrence rate is derived from the Konus-WIND flare catalog <cit.>.

We can identify a solar cycle variation of the radial density gradient at P∼130 GV (solid circles in Figure <ref>b), where the density gradient is suppressed in the solar activity minima around 1997, 2009, and 2020.
The density gradient at P∼30 GV (open circles in Figure <ref>b) also shows minima at least around 2009 and 2020, but its relative variation expressed in the logarithmic scale is smaller than that of P∼130 GV.
We can conclude that the solar cycle variation of the NS anisotropy's rigidity spectrum as seen in Figure <ref> is attributed to this rigidity-dependent solar cycle variation of the density gradient.

Despite the hard spectrum of the NS anisotropy overall in Figure <ref>, the radial density gradient at P∼30 GV is generally larger than P∼130 GV in Figure <ref>b, indicating a soft spectrum of the density gradient.
This soft spectrum is reasonable from a picture of the smaller solar modulation in higher rigidity GCRs.
From a technical viewpoint, it is caused by the factor 1/P_k^m multiplied in equation (<ref>).
While a quantitative reproduction of the observed anisotropy has been rarely performed by heliospheric simulations except for a small number of studies <cit.>, the radial density gradient, or radial distribution of GCR density, is predicted by some simulation works for a realistic heliosphere model including the termination shock and heliosheath <cit.>. 
Our result on the radial density gradient will provide an observational validation of such heliospheric models.

It is also worth focusing on some years deviating from the general picture of the density gradient with a soft spectrum, as seen in 1999-2001, 2003, 2007, and 2013 in Figure <ref>c.
While being insignificant because of the large errors, these years have rapid suppression of the density gradient only in ∼30 GV GCRs from the previous years, resulting in comparable density gradients between ∼30 GV and ∼130 GV GCRs.
Such a hard spectrum of the density gradient causes a noticeably hard spectrum of the NS anisotropy in these years, as displayed in Figure <ref>.
Most of these years, 1999-2001, 2003, and 2013, correspond to the years when the area or number of low-latitude coronal holes was enhanced on Sun, as seen in the results of automatic identification of the coronal holes <cit.>.
The remaining year, 2007, is just after the coronal mass ejection (CME) event on December 13th, 2006, which was the biggest halo CME since the “Halloween storm” in 2003 as of the study by <cit.>.
From black points in Figure <ref>b, it is found that a few X-class flare events, including the CME on December 13th, 2006, successively occurred after a relatively calm period without X-class flares in 2005-2006. 
A similar situation was found in the last half of 2017, where three successive X-class flares were detected after over ∼2 years with no X-class flares.
After this event, the density gradient shows a spectral hardening in 2017-2018, similar to 2006-2007, as displayed by the open and solid circles getting closer to each other in these years in Figure <ref>b.
These results imply that the post-storm interplanetary state after discrete solar eruptions or the low-latitude coronal holes mentioned above possibly caused unusual modulation effects suppressing the density gradient in lower rigidity GCRs, although more detailed analyses are required.

Assuming a radial balance between the diffusion flow and solar wind convection of GCRs further provides an estimation of their diffusion coefficient, as

    κ_rr(P) g(P) ∼2 + Γ/3 V

where Γ∼ 2.7 is a power-law index of the energy spectrum of the GCR density.
The radial diffusion coefficient κ_rr of GCRs at the position of Earth is expressed as

    κ_rr = κ_∥cos^2 ψ + κ_⊥sin^2 ψ

where κ_∥ and κ_⊥ are diffusion coefficients parallel and perpendicular to the IMF.
Therefore we can estimate the radial diffusion coefficient κ_rr,k in each rigidity bin k as

    κ_rr,k∼2 + Γ/3·V/g_k.

The error range of κ_rr,k can no longer be approximated by a symmetric normal distribution because it is dominated by the error of g_k, which is inversely proportional to κ_rr,k.
We approximate the lower and upper limits of a 1σ confidence interval of κ_rr,k by

    κ_rr,k^- = 2 + Γ/3·V/g_k + σ(g_k)andκ_rr,k^+ = 2 + Γ/3·V/g_k - σ(g_k)

respectively.
A standard deviation of the solar wind speed V is ignored because it is negligible compared to σ(g_k).
Figure <ref>d displays temporal variations of the radial diffusion coefficient κ_rr for ∼30 GV and ∼130 GV GCRs.
It demonstrates a solar cycle variation of the diffusion coefficient in both rigidities as well as the radial density gradient but negatively correlates with solar activity.
We have to note that the diffusion and convection balance in equation (<ref>), known as the force-field approximation <cit.>, is based on the diurnal anisotropy's phase observed to be around 18:00 local solar time <cit.>.
However, <cit.> and <cit.> showed a 22-year cycle excursion of the phase to earlier local time, suggesting that the left side of equation (<ref>) can be ∼4/5 of the right side in the case of such an excursion.
Therefore, this estimation of the diffusion coefficient κ_rr can have a systematic error of ∼20% in addition to equation (<ref>).


Based on a comparable sin^2ψ with cos^2ψ in the average IMF and κ_⊥/κ_∥≪ 1 supported by previous studies <cit.>, we obtain the parallel diffusion coefficient as

    κ_∥,k∼κ_rr,k/cos^2ψ

from equation (<ref>) for the rigidity bin k.
Therefore, the parallel mean-free-path length of GCRs and its error range are estimated as

    λ_∥,k = 3κ_∥,k/v∼3κ_rr,k/vcos^2ψ,  λ_∥,k^- ∼3κ_rr,k^-/vcos^2ψ,andλ_∥,k^+ ∼3κ_rr,k^+/vcos^2ψ

where v ∼ c is the particle velocity and approximated by a light speed c for relativistic GCRs with P > 10 GV.
Blue and red curves in Figure <ref> display rigidity spectra of the parallel mean-free-path length in 2009 and 2015, respectively, as sample years in the solar activity minimum and maximum.
The outer rigidity bins with substantial errors are truncated in this figure.
The gray shaded area is “Palmer consensus” region which was proposed by <cit.> and has been validated mainly by solar energetic particle (SEP) observations and numerical simulations <cit.>.
Our result extends the mean-free-path length estimation into higher rigidity region ≤200 GV, and from Figure <ref> it is expected to smoothly connect to the consensus region by extrapolating the rigidity dependence.

Analyses of the GCR density and its rigidity spectrum also provide estimations of the diffusion coefficient or mean-free-path length in >10 GV region <cit.> by surveying optimal heliospheric parameters fitting to observations.
Our result of the mean-free-path length in Figure <ref> is consistent with their results in the order of magnitude, while our result features the unprecedented temporal resolution on a yearly basis and the wide rigidity range of over one order of magnitude.
It is also noted that their approach using the GCR density requires the assumption on the whole heliospheric structure, such as a spatial distribution of the diffusion coefficient, because the GCR density at Earth reflects an integration of the modulation effects from the outer edge of the heliosphere.
On the other hand, the NS anisotropy can provide an estimation of the diffusion coefficient free from such an assumption as well as the SEP observations in <10 GV region.
This is one of the advantages of the anisotropy observation, leading to a relatively reliable estimation of the diffusion coefficient.



 §.§ Future perspective of the NS conjugate observation


Recently, we started a muon observation at Syowa Station, Antarctica, in 2018 in collaboration with the National Institute of Polar Research (NIPR), Japan.
NIPR also has a collaborative research framework with Iceland, located around an Arctic geomagnetic conjugate point with Syowa Station.
The polar conjugate observation in Iceland and Syowa Station will be an expanded concept of the current Nagoya-Hobart observation, and we provide a result of its quantitative simulation in this section.
The black bold line in Figure <ref> displays the differential response to the NS anisotropy for vertical channels of muon detectors at Husafell, Iceland, and Syowa Station.
The virtual detector at Husafell is set to have an equivalent geometry with the Syowa detector having a 1×2 m^2 detection area.
Other lines in Figure <ref> are the responses of channel pairs used in this study, i.e., repeats of those in Figure <ref>.
Notably, only the Husafell-Syowa pair has a substantial response in the lower rigidity (<20 GV) region, without a polarity reversal into a negative response seen in other channel pairs.
The north- and southward inclined channels in Nagoya and Hobart detectors respectively view high-latitude directions, but the geomagnetic field around the mid-latitude stations deflects cosmic-ray trajectories into lower-latitude directions as briefly discussed in Section <ref>.
This effect is more substantial for lower-rigidity cosmic rays and prevents us from gaining a significant response to the NS anisotropy in lower-rigidity regions only by mid-latitude stations.
Arctic and Antarctic regions, such as Iceland and Syowa Station, are relatively free from this effect and are predicted to allow us to observe the NS anisotropy in lower-rigidity GCRs accurately.
While having less angular resolutions than muon detectors, neutron monitors in Arctic and Antarctic regions are also expected to be responsive to the NS anisotropy in further lower rigidity (∼1 GV) regions.
Therefore, the concept of Iceland-Syowa conjugate muon observation can be a bridge to connect the muon and neutron detector networks.



§ CONCLUSION

A new analysis method to derive the rigidity spectrum of GCR anisotropy from ground-based observations has been developed and demonstrated by revealing a yearly variation of the NS anisotropy's spectrum.
In this method, atmospheric temperature effects on cosmic-ray muons are directly corrected by the meteorological reanalysis data, allowing for combining multiple muon detectors in a network observation free from individual local effects.
General graph matching in graph theory is adopted to survey optimal combinations of directional channels, which ensures high sensitivity to the NS anisotropy.
The highlight of our analysis method is Bayesian estimation with the Gaussian process, which has the potential to be applied to any unfolding problem of the ground-based observations to derive GCR properties in space.
The Gaussian process only confines the acceptable ranges of the spectrum value and its smoothness without supposing any analytical function, providing a sufficient tolerance of varying spectral shapes to trace dynamic variations of the GCR spectrum in solar modulation phenomena.

Previous works deriving rigidity spectra of the anisotropy use simultaneous observation data by multiple detector types, i.e., ground-based muon detectors, underground muon detectors, and neutron monitors <cit.>.
On the other hand, this paper succeeded in deriving the NS anisotropy's spectrum only by ground-based muon detectors in Nagoya and Hobart.
This relaxed requirement for observational data will make the anisotropy analysis a more common approach in cosmic-ray research.
The minimized number of observations also ensures a uniform dataset for an extended period, allowing this study to reveal the solar cycle variation of the NS anisotropy's spectrum.
Softening of the spectrum in the solar activity minima was discovered, and it is attributed to the rigidity-dependent variation of the radial density gradient of GCRs.
The diffusion coefficient or mean-free-path length of GCRs is subsequently derived based on the force-field approximation, and it is demonstrated that our analysis expands the mean-free-path length estimation into ≤200 GV region from <10 GV region achieved by SEP observations.
The rigidity-dependent diffusion coefficient has been a critical problem in elucidating the GCR propagation <cit.>, emphasizing the importance of the anisotropy observation and our result.

In addition to expanding the muon detector network into polar regions as described in Section <ref>, applying our analysis scheme to a broader range of cosmic-ray studies is desired.
Reconstruction of the three-dimensional anisotropy or density spectrum of GCRs is one such scientific target, and it can, in principle, be achieved by generalizing the Gaussian process used in this paper.
Short-term disturbance phenomena, including CME events, are also crucial topics along with the solar cycle variation analyzed in this paper, and an application of our analysis method to such an event will be presented in the future.


This research was supported by the “Strategic Research Projects” grant from ROIS (Research Organization of Information and Systems) and by JSPS KAKENHI Grant Number JP22KK0049.
The GMDN project is partially supported by ROIS-DS-JOINT (030RP2023), JARE AJ1007 program, and JSPS KAKENHI Grant Number JP24K07068.
The observations with Nagoya and Hobart muon detectors are supported by Nagoya University and Australian Antarctic Division, respectively.




§ DATASET PREPARATION

In this appendix, we describe the preparation procedure of muon counting-rate data performed before the analysis procedure in the main text, mainly the correction for atmospheric temperature effect on muons.
The temperature correction method by <cit.> defines a mass-weighted temperature for each station as

    H_ st = ∑_h=1^h_max-1 w_h T_h + T_h+1/2

where T_h is the atmospheric temperature at altitude h above the location of each station and h_max indicates the top of the atmosphere.
The weight w_h is derived as w_h = (z_h - z_h+1)/z_0 where z_h is the atmospheric depth at the altitude h and z_0 = z_1 - z_h_max.
Muon counting rates N_i(t) corrected for atmospheric pressure effects on each time t are published by <cit.>.
The atmospheric temperature effect Δ N_i(t) on N_i(t) is expressed as 

    Δ N_i(t) = α_i  H_ st(t) - H_ st^0

for the directional channel i belonging to each station.
The constant H_ st^0 is an average level of the mass-weighted temperature and set at 253 or 250 [K] for Nagoya or Hobart station, respectively.
The coefficient α_i for each directional channel is derived by <cit.> and ranges from -0.28 to -0.24 [%/K] in Nagoya station's channels and from -0.24 to -0.22 [%/K] in Hobart station's channels. 
The temperature-corrected counting rate is derived as

    N_i^ corr(t) = N_i(t) - Δ N_i(t).


The altitude profile of the atmospheric temperature is provided as meteorological reanalysis data by GDAS (Global Data Assimilation System) from 2005 and NCEP/NCAR (National Centers for Environmental Prediction and National Center for Atmospheric Research) until 2004, both published in NOAA's Air Resources Laboratory website <cit.>.
It provides the altitude and temperature T_h at each isobaric surface from 20 hPa to 1000 hPa at a designated location.
We set the ground-level altitude, h=1, at the isobaric surface with the lowest altitude above 500 m.
The top of the atmosphere, h_max, is approximated by the highest isobaric surface in the data, corresponding to the pressure of 20 hPa.
In calculating the weight w_h, the atmospheric pressure of each isobaric surface is used in place of the atmospheric depth z_h.
Temporal resolutions of the meteorological reanalysis data are 3 hours in GDAS data and 6 hours in NCEP/NCAR data, respectively, which are insufficient for hourly muon counting-rate data.
We first calculate the mass-weighted temperature H_ st(t) on a 3-hour or 6-hour basis using the provided data.
Then, it is interpolated by a linear interpolation between timestamps, and hourly H_ st(t) for the temperature correction is derived.

The corrected muon counting rate is converted into a deviation from its Bartels rotation average,

    I_i(t) = N_i^ corr(t)/N_i^ BR

where N_i^BR is N_i^ corr(t) averaged over each Bartels rotation.
Hourly counting rates >2% below or above the Bartels rotation average N_i^ BR are omitted to eliminate disturbance events such as the Forbush decrease.



§ SENSITIVITY OF EACH CHANNEL PAIR TO THE NS ANISOTROPY

Optimization of the channel pairing described in Section <ref> is based on the sensitivity a_ij defined for each channel pair ij.
In this appendix, we describe its derivation procedure.
From the discussion with equation (<ref>) in Section <ref>, the perturbation on η_ij^TA from anisotropy components other than the NS anisotropy is estimated as

    Δη_ij^TA = ∫_P=0^∞( ϵ_0^TAc_ij^0 + ϵ_c^TAc_ij^d + ϵ_s^TAs_ij^d ).

Previous studies <cit.> report nearly flat spectra for the anisotropy components and a soft spectrum proportional to ∼ P^-1 for the short-term density variations.
Based on these suggestions, we temporarily adopt an ad hoc simplification on the rigidity spectra only for the channel pairing as

    ξ_z^TA(P) = ξ_z^TA,  ϵ_0^TA(P) = ϵ_0^TA(P/60[ GV])^-1,  ϵ_c^TA(P) = ϵ_c^TA,andϵ_s^TA(P) = ϵ_s^TA.

The constant 60 GV for the density variation ϵ_0^TA(P) is introduced as a representative rigidity of muon detectors.
Equations (<ref>) and (<ref>) are simplified as

    η_ij^TA∼ξ_z^TA c_ij^zandΔη_ij^TA∼ϵ_0^TA c_ij^0 + ϵ_c^TA c_ij^d + ϵ_s^TA s_ij^d

where c_ij^z, c_ij^0, c_ij^d, and s_ij^d are coupling coefficients for the assumed spectra, defined as

    c_ij^z = ∫_P=0^∞c_ij^z,  
        c_ij^0 = ∫_P=0^∞(P/60[ GV])^-1c_ij^0,  
        c_ij^d = ∫_P=0^∞c_ij^d,and
        s_ij^d = ∫_P=0^∞s_ij^d.

From equation (<ref>), the optimized pairs can be defined as those maximizing the response ∂η_ij^TA/∂ξ_z^TA∼ c_ij^z to the NS anisotropy while minimizing the response to other anisotropy components which is estimated as

    (∂Δη_ij^TA/∂ϵ_0^TA)^2 +
        (∂Δη_ij^TA/∂ϵ_c^TA)^2 +
        (∂Δη_ij^TA/∂ϵ_s^TA)^2
        ∼(c_ij^0)^2 + (c_ij^d)^2 + (s_ij^d)^2.

Therefore, we define the sensitivity of each channel pair ij to the NS anisotropy as

    a_ij = {[ c_ij/√((c_ij^0)^2 + (c_ij^d)^2 + (s_ij^d)^2)              for c_ij^z > 0.7 and a_ij > 1.6;                                            0               for c_ij^z ≤ 0.7 or a_ij≤ 1.6. ].

The thresholds, c_ij^z > 0.7 and a_ij > 1.6, are introduced to truncate channel pairs with worse sensitivities than Nagoya-GG, whose response and sensitivity are c_ij^z=0.63 and a_ij=1.51 respectively.




§ LIKELIHOOD FUNCTION AND GAUSSIAN PROCESS FOR THE SPECTRUM PARAMETERS

In equation (<ref>), the conditional probability distribution 𝒫(η | θ) of observed values η_l^TA's for the parameter θ is equivalent to a likelihood function for θ, which is generally used in the maximum likelihood estimation.
The parameters Θ_L and Σ_L in equations (<ref>) and (<ref>) are respectively identical to the mean vector and covariance matrix of the likelihood function expressed in a multivariate normal distribution.
In this appendix, we prove these equations and briefly inspect the profile of the likelihood function.
The effect of the Gaussian process is also described by visualizing the distributions of the likelihood and posterior probability of Bayesian estimation.
        
The exponent in equation (<ref>) is transformed as

    -1/2∑_l η_l^TA - (Cθ)_l^2/(σ_l^TA)^2   = -1/2∑_l 1/(σ_l^TA)^2∑_k=1^N∑_k'=1^N C_lk C_lk'θ_k θ_k'
               - 2 ∑_k=1^Nη_l^TA C_lkθ_k
               + (η_l^TA)^2
       = -1/2∑_k=1^N∑_k'=1^Nθ_k
               (∑_l C_lk C_lk'/(σ_l^TA)^2) θ_k'
               - 2 ∑_k=1^Nθ_k ∑_l C_lkη_l^TA/(σ_l^TA)^2
               + ∑_l (η_l^TA/σ_l^TA)^2
       = -1/2( θ^⊤Qθ
               - 2 θ^⊤p + u )

where u = ∑_l (η_l^TA / σ_l^TA)^2 and Q and p are defined in equation (<ref>).
On the other hand, a square-completed quadratic form for θ is defined and expanded as

    (θ - Θ_L)^⊤Q(θ - Θ_L) + u'
               = θ^⊤Qθ
                - θ^⊤QΘ_L
                - Θ_L^⊤Qθ
                + Θ_L^⊤QΘ_L + u'
       = θ^⊤Qθ
                - θ^⊤( QΘ_L + Q^⊤Θ_L )
                + Θ_L^⊤QΘ_L + u'
       = θ^⊤Qθ
                - 2 θ^⊤QΘ_L
                + Θ_L^⊤QΘ_L + u'

where Q^⊤ = Q is used.
Comparing the parentheses in equation (<ref>) and equation (<ref>) for θ, we obtain p = QΘ_L, which is also identical to the definition of Θ_L in equation (<ref>).
The parameter u' is also defined from the comparison as u' = u - Θ_L^⊤QΘ_L.
Replacing the parentheses in equation (<ref>) with the square-completed expression in equation (<ref>) and revisiting the conditional probability distribution 𝒫(η | θ) in equation (<ref>), we obtain

    ℒ(θ) ∝exp- 1/2(θ - Θ_L)^⊤Σ_L^-1(θ - Θ_L) + u' ∝𝒩(θ | Θ_L, Σ_L)

where 𝒫(η | θ) is replaced by a symbol ℒ(θ) to express the likelihood function for θ.
The definition Σ_L = Q^-1 in equation (<ref>) is also used.
Therefore, it is proven that the likelihood function is proportional to a multivariate normal distribution for θ where parameters Θ_L and Σ_L in equation (<ref>) are identical to the mean vector and covariance matrix of the distribution, respectively.
The posterior distribution in equation (<ref>) is equivalent to a product of multivariate normal distributions ℒ(θ) ∝𝒩(θ | Θ_L, Σ_L) and 𝒫(θ) ∝𝒩(θ | Θ_G, Σ_G).
From a formula for a product of multivariate normal distributions, the mean vector Θ and covariance matrix Σ of the posterior distribution are derived from Θ_L, Σ_L, Θ_G, and Σ_G, as described in equation (<ref>).


Figure <ref>a shows a conditional distribution of the likelihood ℒ(θ) in the θ_3 - θ_4 parameter space for the sample year 2015, where other parameters, θ_k for k  3  or  4, are fixed at those in the mean spectrum Θ derived by equation (<ref>).
The distribution is normalized so that its integration in the θ_3 - θ_4 space equals 1.
The likelihood distribution is widely extended along a θ_3 + θ_4 = const. line, indicating that these parameters cancel out each other in the likelihood.
This demonstrates the difficulty to uniquely determine these parameters, or the spectrum values in adjacent rigidity bins by using the likelihood function.
Figure <ref>b demonstrates how the Gaussian process works to overcome this problem.
It is a θ_3 - θ_4 parameter space distribution of the posterior probability 𝒫( θ | η ) in 2015, derived by multiplying the Gaussian process 𝒫(θ) to the conditional probability 𝒫( η | θ ), or likelihood function ℒ(θ), in equation (<ref>). 
The high-confidence region is shrunk in Figure <ref>b compared to Figure <ref>a, within a reasonable region where the adjacent rigidity bins θ_3 and θ_4 have comparable values with each other.
This limitation by the Gaussian process prevents a discontinuous spectrum, enabling us to determine the spectrum while keeping a sufficient agreement between the derived spectrum and observed data.
On the other hand, how the likelihood works is visualized by the more significant likelihood in the θ_4>θ_3 region than θ_4<θ_3 region in Figure <ref>a.
It suggests that the spectrum value θ_k increases with the rigidity or its index k, resulting in the hard spectrum in 2015 as displayed in Figure <ref>b.

<cit.> also approximated a rigidity spectrum of the semi-diurnal anisotropy amplitude by a step function, as well as our formulation for the NS anisotropy in equation (<ref>).
Then, they derived the spectrum by a least-square method which is equivalent to the maximum likelihood estimation.
In their analysis, the spectrum was split into only seven rigidity bins in a range from ∼10 GV to ∼2000 GV, and not only ground-based muon detectors but also underground detectors were used.
This rough rigidity resolution and expanded dataset likely allowed them to derive the spectrum only by the likelihood function ℒ(θ) without the Gaussian process used in this study.
However, this approach is a trade-off with some accuracy, such as the rigidity bin and temporal resolution.
They derive only a several-year average of the spectrum, unlike our analysis deriving a solar cycle variation of the NS anisotropy spectrum on a yearly basis.



§ HYPERPARAMETERS OF THE GAUSSIAN PROCESS

In this appendix, we visualize hyperparameter dependences of the Gaussian process and describe a hyperparameter tuning performed for this study.
As described in Section <ref>, the Gaussian process expresses the probability density function of the rigidity spectrum θ by a multivariate normal distribution 𝒩(θ | Θ_G, Σ_G).
The covariance matrix Σ_G defines the correlation coefficients of the spectrum values across individual rigidity bins, confining the acceptable range of the smoothness of the spectrum.
The correlation coefficient is determined by the hyperparameter b in equation (<ref>).
The coefficient is reduced along with the rigidity difference between rigidity bins k and k', |q_k - q_k'| in equation (<ref>), and reaches 1/e ∼ 0.37 where |q_k - q_k'| = b [log_10( GV)].


The upper panel in Figure <ref>a displays the θ_3 - θ_4 parameter space distribution of the Gaussian process 𝒩(θ | Θ_G, Σ_G), where the hyperparameters σ_G and b in equation (<ref>) are set at those used in Section <ref>, σ_G=0.1% and b=0.4 [log_10( GV)].
Components of the mean vector Θ_G are set at a constant value θ^c = 0.15%.
From Table <ref>, the rigidity difference between rigidity bins k=3 and k'=4 is |q_k - q_k'| = 0.2 [log_10( GV)].
This results in the correlation coefficient exp(-|q_k - q_k'|^2/b^2) = 0.78, visualized by the probability density distribution along the x=y line in the upper panel of Figure <ref>a.
The lower panel in Figure <ref>a displays 10 random spectra θ's following the probability distribution of the Gaussian process.
Each vertical solid line denotes the median rigidity P_k^m [GV] in the rigidity bin k=3 or 4.
An average level of the spectra is determined by Θ_G or θ^c at 0.15%, while the deviation of spectrum values is limited from ∼0.0% to ∼0.3% by the hyperparameter σ_G.
In each line, discontinuity of the spectrum values in adjacent rigidity bins, such as θ_3 and θ_4, is confined by the correlation coefficient determined by b.

Figure <ref>b is the same as Figure <ref>a but only the hyperparameter b is set at a different value b=1.0 [log_10( GV)].
The correlation coefficients between spectrum values across rigidity bins become more significant than the case of b=0.4 in Figure <ref>a, where exp(-|q_k - q_k'|^2/b^2) = 0.96 between θ_3 and θ_4.
This high coefficient leads to smoother rigidity spectra than the case of b=0.4, or nearly linear spectra, as demonstrated by the lower panel of Figure <ref>b. 
In this manner, we can confine the smoothness of the spectrum by introducing the Gaussian process as a prior distribution in Bayesian estimation as performed in Section <ref>, without assuming any analytical function.

Hyperparameters of the Gaussian process have to be predefined when we perform Bayesian estimation of the rigidity spectrum in equation (<ref>).
The mean vector Θ_G is determined every year as follows.
In the same way as Appendix <ref> and equation (<ref>), assuming a flat spectrum, θ_1 ∼θ_2 ∼…∼θ_N ∼θ^c = const., simplifies equation (<ref>) as 

    η̃_l^TA∼ c_l^z θ^c

where c_l^z = ∑_k=1^N C_lk.
In this case, the NS anisotropy θ^c is estimated as a weighted mean of the observed value η_l^TA divided by c_l^z, as

    θ^c ∼∑_l w_l η_l^TA / c_l^z/∑_l w_l

where w_l = (c_l^z / σ_l^TA)^2.
We adopt Θ_G,1 = Θ_G,2 = … = Θ_G,N = θ^c for Θ_G using the observed value η_l^TA in each year, as a reference to the average level of the spectrum.


On the other hand, we use common values of σ_G and b for all years to ensure consistency of the precondition in this study. 
Too loose constraints of the Gaussian process, corresponding to larger σ_G and smaller b, cause an over-fitting of the spectrum to observed data.
On the other hand, too strict constraints cause an under-fitting, and the reproducibility of the observed values by the derived spectrum is degraded.
From the mean spectrum Θ of the posterior distribution derived in equation (<ref>), a coefficient of determination is defined as

    R^2 = 1 - ∑_l (η_l^TA - η̃_l^TA)^2/∑_l (η_l^TA - <η_l^TA>)^2

where η̃_l^TA is an expected value of the observable η_l^TA derived by inserting the mean spectrum Θ as θ in equation (<ref>).
An average of η_l^TA for all channel pairs l's is expressed by <η_l^TA>.
Now we focus on observed data η_l^TA's in 2015 as a sample year.
For this year's data, the posterior distribution is calculated by equation (<ref>) for all combinations of σ_G and b in a sufficiently wide range of these parameters.
Then, we derived R^2 for each combination of the hyperparameters.
Figure <ref> displays the result, R^2, as a function of b for each σ_G value.
An increase of R^2 with increasing σ_G saturates around σ_G ∼ 0.1%, indicating that a constraint as strict as σ_G ∼ 0.1% can be imposed without losing R^2.
In overall lines in Figure <ref> for σ_G ≤ 0.1%, R^2 is maximized around b = 0.4 [log_10( GV)].
Based on these inspections, we adopt σ_G=0.1% and b = 0.4 [log_10( GV)] as optimal hyperparameters of the Gaussian process for our dataset.


aasjournal

