Abstract
Introduction:
Future crewed missions to Mars will encounter substantially elevated radiation levels compared to low Earth orbit operations. To address this challenge, we present the Space-Dependent Energetic cosmic ray Modulation using MAgnetic spectrometer (SDEMMA) model, a novel framework for modeling galactic cosmic ray (GCR) dynamics in deep-space environments.
Methods:
The model employs stochastic differential equations with outer boundary conditions derived from contemporary local interstellar spectrum models. Time-dependent diffusion and drift coefficients were optimized through Markov Chain Monte Carlo parameter fitting against 2006-2019 observational data from the space-borne magnetic spectrometers of AMS-02 and PAMELA.
Results:
SDEMMA extends GCR spectral calculations to radial positions beyond 1.0 AU, explicitly resolving radial gradients under diverse heliospheric conditions. The framework provides spatiotemporally resolved GCR spectra for charge numbers Z=1–28 at rigidities >0.2 GV, covering the inner heliosphere between Earth and Mars and currently the 2006-2019 epoch.
Discussion:
Implementation demonstrates the model's operational utility: dose equivalent rates behind 30 g/cm2 polyethylene shielding during a flux minimum range from 14-17 cSv/yr, with variance attributable to quality factor selection.
1 Introduction
While radiation exposure for astronauts in a low Earth orbit or even during a journey to the Moon can now be considered less challenging, this problem is still not well understood in the next natural step for a journey to the Mars. Due to the absence of the geomagnetic field shielding and a longer trip which last for years, the radiation level is much higher. Few in situ measurements have been performed. The Radiation Assessment Detector (RAD) on the Mars Science Laboratory (MSL), also known as the Curiosity rover, made the first-ever measurement during the transit from Earth to Mars in the 2011 Mars mission time window (). The total measured dose equivalent rate is mSv/d. Analysis also shows that behind an average shielding of , the galactic cosmic rays (GCR) made the dominant dose contribution in the energy range of MeV/n to GeV/n, while the solar energetic particles only contributed approximately . Then the dose contribution of GCR is calculated by a rescaling of the MSL/RAD results under different solar modulation conditions in . During the 2016 time window, a Liulin-MO dosimeter on the ExoMars Trace Gas Orbiter conducted a second measurement that was claimed to be consistent with the first MSL/RAD measurement (), given the modulation condition difference.
Such GCRs induced radiation dose in the transit orbit represents the most important radiation exposure during the human exploration of Mars (), since GCRs are a continuous radiation source that is hard to shield against. On the other hand, the SEP events can be shielded in a better shelter part of a spacecraft for a few hours. And once on the planet, Mars’ atmosphere provides some protection, and a sub-surface shelter gives even more protection. Therefore, this important problem deserves further study. Without expensive in situ measurements, our goal is to develop a model-based calculation for the cumulative radiation dose experienced by astronauts on the journey to Mars (; ). This calculation scheme is aimed to (a) cover various possible Mars mission time windows under different solar modulation conditions, (b) provide full radial dependence, and (c) be applied to realistic target astronaut phantoms with flexible shielding. The goal involves researches in several broad field, and we plan to achieve it in a step by step manner. In the present paper we will mainly address the more “academic” GCR related issue (a and b), and keep the discussion of the space dosimetry (c) only at an illustrative level. The more realistic “engineering” dosimetric values will be provided in a forthcoming publication.
GCR spectra are influenced by solar activity, with the solar wind and interplanetary magnetic field playing crucial roles in their modulation. This variability with time is continuously captured by numerous experiments conducted on or around the Earth. Here, we use data in time series from space-borne magnetic spectrometers, specifically the PAMELA and the AMS-02 (; ; ; ; ; ). Furthermore, the encountered GCR spectra of the astronauts also vary at different radial locations. This variability cannot be covered by most previous measurements, including the PAMELA and the AMS-02. Our approach is to explicitly expand the calculation previously limited to the 1.0 AU slice (; ; ) to other radial locations in the inner solar system. With at least the solar modulation effect, the calculated GCR spectra dataset forms a GCR model, which we have named the Space-Dependent Energetic cosmic ray Modulation using MAgnetic spectrometer (SDEMMA) model. Except that its 1.0 AU spectra has been integrated with the ICRP123 fluence-to-dose conversion coefficient () for the unshielded astronaut dose rate (), this new GCR model is comparable to the Badhwar-O’Neill series of models (; ), the HelMod model (), etc. It can be used for general purposes beyond space dosimetry, and will be long-term supported.
As an example for its dosimetric application, we use this model to calculate the astronaut radiation dose rate between the Earth’s and the Mars’ orbit. While the SDEMMA GCR model provides the number of the incident particles for each species at each kinetic energy, the dose calculation needs another factor, which is the expected dose equivalent caused by a single incident particle (the fluence-to-dose-equivalent conversion coefficient). We have used several sets of fluence-to-dose-equivalent conversion coefficient, including the ones calculated by ourselves using the particle physics toolkit GEANT4 (; ; ). The dose equivalent rates are obtained in time series upon integration with the two factors.
The paper is organized as follows. In Section 2 we review the GCR spectra calculation method, including the stochastic differential equation (SDE) approach and the local interstellar spectra as the outer boundary conditions. We pay special attention to the heliospheric environment modeling, as well as the data from the PAMELA and the AMS-02 experiments and the fitting procedure to determine the diffusion and drift coefficients. In Section 3 we present a detailed discussion of the GCR flux, which depends on rigidities, GCR species, time, and radial positions. These two sections actually define the SDEMMA model. Then we introduce the fluence-to-dose-equivalent conversion coefficient and calculate the dose equivalent rate in Section 4, using three shielding settings: the unshielded case for uncertainty demonstration, the MSL/RAD shielding case for validation, and the optimized shielding thickness for the reference values. Finally, we summarize in Section 5.
2 The GCR spectra calculation scheme
The spectra of GCRs are related to their phase space density through , where is the position, is the momentum with as its magnitude, and is the phase space density. The evolution of phase space density is governed by the Boltzmann equation, which states that the total derivative of the phase space density is determined by collision terms. In the environment of the heliosphere, it becomes the Parker’s transport equation, which readsHere, and are the background solar wind and pitch angle-averaged drift velocities, and in the first term of the right hand side they give the convection and the drift effects, respectively. The middle term represents the adiabatic cooling effect caused by expansion of the solar wind. The diffusion coefficient tensor, , is a result of the small-scale turbulence of the heliospheric magnetic field (HMF).
As an application of the Feynman-Kac formula treatment of the stochastic diffusion process, the Parker’s transport Equation 1 can be reformulated into an equivalent set of 3D SDEs for the GCR phase space coordinates ()with , and
Here, three spatial dimensions have been considered. is the backward time. The differential random noises superimposed on the deterministic motion describe the Wiener diffusion process. The phase space coordinates are not for a single GCR particle, but for a macroscopically small but microscopically large phase space region which still contains a large number of GCR particles. The initial-boundary value problem for the phase space density distribution can now be solved through a Monte Carlo simulation of a Markov stochastic process for each small piece of phase space, which avoids the need for numerical solutions of the complicated partial differential equation.
The local interstellar spectra for protons from , for iron from , and for all other elements from are used as the outer boundary conditions (). They are implemented at 120 AU where the heliopause locates. The spectra of the all the GCR elements can be calculated with these boundary conditions.
2.1 Modeling the heliospheric environment
Our modeling of the heliospheric environment is identical to that of . Here, we provide a brief review of the four specific models in the Radial-Tangential-Normal coordinates: the solar wind velocity (Equation 4), the HMF (Equations 5-8), the diffusion tensor (Equations 11-13), and the drift velocity (Equations 9, 10). Simply put, there are four time-dependent parameters in the first two models ( and ) that are determined observationally: , , and . And the last two models contain five time-dependent parameters that are fitted to data using the Markov Chain Monte Carlo (MCMC) processes: in the diffusion tensor , and in the drift velocity . One can refer to Figure 1 for a glance.
FIGURE 1
2.1.1 Solar wind
The solar wind velocity up to the termination shock is given by where is the observed solar wind speed near Earth, for which we use the daily OMNI data. AU is the radius of the Sun, AU is the radial position of the Earth. is the observational tilt angle of the wavy heliospheric current sheet (HCS) from the L model, and we use the Wilcox Solar Observatory data which is given every 10 days. The first bracketed factor describes the radial acceleration behavior of the solar wind. After the acceleration, the speed remains almost constant until reaching the termination shock. The second bracketed factor describes the behavior in the polar direction. At two polar directions ( and ) and/or for a small tilt of the HCS (or during weak solar activity), the factor tends to be large and approaches 1.875. And for locations near the equator and/or for a large (or during strong solar activity), the factor tends to be small and approach 1.075 (). This depiction of the spatial and temporal distribution of the solar wind is valid up to the termination shock.
Outside the termination shock of 90 AU, the solar wind transitions from supersonic to subsonic speeds. Here, the velocity is assumed to decrease to of the value just inside the termination shock (), then follow an inverse square law as further moves out. The heliosheath is acting as a “modulation barrier”, in which our solar wind does not expand, exerting no adiabatic cooling effect on the incoming GCR particles. The enhanced modulation observed in the heliosheath is not accounted for in our calculation, and our modeling is not entirely accurate or comprehensive. However, we expect that this inconsistency will not have a significant numerical impact on our final spectra calculation, as our focus of the inner solar system region of AU is very far from the heliosheath. Moreover, there is some dependence of the inferred local interstellar spectra on the modulation process in the heliosheath, consequently the calculated GCR spectra (; ).
2.1.2 Heliospheric magnetic field
The used 3D HMF is the Parker’s spiral with a polar region enhancement ().Here, is the observed polarization of the HMF, being positive (negative) indicates that the HMF points outward (inward) in the northern hemisphere of the Sun. is the observational HMF strength near the Earth. The sources of the and are also the above Wilcox Solar Observatory data and the daily OMNI satellite data, respectively. The is the signum function. s is the angular rotation speed of the Sun. The colatitude is the separation between the north and south hemispheres caused by the undulating HCS at the phase angle . In addition to the standard Parker spiral, the term acts as a correction for the observed HMF increase at a large radial distance in the polar heliosphere. This term ensures divergence-free behavior, as introduced by . Here, the perturbation parameter in the equatorial plane is set to be (; ). To avoid a singularity, is set equal to the value when or ().
2.1.3 Diffusion and drift coefficients
In contrast to the two aforementioned models which are determined through observation and have become quite distinctive, there are numerous options in the literature for diffusion and drift coefficients. Here, we use a set of parameterized coefficients as provided by . Compared with diffusion coefficients that are more theoretically motivated (e.g., ), our parametrized coefficients usually have the advantage of reproducing observational spectra with greater precision, but at the cost of complexity in making forecast.
The drift velocity is given by withHere, is the rigidity of the particle, and for nonrelativistic case the particle momentum . and are the mass and speed of the particle, respectively, and is the particle’s charge. is the overall dimensionless drift coefficient to be fitted, which ranges from 0 to 1, with describing undisturbed drift without the competition of scattering. is the ratio of the particle’s speed to the speed of light. Note that this ratio is an implicit function of rigidity, as . GV empirically, and the corresponding factor smoothly connects the high rigidity drift region and the low rigidity scattering region, where small-scale turbulence plays a central role. is the magnitude of the HMF.
The diffusion tensor is given by its components ().Here is the overall diffusion coefficient to be fitted in units of , GV empirically. The coefficient is a smooth connection of two asymptotic power laws, with their indices of and to be fitted for the low and high rigidity region, respectively. The last parameter to be fitted is , which determines the transition to the asymptotic power law on the high rigidity side.
Finally we would like to make comment on the validity of the spectra calculation at locations different from the 1.0 AU. Since all of the above models hold for the entire spatial heliosphere region, for given observational and fitted parameter in time series, the calculation of GCR spectra is equally accurate at other locations within the heliosphere, especially at the range of AU between the orbits of Earth and Mars. In fact, the current work is extending the scope of our previous work by generalizing the 1.0 AU slice provided in ; to other radial locations. The same generalization is conceptually claimed by other similar SDE works as well (e.g., ).
2.2 The data source and fitting
As mentioned, the Alpha Magnetic Spectrometer (AMS) 02 detector is one of the most effective tools for measuring GCR spectra. It was launched in May 2011 and installed on the International Space Station at an altitude of about 400 km. The detector can measure all the GCR elements with atomic numbers ranging from 1 to 28, with a strong capability for discriminating between elements and controlling errors. The PAMELA detector is another space-borne magnetic spectrometer before the AMS-02, which was launched in June 2006 into an orbit at an altitude of km. As magnetic spectrometers, both of them utilize the Lorentz force to bend the trajectory of incident charged particles. By measuring the deflection as well as the energy loss in the detector, the charge of GCRs can be precisely determined (), which is crucial for dosimetry. On the other hand, the two spectrometers are still in the geomagnetic field, which will partially shield the incident GCR from being measured, particularly for low rigidities. But such shielding effect depends on geographic location and direction. By selecting the events above the geographic-location-and-direction dependent geomagnetic cutoff rigidity determined in the backtracing process (e.g., see ), the AMS and PAMELA group were able to discriminate and measure the primary GCR down to low rigidities in high latitude region. Both spectrometers have published time-dependent spectral measurements for the GCR proton and helium (the solar energetic particle flux has been manually removed already in the published data), which consist of approximately and of GCR flux respectively. Measurements for heavier elements as well as extended time span are underway at the AMS-02.
We determine the time-series of the aforementioned five parameters by fitting them to the time-dependent spectra measured by AMS-02 and PAMELA in the proton and helium channels (; ; ; ; ; ). In Figure 1 we have plotted the four observables (, , , and ) in the left panel and the five fitted parameters (, and ) in the right panel as time series. For each Bartels rotation, or every 27 days, fitting is performed using a Markov Chain Monte Carlo method (; ). In the relevant SDE simulation, it takes some time for each pseudoparticle to move backward and reach the heliopause. During this time, we use the time-dependent solar wind and HMF. Specifically, we implement a stepwise approach using the observed time series of , while assuming that the fitting parameters remain constant. The 6-month average of the fitted time series in Figure 1 is only used to provide a clear depiction of the trend for each quantity, and we utilize the best fit at each time point. There is an overlap period for the PAMELA and AMS-02 experiments, but in the plot, they are treated as independent fits.
The GCR spectra data are provided in the form of time series only for the proton and helium channels. Therefore, the fitting of diffusion and drift coefficients is performed only for these two species. However, was able to fit both the proton channel and the helium channel with the same parameter set. This indicates the universality of the diffusion and drift coefficients for different GCR species (; ; ; ; ). In fact, the diffusion and drift coefficients are statistical measures of the small and large scale magnetic field irregularities. For GCRs with the same rigidity, the local curvature radius of the trajectory of an incident charged particle is determined by the ratio of the rigidity to the magnetic field strength, regardless of the particle’s type or mass. It implies the same trajectory for particles with the same rigidity, explaining the universality of the drift coefficient for all GCR species. On the other hand, the diffusion coefficient is determined by the interaction between particles and waves/turbulence. According to quasi linear theory, the particle with the same gyro-radius resonates with the same waves, thus the resonant condition determines that particles with the same rigidity have the same diffusion coefficient. So here we use the best fits of the combined proton and helium data as the universal coefficients, and apply them to calculate the other 26 elements.
3 The calculated GCR spectra
The GCR flux is calculated on a four-dimensional grid, which is summarized in Table 1. Figure 2 shows sample calculations for three dimensions at 1.0 AU, namely during solar minimum and maximum as two representative dates, in each panel for five representative GCR species and the whole rigidity range.
TABLE 1
| Sampling variables | Range | Number of samplings | Sampling description |
|---|---|---|---|
| Rigidity | 0.2–100 GV | 91 | Equal spacing in logrithm |
| GCR species | to 28 | 28 | All elements |
| Date | Jun 2006 to Oct 2019 | 166 | Average of every Bartels period if measured |
| Radial Location | 1.0 to 1.6 AU | 7 | Every 0.1 AU |
The list of variables for our current SDEMMA spectra. The major part of the calculation consists of data points. The flux for each point is based on a statistics of 3,000 pseudoparticles.
FIGURE 2
3.1 The rigidity
The calculated range of rigidity is always from 0.2 to 100 GV for all elements. As a result of the MCMC fitting process, the calculated spectra accurately reproduce the measured proton and helium spectra in the rigidity region where PAMELA and AMS-02 have conducted direct measurements (; ), as a validation of the model. However, due to the finite thickness of the apparatus, low-energy/rigidity GCRs are blocked in the detector and cannot be detected, causing the measured rigidity to terminate at 1 GV for a proton and 1.65 GV for helium for the AMS-02. Beneath these minimal rigidities, one must work with the calculated spectra. As seen in Figure 2, the flux asymptotically vanishes as the rigidity decreases on the low-rigidity side (). The GCR spectra also decrease on the high rigidity side, following a power law with a well-known spectral index of approximately .
3.2 The GCR species
Although heavy ions are less abundant, their dose contribution can be enhanced by powers of their nuclear charge , due to the nature of the Bethe-Bloch stopping power as well as the biological effectiveness. As a result, their dose contribution can be comparable to that of protons or helium. In the 1.0 AU calculation of , it is noted that the five elements shown in Figure 2 together with magnesium are the most significant contributors to radiation dose, accounting for over of the unshielded astronaut effective dose equivalent. These elements also represent all fractional contributions greater than from a single element. Due to the smallness of the radial gradient, as demonstrated below, we expect that this behavior will persist across the entire spatial range.
The AMS-02 measurements for heavy elements extend to silicon and iron channels, but they are now in the form of averages over many years, and results with even shorter period have not been published yet. On the other hand, the ACE/CRIS experiment provides time-dependent measurement for heavy elements, but the detector type is calorimeter rather than magnetic spectrometer. Since we define our model using only the time-dependent data from space-borne magnetic spectrometers, the ACE/CRIS data for heavy elements are used only for cross check, but not for fitting. We respect the aforementioned universality of the diffusion and drift coefficients in calculations of other heavy elements for theoretical consistency. In we have compared our calculated spectra to the measurements averaged over many years for the elements with available data, and found good agreement.
3.3 The time dependence
Solar activity is known for its stochastic fluctuations in addition to its well-known 11-year cycle, which consists of a solar maximum and a solar minimum in each cycle. Many issues related to the time dependence have already been addressed in the previous discussion of Section 2.2. One can also refer to for more information about the GCR spectra prior to June 2017, and a comparison in the form of dose rate time series in for the 1.0 AU slice result.
Starting from June 2006 of the beginning of the PAMELA experiment (; ; ), the original fitting up to June 2017 in is further extended to October 2019, in order to match the latest available AMS-02 time-dependent proton and helium data (; ; ). New data after October 2019 will be implemented as soon as they become available.
Currently the model does not include forecast for future GCR spectra, which is crucial for future spacecraft design and space mission planing. For future solar modulation, all the GCR spectra forecast should be based on the forecast of future heliospheric environment parameters. A new GCR spectra forecast approach has been developed with the machine learning technique in , which uses the same previous GCR spectra data but not the Parker’s transport equation. A more traditional approach of forecast based on solving the Parker’s transport equation, more similar to the current forecast of Badhwar-O’Neil model and HelMod model, is also under development.
3.4 The radial dependence
In the radial gradient was measured to be /AU in the MeV/n helium channel and the MeV/n carbon channel during 1997–2006. The two above channels both correspond to a bin of central rigidity of 1.122 GV in our calculated spectra. The fact that higher GCR flux is observed at an outer position can be simply explained by the scattering of the GCR particles in the heliosphere, that part of the GCRs can be scattered back into the interstellar medium as they enter the heliosphere. As they enter deeper, more GCRs will be scattered back, and the flux will be lower.
In Figure 3 we present several flux ratios as functions of the radial position. We can observe a distinct positive radial gradient, at least for the two solar maximum configurations. The radial gradient for the astronaut’s radiation dose should be closer to the radial gradient of the total integrated flux, rather than that of the representative low rigidity bin of 1.122 GV. In our calculation, we can clearly distinguish the radial dependence between a solar minimum and a solar maximum, with that at a solar maximum noticeably larger. The measurement made by over a period of 9 years can be considered as a weighted average, yeilding a radial gradient between our solar minimum case and the maximum case of the corresponding 1.122 GV bin.
FIGURE 3
In Figure 4 we plot the flux ratios as a function of rigidity for two extreme epochs. The previous Figure 3 has already shown that the radial gradient is greater at low rigidities and smaller at high rigidities. Such phenomena are being continuously demonstrated now. Here we do not directly use the calculation results described in Table 1 with the default statistics of 3,000 pseudoparticles each, but increase the statistics to each to reduce the statistical noise. In summary, as we have demonstrated, the radial gradient consistently remains small for any GCR element at any given time. Ignoring this gradient may result in a maximum error of only a few percents.
FIGURE 4
Note that there are various techniques for solving Parker’s transport equation. Beyond the SDE method, there is another class of method, which is the alternating directional implicit (ADI) method (). Different from the SDE method which introduces the stochastic motion for a traced pseudoparticle and solves the problem in a Monte Carlo way, the ADI method focus on directly solving the partial differential equation in a deterministic way. In the latter approach a static equilibrium phase space density flow configuration is solved, then information across the whole spatial region is immediately extracted. On the other hand, in the SDE method the calculation of spectrum can only be done spatial point by spatial point, so if a radial gradient is the primary goal the SDE calculation is relatively less effective. But since the sampled pseudoparticle is simulated in real time, as mentioned before we can adopt a dynamical heliospheric environment with the observed time series. This is an advantage over the ADI method, the assumed static heliospheric environment of which is indeed inconsistent with its nature of time variability.
3.5 Comparison with other models
In Table 2 we present a simple comparison of our SDEMMA model with other recent GCR models, including the latest version of the Badhwar-O’Neill models (; ) and the HelMod model (). All versions of NASA’s Badhwar-O’Neill models are based on a one-dimensional (radial) solution to Parker’s transport equation. Only one time-dependent variable, the modulation potential , is developed to control the modulation effect, but for different element it is fitted to different values. In the latest 2020 version, the data set used for fitting is the ACE/CRIS data, when available, and the SSN is used otherwise. On the other hand, similar to ours, the HelMod model is also based on the SDE solution of Parker’s full (two-dimensional) transport equation. The SSN and CR neutron monitor data are used to fit the modulation. For both the Badhwar-O’Neil 2020 model and the HelMod model, the AMS-02 and PAMELA time-dependent spectra are used for calibration at a level of average over many years, but not for exact direct fitting of the solar modulation time series. It is natural that our GCR model provides better agreement with the AMS-02 and PAMELA data as the baseline, with an average absolute relative difference () of less than one percent in the proton and helium channels. In comparison, the Badhwar-O’Neill 2020 model and the HelMod model show an average absolute relative difference of about .
TABLE 2
| Badhwar-O’Neil (2020) | HelMod | SDEMMA | |
|---|---|---|---|
| Open Access | No | Yes | Yes |
| Transport Equation | 1D, analytical solution | 2D, SDE | 3D, SDE |
| Fitted Variable | and | ||
| Fitted Data | ACE/CRIS, SSN | SSN, NMCR | AMS-02, PAMELA |
| Fitted Way | Time series as a whole | Time series as a whole | Point by point in time series |
| proton helium | |||
| Forecast | Yes | Yes | In development |
Comparison of the Badhwar-O’Neil 2020 model, the HelMod model and our SDEMMA model.
A dedicated and more comprehensive comparison of GCR models is given in , which includes our model, the BON2020 model, the HelMod model, as well as the CREME model (; ) and DLR model ().
4 Induced dose equivalent rate between Earth’s and Mars’ orbit
Based on the SDEMMA model, we calculate the astronaut dose equivalent rate induced by the isotropic GCR flux between the Earth’s and the Mars’ orbit. It is given by
By cancelling the factors such as , and , the first equation reduces eventually to as the definition of dose equivalent rate with all particles’ contribution. And in the second equation, the two factors reduce respectively to the GCR spectra and the fluence-to-dose-equivalent conversion coefficients . Here is number of incident particle of radiation “R”, is a cross section of the incident particle beam, is the solid angle for the incident flux. is the energy loss of radiation “R”, is the corresponding target mass, so is the energy deposited in unit mass, namely the absorbed dose. is the quality factor which converts the physical absorbed dose to the medical dose equivalent, which has two definitions: the ICRP60 definition () and the NASA definition (). is the fluence for radiation “R”.
The fluence-to-dose-equivalent conversion coefficient has the meaning of the expected dose equivalent caused by a single incident particle, if the incident particle is integrated over all contributing area on the normal plane of its “beam” direction. It is calculated by particle physics Monte Carlo code. We have calculated the isotropic fluence-to-dose-equivalent conversion coefficients using the Monte Carlo toolkit GEANT4 (; ; ) using the ICRU sphere1 as the dose counter, for all the GCR elements and energy points ranging from 1 MeV/n to 100 GeV/n for each element. Except for using the simplified target of the ICRU sphere rather than the detailed human voxel phantom (), the calculation procedures are the same as the series of independently developed dose coefficient sets by ourselves (), such as the linear energy transfer and quality factor calculation methods. Compared with the only published dose coefficient set of ICRP123 (2013) for the same elements and energy range and the detailed human voxel phantom, we have the freedom to considered characteristic shielding configurations as an essential improvement of the naked human phantom assumption. This is achieved by setting an extra shielding shell around the target dose counter in the world construction of the simulation, where the shell has a certain dimension, shape, and chemical composition. Then the simulated incident particle will interact with the shielding structure before it hits the dose counter. The shielded dose coefficients are defined as the expected (average) dose equivalent caused by a single incident particle before entering any shielding.
4.1 Unshielded case
The first calculation set is for the unshielded case, and the results using the ICRP60 quality factor are collected in Figure 5. We have tested several combinations of different GCR spectra models (our SDEMMA model and the alternative HelMod model ()) with different fluence-to-dose-equivalent conversion coefficients (the ICRP123 dose coefficient (ICRP123, 2013), our independently calculated dose coefficient set using the ICRP110 human voxel phantom (), and our independently calculated dose coefficient set using the ICRU sphere). Since the SDEMMA model provides GCR spectra at seven different radial locations at per 0.1 AU, the range of the seven values form a small vertical bar for each time point. Moreover, based on the ICRP110 human voxel phantom, the first three calculations yield dose equivalent rates for each of the 15 sensitive organs/tisse “T”, and the shown is the weighted sum of the s by the tissue weight factor () (the effective dose equivalent). We can see that the dose equivalent rate differences between the HelMod and the SDEMMA GCR models are smaller than those between different dose coefficient sets of . for detailed human voxel phantom are lower than those for the ICRU sphere.
FIGURE 5
4.2 MSL/RAD shielding
The MSL/RAD shielding during the cruise stage in the transit orbit is too complicated to simulate exactly. In simulation we simplify the shielding to three aluminum layers: of the detector acceptance is shielded by a mass thickness of 1 g/cm2, is shielded by a mass thickness of 9 g/cm2, and the remainder is shielded by a mass thickness of 56 g/cm2. The average is 16 g/cm2, the same as the averaged MSL/RAD shielding thickness. On the other hand, beneath the shielding the mass thickness of the detector is small, giving negligible self-shielding. Therefore, the target dose counter has been manually divided into the surface part (with 2 mm thickness) and the remainder part. The surface part emulates the detector with negligible self-shielding, and the remainder should give values closer to the effective dose equivalent of real human phantom.
We have performed two round of simple simulations. The first one also using the same isotropic GCR incidence and the ICRU sphere gives a dose equivalent rate of 2.07 mSv/d averaged over the MSL/RAD time window. While the isotropy is the case for the astronaut in deep space, the MSL/RAD detector has a small field of view of (). While an energetic particle can penetrate the dose counter and induce dose on the other side opposite to its incidence, this cannot be measured by the MSL/RAD. The second simulation uses a plane beam and target, and gives a dose equivalent rate of 1.20 mSv/d if the flux is scaled to solid angle. The MSL/RAD measured dose equivalent rate of mSv/d purely from the GCR contribution () is between the two. Note that the remainder part of the dose counter in the two simulations both give smaller dose-equivalent rates , due to the self-shielding.
4.3 Optimized 30 g/cm2 shielding
The last calculation set is for shieldings with mass thickness of 30 g/cm2. In addition to the Bethe-Bloch stopping effect which slows down the incident particle and decreases the dose, shielding materials also produce secondary particles when hit by energetic incident particle, and those secondary particles can contribute dose when hitting the astronaut as well. For light proton and helium the secondary particles will eventually increase the dose. So there is an optimized depth which balances the production of secondary particles (increase with shielding depth) and the Bethe-Bloch stopping, and the depth is determined to be 30 g/cm2 ().
Figure 6 shows our final calculation with the ICRU sphere and optimized shielding mass thickness. As for shielding material, aluminum is currently widely used in space as structural material, but polyethylene is an optimized baseline choice for its high electron-number-to-mass ratio which facilitates the Bethe-Bloch stopping effect (). We have also varied the quality factor definition. As seen, at the flux minimum in the covered period, the dose equivalent rates induced by GCR behind optimized polyethylene shielding are cSv/yr using the ICRP60 quality factor and cSv/yr using the NASA quality factor, while at the flux maximum the rates are cSv/yr and cSv/yr, respectively.
FIGURE 6
The fluence-to-dose-equivalent conversion coefficient and dose equivalent rate calculation using the detailed human voxel phantom will be provided in a forthcoming publication.
5 Summary
We have presented a model for the modulated galactic cosmic ray spectra called SDEMMA. This model is based on data from the space-borne magnetic spectrometers PAMELA and AMS-02. In this model, we use the 3D stochastic differential equation method (Equations 2, 3) to calculate spectra, incorporating the recently developed local interstellar spectra of galactic cosmic rays for all the elements. Heliospheric environment modelings based on observational inputs and MCMC fittings are used. This model enables us to precisely reproduce the time-dependent measurements of PAMELA and AMS-02 (; ; ; ; ; ) with a margin of error of just a few percent, regardless of energy and GCR species. We extend the spectra calculation to other radial locations in the inner solar system beyond the previous focus at 1.0 AU, which is relevant to the Mars mission.
We have also developed a set of fluence-to-dose-equivalent conversion coefficients based on simplified dose counter but several shielding considerations. Combining the SDEMMA GCR model and the dose coefficients, the astronaut radiation dose equivalent rates on the transfer orbit are calculated using Equation 14, for the considered period with an explicit radial dependence.
Statements
Data availability statement
The SDEMMA dataset for GCR spectra time series can be downloaded as a zip file from https://en.iat.cn/resource.
The OMNI solar wind and HMF data are available from the GSFC website https://omniweb.gsfc.nasa.gov. The Wilcox Solar Observatory HCS data is available from the WSO website http://wso.stanford.edu/. The sun spot number data is available from the Solar Influences Data analysis Center website https://www.sidc.be.
The PAMELA proton and helium spectra time series can be extracted from https://tools.ssdc.asi.it/CosmicRays/chargedCosmicRays.jsp. The AMS-02 proton and helium spectra tables in time series are available via the corresponding supplemental material and data links of Phys. Rev. Lett. 121, 051101 (2018), Phys. Rep. 894, 1 (2021), and Phys. Rev. Lett. 128, 231102 (2022) at the AMS group publication page https://ams02.space/publications.
The HelMod GCR spectra time series are available from https://www.helmod.org/index. php?view=article&id=76:transfer-orbit-fluence&catid=14. The fluence-to-dose-equivalent conversion coefficients for the ICRU sphere can be downloaded as a zip file from https://en.iat.cn/resource.
The ICRP123 fluence-to-dose-equivalent conversion coefficients are available via the “Supplemental Material” link at the ICRP123 publication URL https://www.icrp.org/publication.asp?id=ICRP%20Publication%20123.
Author contributions
XS: Data curation, Formal Analysis, Investigation, Methodology, Resources, Software, Validation, Writing–review and editing. RH: Conceptualization, Funding acquisition, Project administration, Supervision, Visualization, Writing–original draft, Writing–review and editing. SX: Data curation, Writing–review and editing. XC: Data curation, Writing–review and editing. XL: Funding acquisition, Methodology, Software, Writing–review and editing.
Funding
The author(s) declare that financial support was received for the research and/or publication of this article. The presented work was supported by the Shandong Institute of Advanced Technology start funding (2020106R01, 2020106R02) and the NSFC grants (U2106201).
Acknowledgments
We are grateful to the useful discussion with Jingnan Guo, Tong Su, Weiwei Xu and Vladimir Mikhailov. The generative AI of Wordvice AI and DeepSeek-R1 have been used for English improvement.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Footnotes
1.^The ICRU sphere is a phantom used in radiation protection. It has a diameter of 30 cm, and a hypothetical tissue equivalent material of density 1 g/cm3 and composition of oxygen , carbon , hydrogen and nitrogen .
References
1
AdamsJ. H.BarghoutyA. F.MendenhallM. H.ReedR. A.SierawskiB. D.WarrenK. M.et al (2012). Creme: the 2011 revision of the cosmic ray effects on micro-electronics code. IEEE Trans. Nucl. Sci.59, 3141–3147. 10.1109/TNS.2012.2218831
2
AdrianiO.BarbarinoG. C.BazilevskayaG. A.BellottiR.BoezioM.BogomolovE. A.et al (2013). Time dependence of the proton flux measured by pamela during the 2006 july-2009 december solar minimum. Astrophys. J.765, 91. 10.1088/0004-637x/765/2/91
3
AgostinelliS.AllisonJ.AmakoK.ApostolakisJ.AraujoH.ArceP.et al (2003). Geant4—a simulation toolkit. Nucl. Instrum. Methods Phys. Res. Sect. A Accel. Spectrom. Detect. Assoc. Equip.506, 250–303. 10.1016/S0168-9002(03)01368-8
4
AguilarM.AisaD.AlpatB.AlvinoA.AmbrosiG.AndeenK.et al (2015). Precision measurement of the proton flux in primary cosmic rays from rigidity 1 GV to 1.8 TV with the alpha magnetic spectrometer on the international space station. Phys. Rev. Lett.114, 171103. 10.1103/PhysRevLett.114.171103
5
AguilarM.Ali CavasonzaL.AlpatB.AmbrosiG.ArrudaL.AttigN.et al (2018). Observation of fine time structures in the cosmic proton and helium fluxes with the alpha magnetic spectrometer on the international space station. Phys. Rev. Lett.121, 051101. 10.1103/PhysRevLett.121.051101
6
AguilarM.Ali CavasonzaL.AmbrosiG.ArrudaL.AttigN.BaraoF.et al (2021). The alpha magnetic spectrometer (AMS) on the international space station: Part II - results from the first seven years. Phys. Rep.894, 1–116. 10.1016/j.physrep.2020.09.003
7
AguilarM.Ali CavasonzaL.AmbrosiG.ArrudaL.AttigN.BaraoF.et al (2022). Properties of daily helium fluxes. Phys. Rev. Lett.128, 231102. 10.1103/PhysRevLett.128.231102
8
AllisonJ.AmakoK.ApostolakisJ.AraujoH.Arce DuboisP.AsaiM.et al (2006). Geant4 developments and applications. IEEE Trans. Nucl. Sci.53, 270–278. 10.1109/TNS.2006.869826
9
AllisonJ.AmakoK.ApostolakisJ.ArceP.AsaiM.AsoT.et al (2016). Recent developments in Geant4. Nucl. Instrum. Methods Phys. Res. Sect. A Accel. Spectrom. Detect. Assoc. Equip.835, 186–225. 10.1016/j.nima.2016.06.125
10
BoschiniM.Della TorreS.GervasiM.La VaccaG.RancoitaP. (2022). The transport of galactic cosmic rays in heliosphere: the HelMod model compared with other commonly employed solar modulation models. Adv. Space Res.854, 2636–2648. 10.1016/j.asr.2022.03.026
11
BoschiniM. J.Della TorreS.GervasiM.GrandiD.JóhannessonG.La VaccaG.et al (2021). The discovery of a low-energy excess in cosmic-ray iron: evidence of the past supernova activity in the local bubble. Astrophys. J.913, 5. 10.3847/1538-4357/abf11c
12
BoschiniM. J.Della TorreS.GervasiM.La VaccaG.RancoitaP. G. (2016). Propagation of cosmic rays in heliosphere: the HELMOD model. Adv. Space Res.207, 2859–2879. 10.1016/j.asr.2017.04.017
13
BoschiniM. J.TorreS. D.GervasiM.GrandiD.JóhannessonG.VaccaG. L.et al (2020). Inference of the local interstellar spectra of cosmic-ray nuclei Z ≤ 28 with the GalProp–HelMod framework. Astrophys. J. Suppl.250, 27. 10.3847/1538-4365/aba901
14
Badhwar-O’Neill(2014). Galactic cosmic ray flux model description. NASA Technical Reports-2015-218569
15
ChenL.ChenX.HuoR.XuS.XuW. (2025). Astronaut dose coefficients calculated using GEANT4 and comparison with icrp123. Radiat. Environ. BiophysicsXX, XXX.
16
ChenX.XuS.SongX.HuoR.LuoX. (2023). Astronaut radiation dose calculation with a new galactic cosmic ray model and the AMS-02 data. Space weather.21, e2022SW003285. 10.1029/2022SW003285
17
CortiC.PotgieterM. S.BindiV.ConsolandiC.LightC.PalermoM.et al (2019). Numerical modeling of galactic cosmic-ray proton and helium observed by ams-02 during the solar maximum of solar cycle 24. Astrophysical J.871, 253. 10.3847/1538-4357/aafac4
18
CucinottaF. A.KimM.-H. Y.ChappellL. J. (2011). Space radiation cancer risk projections and uncertainties - 2010. NASA Tech. Reports-2011-216155.
19
DuY.-L.SongX.LuoX. (2025). Deep learning the forecast of galactic cosmic-ray spectra. Astrophysical J. Lett.978, L36. 10.3847/2041-8213/ada427
20
FiandriniE.TomassettiN.BertucciB.DonniniF.GrazianiM.KhialiB.et al (2021). Numerical modeling of cosmic rays in the heliosphere: analysis of proton data from ams-02 and pamela. Phys. Rev. D.104, 023012. 10.1103/PhysRevD.104.023012
21
GieselerJ.HeberB.DunzlaffP.Müller-MellinR.KlassenA.Gomez-HerreroR. (2008). “The radial gradient of galactic cosmic rays: ulysses KET and ACE CRIS Measurements,” in International cosmic ray conference. Vol. 1 of international cosmic ray conference, 571–574.1. Available online at https://ui.adsabs.harvard.edu/abs/2008ICRC
22
GuoJ.SlabaT. C.ZeitlinC.Wimmer-SchweingruberR. F.BadaviF. F.BöhmE.et al (2017). Dependence of the martian radiation environment on atmospheric depth: modeling and measurement. J. Geophys. Res. Planets122, 329–341. 10.1002/2016JE005206
23
GuoJ.WangB.WhitmanK.PlainakiC.ZhaoL.BainH. M.et al (2024). Particle radiation environment in the heliosphere: status, limitations, and recommendations. Adv. Space Res.10.1016/j.asr.2024.03.070
24
GuoJ.ZeitlinC.Wimmer-SchweingruberR. F.HasslerD. M.PosnerA.HeberB.et al (2015). Variations of dose rate observed by MSL/RAD in transit to Mars. Astronomy and Astrophysics577, A58. 10.1051/0004-6361/201525680
25
HelMod (2024). Available online at: https://www.helmod.org/index.php?view=article&id=76:transfer-orbit-fluence&catid=14.
26
ICRP103 (2007). 2007 recommendations of the international commission on radiological protection. ICRP publication 103. Ann. ICRP37.
27
ICRP110 (2009). Adult reference computational phantoms. ICRP publication 110. Ann. ICRP39.
28
ICRP123DietzeG.BartlettD. T.CoolD. A.CucinottaF. A.JiaX.et al (2013). ICRP, 123. Assessment of radiation exposure of astronauts in space. ICRP Publication 123. Ann. ICRP42, 1–339. 10.1016/j.icrp.2013.05.004
29
ICRP60 (1991). 1990 recommendations of the international commission on radiological protection. ICRP publication 60. Ann. ICRP21.
30
JokipiiJ. R.KotaJ. (1989). The polar heliospheric magnetic field. Geophys. Res. Lett.16, 1–4. 10.1029/GL016i001p00001
31
LangnerU. W.PotgieterM. S. (2004). Solar wind termination shock and heliosheath effects on the modulation of protons and antiprotons. J. Geophys. Res. Space Phys.109. 10.1029/2003JA010158
32
LangnerU. W.PotgieterM. S.WebberW. R. (2003). Modulation of cosmic ray protons in the heliosheath. J. Geophys. Res. Space Phys.108. 10.1029/2003JA009934
33
LiH.WangC.RichardsonJ. D. (2008). Properties of the termination shock observed by voyager 2. Geophys. Res. Lett.35. 10.1029/2008GL034869
34
LiuW.GuoJ.WangY.SlabaT. C. (2024). A comprehensive comparison of various galactic cosmic-ray models to the state-of-the-art particle and radiation measurements. Astrophysical J. Suppl. Ser.271, 18. 10.3847/1538-4365/ad18ad
35
MarcelliN.BoezioM.LenniA.MennW.MuniniR.AslamO. P. M.et al (2020). Time dependence of the flux of helium nuclei in cosmic rays measured by the pamela experiment between 2006 july and 2009 december. Astrophys. J.893, 145. 10.3847/1538-4357/ab80c2
36
MartucciM.MuniniR.BoezioM.FeliceV. D.AdrianiO.BarbarinoG. C.et al (2018). Proton fluxes measured by the PAMELA experiment from the minimum to the maximum solar activity for solar cycle 24. Astrophys. J.854, L2. 10.3847/2041-8213/aaa9b2
37
MatthiäD.BergerT.MrigakshiA. I.ReitzG. (2013). A ready-to-use galactic cosmic ray model. Adv. Space Res.51, 328. 10.1016/j.asr.2012.09.022
38
MoraalH.PotgieterM. S. (1982). Solutions of the spherically-symmetric cosmic-ray transport equation in interplanetary space. Astrophysics Space Sci.84, 519–533. 10.1007/BF00651330
39
NaitoM.KodairaS.OgawaraR.TobitaK.SomeyaY.KusumotoT.et al (2020). Investigation of shielding material properties for effective space radiation protection. Life Sci. Space Res.26, 69–76. 10.1016/j.lssr.2020.05.001
40
NgobeniM. D.AslamO. P. M.BisschoffD.PotgieterM. S.NdiitwaniD. C.BoezioM.et al (2020). The 3d numerical modeling of the solar modulation of galactic protons and helium nuclei related to observations by pamela between 2006 and 2009. Astrophysics Space Sci.365, 182. 10.1007/s10509-020-03896-1
41
O’NeillP.GolgeS.SlabaT. C. (2015). Galactic cosmic ray flux model description. NASA Technical Reports-2015-218569
42
OpherM.RichardsonJ.ZankG.FlorinskiV.GiacaloneJ.SokolJ. M.et al (2023). Solar wind with hydrogen ion charge exchange and large-scale dynamics (shield) drive science center. Front. Astronomy Space Sci.10. 10.3389/fspas.2023.1143909
43
PotgieterM. S. (2013). Solar modulation of cosmic rays. Living Rev. Sol. Phys.10, 3. 10.12942/lrsp-2013-3
44
PotgieterM. S.MoraalH. (1985). A drift model for the modulation of galactic cosmic rays. Astrophys. J.294, 425–440. 10.1086/163309
45
PotgieterM. S.VosE. E.BoezioM.De SimoneN.Di FeliceV.FormatoV. (2014). Modulation of galactic protons in the heliosphere during the unusual solar minimum of 2006 to 2009. Sol. Phys.289, 391–406. 10.1007/s11207-013-0324-6
46
QinG.ShenZ.-N. (2017). Modulation of galactic cosmic rays in the inner heliosphere, comparing with pamela measurements. Astrophysical J.846, 56. 10.3847/1538-4357/aa83ad
47
QinG.ZhangL.-H. (2014). The modification of the nonlinear guiding center theory. Astrophysical J.787, 12. 10.1088/0004-637X/787/1/12
48
RaathJ. L.PotgieterM. S.StraussR. D.KoppA. (2016). The effects of magnetic field modifications on the solar modulation of cosmic rays with a SDE-based model. Adv. Space Res.57, 1965–1977. 10.1016/j.asr.2016.01.017
49
SemkovaJ.KolevaR.BenghinV.DachevT.MatviichukY.TomovB.et al (2018). Charged particles radiation measurements with liulin-mo dosimeter of frend instrument aboard exomars trace gas orbiter during the transit and in high elliptic mars orbit. Icarus303, 53–66. 10.1016/j.icarus.2017.12.034
50
SlabaT. C.WhitmanK. (2019). The badhwar - O’neill 2020 model. NASA Technical Reports-2019-220419
51
SongX.LuoX.PotgieterM. S.LiuX.GengZ. (2021). A numerical study of the solar modulation of galactic protons and helium from 2006 to 2017. Astrophys. J. Suppl.257, 48. 10.3847/1538-4365/ac281c
52
TomassettiN.BarãoF.BertucciB.FiandriniE.FigueiredoJ.LousadaJ.et al (2018). Testing diffusion of cosmic rays in the heliosphere with proton and helium data from ams. Phys. Rev. Lett.121, 251104. 10.1103/PhysRevLett.121.251104
53
TylkaA.AdamsJ.BobergP.BrownsteinB.DietrichW.FlueckigerE.et al (1997). Creme96: a revision of the cosmic ray effects on micro-electronics code. IEEE Trans. Nucl. Sci.44, 2150–2160. 10.1109/23.659030
54
WangB.-B.BiX.-J.FangK.LinS.-J.YinP.-F. (2019). Time-dependent solar modulation of cosmic rays from solar minimum to solar maximum. Phys. Rev. D.100, 063006. 10.1103/PhysRevD.100.063006
55
ZeitlinC.HasslerD. M.CucinottaF. A.EhresmannB.Wimmer-SchweingruberR. F.BrinzaD. E.et al (2013). Measurements of energetic particle radiation in transit to mars on the mars science laboratory. Science340, 1080–1084. 10.1126/science.1235989
56
ZhangM. (1999). A markov stochastic process theory of cosmic-ray modulation. Astrophys. J.513, 409–420. 10.1086/306857
Summary
Keywords
galactic cosmic ray, solar modulation, stochastic differential equation, spatial dependence, fluence-to-dose conversion coefficient, dose equivalent rate
Citation
Song X, Huo R, Xu S, Chen X and Luo X (2025) The SDEMMA model for galactic cosmic ray and its dosimetric application. Front. Astron. Space Sci. 12:1383946. doi: 10.3389/fspas.2025.1383946
Received
08 February 2024
Accepted
11 March 2025
Published
04 April 2025
Volume
12 - 2025
Edited by
Yulia Bogdanova, Rutherford Appleton Laboratory, United Kingdom
Reviewed by
Mike Hapgood, Rutherford Appleton Laboratory, United Kingdom
Bingbing Wang, University of Alabama in Huntsville, United States
Updates
Copyright
© 2025 Song, Huo, Xu, Chen and Luo.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Ran Huo, huor@iat.cn
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.