Constraining the Thickness of the Galactic Halo through Cosmic-Ray Anisotropy Using the Spatial-Dependent-Propagation Model

Abstract

The spatial-dependent-propagation (SDP) model with a nearby source works well to reproduce the coevolving features of both cosmic-ray (CR)-nuclei spectra and anisotropy. However, it is well known that the Sun is actually deviating from the galactic disk. This will lead to a dominating anisotropy in the direction perpendicular to the galactic disk, which is discrepant with current observations. Thus, it is necessary to further investigate the effect of the solar offset on anisotropy. In this work, to the best of our knowledge, this is the first time that the combined studies of the solar offset, nuclei spectra, and anisotropy have been performed based on the SDP model. As a result, to reproduce CR spectra and anisotropy, the thickness of the inner halo $\left(Z_{IH}\right)$ needs to increase linearly with the displacement of the Sun. We also know that the PeV anisotropy could be used to estimate the value of the diffusion coefficient, thus breaking the degeneracy between the diffusion coefficient and halo thickness. Therefore, it is a good approach to constrain the halo thickness. Moreover, the anisotropy in the PeV energy region, as a new probe, might also shed new light on constraining the solar offset. It is hoped that the anisotropy of the energies from ∼TeV to $\mathrm{PeV}$ can be finely measured with the LHAASO experiment, leading to a better understanding of the thick halo.

1. Introduction

In recent years, great progress has been made on the spectral measurement of cosmic rays (CRs) with balloon-borne and space-borne experiments. The fine structure of spectral hardening of nuclei at 200 GV was observed by ATIC-2 [1], CREAM [2], and PAMELA [3]. AMS-02 also confirmed the hardening with unprecedented precision, although the spectral shapes have slight discrepancies from each other [4]. At higher energies, the DAMPE observation clearly revealed that the proton spectrum further experiences a spectral softening at ∼13.6 TeV, with the spectral index changing from ∼2.60 to ∼2.85 [5]. Hints of such spectral features were also found previously in CREAM [6] and NUCLEON measurements [7]. These observations have stimulated significant concern among theoretical physical scientists. Some solutions that have been proposed to explain the spectral features include the contribution from nearby SNRs [8], the reacceleration mechanism of old SNR sources [9,10], the combination effects from different group sources [11,12], and the spatial-dependent propagation (SDP) of CRs [13,14,15].

Although the structures of the spectra and anisotropy could be explained, one critical factor has been ignored in previous works, i.e., the Sun’s offset. Usually, the solar system is assumed to be located at the midplane of the galactic disk, and the source distribution is symmetric above and below the disk. Yet, it has long been known that the Sun is located slightly above the galactic plane (toward the north galactic pole). The inferred distance above the midplane is from several to ∼30 pc [16,17,18,19,20,21,22,23,24,25,26,27]. The offset would induce a net vertical flow outward from the galactic plane, which generates a corresponding anisotropy component. This is totally inconsistent with experimental observations. However, the estimated solar offset has a significant spread of values when employing a variety of different methods, although recent analyses have revealed smaller uncertainties. The combined study of nuclei spectra and anisotropy may shed light on this topic. In this work, further investigation aimed at the anisotropy problem is performed.

The rest of this paper is organized as follows. In Section 2, the model is described, and the calculated results compared with observations are presented in Section 3. Finally, conclusions are given in Section 4.

2. Model Description

Supernova remnants (SNRs) are recognized as significant accelerators of galactic CRs, capable of reaching energies up to a few PeV [28,29,30]. Before arriving at Earth, those CRs have traveled in the galaxy for ∼${10}^7$ years after diffusing away from the acceleration sites [31]. During the journey, the impacts due to fragmentation and radioactive decay in the interstellar medium (ISM) result in the production of secondary particles. Additionally, the electrons suffer energy loss in the interstellar radiation field (ISRF) and magnetic field. This journey can be described by the propagation equation as

$$\begin{array}{ccc}\frac{\partial \psi (\overrightarrow{r},p,t)}{\partial t}\hfill & =& \mathrm{q}(\overrightarrow{r},p,t)+\overrightarrow{\nabla }·\left(D_{xx}\overrightarrow{\nabla }\psi -{\overrightarrow{V}}_c\psi \right)\hfill \\ \hfill & +& \frac{\partial }{\partial \mathrm{p}}\left[p^2{D}_{pp}\frac{\partial }{\partial p}\frac{1}{p^2}\psi \right]-\frac{\partial }{\partial p}\left[\dot{p}\psi -\frac{p}{3}\left(\overrightarrow{\nabla }·{\overrightarrow{V}}_c\psi \right)\right]\hfill \\ \hfill & -& \frac{\psi }{{\tau }_f}-\frac{\psi }{{\tau }_r},\hfill \end{array}$$

(1)

where $\mathrm{q}(\overrightarrow{r},p,t)$ is the acceleration source, $\psi (\overrightarrow{r},p,t)$ is the density of CR particles per unit momentum p at position $\overrightarrow{r}$, ${\overrightarrow{V}}_c$ is the convection velocity, $\dot{p}\equiv dp/dt$ is the momentum loss rate, and ${\tau }_f$ and ${\tau }_r$ are the characteristic timescales for fragmentation and radioactive decay, respectively. $D_{xx}$ and $D_{pp}$ are the diffusion coefficients in coordinate and momentum space, respectively. For simplicity, the process of convection are not considered in this work.

2.1. Spatial-Dependent Diffusion

Located almost in the galactic disk, the galactic astrophysical objects would generate large irregularities in the form of turbulence and affect the properties of particle diffusion. In the SDP model, the diffusion coefficient is assumed to be anticorrelated with the CR source distribution [13,32]; the diffusive volume is divided into two regions as an inner halo (IH, $\left|z\right|\le Z_{\mathrm{IH}}=\xi Z_{\mathrm{H}}$) and an outer halo (OH, $\left|z\right|>Z_{\mathrm{IH}}$). Here, $Z_{\mathrm{H}}$ is the half-thickness of the diffusive halo, and $\xi$ is the ratio of the inner halo to the entire halo with a typical value of ∼0.1. In the IH region, the level of turbulence is expected to be high due to supernova-explosion activities; hence, the diffusion coefficient is relatively small. The OH region contains scarce active sources; thus, particles diffuse much faster. The diffusion coefficient is parameterized as [32,33]

$$D_{xx}(r,z,\mathcal{R})=D_{0}F(r,z){\beta }^{\eta }{\left(\frac{\mathcal{R}}{{\mathcal{R}}_0}\right)}^{{\delta }_{0}F(r,z)},$$

(2)

where

$$F(r,z)=\left\{\begin{array}{cc}g(r,z)+\left[1-g(r,z)\right]{\left({\displaystyle \frac{z}{Z_{\mathrm{IH}}}}\right)}^n,\hfill & \left|z\right|\le Z_{\mathrm{IH}}\hfill \\ 1,\hfill & \left|z\right|>Z_{\mathrm{IH}},\hfill \end{array}\right.$$

(3)

$$g(r,z)=\frac{N_m}{1+f(r,z)},$$

(4)

where $\beta$ is the particle’s velocity in units of light speed, $D_0$ and ${\delta }_0$ are constants representing the diffusion coefficient and its high-energy rigidity dependence in the outer halo, and $\eta =0.05$ is a phenomenological constant in order to fit the low-energy data. The $N_m$ is a constant factor reflecting the property of irregular turbulence in the halo. The constant n describes the smoothness of the parameters at the transition between the two halos. $f(r,z)$ is the source-density distribution. The numerical package DRAGON [34] was used to solve the transport equation. In this work, the diffusion–reacceleration model has been adopted. The momentum-diffusive coefficient, $D_{pp}$, correlates with $D_{xx}$ via $D_{pp}{D}_{xx}=\frac{4p^2{v}_A^2}{3\delta (4-{\delta }^2)(4-\delta )}$, where $v_A$ is the Alfvén velocity, p is the momentum, and $\delta$ is the rigidity-dependence slope of the diffusion coefficient. The injection spectrum of background sources is assumed to be a power law of rigidity with a high-energy exponential cutoff, i.e., $\mathrm{q}\left(\mathcal{R}\right)\propto {\mathcal{R}}^{-\nu }{e}^{(-\mathcal{R}/{\mathcal{R}}_{\mathrm{c}})}$. The cutoff rigidity of each element could be either Z- or A-dependent. The spatial distribution of sources takes the form of the SNR distribution [35], $f(r,z)\propto {(r/r_{\odot })}^{1.69}{e}^{[-3.33(r-r_{\odot })/r_{\odot }]}{e}^{(-|z|/z_s)}$, where $r_{\odot }=8.5$ kpc and $z_s=0.2$ kpc.

2.2. Local Source

The time-dependent propagation of CRs from the local source is obtained using the Green’s function method, assuming a spherical geometry with infinite boundary conditions. The solution is

$$\varphi (r,\mathcal{R},t)={\displaystyle \frac{q_{\mathrm{inj}}\left(\mathcal{R}\right)}{{\left(\sqrt{2\pi }\sigma \right)}^3}}exp\left(-{\displaystyle \frac{r^2}{2{\sigma }^2}}\right),$$

(5)

where $q_{\mathrm{inj}}\left(\mathcal{R}\right)\delta \left(t\right)\delta \left(\mathbf{r}\right)$ is the instantaneous injection spectrum of a point source and $\sigma (\mathcal{R},t)=\sqrt{2D\left(\mathcal{R}\right)t}$ is the effective diffusion length within time t. The diffusion coefficient $\mathrm{D}\left(\mathcal{R}\right)$ is taken as the value near to the solar system. The injection spectrum is again parameterized as a cutoff power-law form, $q_{\mathrm{inj}}\left(\mathcal{R}\right)=q_0{\mathcal{R}}^{-\alpha }{e}^{(-\mathcal{R}/{\mathcal{R}}_{\mathrm{c}}^{\prime })}$. The normalization ${\mathrm{q}}_0$ and spectral index $\alpha$ are determined through fitting to the galactic-cosmic-ray (GCR) energy spectra and anisotropy data simultaneously.

2.3. Solar Offset

In the literature, several methods have been adopted to measure the solar offset from the galactic plane, which can be roughly classified into two categories: the matter-based (${\mathrm{H}}_{\mathrm{II}}$, molecular and methanol masers) and star-based (various stars, pulsars, optical stars, Wolf–Rayet stars, cepheid variables, magnetars, and open clusters) methods. Figure 1 presents estimated offset distances, which are summarized in Table 1 of ref. [27], from which it can be seen that the solar offset spans a large range, i.e., from $5\mathrm{pc}$ to ∼30 pc, although some analyses have small uncertainties. In the following, the effect of solar offset to the anisotropy is considered in the calculation.

3. Results

There are several components that might affect CR anisotropy that might remedy the nonconformity of anisotropy because of the above-inferred solar offset, namely, the source distribution, inner-halo thickness, and outer-halo thickness. At the beginning of this work, the effect of each component on the anisotropy was studied by fixing other components and taking the solar offset into account. We found that, except for the inner-halo thickness, the other two components have little influence on the anisotropy. In the following, the model results with a different ${\mathrm{Z}}_{\mathrm{IH}}$ are mainly presented.

First, the influence of halo thickness on the CR flux density with the different solar offset is discussed. Then, the propagation parameters under various ${\mathrm{Z}}_{\mathrm{IH}}$ values are determined by fitting the ratio B/C, and the background component of proton spectra are obtained with the DRAGON package. Note that the contributions from the local source mentioned above are also included. Finally, the anisotropy amplitudes and their correlation with the vertical solar displacement from the galactic disk are given.

3.1. Effect of Halo Thickness

As shown in Figure 2, compared with the consistent radial anisotropy (slope of the CR flux), the anisotropy in the vertical direction changes significantly with an increasing vertical distance from the point $({\mathrm{R}}_0,{\mathrm{Z}}_0)$. Clearly, the thickness of the halo has a strong (little) influence on the vertical (radial) anisotropy. Therefore, a thicker halo would counteract the vertical CR stream caused by the solar offset, thereby reducing the total anisotropy.

3.2. B/C and CR Spectra

The transport parameters for the SDP model, with various values of $Z_{\mathrm{IH}}$, are determined based on the B/C ratio, as illustrated in Figure 3. This figure compares the model intensities with data points from AMS-02 [36], CALET [37], and the latest data released from DAMPE [38]. The results show that the model fluxes, for different galactic-halo thicknesses, align well with the observations. Detailed parameters can be found in

Table 1. Propagation parameters ${}^{†}$.

${\mathit{Z}}_{\mathbf{IH}}$ [kpc]	${\mathit{Z}}_{\odot }$ [pc]	${\mathit{D}}_0$[$\times {10}^{28}{\mathbf{cm}}^2{\mathbf{s}}^{-1}$]	${\mathit{\delta }}_0$	${\mathit{N}}_{\mathbf{m}}$
0.5	0	0.57	0.55	0.62
0.75	5	6.9	0.65	0.39
1.0	10	8.75	0.65	0.39
1.15	13	9.24	0.69	0.39
1.3	16	10.07	0.69	0.39
1.45	19	11.05	0.7	0.39

${}^{†}$ $\mathrm{n}$ and ${\mathrm{v}}_{\mathrm{A}}$ are adopted as 3.5 and $6\mathrm{km}{\mathrm{s}}^{-1}$, respectively.

. The quantitative chi-square values are listed in

Table 2. Chi2/ndf for each distribution.

${\mathit{Z}}_{\mathbf{IH}}$ (kpc)	B/C	Proton	Helium	All-Particle	Aniso	Combined
0.5	255.84/99	203.17/112	238.95/105	13.63/33	3.42/3	715.01/352
0.75	91.96/99	406.94/112	134.08/105	19.76/33	2.43/3	655.17/352
1.0	106.25/99	295.72/112	134.78/105	24.2/33	1.29/3	562.24/352
1.15	94.12/99	338.66/112	147.66/105	23.81/33	14.53/3	618.78/352
1.3	149.05/99	515.05/112	157.51/105	26.92/33	30.47/3	879/352
1.45	233.72/99	610.67/112	157.51/105	22.87/33	50.79/3	1069.94/352
data ${}^{†}$	A(>15GeV) +C1+D	A(>30GeV)+ C1+C2+D+K	A+C1 +C2+D+K	$\mathrm{H}\ddot{\mathrm{o}}\mathrm{randel}$	>PeV

${}^{†}$ The abbreviations A, C1, C2, D, and K correspond to AMS-02, CALET, CREAM, DAMPE and KASCADE-Grande, respectively.

.

After determining these fundamental parameters, we calculated and compared the propagated spectra of primary CR elements with experimental observations. For clarity and conciseness, we focused on presenting the proton and helium spectra, as depicted in Figure 4. We utilized measurements from AMS-02 [36,37], CALET [39], CREAM [2,40], DAMPE [5], KASCADE [41], KASCADE-Grande [42], and IceTop [43] for comparison. Dotted-dashed and dashed lines are fluxes from background and local sources, respectively, and solid lines are their sums. The results indicate that the inclusion of a source component can simultaneously account for the spectral-hardening features at ∼200 GV and the softening features at ∼10 TV. It is worth noting that there are slight variations in the local-source component under different ${\mathrm{Z}}_{\mathrm{IH}}$ scenarios to better fit the data. Furthermore, Figure 5 presents the all-particle CR-energy spectrum predicted by the SDP model, alongside the $\mathrm{H}\ddot{\mathrm{o}}\mathrm{randel}$ all-particle spectrum derived from extensive experimental data Hörandel [44]. Table 2 displays the chi-square values obtained from comparing model calculations with the observations. In the chi-square calculation process, the energy spectra of protons below 15 GeV and helium nuclei below 30 GeV are excluded due to solar modulation.

3.3. Anisotropy and Solar Offset

It has been proposed in the literature that the energy-dependent anisotropy and softening features in the energy spectra may share a common origin [45]. Figure 6 demonstrates the successful reproduction of energy-dependent CR anisotropy amplitudes in this model, encompassing all major CR elements. By considering the contribution of the local source, the anisotropy of all components can be divided into three energy regions. In the lowest-energy region (region I), the anisotropy is primarily influenced by the propagation coefficient and local sources. At this stage, the former dominates, although the latter becomes increasingly important as the energy increases. The contribution of the local source to the anisotropy reaches its maximum value at tens of TeVs and then diminishes, playing a crucial role in reproducing the dip features observed in region II. In the highest-energy range (region III), the anisotropy increases as the galactic halo thickens. In this region, the effect of the local source becomes limited, and the propagation effect once again dominates the CR anisotropy.

Owing to the limited contribution of the nearby source in the region III, the PeV anisotropy is a good approach to constrain both the halo thickness and the solar offset. As is known, the dipole anisotropy estimation is proportional to ${\mathrm{D}}_{\mathrm{xx}}$[75], while the ratio of the parameterized diffusion coefficient to the halo thickness, i.e., ${\mathrm{D}}_{\mathrm{xx}}/{\mathrm{Z}}_{\mathrm{H}}$($Z_H=Z_{\mathrm{IH}}+Z_{\mathrm{OH}}$), can be fixed by observing of the $\mathrm{B}/\mathrm{C}$ ratio. Therefore, the thickness of the halo also affects the amplitudes of the anisotropy. That is to say, CR anisotropy is not only related to the diffusion coefficient ${\mathrm{D}}_{\mathrm{xx}}$ but also to the thickness of the halo and the distance of the Sun deviating from the galactic disk. As previously mentioned, we have found through detailed analysis that only the inner halo has a significant impact on anisotropy. To reproduce reversal structures of the dipole anisotropy, the thickness of the halo is required to increase linearly with the solar offset, as shown as in the left of Figure 7. This figure also shows the values of halo thickness with the solar displacement measured with matter- and star-borne methods in Figure 1, under the relationship of halo thickness and solar offset estimated with PeV CR anisotropy.

To quantitatively assess the model’s performance, we compared the calculated amplitudes of dipole anisotropy with observed values in region III. We utilized observational data points above PeV. The chi-squared values for different values of ${\mathrm{Z}}_{\mathrm{IH}}$ (0.75, 1.0, 1.15, 1.3, and 1.45 kpc) were 2.4, 1.3, 14.5, 30.5, and 50.8, respectively, with a number of degrees of freedom equal to three. For the sake of simplicity, we temporarily disregarded the correlations between different energy spectra. Consequently, we aggregated the chi-square values of the B/C ratio, proton spectrum, helium spectrum, and all-particle spectrum, and the 4 PeV anisotropy points for the same set of model parameters, employing the additivity property of the chi-square distribution. The resulting total chi-square value is provided in the last column of Table 2. On the right side of Figure 7, the combined chi-square value is depicted as it varies with different values of halo thickness. It is evident that the reduced chi-square value reaches its minimum when the halo thickness, $Z_{\mathrm{IH}}$, is approximately 1 kpc, which corresponds to a solar offset of 11 pc. When the halo thickness exceeds 1.2 kpc, the chi-square value increases sharply. Similarly, when the halo thickness is less than 0.8 kpc, the corresponding chi-square value also increases.

4. Conclusions

Measurements of CRs have entered a precise era, thanks to the fast development of space-borne and ground-based experiments in recent years. Based on the new features of CR spectra, including the spectral hardening at ∼200 GV and softening at ∼10 TV, together with the long-time enigma of the energy-dependent evolution of the dipole anisotropy features, an SDP framework with contributions from a local CR source was established and could explain most of these new observational facts [45].

To the best of our knowledge, this work is the first time that the SDP model has been extended to simultaneously study the CR spectra, anisotropy of CRs, and the solar offset from the galactic disk. It was found that the primary CR protons, all particle spectra, and the dip structures of the CR anisotropy amplitude could still be well-reproduced after considering the actual solar location. At this moment, the thickness of the halo needs to increase linearly with the displacement of the Sun from the galactic disk because a thicker halo is required to counteract the vertical anisotropy due to the solar offset. Recent measurements, i.e., Fermi bubbles [76] and large-scale X-ray bubbles [77], also support the point of a thick halo. Moreover, owing to the limited contribution of nearby sources in the PeV energy region, the anisotropy can be used to estimate the value of $D_{xx}$, thus breaking the degeneracy of ${\mathrm{D}}_{\mathrm{xx}}/{\mathrm{Z}}_{\mathrm{H}}$, which is tuned using the B/C ratio. Therefore, it is also a good approach to constrain the halo thickness. However, with present observations of CR anisotropy, our model shows that the solar offset prefers to be around 11 pc. It is worth noting that this might be a new (model-dependent) method of constraining the solar location. It is expected that future precise measurements of the anisotropy in the energy range from ∼TeV to PeV, e.g., using LHAASO [78], could provide a more accurate determination of the thickness of the halo, and, in addition, be used to test the proposed model and constrain the Sun’s vertical location.