Astronaut Radiation Dose Calculation With a New Galactic Cosmic Ray Model and the AMS‐02 Data

We present a new calculation of the astronaut dose rate from the galactic cosmic rays in free space at 1 AU. We use the unshielded isotropic fluence‐to‐dose conversion coefficients given in the International Commission on Radiological Protection publication 123. A new 3D and time‐dependent solar modulation model based on Parker's transport equation as originally developed in Song et al. (2021, https://doi.org/10.3847/1538-4365/ac281c) is used to calculate the galactic cosmic ray spectra at 1 AU. This model uses the recent local interstellar spectra of Corti et al. (2019, https://doi.org/10.3847/1538-4357/aafac4), M. J. Boschini et al. (2020, https://doi.org/10.3847/1538-4365/aba901, 2021a, https://doi.org/10.3847/1538-4357/abf11c) to reproduce the PAMELA and AMS‐02 observations between 2006 and 2019. The radiation dose calculated from our model and from the AMS‐02 spectra in the same rigidity region agrees better than 1% for proton and helium in a time‐dependent way, and at 2% level for six most contributing cosmic ray elements averaged over 7 or 8.5 years. The time‐dependent dose rate analysis over 13 years shows an effective dose equivalent rate of 55–58 cSv/yr at solar minimum (January 2010) and 26 cSv/yr at solar maximum (February 2014).

Here  = ∕( ) is the rigidity, N R is the number of incident particle of radiation type "R," A is the cross section of the flux. R() is the isotropic flux measured in ( m 2 × s × sr × GV ) −1 , and integration with time is also implied to return to 2 R∕( ) and obtain dose, otherwise we will get dose rate (in this paper we use the same symbol for dose and its rate in all respective definitions). The absorbed dose D T,R = dE T,R /dm is the energy E T,R deposited per unit mass into tissue/organ "T" induced by such N R particles, Φ is the fluence, and Q is the quality factor. A summation over the human sensitive tissue/organ with the tissue weighting factor w T gives the effective dose equivalent from a specific radiation species, and a summation over the GCR radiation species gives the total effective dose equivalent The fluence-to-dose conversion coefficient (D T,R Q)/(dN R /dA) with dimension Sv × m 2 has the physical meaning of the expected dose equivalent caused by a single particle of certain type and energy/rigidity to certain human tissue/organ, if the incident particle is integrated over all contributing cross section and averaged over directions assuming isotropy. It can be calculated by Monte Carlo simulations with particle physics code such as Geant4, PHITS, FLUKA, etc. This calculation involves definition or knowledge from many fields such as particle physics, atomic physics, molecular biology, human anatomy and so on. The International Commission on Radiological Protection (ICRP) has made regulating definitions in the relevant "biological" or "medical" aspect. The Adult Reference Computational Phantoms. ICRP Publication 110 (2009) 60 (1991) and a more recent NASA definition  give the radiation quality factor Qs, as measures of the relative biological effectiveness, to change the physical-based D T,R into the biological-based operational dose equivalent. The quality factor Q ICRP60 is purely based on the linear energy transfer (LET, or dE/ dx), and the Q NASA is based on the effective charge as well as the LET. Eventually all the efforts are combined in the ICRP (2013), giving the fluence-to-dose conversion coefficient for direct utilization. The standard deviation of these conversion coefficients were estimated to be 5% (Sato et al., 2009(Sato et al., , 2010.

The GCR Model
Back to the GCR spectra side, we will soon see that the AMS-02 measurements cannot be directly used as the input.
Here we use a recent 3D and time-dependent GCR numerical model (Song et al., 2021) as the spectrum generator for our dosimetric calculation, together with the local interstellar spectra (M. J. Boschini et al., 2020Boschini et al., , 2021Corti et al., 2019) as boundary conditions.
The original Parker's transport equation for the phase space density where sw and d are the background solar wind and pitch angle averaged drift velocities respectively, and ⃖⃖ ⃗ is the diffusion coefficient tensor. In principle d and ⃖⃖ ⃗ are determined by time-dependent observational heliospheric parameters only (e.g., the solar wind speed, the heliospheric current sheet tilt angle, the heliospheric magnetic field strength and its polarity), but due to the difficulty in establishing direct relations, in practice for every measurement in each Bartel rotation we use empirical parametric formulas with extra parameters fitted to the PAMELA and AMS-02 data (Song et al., 2021).

of 14
In practical calculation Equation 4 is reformulated into an equivalent stochastic differential equation (SDE) set (M. Zhang, 1999) for the GCR phase space coordinates Here s = −t is the backward time. The differential random noises ⃗ superimposed on the deterministic motion describe the Wiener diffusion process, which are given by the small scale heliospheric magnetic field irregularity.
The local interstellar spectra from Corti et al. (2019) for proton, from M. J. Boschini et al. (2021) for iron and from M. J. Boschini et al. (2020) for all the other elements are used as the boundary conditions. Compared with Song et al. (2021), the local interstellar spectrum for helium has changed from the one of Corti et al. (2019) to the one from M. J. Boschini et al. (2020). In the rigidity region  ≳ 1 GV which gives the dominant dose contribution the two have little difference, but for  ≲ 1 GV the latter is smaller, by up to 17% at 0.6 GV. Actually in this region the former is slightly above the Voyager direct measurement and the latter is slightly below the Voyager direct measurement.
Our simulation covers the original 11-year solar cycle studied in Song et al. (2021) which starts at July 2006, and extends further to October 2019 to match the most recent published AMS-02 measurements (Aguilar, Cavasonza, Allen, et al., 2021;Aguilar et al., 2021b). The simulation expands more than 13 years, or almost 180 Bartel rotations. When time-dependent measurement is available, calculation of spectra is made for a Bartel rotation with the fitted parameters given either directly by Song et al. (2021) or by similar new fitting during the year 2017-2019.

The Extension Beyond the AMS-02 Measured Rigidity Minimum
One important limitation of the AMS-02 results is that at the low rigidity side the measurement has a finite rigidity minimum (e.g., see Aguilar et al., 2015), which is mainly due to the energy loss process in the detector and should not be confused with the geomagnetic "cutoff." In Figure 1 we show the AMS-02 measured spectra as well as our calculated spectra for proton and iron. The measured spectra have minimums at 1 GV for proton and 2.65 GV for iron respectively. For such spectrum the effective dose equivalent (rate) from a specific radiation species can only be calculated above the minimum M as On the other hand, the modulation model calculated spectra do not suffer such problem, and extend down to low rigidities. We can see in Figure 1 that below the rigidity minimums the spectra deviate sharply from the power law behavior with the famous spectral index of −2.7. They develop local peaks not far below the AMS-02 minimums, then the spectral index turns to positive values, avoiding large fluence as rigidity vanishes. In fact Moraal and Potgieter (1982) has shown a spectral index of +3 or so for small rigidities. In table 1 (see below) we give the fractions of the dose contributions below the minimums in the rows of "1 − E,R,Model(M)∕ E,R,Model ", where E,R,Model(M) is given by Equation 7 and H E,R,Model is given by Equation 2, with the spectrum both given by model calculation.
Starting from May 2011, the published GCR spectra of the AMS-02 experiment cover the solar maximum, but not the solar minimum. As shown in Figure 1, the last available proton data in Aguilar et al. (2021b) is October 2019, which is before the true next solar minimum in the year 2020. We show both the calculated October 2019 spectrum for direct comparison to the AMS-02 data, as well as the solar minimum spectrum of January 2010. The total fluence of January 2010 is larger than that of October 2019 by 2.8% from model calculation. The cross of the two spectra exhibits the rigidity dependence of the solar modulation, which had been seen in the comparison of the PAMELA 2006 proton data (Adriani et al., 2013) and its 2011 proton data (Martucci et al., 2018). Therefore, throughout the paper we define the solar maximum and minimum as simply giving the minimal and maximal dose rate, at February 2014 and January 2010 respectively. 4 of 14

The Extension to the Elements Absent in AMS-02 Measurements
Currently the AMS-02 has published measurements for the Z = 1-14 and Z = 26 elements. For the elements which have not yet been measured by AMS-02, we use our model to calculate their spectra. Including all the high Z and energy GCR species, all the 28 GCR spectra are calculated at 1 AU point by point in a time series. The diffusion and drift parameter sets for each time point (Song et al., 2021) have been shown to be universal for both the proton and helium with the help of data. Here we assume that such theoretically motivated universality also holds for the other heavy elements (M. J. Boschini et al., 2021;Corti et al., 2019;Fiandrini et al., 2021;Ngobeni et al., 2020;Tomassetti et al., 2018). Further detailed measurements of the time-dependent fluxes of these heavy elements from AMS-02 and other experiments are important to verify this assumption.
Note that the contribution from electrons is not accounted in this paper, as the electron fluence-to-dose conversion coefficient is not given in ICRP123. In fact the measured GCR electron flux is a factor of ∼50 (at 1 GV) to ∼1,000 (at 1,000 GV) below that of proton (Aguilar et al., 2021a), its contribution is expected to be below ∼1%.  (Aguilar, Cavasonza, Allen, et al., 2021;Aguilar et al., 2020Aguilar et al., , 2021a  spectra have minimums at low rigidities, so we make the same cuts on the model side. However still some part of dose are from the GCR below M , so the fraction below the cut is also shown. Here the effective dose equivalent rates are averaged over genders. Bold values represent that they are more important than the surrounding values and deserve attention in a first glance.

of 14
In Figure 2 we show the unshielded gender averaged effective dose equivalent fractional contributions for all the 28 elements, respectively at the solar minimum and maximum, and respectively with the quality factor definition of ICRP60 and NASA. We can see that the top six contributing elements are always proton (p), helium (He), oxygen (O), magnesium (Mg), silicon (Si) and iron (Fe), and the six elements add to more than 74% of the total contribution. The individual contribution from the other elements is always below 5%.
In Figure 3 we show the fractional contributions for all the 28 elements to the unshielded gender averaged skin dose equivalent. As seen, the top six contributing elements to the effective dose equivalent also give leading contributions. The literature value of Durante and Cucinotta (2011) is also shown, and we reserve further discussion of the comparison in Section 3.

The Model Precision Checking
Detailed comparisons between the AMS-02 (and PAMELA) measurements and the model calculation have been performed, in the form of time series for the proton and helium channels (see Figures 5 and 6 of Song et al., 2021). As we add new data from 2017 to 2019 here, we have checked that similar fitting gives comparable precision of ∼2%. In order to explicitly show the agreement level, in Figure 4 we plot the E,p(M) and E,He(M) rates  Boschini et al., 2016Boschini et al., , 2022 for comparison. Using the calculated effective dose equivalent rates on a series of 113 time points, we find that for our model the absolute relative difference |Rd| in the definition of for example, O'Neill et al. (2015) is 0.64% for proton and 0.98% for helium, and for the HelMod model they are 2.9% and 3.8%, respectively.
Here we give some detail of our dose rate calculation, which will be used throughout this paper. Different from the model calculated spectrum, for the AMS-02 measured spectra each flux value is assigned to a finite rigidity bin interval but not to a single rigidity value. According to its underlying universal assumption of the ∝  −2.7 flux behavior even in a bin, for a bin with rigidity interval min to max the single reference rigidity is (Lafferty & Wyatt, 1995) So we set the AMS-02 measured differential flux for each bin to such reference rigidities, then the AMS-02 measured spectra take the same form as the theoretically calculated ones. In using Equation 7 the spectrum, linearly interpolated in log-log scale, is multiplied by the ICRP123 fluence-to-dose conversion coefficient, linearly interpolated in log-linear scale, then numerically integrated for the  > M region. Still for the AMS-02 measured spectra, between M (the min of the lowest bin) and the lowest ref (e.g., for proton it is between 1 GV and about 1.08 GV since the AMS-02 lowest bin is 1 − 1.16 GV) constant differential flux of the lowest bin is used. And we have checked that, for example, removing such small interval in dosimetric integration does not change the comparison between different spectra sources visibly. If there is no minimum M specified as in Equation 3 the lowest rigidity used is 0.2 GV, and the highest rigidity used is always 100 GV, as the rigidity range of our model calculation. Variation of such practical bounds does not affect the dose results visibly since the GCR spectrum is suppressed at the two ends, we have checked this point with the HelMod spectra the two bounds of which can be extended to 0.14 GV and about 1500 GV respectively.
Except for the proton and helium which have data in time series, the other four top six contributing GCR species of oxygen, magnesium, silicon and iron only have AMS-02 direct measurements available in time averages over 7 (or 8.5 for iron) years (Aguilar, Cavasonza, Allen, et al., 2021;Aguilar et al., 2020Aguilar et al., , 2021a. Therefore we compare the averaged effective dose equivalent rates from these elements in the corresponding years. Using Equation 7 we calculate the E,R,Model(M) rates point by point in the time series, and average during exactly the same 7 (or 8.5) years. Comparing them to the E,R,AMS-02(M) rates directly calculated using the AMS-02 spectra, as seen in Table 1 the agreement is better than or equal to 2%. On the other hand, we can see that for some elements the below-minimum fractions can be relatively high, such as the oxygen and magnesium. In Figure 5 we show both the GCR spectra and the ICRP123 fluence-to-dose conversion coefficients as functions of rigidity, for proton, oxygen and iron respectively. Compared with iron, oxygen does have a lower rigidity minimum from measurement, but its fluence-to-dose conversion coefficient peaks (corresponding to the Bragg peak) at an even lower rigidity which is not covered by the rigidity minimum, so its below-minimum fraction is higher. But since these more uncertain species make relatively subdominant contribution to the total dose (see Figure 2), even if the ∼2% error gets doubled and conspires in the same direction below the minimum, the resultant final total dose error is less than ∼0.4%.

Results and Discussion
The (effective) dose equivalent rates as a time series from 2006 to 2019 are presented in Figure 6. The detailed results of the (effective) dose equivalent rates from model calculation are shown in Tables 2 and 3, for solar . Comparison of the proton and helium channel effective dose equivalent rate contributions using the Q ICRP60 definition in the form of interpolated time series, above their specific M . For the AMS-02 proton and helium measurements we use the Bartel rotation averaged spectra (Aguilar et al., 2018) before May 2017, after that and till October 2019 we average the daily spectra (Aguilar et al., 2021b(Aguilar et al., , 2022 in every Bartel rotation, and if the daily spectrum misses the lowest bin then it is not counted in the time average. For the HelMod model (M. Boschini et al., 2022) the galactic cosmic ray data is taken from https://www.helmod.org/index.php/helmod-calculators?view=article&id=60:heliospheric-modulator-v2018&catid=9:online-resources. In addition to the same M s, we have also used the same upper bounds for the three spectra, which are determined by the AMS-02 proton and helium spectra as the ref s (see Equation 8) of their last bins. In the upper panel the three are directly compared, and in the lower panel the two ratios of the rates from the models normalized with those from the AMS-02 measurement are shown. The ratios using Q NASA are almost identical respectively. 8 of 14 minimum and maximum respectively. At solar minimum the effective dose equivalent rate is 55-58 cSv/yr, and at solar maximum it is 26 cSv/yr. The exact number depends on the choice of the quality factor, and on the gender to a less extent.

Comparison Between the Quality Factor Choices
While the choice of the quality factor definition does not affect the physical processes or the medical facts, the new NASA quality factor generally assigns larger effective dose equivalent than the ICRP60 quality factor. When comparing the two definitions in detail in Tables 2 and 3, Q NASA always assigns a smaller dose equivalent to the red bone-marrow (for leukemia). For the other 14 organs or tissues (for solid cancer), Q NASA assigns a larger dose equivalent for light GCR species and a smaller dose equivalent for heavy species, with the transition between oxygen and magnesium.
The relative difference between the two effective dose equivalent rate at solar minimum is about 4%, and at organ dose equivalent level the difference is even larger. This difference is larger than the uncertainty from the GCR spectra, increasing the significance for future relevant work to further clarify or improve the health risk evaluation.
Please note that throughout the paper we use the ICRP103 tissue weighting factor, even in combination with the NASA quality factor. Cucinotta et al. (2011) has also provided several different sets of tissue weighting factors, for US male or female, in combination with average population or never-smokers. The latter is "intended strictly for internal NASA usage" so we do not use, but just point out that for certain organ such as lung the w T difference can be a factor of 4 times large.
On the other hand, the aforementioned standard deviation of the fluence-to-dose conversion coefficient of 5% (Sato et al., 2009(Sato et al., , 2010 also exceeds the uncertainty of the GCR spectra.

Comparison Between the Effective Dose Equivalent and Dose Equivalent
In Figure 6 we have shown the results for both the effective dose equivalent rate which is the w T weighted sum of the dose equivalent rates for all the 15 sensitive tissues/organs, and the skin dose equivalent rate as one instance.  The former has better medical meaning according to ICRP103, while the latter is conceptually more appropriate to be compared to most previous calculations and measurements.
According to the approximation which is widely made in most previous dose equivalent calculations or measurements in literature, the absorbed dose (mainly due to the absorption of the low energy secondary particles) is given by the energy loss of the high energy primary particle as Note. We expand both genders (Male vs. Female) and the quality factor definitions (Q ICRRP60 vs. Q NASA ). For each case we show the specific contributions from the top six contributing elements, together with the sum of all the 28 elements. We also show the 15 sensitive organ dose equivalent as defined in the ICRP103 with the corresponding w T factor themselves, and finally the weighted total effective dose equivalent. Bold values represent that they are more important than the surrounding values and deserve attention in a first glance.
CHEN ET AL.

10.1029/2022SW003285
10 of 14 then the quality factor Q is multiplied to obtain the dose equivalent. In further simplified calculations (e.g., Mewaldt et al., 2005) or the LET based measurements (e.g., Zeitlin et al., 2013;S. Zhang et al., 2020), a single LET value (the Bethe-Bloch stopping power) is determined by the incident particle's energy or measured, then the dose determination reduces to the particle counting. Such simplified dose determination only holds for an infinitely thin layer at surface, but at any finite depth after energy loss the LET value will change. For realistic astronaut phantom although the skin still has some finite thickness, it is the closest organ to an infinitely thin layer at the surface. So its dose equivalent best matches the literature dose equivalent in the above approximation, when using the ICRP123 fluence-to-dose conversion coefficients.
In Tables 2 and 3 one can directly read that the corresponding skin dose equivalent rates are 77-83 cSv/yr at solar minimum and about 33 cSv/yr at solar maximum, larger than the effective dose equivalent rate by over Table 3 The Same as  1 3.0 1.6 1.6 1.2 1.3 1.6 1.7  Note. Bold values represent that they are more important than the surrounding values and deserve attention in a first glance.

of 14
40% to about 25%, respectively. The skin dose equivalent (or the literature dose equivalent) being larger than the effective dose equivalent is a natural consequence of the self-shielding effect, that unlike the skin, most medically important organs/tissue (with large w T s) are buried deep in human body, so the radiation flux there has been significantly reduced or changed.
Another observation from the comparison between the effective dose equivalent and the skin dose equivalent is the difference of the fractional contribution. Comparing Figures 2 and 3, from the effective dose equivalent to the skin dose equivalent the proton fractional contribution reduces from about 30% to about 20%, and the iron fractional contribution increases from about 13% to about 19%. This can also be attributed to the self-shielding effects. Because the nucleus' LET is proportional to Z 2 but the kinetic energy is generally proportional to its nucleon number A, for heavy ions the former wins and they tend to be stopped at a shallower position, enhancing the surface energy deposit.

Comparison With Literature Results
In literature many GCR spectra level works do not explicitly report dosimetric values, but they are still considered as equivalent to dosimetric studies since the integration with Equation 3 is straightforward especially with the help of ICRP (2013). The aforementioned HelMod model (M. Boschini et al., 2022) is such a recent example, which also uses SDE code (M. Boschini et al., 2016) and recent local interstellar spectra similar to ours, as well as the AMS-02 experimental results for calibration. In Figure 4 we can see that our model traces much more closely the time variation of the GCR flux than the HelMod model, and gives better time-dependent dose rate. This is because our model do a diffusion and drift coefficients fitting for every Bartal rotation, as a tradeoff for the forecast capability.
The most recent 2020 version of the Badhwar-O'Neill models is also calibrated to the AMS-02 data, but more driven by the ACE/CRIS data. Compared to it our model has a similar advantage in the time-dependent performance if compared with the AMS-02 data. Without access to their GCR data we would simply refer readers to the comparison in Figure. 8 of Slaba and Whitman (2019), and quote a few numbers from their Table 2 for quick comparison: the BON2020 model has |Rd| = 3% for both proton and helium, while the early BON2014 model has |Rd| = 7% for proton and |Rd| = 10% for helium. Please note that for |Rd| they are using the integral flux rather than the effective dose equivalent rate.
Since the AMS-02 and PAMELA have the highest precision measurements our model are tuned to reproduce them, but we have crosschecked with other available data, in particular the ACE/CRIS data. We have made comparison figures for our model results and the ACE/CRIS data as time series, for elements of carbon, oxygen, neon, magnesium, silicon, iron, and nickel. For the overlapping time it also has time-dependent flux measurements (we take 189 time points and average the ACE/CRIS spectra over 27 days each), and for a total of 24 ions from boron nucleus to nickel nucleus it has data for seven energy bins each. The estimated relative error is 9%. Simply put, using a total of 31,632 ≈ 189 × 24 × 7 (for cobalt a few data points are missing) we find |Rd| = 14.8%, and for the relative differences Rd which allows cancellation between positive and negative errors we find Rd = −0.73%.
As dosimetric values, our results in Figure 6 can be compared directly to available literature results, both for the effective dose equivalent if realistic human phantom is considered, and for the skin dose equivalent to other dose equivalent calculation or measurement. The effective dose equivalent rates are within the range between differ ent models shown in Figure 1 of Slaba and Blattnig (2014a). At the dose equivalent level, our results are larger than the literature results of Schwadron et al. (2010) at the 2009-2010 solar minimum, and smaller than those of Mewaldt et al. (2005). Note that in the perhaps simplified work of Mewaldt et al. (2005)  impossible without this NASA code but only the public unshielded ICRP123 fluence-to-dose conversion coefficient for us. The generated secondary particles in the thin (∼1 g/cm 2 ) layer shielding are complicated, and even controversial in whether they enhance or suppress the dose (e.g., Mewaldt et al., 2005;Schwadron et al., 2010;Zeitlin et al., 2013).
For the fractional contribution from different elements, at the effective dose equivalent level our results (see Figure 2) are quite similar to Naito et al. (2020) which also uses ICRP123. At the dose equivalent level our results are compared to Durante and Cucinotta (2011) in Figure 3. Only in the combination of solar maximum with the Q ICRP60 definition, iron gives slightly larger fractional contribution than proton, which is consistent with the literature result of Durante and Cucinotta (2011). For other combinations in Figure 3 proton still gives the largest fractional contribution. Since we cannot find the original description for Figure 5 in Durante and Cucinotta (2011), the source of such discrepancies cannot be completely identified. But at least part of the differences can be due to the difference in target. Our result is based on the realistic human phantom, and the dose from certain part of the skin can be contributed by a GCR particle from the opposite direction passing the body. On the other hand, Durante and Cucinotta (2011) probably implicitly assumes a simplified geometry which does not have such contribution from the back, while its assumed 5 g/cm 2 aluminum shielding contributes in the opposite direction to this discrepancies. The target material chemical composition difference between skin and water should contribute to this difference as well. In addition to target, the difference in the flux and the underlying physics models may also be relevant. Therefore for human GCR dose we suggest the new reference results in Figure 2, rather than the update of the historic Figure 3.
Most direct dose measurements correspond to different GCR spectra more specific to the environment, so comparison is not straightforward. For example, the detector on the low earth orbit experiences smaller flux than our "free space" case due to both the Earth's geometric shielding and the magnetic shielding of the magnetosphere. By modeling the GCR propagation in the magnetosphere and its interaction with the top part of the atmosphere, in future we would be able to reproduce the GCR spectra and the radiation dose on such orbit.
In conclusion, the recent GCR model with the AMS-02 data have significantly reduced the radiation dose uncertainty from the GCR spectra. At ≲2%, its uncertainty is subdominant to other uncertainty sources in the computation of radiation dose. The new (effective) dose equivalent rates as a time series from 2006 to 2019 and the rate details at the solar minimum (January 2010) and maximum (February 2014) are provided for direct practical use.