On the Comprehensive 3D Modeling of the Radiation Environment of Proxima Centauri b: A New Constraint on Habitability?

The combined influence of stellar energetic particles and galactic cosmic rays (GCRs) on the radiation environment, and hence potential habitability, of Earth-like exoplanets is relatively unknown. The present study, for the first time, comprehensively models the transport of these particles in a physics-first manner, using a unique suite of numerical models applied to the astrosphere of Proxima Centauri. The astrospheric plasma environment is modeled magnetohydrodynamically, while particle transport is modeled using a 3D ab initio GCR modulation code, as opposed to previous 1D approaches to this problem. StEP intensities are also calculated using observed stellar event profiles for Proxima Centauri as inputs. Computed intensities are then used to calculate possible atmospheric ionization effects and dose rates. We demonstrate that the contribution of GCRs to these quantities is indeed significant, contrary to the conclusions of previous studies. Furthermore, we propose a novel potential constraint on exoplanetary habitability, namely the rotational period of the host star, based on the unique 3D modeling approach presented here.


Introduction
One of the most interesting stellar systems to study in the context of exoplanetary habitability is the system of our nearest neighbor, Proxima Centauri, an M5.5Ve flaring star, in which a rocky, presumably Earth-like exoplanet (Proxima Centauri b) is orbiting the star at a distance of only 0.0485 au.Despite its close proximity, Proxima Centauri b (hereafter Prox Cen b) is located well within the habitable zone of the stellar system (Hill et al. 2023) and, with a minimum mass of 1.27 Earth masses, may have bulk characteristics similar to Earth (Anglada-Escudé et al. 2016).Despite receiving only 65% of Earth's mean total stellar irradiation, an Earth-like atmosphere with, for example, a surface carbon dioxide content of a few percent might lead to habitable conditions (Meadows et al. 2018).However, stellar space weather (i.e., (E)UV, stellar flares, coronal mass ejections, and energetic particles) might play a crucial role in the context of habitability (Tabataba-Vakili et al. 2016;Fraschetti et al. 2019;Herbst et al. 2019a;Scheucher et al. 2020).Galactic cosmic rays (GCRs), consisting of charged particles with considerable energies, can also reasonably be expected to play some role in the potential habitability of exoplanets in the astrospheres such as that considered here, whether in terms of their influence on the possible atmospheres of these planets or in terms of their contribution to the planetary radiation environment.Modeling these influences, however, requires a thorough understanding of not only astrospheric and stellar but also planetary magnetospheric and atmospheric processes.Both GCRs and stellar energetic particles (StEPs) must be included when modeling the influence of energetic particles on (exo)planetary atmospheres.However, their relative importance is the subject of some debate: while StEPs could lead to significant changes in atmospheric climate and chemistry-and with that on atmospheric biosignatures like ozone and methane (Herbst et al. 2019a;Scheucher et al. 2020)-GCRs are often argued to be negligible in the context of exoplanetary habitability (Sadovski et al. 2018;Mesquita et al. 2021Mesquita et al. , 2022)).
Several studies have investigated the modulation of GCRs in different astrospheres, drawing varying conclusions on how significant such modulation would be (Herbst et al. 2022).Although one study suggests that GCR intensities at exoplanetary orbits may be significantly larger than intensities at Earth (Herbst et al. 2020b), the majority of studies report the opposite, finding intensities close to, or several orders of magnitude below Earth levels (Sadovski et al. 2018;Struminsky et al. 2018;Mesquita et al. 2021Mesquita et al. , 2022;;Rodgers-Lee et al. 2021a, 2021b).However, all of these studies suffer from the significant limitation that they either solve a 1D (Parker 1965) GCR transport equation (TPE) or rely on the force-field solution (Gleeson & Axford 1968) for that equation.Given the observed and modeled three-dimensional structures of the astrospheres in question, this is a severe limitation, as it would not be possible to realistically model the complex physical mechanisms governing the transport of GCRs (e.g., Engelbrecht et al. 2022).For example, neither approach can be employed to model the influence of drifts due to gradients in, and the curvature of, the 3D astrospheric magnetic field (AMF), long known to be important in the heliospheric transport of GCRs (Jokipii et al. 1977;Jokipii & Thomas 1981;Kota & Jokipii 1983;Potgieter & Burger 1990), and an inherently 3D process.Therefore, modeling this process in less than 3D is unphysical (Engelbrecht & Di Felice 2020).The force-field approach is particularly restrictive in terms of the unphysical assumptions made in its derivation, i.e., a one-dimensional, spherically symmetric system with no adiabatic energy changes (Caballero-Lopez & Moraal 2004), which severely limits its predictive capacity (Engelbrecht & Di Felice 2020).The 1D approaches to this problem have similar limitations (Savopulos & Quenby 1997, 1998;Light et al. 2022).At the least, a fully 3D approach to modeling GCR modulation would allow for the realistic modeling of the inherently 3D mechanisms governing GCR transport, leading to more realistic estimations of the role that these particles play in the potential habitability of exoplanets in various astrospheres.The present study, for the first time, studies the modulation of GCRs using a fully 3D stochastic solver of the Parker TPE that has been successfully used in heliospheric GCR modulation studies (see, e.g., Engelbrecht & Burger 2015;Moloto et al. 2018;Engelbrecht & Moloto 2021), thereby properly taking into account, also for the first time, the influence of GCR drifts and diffusion both parallel and perpendicular to the AMF.These mechanisms alone have long been known to influence the heliospheric modulation of GCRs significantly (see, e.g., Engelbrecht et al. 2022, and references therein).We demonstrate that this is also true for Proxima Centauri and could be reasonably expected to be so for other cool star astrospheres.We also demonstrate that the stellar rotation rate, through its influence on the AMF geometry, and hence on the 3D diffusion of energetic charged particles, can have a hitherto unconsidered yet significant influence on the radiation environment of exoplanets and thus affect their potential habitability, providing a further potential constraint for exoplanet habitability.This is done by employing a suite of models pertaining to various aspects of Prox Cen b's environment that are relevant to energetically charged particle transport in its astrosphere.The astrosphere itself, including large-scale plasma quantities such as the AMF and stellar wind, is modeled using the CRONOS magnetohydrodynamic (MHD) code (Kissmann et al. 2018).These outputs are, in combination with analytical models for turbulence quantities constructed using the observed behavior of their heliospheric counterparts in mind (see, e.g., Herbst et al. 2022), used as inputs for a state-of-the-art 3D stochastic solver for the Parker (1965) GCR TPR introduced by Engelbrecht & Burger (2015) in order to calculate GCR differential intensities at the location of Proxima Cen b.
Solar-terrestrial studies have shown that the solar space weather strongly impacts planetary atmospheres by inducing the production of secondary particle cascades that lead to changes in atmospheric electric properties (Mallios et al. 2022), climate and chemistry changes (Grenfell et al. 2022;Reddmann et al. 2023), and enhancements of, e.g., radiation exposure (Guo et al. 2019).The impact of stellar space weather, particularly cosmic rays (CRs), may play an even more crucial role in determining atmospheric climate and chemistry changes, and hence the habitability of Earthlike exoplanets (Herbst et al. 2019b;Scheucher et al. 2020).While GCRs vary on longer timescales related to possible stellar cycles (Engelbrecht et al. 2022), StEPs are accelerated during transient stellar phenomena (e.g., stellar flares and coronal mass ejections), displaying short-term variability.Here, we investigate, using a numerical transport model (Strauss & Fichtner 2015), the acceleration and transport of 72 stellar flares occurring over a ∼52 day interval using Transiting Exoplanet Survey Satellite (TESS) measurements (Vida et al. 2019), assuming identical transport conditions to what was assumed for GCR transport.The corresponding series of flare-induced StEP intensities in the Proxima Centauri system is simulated and found to contain several events comparable to large ground-level enhancement (GLE) events observed at Earth over a relatively short period.
With that, the question arises as to whether increases in GCR and StEP intensities would correspond to a significant effect in the radiation environment of Prox Cen b.We demonstrate the impact of the modeled particle spectra, assuming a terrestrial atmosphere consisting of 50.42% N 2 , 27.58% CO 2 , and 21.06% O 2 with a relative surface humidity of 50%, by investigating changes in atmospheric ionization and dosimetric quantities (i.e., absorbed dose rates and dose equivalent).The next section provides details of the suite of models employed here.Subsequently, the results so calculated are presented, followed by a section dedicated to a summary of said results, and the conclusions drawn therefrom.

The Model Suite
Given the comprehensive, multi-tiered modeling approach taken here, this section provides additional motivation and details pertaining to each step taken in this study.

Modeling the Astrospheric Plasma Environment of Proxima Centauri
In the present study, a dual approach is taken to modeling the plasma environment of Proxima Centauri.Broadly, large-scale plasma quantities are modeled using a 3D MHD code.Where necessary, as discussed below, analytical models for these quantities are employed.

MHD Modeling
The astrosphere of Proxima Centauri is modeled, as a first approach, using the single-fluid 3D finite-volume MHD code CRONOS (Kissmann et al. 2018), motivated by the relative simplicity of such an approach, and the ease of comparison with the results of previous single-fluid MHD studies of astrospheres (e.g., Herbst et al. 2020b).Future work will take advantage of the multi-fluid prescription available for this code, in order to more fully investigate, for example, the influence of interstellar neutrals, which has been noted from heliospheric studies to be significant (see, e.g., Alouani-Bibi et al. 2011;Pogorelov et al. 2017;Kleimann et al. 2022;Michael et al. 2022;Sokółet al. 2022), and to model the large-scale transport of turbulence quantities (for heliospheric examples, see, e.g., Wiengarten et al. 2016;Usmanov et al. 2018;Fichtner et al. 2020;Kleimann et al. 2023).The code is based on a Riemann solver performing simulations on a star-centered spherical grid with a resolution of N r × N θ × N f = 1024 × 64 × 128 cells following a similar approach to that in Scherer et al. (2020).The single-fluid MHD results are based on solving the 3D ideal MHD equations, which are given by with ρ representing the mass density, v the velocity, p the pressure and the total energy by . The unit tensor is represented by Ĩ and the dyadic/tensorial product by ⊗.The cooling and heating functions, −R L , are treated as is done in Scherer et al. (2020).The stellar wind velocity of Proxima Centauri is assumed to be 1500 km s −1 following the MHD simulation results of Garraffo et al. (2016), and a stellar wind density of 0.25 cm −3 at 1 au, which is calculated from the mass-loss rate of 2 × 10 −15 M e yr −1 (see Wood et al. 2001) and stellar wind velocity at the inner boundary of the model located at 30 au (see, e.g., Wilkin 2000;Herbst et al. 2020b).
The stellar magnetic field of Proxima Centauri is scaled, following the approach of Herbst et al. (2020b), to be 2.4 nT at 1 au, based on observations of the surface AMF of Proxima reported by Reiners & Basri (2008).The interstellar medium (ISM) density surrounding the astrosphere is assumed to be 0.06 cm −3 with a velocity of 30 km s −1 .Along the inflow line, that is a line parallel to the inflow vector passing through the star, the termination shock (TS) distance (r ts is 55.05 au and the astropause distance is 80.65 au (see Figure 1 and its accompanying discussion below).There is no bow shock because the fast-magnetosonic Mach number in the ISM is less than 1.It is of interest to compare this result with that calculated using the analytic expression for the TS distance from Parker (1963) or Wilkin (2000), given by , 5 ts 0 0,sw sw 2 ISM ISM 2 ISM ISM 2 using either the stellar mass-loss rate  M or the stellar wind number density n 0,sw at a reference distance r 0 = 1 au with the ISM speed v ISM and the stellar wind speed v sw , as well as the ISM number density n ISM .With the parameters used here, i.e., stellar input parameters n 0,sw = 0.25 cm −3 , v sw = 1500 km s −1 and ISM parameters n ISM = 0.06 cm −3 ,v ISM = 30 km s −1 one finds from Equation (5) that r ts = 102.06au, a significant overestimation of the MHD result found in this study.The reason for this is that Equation (5) is only valid for a pure hydrodynamic approach, where no magnetic field is present.In this case, the inflow line is a streamline, passing through the star, the TS, and the stagnation point at the astropause.In the case of MHD simulation, the inflow line must not be a streamline and the above approach fails.For a detailed discussion of this see Scherer et al. (2020).

Analytical Modeling
For regions within the inner boundary of the MHD solver described above, analytical models for quantities such as the AMF and stellar wind speed, normalized so as to assume the same values in the MHD code at 30 au, need to be employed.Several assumptions must be made to analytically model the large-scale plasma quantities required to solve the Parker TPE within the astrospheric TS.First, it is assumed, based on the MHD simulation results, that the AMF of Proxima Centauri can be described using a Parker (1958) magnetic field model, given by where |A| can be chosen so as to ensure an AMF magnitude of 2.4 nT at r o = 1 au.Here, we assume a similar magnetic polarity cycle to that observed in the heliosphere so that the sign of A indicates the polarity of the field, such that when A > 0, the field in the northern hemisphere points away from the star, and toward it in the southern hemisphere.These directions are reversed for A < 0. The Parker field lines have a spatial structure on cones of constant colatitude θ, wound at an angle given by which in turn is a function of the stellar wind speed V sw , assumed here to be a constant 1500 km s −1 following the MHD simulation results, and stellar rotation rate Ω, chosen here to correspond to the observed (Collins et al. 2017) rotation period of 82.6 days for Proxima Centauri.As for the heliospheric magnetic field (HMF), we assume that the location for the field source surface is at r s = 0.005 au.This explicitly assumes that Prox Cen b lies outside the star's Alfvén surface (Klein et al. 2021).Although this is not yet resolved in the MHD simulations, we also include an astrospheric current sheet (ACS) in a manner standard to heliospheric GCR modulation studies.The ACS angle is given by Kota & Jokipii (1983) as ns s 1 with α the stellar tilt angle, the angle between the star's magnetic and rotational axes.As a first approach, this angle is assumed to be zero, implying a flat rather than wavy current sheet, although observations indicate larger tilt angles are possible for Proxima Centauri (Klein et al. 2021).
Although limited observations and MHD simulation results assist in modeling large-scale plasma quantities, no such observations or simulations exist for small-scale plasma turbulence quantities needed to self-consistently model the diffusion coefficients for GCRs and StEPs (Herbst et al. 2022).However, given the relative similarity between the astrosphere of Proxima Centauri and the heliosphere, estimates can be made based on theory and observation in the heliosphere (Matthaeus & Velli 2011;Bruno & Carbone 2016;Oughton & Engelbrecht 2021).As a point of departure, the slab/2D composite model of turbulence Bieber et al. (1994) is assumed.In the heliosphere, magnetic variances and correlation scales have been shown to scale with the HMF magnitude (Matthaeus et al. 1986).As such, the heliospheric values of these quantities at 1 au (Smith et al. 2006;Weygand et al. 2011) are scaled accordingly with the ratio at 1 au of the AMF and HMF magnitudes B A and B H .The present study, as a first approach, assumes power-law radial dependences for these quantities, following Voyager observations in the heliosphere (Zank et al. 1996;Smith et al. 2001;Cuesta et al. 2022;Burger & McKee 2023).Such an approach to modeling turbulence quantities has been previously employed with success in heliospheric modulation studies (e.g., Engelbrecht & Moloto 2021).It should also be noted that the results of various turbulence transport models (e.g., Oughton et al. 2011;Zank et al. 2012Zank et al. , 2018) ) in the heliosphere support such power-law scalings for both magnetic variances and correlation scales, if pickup ion-driven wave formation (see Zank 1999;Isenberg et al. 2003) is neglected (Breech et al. 2008;Engelbrecht & Burger 2015).The scalings for the total magnetic variance and slab/2D correlation scales are with the slab correlation scale calculated using a constant factor observed in the heliosphere (Weygand et al. 2011).Slab and 2D variances are calculated from the total magnetic variance assuming an 80/20 2D/slab variance anisotropy (Bieber et al. 1994), although this has been known to vary in the heliosphere (e.g., Oughton et al. 2015).Furthermore, as a first approach, turbulence power spectra are assumed to consist of a wavenumber-independent energy-containing range and a Kolmogorov inertial range.These quantities are then employed in modeling diffusion coefficients as discussed in Section 2.2.1.It should be noted, however, that the above arguments are based on the behavior of transverse turbulence as observed and modeled for the supersonic solar wind.As such, these scalings do not take into account the observed predominance of compressive turbulence in the heliosheath (e.g., Burlaga et al. 2018;Fraternale et al. 2022).It would not be unreasonable to assume similar turbulence in the astrosheath of Proxima Centauri.However, it is still somewhat unclear how to model the spatial behavior of such turbulence, even in the heliosheath (see, e.g., Fichtner et al. 2020;Kleimann et al. 2023), and as such the approach outlined above is taken.

Modulation of GCRs within the Astrosphere of Proxima Centauri
The transport of CRs in the heliosphere is described by the Parker TPR, given by where f (r, p, t) = p −2 j T is the omnidirectional cosmic-ray (CR) distribution function, related to the CR differential intensity j T (Moraal 2013), and a function of position r, momentum p and time t.Diffusion and drifts due to gradients in, and the curvature of, the AMF, as well as due to the presence of a possible ACS, are modeled via the 3D CR diffusion tensor K, the outward convection of GCRs with the stellar wind by the term V sw • ∇f, and adiabatic energy changes with , with V sw the stellar wind speed (see, e.g., Quenby 1984;Jokipii 1989;McDonald 1998;Potgieter et al. 2001;Engelbrecht et al. 2022).Lastly, the term S represents any potential CR sources or sinks in the astrosphere.In heliospheric CR transport studies, this latter term is often used to model the Jovian source of low-energy electrons (e.g., Ferreira et al. 2001;Vogt et al. 2020).Presently, this term is neglected due to a lack of information as to whether such sources do, in fact, exist within Proxima Centauri's astrosphere.It is insightful to write the diffusion tensor K in AMF-aligned coordinates, such that (Burger et al. 2008) ,2 in terms of diffusion coefficients parallel and perpendicular to the AMF, which in turn can be written in terms of mean free paths (MFPs) via, e.g., κ ⊥ = vλ ⊥ /3, with v is the particle speed (Shalchi 2009).These quantities are discussed in more detail in Section 2.2.1.Off-diagonal elements represent drift coefficients, discussed in Section 2.2.2.The present study solves Equation (12) following the stochastic approach presented by Engelbrecht & Burger (2015), outlined in what follows.The Parker TPE can be written as a set of equivalent Itō stochastic differential equations (see, e.g., Zhang 1999aZhang , 1999b;;Strauss & Effenberger 2017;Moloto et al. 2019) satisfy a Weiner process, where η (t) ä (0, 1) represents a pseudo-random, Gaussian distributed number here generated using the Mersenne Twister algorithm (Matsumoto & Nishimura 1998).The above equations are solved time backward with an absorbing inner boundary condition.The physical motivation for this is that, at some point, particles at all pitch angles will be mirrored and no particles can actually reach the inner boundary r in = 0.005 au = 1, so that f (r in ) = 0.The way f changes from the boundary condition is determined by the physics in the TPE, and as long as the choice of boundary condition and transport coefficients are fairly realistic, and the correct TPE is used, the result should be reasonable.With only mirroring and no scattering, a highly anisotropic distribution will develop.However, even when a low level of scattering is included in a model the distribution quickly isotropizes.This has already been studied in the heliospheric context by Strauss et al. (2022), and for most cases, the isotropic (Parker) formalism (without mirroring but with absorbing boundary) compares reasonably well with the focused equation (with mirroring).It does not agree as well as the modified Parker equation with weak focusing (or mirroring), but does not completely deviate from the numerical solution.The distributions are also shown in that paper to remain fairly isotropic even for weak scattering conditions.Furthermore, similar to the heliosphere, the Parker TPE should be valid for most of Proximaʼs astrosphere, where most of the modulation of GCRs occurs.Mirroring is not included, as the assumption is made, again based on what is known from the heliosphere, that drift and diffusion would more strongly influence GCR transport than mirroring.Equations (14) are solved using the Euler-Maruyama scheme (see, e.g., Pei et al. 2010;Strauss & Effenberger 2017;Moloto et al. 2019).Here, the evolution of N = 10 4 pseudo-particles, each with a specified energy at an initial phase-space point in time ( ) x t , i o o , are followed iteratively until they reach their exit points and times, along with appropriately modified energies ( ) x t , i e e at some prespecified boundary.Each pseudo-particle represents a realization of the GCR distribution function, and therefore, a phase-space density element, and thus, should not be confused with actual particles (Strauss & Effenberger 2017).When this has occurred for all N pseudo-particles, an average omnidirectional GCR distribution function, and hence, the GCR differential intensity, can be calculated at the initial point using (e.g., Strauss et al. 2011) where f B the omnidirectional GCR distribution function at the boundary, which can be written in terms of the specified intensity of these particles at the boundary j B , whether this be a local interstellar spectrum (LIS) or a spectrum at a particular location in the astrosphere (e.g., the TS).Assuming a generalized, three-component AMF, it can be shown that tensor B i,j can be written in terms of the GCR diffusion coefficients (see, e.g., Pei et al. 2010;Engelbrecht & Burger 2015) which for a generalized 3D AMF can be calculated as (Burger et al. 2008) ,3 2 which depend on various angles pertinent to the geometry of the AMF, such as the AMF winding angle ψ, which can be calculated using ( ) The components for vector A are given, in terms of GCR kinetic energy E (Engelbrecht & Burger 2015) as Note that the signs of the drift V d and stellar wind velocities V sw are reversed, due to the time-backward nature of the numerical technique, and that E o denotes the rest-mass energy of the GCR species.For more information on the stochastic approach to solving the Parker TPR, the interested reader is invited to consult Kloeden & Platen (1999), Gardner (2009), and Strauss & Effenberger (2017).

Modeling Diffusion Coefficients
Currently, many scattering theories have been proposed to model the diffusion coefficients of GCRs parallel and perpendicular to some uniform background magnetic field, often with very different predictions for these quantities in terms of spatial and energy dependences, and are often mathematically intractable (Shalchi 2020;Engelbrecht et al. 2022).Given this uncertainty, the approach taken in the present study is to employ expressions for GCR parallel and perpendicular MFPs that have been used successfully in heliospheric GCR modulation studies (Engelbrecht et al. 2022).The parallel MFP expression employed here is derived from magnetostatic quasilinear theory and is given by Teufel & Schlickeiser (2003) as assuming a slab turbulence spectral form with a wavenumberindependent energy-containing range and an inertial range with spectral index −s.The quantity = R R k L min is the product of the wavenumber k min at which the inertial range commences and the maximal Larmor radius R L .Furthermore, dB sl 2 denotes the total slab variance, and B 0 the AMF magnitude.The perpendicular MFP is modeled using the nonlinear guiding center (NLGC) theory (Matthaeus et al. 2003), assuming a similar 2D turbulence spectral form as was the case for the parallel MFP, and is given by Shalchi et al. (2004) as with dB 2D 2 the total 2D variance, λ 2D the length scale at which the 2D spectral inertial range commences, and ν half the 2D inertial range spectral index.Based on a comparison of theory with numerical test-particle simulations, it is assumed that α 2 = 1/3 (Matthaeus et al. 2003).
The above MFPs depend on various turbulence quantities that cannot yet be ascertained for an astrosphere such as that of Proxima Centauri.However, these quantities can be modeled, and to some degree, estimated based on extensive observational evidence and modeling in the heliosphere.The present study employs simple power-law radial scalings based on heliospheric observations, as discussed in Section 2.1.2.

Modeling GCR Drift Effects
The components of the drift velocity are calculated following the approach of Burger (2012), which allows for the incorporation of drift effects due to gradients and curvatures of the AMF and due to the presence of an ACS.Although other techniques for implementing drift effects in GCR modulation codes have been proposed (see, e.g., Burger et al. 1985;Engelbrecht et al. 2019), this method has been demonstrated to work well for solar minimum conditions (e.g., Mohlolo et al. 2022).The drift velocity is calculated using , with e B a unit vector in the direction of the AMF.The components are given by Burger (2012) as where the change of sign of the AMF over the ACS is modeled using a hyperbolic tangent function, which in itself is a function of the assumed current sheet angle θ ns (Equation ( 8)).Although it has long been known that turbulence acts to reduce particle drift effects (Minnie et al. 2007;Engelbrecht et al. 2017), it is assumed that turbulence levels in the astrosphere of Proxima Centauri will be too low to significantly influence particle drifts.Therefore, the drift coefficient contained in the Parker TPR is modeled assuming weak scattering (see Forman et al. 1974), such that κ A = vR L /3.

Boundary Differential Intensity Spectra
The present study employs expressions for the local interstellar intensity spectra of GCR protons and helium, both applicable to the nearby region of the heliosphere, and thus, given the proximity of our subject, likely estimates for the intensities of these particles in Proxima Centauri's vicinity, given that both stars encounter local ISM conditions that are not expected to be vastly different due to their location in the Local Bubble (see, e.g., Frisch 1998;Maíz-Apellániz 2001;Linsky & Moebius 2023), implying a relative isotropization of GCR intensities (e.g., Gebauer et al. 2015).It should, however, be noted that such an assumption cannot necessarily be made for stars further afield (Jasinski et al. 2020) in the same units as the proton boundary spectrum.

Stellar Energetic Particle Modeling
For the stellar energetic particle (StEP) simulations presented in this study, a numerical transport model (Strauss & Fichtner 2015;Strauss et al. 2020) is used that was initially developed to simulate the transport of solar energetic particles close to the Sun.A spatially 2D version of the focused TPR is solved, where diffusion perpendicular to the mean magnetic field is included, • ( ) is the parallel streaming motion along the mean magnetic field, L is the magnetic focusing length, , D μμ is the pitch-angle diffusion coefficient with μ the cosine of the pitch angle, and ( ) D x includes spatial diffusion perpendicular to the mean magnetic field (Strauss et al. 2017).Note that Equation (27) does include mirroring in the focusing term.
The pitch-angle diffusion coefficient is related to the parallel MFP via (e.g., Bieber et al. 1994) where v is particle speed.Here, we calculate λ || from first principles and then use Equation (28) to calculate the corresponding D μμ, which is needed in the numerical model.Specifically, we assume the form (Dröge et al. 2010) of and use Equation (28) to calculate the value of D 0 , using the quasilinear theory expression for λ ∥ (see Section 2.2.1) employed in the GCR modulation model.Similarly, for the perpendicular diffusion coefficient, we assume the form of ,0 from which an isotropic perpendicular diffusion coefficient can be calculated as which is again related to the perpendicular MFP by λ ⊥ = 3κ ⊥ /v.We now specify λ ⊥ from first principles, using the NLGC expression for this quantity employed in the GCR modulation code (see Section 2.2.1) in order to calculate the appropriate value of D ⊥,0 following the approach outlined by Engelbrecht (2019), which is then included into the model.This ensures a self-consistent treatment of diffusion, depending on basic turbulence quantities, across the GCR and StEP transport modeling approaches presented here.The only free parameter, and the most difficult to determine from first principles, is the inner boundary condition that must be specified near the Sun or star.A proxy derived from remotesensing observations is necessary to specify the flux of StEPs injected from the acceleration region into the interplanetary medium.A pragmatic approach is implemented for this work, guided by the relevant physics.At the model's inner boundary r = r 0 , we specify the following boundary condition: which contains a Gaussian distribution in azimuthal angle, a pulse-like temporal profile for each StEP injection, δ(...), and C is an elusive normalization constant that needs to be computed post facto (Strauss & le Roux 2019).In this work, we inject a time series of StEP events based on TESS measurements (Vida et al. 2019).The model is solved for 10 MeV protons from which a spectrum is constructed (see details later in this section).These flare parameters are used to determine the time of each flare, τ i , and the flare energy for each event I i , where it is assumed that the StEP flux at the boundary is proportional to I i .The flare energy is therefore assumed to be a proxy for the StEP flux in the acceleration region (Steyn et al. 2020).We, therefore, end up with i = 1, 2, 3,K,72 StEP events occurring over ∼50 days.However, we do not have information about the position of the flares on the stellar disk, f i , or the width (size or broadness) of each event, σ i .We opted to generate these parameters by drawing a random number from a uniform distribution for f i and a Gaussian distribution for σ i , shown in Figure 3. Here, the blue data points show these generated quantities along with the injection amplitude of each event, I i , in the top panel.Lastly, we note that TESS observations only indicate flaring observed on the half of the stellar disk pointed to Earth.Statistically, if a full stellar disk is considered, the TESS flares should be double to account for occulted stellar activity.We do so by mirroring the original TESS time series in time and adding this second time series to the original to obtain 144 events over the ∼50 day period.These additional mirrored events are indicated in Figure 3 by the red data points.
The last step is to calculate the omnidirectional particle intensity and normalizing the intensity to a known value.We perform the normalization by running the model for the biggest TESS event, assuming a circumstellar source, and calculating the modeled j for a magnetically well-connected Prox Cen b.The strongest flare observed in the TESS measurements (Vida et al. 2019) utilized in this study had a bolometric energy of 1.74 × 10 32 erg.Assuming that a bolometric energy of 10 33 erg corresponds to a soft X-ray (SXR) energy of 1 × 10 −2 W m −2 (X100) (Benz 2017), this event had an SXR energy of 1.74 × 10 −3 W m −2 (X17.4), which is comparable to the flare that led to one of the strongest solar particle events ever measured as GLE event at the Earth's surface (i.e., GLE05 on 1956 February 23, X28 ± 14, see Papaioannou et al. 2023).Assuming that the physical principles for these comparably low-energy flares are the same for all cool stars, known solar relations, describing the currently known upper solar limits, can be utilized to derive stellar particle intensities at the exoplanetary orbit.Based on the findings in Papaioannou et al. (2023), a flare intensity of X17.4 corresponds to an event-integrated E > 10 MeV flux of 2.81 × 10 7 • (1.74 × 10 −3 W m −2 ) 1.40 = 3.85 × 10 3 pfu (#/(cm 2 s sr)) at a distance of 1 au, which is less than the actual GLE05 event (∼10 4 pfu).Scaling this value with 1/r 2 to the location of Prox Cen b (i.e., 0.04857 au) thus results in an upper limit of the eventintegrated E > 10 MeV proton flux of 1.632 × 10 6 pfu.
In the final step, the omnidirectional particle intensities are scaled according to 1.632 10 pfu , 3 4 scaled 6 with ( ) j max = 111.056(arb.units).Since the strongest observed flare intensity within the 52 day observation period is comparable to the one that led to the wellstudied GLE05, in a first approximation, we further utilize the GLE05 particle spectrum to derive plausible time-dependent stellar particle spectra at the location of Prox Cen b.Therefore, the spectrum is scaled according to j scaled /10 4 pfu.

Modeling the CR-induced Atmospheric Ionization and
Radiation Exposure of Prox Cen b Studying the impact of CRs on exoplanetary atmospheres is of utmost importance to understand chemistry and climate changes -and with that, changes in biosignatures like ozone and methane -and to get information on radiation effects on the surface of an exoplanet.Thereby, an enhancement of atmospheric ionization causes chemistry and climate changes due to the build-up of electromagnetic branches (i.e., electrons, positrons, etc.) induced by CRs, whereas an enhancement of the radiation exposure is correlated with the evolution of hadronic secondary particle branches (i.e., neutrons, protons, etc.).Both quantities can be modeled by utilizing the Atmospheric Radiation Interaction Simulator (AtRIS; Banjac et al. 2019).
Thereby, the CR-induced atmospheric ionization Q is given by with i representing the investigated primary CR particle type (here, protons and alpha particles) and J i (E) as the corresponding modulated differential primary particle flux at Prox Cen b.Y i (E, x) represents the ionization yield given by , with E ion as the mean specific energy loss at a certain atmospheric depth x per simulated primary particle of type i (here, protons and alpha particles) of energy E, while E ion , representing the average atmospheric ionization energy of an Earth-like atmosphere, is assumed to be 32 eV.The quantity E c reflects the planetary magnetic fielddependent cutoff energy, a measure of the ability of a particle of a certain energy to enter the planetary magnetic field at a given location.Since information on exoplanetary magnetic fields is currently missing in the literature, in this study, we assume a weak magnetic field, allowing primary particles with energies as low as 100 MeV to pass the exoplanetary magnetic field.
The precalculated relative ionization efficiency ( ) E R j i ,  , defined as the ratio between the average ionization energy a particle of type j is causing in a well-defined phantom like the ICRU phantom mimicking the human body (E d , see, e.g., McNair 1981) and E i , can further be used to model the average absorbed dose Dj given by


with m ph as the mass of the phantom given by • • r p r 4 3 ph 3 (see, e.g., Herbst et al. 2020a).Finally, the results are convoluted with the primary particle spectrum and summed up over all energy bins and particle types.

Results
The top panel of Figure 1 shows the MHD-modeled largescale plasma quantities used as inputs for the GCR modulation code (see Section 2.1.1 for details).Structurally, Proxima Centauri's astrosphere resembles the heliosphere, similar to previous studies (Herbst et al. 2020b;Mesquita et al. 2022).Within the TS, the stellar wind velocity has a maximum value of 1500 km s −1 .The stellar wind speed reduces drastically in the tail direction and the inner astrosheath past the TS.The TS distance in the nose direction is located at ∼55 au and varies as a function of latitude and longitude, such that in the tail direction, for example, it extends to ∼105 au.Proxima Centauri does not appear to generate a bow shock.This may be due to the single-fluid approach taken here and could change when a multi-fluid approach is applied (Herbst et al. 2020b).The AMF magnitude shown in the bottom half of the top panel of Figure 1 shows the AMF's 1/r 2 dependence, similar to that of a Parker HMF.The bottom left panel of Figure 1 shows GCR proton parallel and perpendicular MFPs, as well as the proton Larmor radius, as a function of rigidity at 0.0485 au, alongside the 1 au heliospheric Palmer (1982) consensus range values (see Section 2.2.1).MFPs are calculated using MHD model outputs for the large-scale AMF magnitude and power-law scalings for turbulence quantities as outlined in Section 2.1.2.The parallel MFP is considerably larger than is expected at Earth, the opposite being true of the perpendicular MFP due to the smaller AMF magnitude and the lower turbulence levels assumed here.The Larmor radius at this radial distance is considerably smaller than the MFPs, implying that drift effects may be less significant.However, given the decrease of the AMF magnitude at larger radial distances, this quantity will assume larger values.The bottom right panel of Figure 1 displays GCR proton intensities calculated at Prox Cen b for periods of positive (A > 0) and negative (A < 0) magnetic polarity, as well as at various radial distances, using a 3D stochastic CR modulation model (see Section 2.2), as a function of kinetic energy.Also shown are the assumed LIS (Burger et al. 2008) (LIS, green line), as well as typical observations of GCR intensities at Earth (McDonald et al. 1992).Direct comparison with the latter data shows that intensities calculated at the location of Prox Cen b greatly exceed intensities at Earth by almost an order of magnitude, in direct contrast to the results of previous studies performed using simpler GCR modulation models (Sadovski et al. 2018;Mesquita et al. 2021Mesquita et al. , 2022;;Rodgers-Lee et al. 2021a, 2021b).Drift effects play a relatively modest but significant effect, increasing intensities at lower energies by a factor of ∼1.8 from the no drift solution.Intensities for A > 0 (when positively charged particles drift inwards over the poles of the astrosphere) are marginally smaller than those for A < 0 (where particles drift inward along the ACS), in contrast with what is seen at Earth, due to the underwound AMF resulting from Proxima Centauri's slow rotation rate.This implies diffusiondominated transport for this astrosphere, similar to heliospheric periods of higher solar activity (Engelbrecht et al. 2022).At 30 and 50 au, intensities rapidly approach the LIS, implying that for this astrosphere, as opposed to the heliosphere, most GCR modulation occurs within the TS.This is due to a combination of factors: the considerably larger parallel MFP throughout this astrosphere due to the lower AMF magnitudes and turbulence levels, coupled with the highly underwound AMF relative to the HMF resulting from a combination of a considerably slower rotation rate for Proxima Centauri and a larger stellar wind speed, allowing for significantly easier access of GCRs into this astrosphere, which results in correspondingly large intensities.
To test this, and motivated by the fact that the majority of GCR modulation here appears to occur within the TS, analytical models for large-and small-scale plasma quantities (see Section 2.1.2) are employed as inputs for the 3D modulation code, with an LIS for both GCR protons and helium (Equations ( 25) and (26), respectively) placed at the TS boundary, the location of which is now modeled using the 3D location of the TS extracted from the MHD simulations shown in Figure 1.This means that the TS location as implemented here does not represent a spherical boundary, but rather one that varies with latitude and azimuth, such that its radial distance in the nose direction is ∼55 au, and in the tail direction it extends to ∼105 au.It should be noted that such an approach effectively neglects astrosheath modulation, even though it is known that modulation in the heliospheric heliosheath is considerable (see, e.g., Stone et al. 2013).Note also that changes in the radial location of the TS can have a strong influence on GCR intensities, as has been demonstrated in heliospheric studies (e.g., Moloto & Engelbrecht 2020).GCR intensities for both species calculated assuming different stellar rotation periods are shown in Figure 2, alongside Earth observations.GCR proton intensities calculated thus for Proxima Centauri's rotation rate closely resemble those calculated using the full MHD inputs, confirming that most modulation occurs within the TS.For both GCR species, intensities decrease due to decreasing rotation period for both full and no drift solutions.Full drift solutions remain larger than intensities at Earth, while no drift solutions drop below Earth intensities for shorter rotation rates, indicating an increase in the significance of this mechanism as stellar rotation periods decrease.Both these phenomena relate to the winding of the 3D AMF, as shorter rotation periods lead to a tighter winding of the Parker spiral field lines (Equation ( 7)).Correspondingly, this results in an increasingly azimuthal field, which influences the inward diffusion of GCRs, as these particles have to increasingly rely on perpendicular diffusion to reach Prox Cen b (Equation ( 17)).Given that these coefficients are significantly smaller than the parallel coefficient (see Figure 1), GCR intensities will decrease with decreasing rotation period accordingly.A strongly wound AMF will also increase the influence of gradient and curvature drifts relative to current sheet drift, as evinced by the fact that for the lower rotation periods, intensities calculated for A > 0 become larger than A < 0. As such, it is apparent that diffusion parallel to the AMF dominates transport processes in this astrosphere, as opposed to the heliosphere, with its more tightly wound magnetic field due to the Sun's faster rotation, where inward transport is dominated by perpendicular diffusion (Engelbrecht et al. 2022).
Figure 3 shows the input parameters to the StEPs model (see Section 2.3).The top panel shows the size (or magnitude) of the StEP event derived from TESS observations (Vida et al. 2019), while the scatter plots give the assumed position and size of the StEP source.The distribution of these input parameters is given in the bottom panel.The position of the StEP source on the stellar disk is assumed to be uniformly distributed, while the size of the source is assumed to follow a Gaussian distribution.Following the approach outlined above, an analytical AMF structure is assumed, while simulations are performed for 10 MeV protons using the same parallel and perpendicular MFPs described above.
The top panel of Figure 4 shows the resulting StEP spatial distribution at two different times after a simulated event.In both panels, colors show the omnidirectional intensity on a logarithmic scale, while white lines show examples of AMF lines, with Prox Cen b at the red circle, which moves through the computational domain with a period of 11.1 days.The event on the left was assumed to be much narrower than that on the right, indicating the large effect of the StEP source size on the resulting particle intensity observed at the planet's location.For the case on the left, the event is not seen on the planet.The relative position of the StEP source with respect to the planet, along with the size of the source itself, are major factors determining the resulting StEP intensity.The lower panels of Figure 4  The upper panels of Figure 5 show the altitude-dependent atmospheric ion-pair production rates (A.1) and absorbed dose rates (A.2) of Prox Cen b (see Section 2.4).The shown results are based on the modeled GCR intensities shown in Figure 2 assuming a stellar rotation period of 82.6 days (colored solid lines).Additionally, the results based on the widely used analytically derived approximation (Rodgers-Lee et al. 2020) (dashed lines), which are marked by a strong GCR modulation in the low-energy regime, are shown.The differences for the solutions computed for positive (A > 0, in red) and negative (A < 0, in blue) AMF polarities, as well as the no drift solutions (in black), manifest only at altitudes above 25 km.Thereby, the atmospheric quantities are highest during A < 0 and lowest for no drift conditions, showing negligibly small differences of up to 5%.However, the discrepancy between the results based on the often used analytic approach and the modeled values is much more severe.Switching from the nonphysical analytic solution to the more advanced numerical solution leads to differences of up to 60% in both quantities at altitudes above 10 km.
The event-induced impact on the abovementioned quantities is studied using the normalized time-dependent energy-integrated intensities shown in

Discussion
Due to the inherently 3D nature of the transport mechanisms influencing the modulation of charged, energetic particles in any astrosphere with a 3D AMF, it is necessary to model these processes three dimensionally.Although there are many uncertainties in such an endeavor, this study demonstrates that a careful combination of MHD modeling and extrapolations based on the heliospheric behavior of various plasma quantities can provide useful insights not possible when modeling in less than 3D.As such, we compute unexpectedly large GCR intensities at Prox Cen b, greater even than observations at Earth, implying that these particles cannot be ignored when considering the radiation environment of this exoplanet.Furthermore, the need for 3D transport modeling is emphasized by the importance of GCR drifts in their astrospheric transport and the unexpectedly important role of parallel diffusion reported here.This latter phenomenon is the reason why GCR intensities at Prox Cen b are so large, as the slow (relative to the Sun) rotation rate of Proxima Centauri, coupled with a relatively high stellar wind speed, leads to an underwound AMF that greatly facilitates the inward transport of GCRs.This The position of the source on the stellar disk is assumed to be uniformly distributed, while a Gaussian distribution for the broadness of the StEP sources is assumed.
is significant because such an underwound field would essentially also provide a highway for stellar energetic particles onto Prox Cen b.This raises the intriguing possibility that stellar rotation periods could provide a relatively simple-toobserve additional constraint on exoplanetary habitability.
This can be seen in the simulated StEP intensities arising from the TESS series.Although these display moderate-Sunlike-bolometric flare energies, within just 52 days, three separate GLE05-type events (comparable to the strongest event directly observed at Earth) are noted, with a correspondingly significant effect on atmospheric ionization and radiation dose.Note that Evryscope observations (Howard et al. 2018) have shown flares with bolometric energies of up to 10 33.5 erg (10 33.5 erg (∼3.16 × 10 −2 W m −2 , comparable to the AD774/775 event (Papaioannou et al. 2023), which would have resulted in an E > 10 MeV intensity of 1.41 × 10 8 pfu at Earth, and thus, an intensity of 5.9 × 10 10 pfu at Prox Cen b, more than 4 orders of magnitudes higher than the highest intensity used in this study.However, it is currently unclear whether known solar relations hold for SXR energies above 10 −2 W m −2 ; thus, the present study confines itself to TESS-like flares.
Accordingly, the results presented in Figures 5(A.1) and (A.2) allow for two conclusions: (a) GCRs cannot be neglected in chemistry/climate as well as habitability studies, and (b)-as suggested by previous 1D modeling efforts (Herbst et al. 2020b) -dedicated 3D modeling efforts are mandatory to derive a more realistic GCR-induced atmospheric production background.We find that previous studies employing the often used analytically derived GCR-induced background to derive, e.g., the radiation exposure at the exoplanetary surface overestimate the corresponding quantities by up to 50% (see panel (C.2)).Panel (B.2) further shows that the CR-induced ion-pair production rate values are orders of magnitudes higher than those induced by the same event at Earth (green band), which is also true for the induced absorbed dose rate values.Thus, careful modeling is needed to reliably reflect the CR-induced atmospheric radiation environments, particularly within exoplanetary atmospheres profoundly different from those we know from the solar system.
Future work will focus on model refinements for Proxima Centauri, and extension of the model to investigate other star systems.A multi-fluid approach to the MHD modeling of the large-scale astrospheric plasma quantities will be implemented, as this will allow for a more nuanced view of the physics influencing GCR transport (see, e.g., Kleimann et al. 2022;Sokółet al. 2022).Turbulence quantities, a key input for particle diffusion coefficients, can be modeled more self-consistently using outputs from turbulence transport models directly coupled to large-scale MHD models such as that used in this study.This has already been done with some success for the heliosphere (see, e.g., Usmanov et al. 2016Usmanov et al. , 2018;;Wiengarten et al. 2016;Kleimann et al. 2023), and can be reasonably expected to significantly influence the behavior of particle transport coefficients.Refinements can also be made in terms of modeling additional factors that are known from heliospheric studies to influence charged particle transport, such as the observed ∼7 yr stellar cycle for Proxima Centauri (Cincunegui et al. 2007;Wargelin et al. 2017), or the observed AMF tilt angles (Klein et al. 2021), which would significantly affect drift effects (see, e.g., Moloto & Engelbrecht 2020, for an analogous heliospheric study).Furthermore, taking into account coronal mass ejectioninduced StEP particle propagation would further refine estimates of the influence of these particles on Prox Cen b, as flare-induced events may only provide a lower limit for expected intensities.Lastly, the model suite presented here can also be applied to other exoplanet-hosting astrospheres, providing unique and potentially valuable insights as to their habitability.

Figure 1 .
Figure 1.The top panel shows the magnetohydrodynamically modeled stellar wind velocity (top half) and magnetic field magnitude (bottom half) within the equatorial plane.The ISM is flowing in from the right-hand side.The bottom left panel shows GCR parallel and perpendicular MFPs and Larmor radii as a function of rigidity at 0.05 au.Heliospheric Palmer (1982) consensus ranges are also shown to guide the eye.The bottom right panel shows GCR proton intensities computed at 0.05 au, the location of Prox Cen b, using inputs from the full MHD model.Red, blue, and black solid lines denote solutions computed for positive (A > 0) and negative (A < 0) AMF polarities, as well as no drift solutions, respectively.Also shown are the assumed local interstellar spectrum (Burger et al. 2008) (green line), as well as differential intensity spectra computed at 30 au (black dotted-dashed line) and 50 au (black dashed line).GCR intensities observed at Earth for both HMF polarities (McDonald et al. 1992) are shown as black and white circles to guide the eye.
again show the time-dependent StEP injection profile assumed in the model, while the bottom panel now also shows the normalized energy-integrated intensities (see Section 2.3) as a function of time, which are used to compute the impact of StEPS on the atmosphere of Prox Cen b.
Figure 4.As a first step, the time-dependent impact of the strongest flare-induced event (see panel (B.1)), which occurred on day 22.79, is computed.Panel (B.2) shows the event-induced atmospheric ionization variations that have been summed up with the modeled GCR-induced A < 0 background.

Figure 2 .
Figure 2. GCR proton and helium intensities computed at 0.05 au, the location of Prox Cen b, using analytical models for the large-scale astrospheric plasma quantities, and a modulation boundary set at the location of the astrospheric TS as calculated from MHD, for varying stellar rotation periods as indicated in the legend.The top panels display proton intensities, while the bottom panels display helium intensities.The left panels show solutions computed with drift effects, and the right panels show solutions where drift effects are switched off.Also shown are the assumed local interstellar spectra (green lines) for protons (Burger et al. 2008) and helium (Boschini et al. 2020).GCR proton and helium intensities observed at Earth (McDonald et al. 1992; Marcelli et al. 2020) are shown as black and white circles and blue triangles, respectively, to guide the eye.

Figure 3 .
Figure 3.The top panel shows the assumed injections of StEPs into the numerical model as derived from TESS observations.Blue data represent the original data, while red is the mirrored data set.The position and broadness of the injected sources are shown, along with the distribution of all parameters at the bottom in green.The position of the source on the stellar disk is assumed to be uniformly distributed, while a Gaussian distribution for the broadness of the StEP sources is assumed.

Figure 4 .
Figure 4.The top panel shows two snapshots of the simulation roughly 74 hr (left) and 430 hr (right) from the start of the simulation.The color scale shows the normalized intensity on a logarithmic scale.The red circle indicates the position of Prox Cen b, while the white lines show example magnetic field lines.The middle panel again shows the assumed injection profile in the model, while the bottom panel shows the simulated intensity at Prox Cen b.

Figure 5 .
Figure5.Upper panels: background production of the altitude-dependent ion-pair production rates (A.1) and absorbed dose rates (A.2).Solid lines correspond to modeled GCR spectra, dashed ones to the analytically derived spectra.Middle panels: time profile of the modeled event corresponding to the strongest flare of the TESS sample occurring on day 22.79 (B.1) and the corresponding induced altitude-dependent ion-pair production rates (B.2) compared to the impact of GLE05 at Earth (in green).Lower panels: summed up GCR-and event-induced ion-pair production rates (C.1) and absorbed dose rates (C.2) at the surface of Prox Cen b.Here, purple lines represent values based on the modeled GCR background, while black lines are those based on the analytically derived GCR background.
and x i (t) representing Itō processes.