Galactic Cosmic Ray Simulator

• GCR simulators are computational and experimental tools that convert input cosmic ray spectra into realistic modulated radiation fields and dosimetric outputs.
• They employ methods such as transport solvers, Monte Carlo codes (e.g., GEANT4, PHITS), and accelerator-based surrogate fields to replicate deep-space conditions.
• They support applications ranging from space mission shielding to biological effect studies, with validation via benchmark comparisons and response-function analysis.

A Galactic Cosmic Ray (GCR) simulator is a computational or experimental system that represents some part of the GCR chain from source spectrum to propagated particle field, secondary production, dosimetry, microdosimetry, detector response, or biological surrogate exposure. In the literature covered here, the term spans several distinct but related classes of tools: heliospheric modulation solvers based on the Parker transport equation, Galactic propagation frameworks, GEANT4- or PHITS-based Monte Carlo transport codes for atmospheres, planetary bodies, and instruments, and accelerator-based mixed-field facilities intended to approximate deep-space GCR exposure on the ground (Cavallotto et al., 10 Feb 2025, Shestakov et al., 3 Sep 2025, Guo et al., 2019, Peter et al., 2022, Lunati et al., 16 Sep 2025). The unifying feature is not a single implementation but a common operational goal: converting a specified GCR input—such as a local interstellar spectrum, a monthly sunspot number, a top-of-atmosphere spectrum, or a weighted beam sequence—into a physically or operationally relevant output such as a modulated spectrum, absorbed dose, ionization profile, albedo flux, or microdosimetric field.

1. Taxonomy and scope

Within the cited literature, GCR simulators fall into four broad categories. First are transport solvers for Galactic or heliospheric propagation, where the primary unknown is a distribution function or spectrum and the governing equation is typically a diffusion-advection-drift form of the Parker transport equation or a related transport equation (Shestakov et al., 3 Sep 2025, Batista et al., 2019). Second are Monte Carlo environment simulators, which inject GCR particles into atmospheres, regolith, shielding, detectors, or biological geometries and explicitly track secondaries and energy deposition with Geant4, PHITS, or related codes (Guo et al., 2019, Mesick et al., 2018, Pereia et al., 26 Feb 2026, Karafasoulis et al., 2022). Third are ground-based surrogate fields, notably NASA’s GCRSim and the later GSI hybrid active-passive simulator, which use accelerator beams and passive modulators to reproduce a mixed-field proxy of a shielded space environment (Peter et al., 2022, Lunati et al., 16 Sep 2025, Pierobon et al., 16 Sep 2025). Fourth are post-processing and response-function tools that replace full transport with precomputed dose or detector-response functions, as in RADMAREE for Mars (Guo et al., 2019).

Class	Representative systems	Typical output
Transport solvers	COSMICA, CRPropa, finite-difference/SDE heliospheric models	Modulated spectra, distribution functions
Monte Carlo environment models	GEANT4, PHITS, DYASTIMA-R, ECRS	Dose, ionization, albedo, secondary fields
Ground-based surrogate fields	GCRSim, SimGCRSim, GSI hybrid simulator	Mixed irradiation fields, LET or microdosimetric spectra
Response-function tools	RADMAREE ready functions	Spectrum-to-dose conversion

This classification indicates that “GCR simulator” is not restricted to one scale or methodology. In some papers it denotes a spectral or transport model near Earth; in others it denotes a beam facility intended to emulate deep-space exposure; in still others it denotes a Monte Carlo pipeline that transforms a prescribed GCR composition into a local radiation observable. A plausible implication is that comparisons between simulators are only meaningful when the simulated level of description—spectrum, field, dose, or biology surrogate—is specified explicitly.

2. Transport-equation simulators for modulation and propagation

A major class of GCR simulators solves the Parker transport equation or a related transport equation for the isotropic distribution function. In the heliospheric modulation paper, the governing equation includes convection by the solar wind, anisotropic spatial diffusion, adiabatic cooling, and drift, with the diffusion tensor written in field-aligned coordinates as

$$\boldsymbol{\hat{\kappa}}= \begin{pmatrix} \kappa_{\parallel} & 0 & 0\ 0 & \kappa_{\perp} & 0\ 0 & 0 & \kappa_{\perp} \end{pmatrix},$$

and the full three-dimensional implementation uses the Parker spiral interplanetary magnetic field, current-sheet drift regularized with a $\tanh$ replacement for the field reversal, and a backward-time stochastic differential equation (SDE) formulation (Shestakov et al., 3 Sep 2025). That work also develops a Crank–Nicolson finite-difference realization for simplified spherical problems and validates both numerical realizations against an analytical Toptygin solution, then benchmarks the SDE implementation against Kota & Jokipii and Burger (Shestakov et al., 3 Sep 2025).

COSMICA is a GPU-native Monte Carlo realization of the same general heliospheric modulation problem. It solves backward-in-time the system of SDEs equivalent to the Parker Transport Equation, is based on the HelMod physical model, and uses a momentum/rigidity-based formulation in which the loss term becomes $L=0$ and one fewer SDE must be solved (Cavallotto et al., 10 Feb 2025). The code implements NVIDIA CUDA parallelization, shared-memory and constant-memory optimizations, rounding of quasi-particle counts to multiples of 32, and distribution of independent energy bins across multiple GPUs and clusters. The paper reports about a 40X speed-up from the initial CPU-to-GPU port and roughly an 80X speed-up compared with the CPU-only implementation after further optimization (Cavallotto et al., 10 Feb 2025).

At Galactic scales, CRPropa extends a framework originally associated with ballistic extragalactic propagation by adding a stochastic solver for the modified Parker transport equation,

$$\frac{\partial n}{\partial t} + \vec{u}\cdot\nabla n = \nabla\cdot(\hat{\kappa}\nabla n) + \frac{1}{p^2}\frac{\partial}{\partial p}\left(p^2 \kappa_{pp} \frac{\partial n}{\partial p} \right) + \frac{1}{3}(\nabla \vec{u})\frac{\partial n}{\partial \ln p} + S(\vec{x},p,t),$$

thereby combining diffusion, advection, momentum diffusion, and source injection with arbitrary three-dimensional magnetic fields such as JF12 (Batista et al., 2019). The code propagates pseudo-particles whose stochastic steps reproduce the transport-equation solution in the ensemble limit. The paper emphasizes that the low-energy Galactic mode is only as predictive as the chosen transport coefficients, which is a methodological constraint rather than an implementation defect (Batista et al., 2019).

A more empirical spectral simulator appears in the ACE/CRIS-based heavy-nuclei model. Rather than solving the Parker equation directly, it maps monthly averaged international sunspot number (SSN, V2.0) to monthly averaged near-Earth differential spectra for $5 \le z \le 28$ over roughly $30$–$500$ MeV/nuc (Fu et al., 2020). Its key parameters $\alpha(t)$, $\bar p(t)$, and $\eta(t)$ are correlated with lagged SSN, with distinct lags for odd and even solar cycles; for example, $\tanh$0 uses 14 months in odd cycles and 6 months in even cycles (Fu et al., 2020). This is still a GCR simulator, but one whose input variable is solar activity rather than a particle source spectrum.

The discrete-supernova stochastic model addresses a different transport level: the age and path-length distribution of GCR protons arriving at the solar system. It solves the stochastic differential equations equivalent to the Parker diffusion-convection equation under simplified assumptions—three-dimensional isotropic diffusion, no Galactic wind, no energy change, and no free-escape boundaries—and propagates pseudo-particles backward until they intersect an active supernova remnant (Miyake et al., 2014). The resulting age and path-length distributions show low-end cut-offs that the authors attribute to the discreteness of supernova occurrences. This suggests that, in some contexts, the defining function of a GCR simulator is not spectral prediction but reconstruction of hidden transport statistics such as grammage and residence time.

3. Monte Carlo simulators for atmospheres, planetary bodies, and radiation fields

A second major class of GCR simulator is the explicit Monte Carlo transport code. In the Martian environment study, GEANT4 version 4.10 and PLANETOCOSMICS are used with Mars Climate Database atmospheric profiles and a 100 m $\tanh$1 regolith layer at Gale Crater to construct response matrices $\tanh$2 for protons and $\tanh$3He ions (Guo et al., 2019). Those matrices are then converted into ready-to-use dose functions so that an arbitrary top-of-atmosphere SEP or GCR spectrum $\tanh$4 can be folded into absorbed dose or dose equivalent without rerunning the full Monte Carlo: $\tanh$5 The functions cover primary energies roughly from 20 MeV to $\tanh$6 MeV, with stored output histograms spanning 1 MeV to $\tanh$7 MeV, and are implemented in RADMAREE as a public spectrum-to-dose converter for Mars (Guo et al., 2019).

For protoplanetary disks, PHITS v2.82 is used as a GCR ionization simulator in a homogeneous cylindrical gas target with height and radius each fixed to $\tanh$8 in column depth (Kataoka et al., 2016). The work explicitly adds neutron decay with a mean lifetime of 887 s, because deep-column ionization depends sensitively on long-lived neutrons. The ionization attenuation length for standard GCRs is updated to $\tanh$9 in pure hydrogen and $L=0$0 for a cosmic-abundance H/He mixture, compared with the canonical $L=0$1 of UN81 (Kataoka et al., 2016). Here the simulator output is not dose but a depth-dependent ionization profile and its asymptotic attenuation length.

Atmospheric shower simulators instantiate the same pattern on Earth. DYASTIMA-R extends DYASTIMA by taking Geant4 atmospheric cascade outputs and converting them into absorbed and equivalent dose in a cylindrical water phantom of height 1.75 m and radius 0.25 m (Paschalis et al., 2016). ECRS, also GEANT4-based, models a full Earth-centered environment with 100 atmospheric layers, IGRF plus Tsyganenko geomagnetic fields, and, in some runs, an 11-km water shell to represent the ocean (Sarajlic et al., 2018). It launches primary cosmic rays from 1.2 Earth radii and records global latitude-longitude maps, altitude distributions, and lateral distributions of secondaries. A concrete runtime example reported in that work is that a 60 GeV proton takes about 9 s without geomagnetic field and about 700 s with geomagnetic field enabled (Sarajlic et al., 2018).

Planetary-body Monte Carlo simulators serve a related but distinct function. The lunar benchmark study uses Geant4 to model GCR protons and alpha particles interacting with a four-layer Moon of radius 1738.1 km, then compares neutron density versus depth and leakage flux to Apollo 17 Lunar Neutron Probe Experiment data (Mesick et al., 2018). The dominant sensitivity is not the high-energy string model but the intra-nuclear cascade treatment: relative to INCL, the Binary cascade produces a neutron density peak/integral about 12% higher, whereas Bertini is about 34% higher (Mesick et al., 2018). The paper concludes that Geant4 reproduces the shape of the neutron-depth profile well but systematically overpredicts the absolute normalization, so it is reliable as a trend simulator but not fully self-calibrating (Mesick et al., 2018).

The 2026 stratospheric chemistry study uses Geant4 in a four-shell spherical stratospheric model to propagate protons, alpha particles, CNO, and Si under Polar and Equatorial cutoff conditions. It finds that 95.6% of the Polar and 96.4% of the Equatorial 15–35 km energy deposition budget comes from cascade secondaries, and converts layer-resolved energy deposition to ion-pair, NO$L=0$2, and HO$L=0$3 production using $L=0$4 eV per ion pair (Pereia et al., 26 Feb 2026). This is again a GCR simulator, but one whose terminal observable is atmospheric chemistry rather than radiation protection.

4. Ground-based mixed-field simulators and surrogate GCR environments

In accelerator-based usage, a GCR simulator is not primarily a code but a radiation field generator. NASA Space Radiation Laboratory’s GCRSim is represented in the neuronal damage study as a 33 ion-energy beam fluence distribution normalized to a total dose of 0.5 Gy, intended to approximate the dose distribution behind 20 g cm$L=0$5 of aluminum shielding (Peter et al., 2022). The composition includes 14 proton energies, 14 alpha-particle energies, one carbon beam, one oxygen beam, one silicon beam, one iron beam at 600 MeV/n, and one substituted iron beam at 1000 MeV/n in place of titanium because Ti cross sections are not supported by Geant4-DNA (Peter et al., 2022). The reduced SimGCRSim field uses only six entries—p, p, $L=0$6, O, Si, Fe—yet the paper reports no statistically significant differences in the average nanoscale physical responses examined (Peter et al., 2022).

The murine dosimetry study treats the NSRL GCRsim facility as a controlled small-animal irradiation system designed to approximate the deep-space GCR environment for mixed-field radiobiology (Hosseini et al., 6 Sep 2025). In the simplified configuration used there, the field includes H, He, O, Si, and Fe, with hydrogen delivered at two different energies (Hosseini et al., 6 Sep 2025). Beam uniformity is checked with a CCD camera or a large pixelated ion chamber and is typically within 3% RMS over the target region, while nominal and delivered dose agree within 10% (Hosseini et al., 6 Sep 2025). The PHITS-based dose-conversion-factor library derived in that paper is therefore an enabling layer for quantitative use of the facility rather than the simulator field definition itself (Hosseini et al., 6 Sep 2025).

The GSI hybrid active-passive simulator is a later European development explicitly optimized to reproduce the 2010 solar-minimum GCR spectrum behind 10 g/cm$L=0$7 Al at 1 AU using only $L=0$8 beams at 1 GeV/u, 0.7 GeV/u, and 0.35 GeV/u combined with passive modulators (Lunati et al., 16 Sep 2025). The system superposes six irradiation configurations: 3 complex modulators and 3 slab-modulator configurations, all preceded by a 32-layer stainless-steel mesh modulator with 60 µm wire diameter and 75 µm pitch (Lunati et al., 16 Sep 2025). The optimization minimizes a weighted mismatch between the six-configuration spectral basis and target elemental spectra $L=0$9, excluding neutrons $$\frac{\partial n}{\partial t} + \vec{u}\cdot\nabla n = \nabla\cdot(\hat{\kappa}\nabla n) + \frac{1}{p^2}\frac{\partial}{\partial p}\left(p^2 \kappa_{pp} \frac{\partial n}{\partial p} \right) + \frac{1}{3}(\nabla \vec{u})\frac{\partial n}{\partial \ln p} + S(\vec{x},p,t),$$0 from the cost function (Lunati et al., 16 Sep 2025).

Its experimental implementation at GSI Cave A uses remote-controllable modulator exchangers, raster-scanned $$\frac{\partial n}{\partial t} + \vec{u}\cdot\nabla n = \nabla\cdot(\hat{\kappa}\nabla n) + \frac{1}{p^2}\frac{\partial}{\partial p}\left(p^2 \kappa_{pp} \frac{\partial n}{\partial p} \right) + \frac{1}{3}(\nabla \vec{u})\frac{\partial n}{\partial \ln p} + S(\vec{x},p,t),$$1 beams, a large-area parallel-plate ionization chamber for monitoring, and target handling with a Universal Robots UR10e robotic arm (Pierobon et al., 16 Sep 2025). The final reconstructed field, obtained by weighting and summing the six measured microdosimetric spectra, yields $$\frac{\partial n}{\partial t} + \vec{u}\cdot\nabla n = \nabla\cdot(\hat{\kappa}\nabla n) + \frac{1}{p^2}\frac{\partial}{\partial p}\left(p^2 \kappa_{pp} \frac{\partial n}{\partial p} \right) + \frac{1}{3}(\nabla \vec{u})\frac{\partial n}{\partial \ln p} + S(\vec{x},p,t),$$2 experimentally and $$\frac{\partial n}{\partial t} + \vec{u}\cdot\nabla n = \nabla\cdot(\hat{\kappa}\nabla n) + \frac{1}{p^2}\frac{\partial}{\partial p}\left(p^2 \kappa_{pp} \frac{\partial n}{\partial p} \right) + \frac{1}{3}(\nabla \vec{u})\frac{\partial n}{\partial \ln p} + S(\vec{x},p,t),$$3 in simulation, with $$\frac{\partial n}{\partial t} + \vec{u}\cdot\nabla n = \nabla\cdot(\hat{\kappa}\nabla n) + \frac{1}{p^2}\frac{\partial}{\partial p}\left(p^2 \kappa_{pp} \frac{\partial n}{\partial p} \right) + \frac{1}{3}(\nabla \vec{u})\frac{\partial n}{\partial \ln p} + S(\vec{x},p,t),$$4 experimentally and $$\frac{\partial n}{\partial t} + \vec{u}\cdot\nabla n = \nabla\cdot(\hat{\kappa}\nabla n) + \frac{1}{p^2}\frac{\partial}{\partial p}\left(p^2 \kappa_{pp} \frac{\partial n}{\partial p} \right) + \frac{1}{3}(\nabla \vec{u})\frac{\partial n}{\partial \ln p} + S(\vec{x},p,t),$$5 in simulation (Pierobon et al., 16 Sep 2025). The paper presents this as the first European ground-based hybrid active-passive GCR simulator and shows reasonable agreement with a microdosimetric spectrum measured during STS-102, while also noting the truncation of the high-energy tail because SIS-18 cannot exceed 1 GeV/u (Pierobon et al., 16 Sep 2025).

A recurrent theme across these facilities is that the surrogate field is intentionally shielded-field-specific. NSRL’s 33-beam GCRSim targets a field behind 20 g cm$$\frac{\partial n}{\partial t} + \vec{u}\cdot\nabla n = \nabla\cdot(\hat{\kappa}\nabla n) + \frac{1}{p^2}\frac{\partial}{\partial p}\left(p^2 \kappa_{pp} \frac{\partial n}{\partial p} \right) + \frac{1}{3}(\nabla \vec{u})\frac{\partial n}{\partial \ln p} + S(\vec{x},p,t),$$6 Al (Peter et al., 2022), whereas the GSI simulator targets 10 g/cm$$\frac{\partial n}{\partial t} + \vec{u}\cdot\nabla n = \nabla\cdot(\hat{\kappa}\nabla n) + \frac{1}{p^2}\frac{\partial}{\partial p}\left(p^2 \kappa_{pp} \frac{\partial n}{\partial p} \right) + \frac{1}{3}(\nabla \vec{u})\frac{\partial n}{\partial \ln p} + S(\vec{x},p,t),$$7 Al at 1 AU (Lunati et al., 16 Sep 2025). This suggests that a ground-based GCR simulator is best understood as an approximation to a mission-relevant transport state rather than a literal replica of the free-space Galactic spectrum.

5. Response functions, dosimetry, microdosimetry, and detector-response layers

Many GCR simulators are operationally useful only when coupled to a response layer. RADMAREE is a clear example: its precomputed dose-per-unit-flux functions convert arbitrary proton or $$\frac{\partial n}{\partial t} + \vec{u}\cdot\nabla n = \nabla\cdot(\hat{\kappa}\nabla n) + \frac{1}{p^2}\frac{\partial}{\partial p}\left(p^2 \kappa_{pp} \frac{\partial n}{\partial p} \right) + \frac{1}{3}(\nabla \vec{u})\frac{\partial n}{\partial \ln p} + S(\vec{x},p,t),$$8He spectra at the top of the Martian atmosphere into absorbed dose rate in $$\frac{\partial n}{\partial t} + \vec{u}\cdot\nabla n = \nabla\cdot(\hat{\kappa}\nabla n) + \frac{1}{p^2}\frac{\partial}{\partial p}\left(p^2 \kappa_{pp} \frac{\partial n}{\partial p} \right) + \frac{1}{3}(\nabla \vec{u})\frac{\partial n}{\partial \ln p} + S(\vec{x},p,t),$$9 and dose equivalent rate in $5 \le z \le 28$0, separately for downward, upward/albedo, and total contributions (Guo et al., 2019). The upward component is generally smaller than the downward dose and is often < 10% near cutoff energies, but becomes relevant where regolith-generated secondaries matter (Guo et al., 2019).

The murine GCRsim dosimetry study formalizes this response concept as a library of fluence-to-absorbed-dose conversion factors,

$5 \le z \le 28$1

with $5 \le z \le 28$2 values in Gy·m$5 \le z \le 28$3·ion$5 \le z \le 28$4 or Gy·cm$5 \le z \le 28$5·ion$5 \le z \le 28$6 depending on tabulation (Hosseini et al., 6 Sep 2025). PHITS version 3.32, the tetrahedral 25 g MOBY phantom, and at least 10$5 \le z \le 28$7 histories per source particle are used to derive organ- and voxel-level DCFs for six irradiation orientations (Hosseini et al., 6 Sep 2025). For the whole body, the delivered-dose-to-organ-dose conversion factor is 1.10 for all six orientations, and calculated absorbed doses for mice irradiated at nominal 0.20, 0.40, and 0.75 Gy agree with ion-chamber measurements within 10% (Hosseini et al., 6 Sep 2025).

At the detector level, the MIDAS study uses GEANT4 to simulate an external GCR field from OMERE/ISO 15390 for solar minimum 1996.4, with nuclei $5 \le z \le 28$8 to $5 \le z \le 28$9 over 0.1 MeV/n to 20 GeV/n in 835 logarithmic bins (Karafasoulis et al., 2022). Tracks are reconstructed from clusters of hit pixels in 50 $30$0m silicon layers, and the measured silicon LET is converted to $30$1 in water via a fitted relation

$30$2

with $30$3 and $30$4 (Karafasoulis et al., 2022). The inferred free-space dose rate is 0.506 mGy/day up to 400 keV/µm and 0.526 mGy/day up to 1000 keV/µm, with corresponding dose equivalent rates of 2.32 mSv/day and 2.57 mSv/day (Karafasoulis et al., 2022). Inside a cylindrical 20 g/cm$30$5 Al spacecraft model, the same pipeline yields 0.2486 mGy/day and 0.634 mSv/day up to 400 keV/µm, showing how shielding alters both composition and quality factor (Karafasoulis et al., 2022).

Microdosimetry provides a finer-grained response descriptor in the GSI experimental implementation. A Far West Technology LET-1/2 tissue-equivalent proportional counter filled with propane at $30$6, equivalent to a 2 µm tissue sphere and located 2 m downstream of the last beamline element, measures lineal-energy spectra $30$7 and dose distributions $30$8 for each of the six configurations (Pierobon et al., 16 Sep 2025). The final weighted reconstruction is compared with Geant4 and with quality-factor prescriptions $30$9 and $500$0. The 0.35 GeV/u complex-modulator configuration shows the largest discrepancy, whereas the slab configurations show satisfactory agreement (Pierobon et al., 16 Sep 2025). In this context the simulator is inseparable from its field-characterization apparatus.

The neuronal GCRSim study pushes the response layer to nanometer scale. Using TOPAS 3.4, TOPAS-nBio 1.0-beta, and Geant4-DNA physics, it scores dose to soma, dendrites, and spines; energy deposition by mediating species; and event frequencies in a realistic CA1 pyramidal neuron (Peter et al., 2022). Ionizations mediate 68 ± 1% of total energy deposition, secondary electrons account for 77 ± 1%, and vibrational excitations make up 70 ± 2% of all energy deposition events by count (Peter et al., 2022). That result does not redefine GCRSim, but it demonstrates how a ground-based mixed-field simulator can be coupled to a track-structure response model that experiments cannot access directly.

6. Validation, limitations, and scientific use

Across these works, validation usually proceeds by benchmark comparison rather than by formal proof. The heliospheric SDE transport code is validated against an analytical steady-state solution and compared successfully with Kota & Jokipii and Burger (Shestakov et al., 3 Sep 2025). COSMICA is presented as reproducing AMS-02 proton spectra within experimental uncertainty while reducing runtime by orders of magnitude (Cavallotto et al., 10 Feb 2025). The ACE/CRIS empirical heavy-ion simulator reports $500$1 values of 3.5% for carbon, 3.9% for oxygen, 4.6% for silicon, and 5.3% for iron in monthly integral intensities, with more than 80% of $500$2 values within $500$3 (Fu et al., 2020).

Monte Carlo simulators are often validated against measured albedo, ionization, or dosimetry. The lunar Geant4 benchmark shows that all tested physics lists overpredict Apollo 17 neutron-depth data, although INCLXX is both the closest to the measurements and about 20–40% faster than the other tested models (Mesick et al., 2018). The murine GCRsim DCF library is validated against ion-chamber measurements within 10% (Hosseini et al., 6 Sep 2025). The GSI hybrid simulator compares experimental microdosimetric spectra with Geant4 and obtains a final quality factor close to the design target, but explicitly notes limitations from the 1 GeV/u SIS-18 beam ceiling, possible modulator imperfections, and nuclear-transport inaccuracies in heavily modulated cases (Pierobon et al., 16 Sep 2025).

Several recurring limitations are methodological rather than incidental. Transport solvers depend on assumed diffusion coefficients, drift-reduction prescriptions, source spectra, and field models (Batista et al., 2019, Shestakov et al., 3 Sep 2025). Empirical spectral simulators inherit the validity range of their calibration data; the ACE/CRIS heavy-ion model is constrained to about 30–500 MeV/nuc and to monthly averaged near-Earth conditions (Fu et al., 2020). Monte Carlo environment simulators can be limited by hadronic physics lists, neutron treatment, unsupported cross sections, or the need for extrapolation. In the neuronal study, unsupported TOPAS-nBio cross sections above 500 MeV/n for protons and above 100 MeV/n for alpha particles require curve fitting and extrapolation (Peter et al., 2022). In the Martian dose tool, the ready functions are built only for protons and $500$4He, albeit with the paper noting that those species contribute about 90% of total dose and about 80% of total dose equivalent when extending to the full GCR population (Guo et al., 2019).

The applications are correspondingly diverse. GCR simulators support Mars mission planning and habitat shielding assessment (Guo et al., 2019), protoplanetary-disk ionization and dead-zone estimates (Kataoka et al., 2016), aviation and spaceflight dosimetry (Paschalis et al., 2016), detector design and onboard monitoring (Karafasoulis et al., 2022), planetary neutron and gamma-ray interpretation (Mesick et al., 2018), atmospheric chemistry coupling (Pereia et al., 26 Feb 2026), and mixed-field radiobiology in mice and neurons (Hosseini et al., 6 Sep 2025, Peter et al., 2022). This suggests that the encyclopedic meaning of the term is functional rather than architectural: a GCR simulator is any rigorously specified system that maps a Galactic-cosmic-ray description into a downstream physical, dosimetric, chemical, or biological quantity under stated transport assumptions and geometry.