Stellar Energetic Particle Transport in the Turbulent and CME-disrupted Stellar Wind of AU~Microscopii

Energetic particles emitted by active stars are likely to propagate in astrospheric magnetized plasma turbulent and disrupted by the prior passage of energetic Coronal Mass Ejections (CMEs). We carried out test-particle simulations of $\sim$ GeV protons produced at a variety of distances from the M1Ve star AU~Microscopii by coronal flares or travelling shocks. Particles are propagated within the large-scale quiescent three-dimensional magnetic field and stellar wind reconstructed from measured magnetograms, and { within the same stellar environment following passage of a $10^{36}$~erg kinetic energy CME}. In both cases, magnetic fluctuations with an isotropic power spectrum are overlayed onto the large scale stellar magnetic field and particle propagation out to the two innnermost confirmed planets is examined. In the quiescent case, the magnetic field concentrates the particles onto two regions near the ecliptic plane. After the passage of the CME, the closed field lines remain inflated and the re-shuffled magnetic field remains highly compressed, shrinking the scattering mean free path of the particles. In the direction of propagation of the CME-lobes the subsequent EP flux is suppressed. Even for a CME front propagating out of the ecliptic plane, the EP flux along the planetary orbits highly fluctuates and peaks at $\sim 2 -3$ orders of magnitude higher than the average solar value at Earth, both in the quiescent and the post-CME cases.


INTRODUCTION
Active low mass K-and M-star circumstellar environments are affected with an occurrence rate much higher than Solar by violent eruptions producing either very energetic coronal flares (Youngblood et al. 2017;Jackman et al. 2020) or possibly escaping Coronal Mass Ejections (CMEs, > 10 31 erg kinetic energy) detected, e.g., via X-ray spectroscopy (Argiroffi et al. 2019); CME candidates are also traced via Doppler shift in Balmer lines (Houdebine et al. 1990) or asymmetries therein (Vida et al. 2019), continuous X-ray absorption during the flare  or dimming in the extreme ultraviolet and X-ray due to the CME mass loss (Veronig et al. 2021). The broadband flare emission (from radio to γ-rays), hence the bolometric detectable energy output, from such stars is routinely investigated (e.g., Paudel et al. 2021) whereas the kinetic energy of the associated CMEs has been estimated only in a handful of cases (e.g., Moschou et al. 2019). Within the heliosphere, the charged particle acceleration at CME-driven shocks has been accurately determined via in-situ measurements to drain ∼ 10% of the total CME energy, regardless of the magnetic obliquity at the shock (David et al. 2022). Comparable energy fractions might be expected at active stars.
The passage of a CME compresses and breaks magnetic field lines leading to a re-arrangement of the large-scale magnetic field topology throughout the astrosphere, from the corona to the interplanetary region, that is traversed by charged particles energized close to the star. Such a disrupted configuration of the stellar wind is more likely to be encountered by outward propagating energetic particles (hereafter EPs) from active stars due to a flaring rate much higher than solar (Youngblood et al. 2017). The flux of EPs onto habitable zone (hereafter HZ) planets in the quiescent winds of active stars was first determined numerically for the case of TRAPPIST-1 (Fraschetti et al. 2019). The flux exceeded the solar value by ∼ 4 orders of magnitude. However, the passage of a very energetic CME is expected to re-shuffle the wind magnetic field over large angular regions out to large distances. To our knowledge, the effect of such a phenomenon has not yet been investigated.
The James Webb Space Telescope (JWST) is expected to open up new pathways toward the observational studies of exoplanet habitability and atmosphere composition and evolution. In particular, exoplanets with radii between 1.7 and 3.5 times the Earth radius (i.e., sub-Neptunes) are favorable targets for atmospheric observations instead of smaller planets as the larger amount of atmospheric H 2 acts as greenhouse gas allowing for stable liquid-water (Pierrehumbert & Gaidos 2011;Hu et al. 2021).
We focus here on AU Microscopii, an M dwarf with flaring activity observed by, e.g., the Hubble Space Telescope (HST) in far ultraviolet (Redfield et al. 2002) or XMM Newton in X-rays (Magee et al. 2003) and a modelled connection between flares (Extreme Ultraviolet Explorer, Cully et al. 1994) and ejected plasmoids self-similarly expanding in a CME fashion. The confirmation of two sub-Neptunian planets orbiting AU Mic (Martioli et al. 2021) makes the system particularly attractive for investigating the effects of the CME passage on the EP propagation from star to planet due to their impact on planet atmosphere and its evaporation (Fulton et al. 2017).
The diffusive transport of EPs originating from solar eruptions is known to be governed by the unperturbed large-scale magnetic field and by its small-scale fluctuations (Jokipii 1966). Existing numerical analyses of the propagation of EPs from young stars surrounded by proto-planetary disks (Rodgers-Lee et al. 2017;Rab et al. 2017;Fraschetti et al. 2018;Padovani et al. 2018;Gaches & Offner 2018) or in exoplanetary environments (Fraschetti et al. 2019) have focused on quiescent stellar conditions. Particle transport is determined by integrating EP trajectories in synthetic 3-dimensional turbulence (Fraschetti et al. 2019) or by solving a suitable transport equation, as done recently in Hu et al. (2022). However, as mentioned above, EP propagation into a realistic astrosphere disrupted by a recent (within 1 − 2 hours) CME passage does not appear to have been discussed previously.
In this paper, we perform a detailed analysis of the propagation of charged particles energized in proximity of AU Mic, i.e., by flares or CME-shocks, through the magnetized stellar wind calculated via the Space Weather Modeling Framework (SWMF) codes, in particular the Alfvén Wave Solar Model (AWSoM, van der Holst et al. 2014), out to the second confirmed planet. Synthetic turbulent magnetic field is added to the large scale unperturbed component (Fraschetti et al. 2019). The propagation within the quiescent state astrosphere is compared with the propagation 90 minutes after the passage of a very energetic CME; a kinetic energy consistent with the best candidate event observed in this star so far (∼ 10 36 erg) is adopted (Katsova et al. 1999;Alvarado-Gómez et al. 2022)).
The outline of this paper follows: in Sect. 2 the observational properties of the AM Mic planetary system are summarized; in Sect.3 the assumptions on the stellar EP origin and propagation properties are emphasized along with the generated magnetic turbulence with intensity and injection scale as parameters. In Sect.4 the main results are presented for the cases of winds in quiescent state (with particles injected as close as the lower corona) and post-CME state. In Sect.5 the EP fluxes impinging on the planets AU Mic b and c for the quiescent and post-CME case are compared; the EPs transport properties for AU Mic and TRAPPIST-1 are also compared; Sect.6 draws the conclusions of this work.

THE LARGE-SCALE MAGNETIZED WIND OF AU MIC
AU Microscopii is a bright, nearby (magnitude 1 = 8.6, d=9.72 ± 0.04 pc, Gaia Collaboration et al. 2018) M dwarf with mass M = 0.5M , radius R = 0.75 R = 5.18 × 10 10 cm (Plavchan et al. 2020) and a rotation period of 4.85 days (Torres et al. 1972). A spatiallyresolved edge-on debris disk surrounds the star with a ∼ 50 au inner radius (Kalas et al. 2004). Transiting Exoplanet Survey Satellite (TESS) light-curves confirmed that AU Mic hosts a Neptune-size planet (AU Mic b, Plavchan et al. 2020) with radius 1.05 R N = 2.6 × 10 9 cm (where R N is the Neptune radius), mass 1.00 M N (with M N is the Neptune mass), orbital distance R b = 0.065 au = 19.1 R = 9.89 × 10 11 cm (orbits are assumed to be circular Martioli et al. 2021), orbital period 8.46 days (Martioli et al. 2021;Klein et al. 2021b), an inclination angle of magnetic field/stellar rotation axis ∼ 19 • (Klein et al. 2021b), and an uncertain alignment between the spin axis of the host star and the orbital vector axis of the planet (Addison et al. 2021). TESS revealed also a second orbiting Neptune-size planet (AU Mic c) with radius 0.84 R N = 2.07 × 10 9 cm, mass 0.13 M N < M c < 1.46 M N , orbital distance R c = 0.11 au = 29. R = 1.5 × 10 12 cm, orbital period 18.8 days (Martioli et al. 2021). The planetary orbital plane for both planets is located within 1 degree from the stellar equator (Martioli et al. 2021).
The 3D magnetized quiescent Stellar Wind (hereafter SW) was computed using the Space Weather Modeling Framework codes (Tóth et al. 2005;van der Holst et al. 2014;Gombosi et al. 2018), evolved from the BATS-R-US MHD code (Powell et al. 1999) that was originally developed for the solar corona. The code uses as inner boundary condition a magnetogram (details in Alvarado-Gómez et al. 2022) describing the surface distribution of the radial magnetic field in the quiescent state and calculates the coronal heating and SW acceleration due to Alfvén wave turbulence dissipation, taking into account radiative cooling and electron heat conduction. The model has been validated with solar wind observations (Cohen et al. 2008); improvement of wave dissipation to electron and anisotropic proton heating lead to good agreement with magnetic field measured by Parker Solar Probe (van der Holst et al. 2022). Further validations have been based on remote observations in solar minimum (Sachdeva et al. 2019) and maximum (Sachdeva et al. 2021) conditions. The code has been adapted and used to simulate SWs and the space weather environments of exoplanets (e.g., Vidotto et al. 2015;Garraffo et al. 2017;Cohen et al. 2020;Evensberget et al. 2021), and has also been used to simulate CME eruptions from highly-magnetized stars (e.g., Cohen et al. 2011;Alvarado-Gómez et al. 2018 or the associated radio emission in Eridani (Ó Fionnagáin et al. 2022).
Stellar surface magnetic field distributions needed to drive stellar wind simulations have generally been based on Zeeman-Doppler Imaging (ZDI) observations. The AU Mic large-scale B-field derived in such a way suffers some uncertainties. Kochukhov & Reiners (2020) have shown that the use of distinct polarizations (circular and linear) from ESPaDOnS and HARPSpol instruments in the ZDI map leads to magnetic field strengths differing by a factor 10 (184 G and 2 kG, respectively). Using SPIRou polarization observations, Klein et al. (2021b) found a 450 G large-scale dipole field inclined by ∼ 20 • with respect to the rotation axis. Both maps were used in a companion paper (Alvarado-Gómez et al. 2022) as an inner boundary to generate the 3D magnetized SW of AU Mic, both in quiescent and CME-disrupted phase; in the present work, we have used the B-field produced in the cases 1 and 3 therein. Cohen et al. (2022), using the same wind reconstruction, have determined the variations in Lyα absorption signatures during transits of AU Mic b due to the passage of a very energetic CME. The Klein et al. (2021b) maps are also implemented in Kavanagh et al. (2021) to generate via AWSOM the AU Mic wind with 2 distinct mass loss rates, that lead to different radio emission.
In this paper, we calculate the propagation of stellar EPs within the turbulent magnetized SW of AU Mic driven via the ZDI maps from Klein et al. (2021b) and Kochukhov & Reiners (2020). In addition to the quiescent SW, we have produced a number of SW configurations disrupted by very energetic CMEs ( kinetic energy ∼ 10 36 erg, Alvarado-Gómez et al. (2022) propagating throughout the entire simulation box. The CME structure is initialized by using the Titov & Démoulin (1999) flux-rope eruption model over the AWSoM field background. We consider herein one such configuration in detail: stellar EPs are propagated within SW snapshots 90 minutes after the CME initialization, i.e., when the CME front has reached > 100 R region . EPs travel much faster than the CME front and the choice of 90 minutes ensures that the astrosphere has been disrupted by the CME as far as the numerical calculation allows.

STELLAR ENERGETIC PARTICLES IN THE
TURBULENT ENVIRONMENT OF AU MIC 3.1. Assumptions for EPs: origin, propagation and abundance The goal of this work is to compare the transport of EPs in a quiescent SW environment of a young active star with the transport within the same environment 90 minutes after the passage of an extraordinarily energetic CME (compared with heliospheric scales).
In the solar context, GOES measurements of proton enhancements at 1 AU confirm that EPs might originate both from SXR flares and CME-driven shocks (Belov et al. 2007). Guided by the heliospheric observations, we assume that young active stars produce EPs via two distinct processes: 1) CME-driven shocks, travelling through the interplanetary medium and therein accelerating EPs, a certain fraction of which are likely to escape the shock at various distances from the host star; 2) coronal flares that release EPs locally accelerated close to the stellar surface. Presumably via different mechanisms, namely diffusive shock acceleration for the former and magnetic reconnection for the latter, both processes contribute multi-MeV to ∼ GeV kinetic energy protons in the heliosphere.
We inject EPs on spherical surfaces with radius R s = 2 R and 5 R concentric with the star, to compare the effect of the injection at two locations where the relative amount of closed-to-open field lines changes with a large spatial gradient; the diffusion in that region determines the EP flux at the planets. Following the approach in Fraschetti et al. (2019), we calculate the time-forward propagation of test particles by using two distinct magnetostatic SW configurations: 1) quiescent interplanetary magnetic field; 2) stellar magnetic field 90 minutes after the initiation of a very energetic CME; in both cases we include the same overlapping small scale turbulence (see Sect. 4.4).
The modelling of TESS flaring rate for AU Mic (Gilbert et al. 2021) is consistent with 1 flare every 3.8 hours for a flux between 0.06% and 1.5% of the stellar flux (Martioli et al. 2021). For the AU Mic bolometric luminosity of 0.09 L (Plavchan et al. 2009), these correspond to fluxes at AU Mic b of 1.8×10 4 erg cm −2 s −1 sr −1 and 4.4 × 10 5 erg cm −2 s −1 sr −1 and, possibly, to very energetic associated CMEs. Thus, the 90-minutes interval elapsed since the CME initialization allows the CME front to reach R > 100 R before the eruption of a subsequent energetic CME, consistently with observations (Gilbert et al. 2021). On the other hand, due to the high flaring occurrence in active stars, i.e., a likely high rate of associated CMEs, a wind configuration disrupted by a CME has a higher filling factor than for the Sun; therefore a significant fraction of EPs originating from the host star encounters typically a non-quiescent wind. In contrast, due to the lower level of solar activity, in the heliosphere most of the EP transport occurs within a quiescent, rather than a CME-disrupted, wind.
The magnetostatic approximation adopted herein is justified as follows. The MHD wind solution and the magnetic turbulence are stationary on the time-scale of EP propagation to a good approximation. The EPs travel close to the speed of light whereas the stellar rotation period of 4.86 d and radius 0.75 R (Klein et al. 2021b) imply a surface stellar rotation speed of a few km s −1 ; the Alfvén speed in the circumstellar is at most a few thousand km s −1 throughout the simulation box, which is much smaller than the EPs speed.
The EP abundance in the circumstellar medium at a given distance from the host star cannot be constrained through direct observation; instead, we use the estimate based on solar scaling relations between EP fluence and far-UV and SXR fluence during flares by Youngblood et al. (2017) ; see Sect.5 below. This scaling provides a time-averaged EP enrichment for time scales comparable with a statistically typical flare duration (Vida et al. 2017). We have verified that a total number of EPs N inj = 10, 240 yields numerical convergence in all cases presented herein.
Finally, we note that the AU Mic detected debris disk is not expected to impact the EP transport as the inner radius of the disk is measured to be 50 au (Kalas et al. 2004), which is much greater than R c .

Turbulent stellar magnetic field
Leveraging the universality of the Kolmogorov scaling within the turbulent inertial range (Armstrong et al. 1995), we assume that the magnetic fluctuations around AU Mic support a 3D Kolmogorov isotropic power spectrum (see also Fraschetti et al. 2019). Spectral analysis of Parker Solar Probe (PSP) measurements near the minimum of the solar cycle 25 in the quiescent inner heliosphere (as close as 0.2 AU to the Sun) have shown (Zhao et al. 2020) that the power spectrum of the magnetic turbulence in the direction aligned with the magnetic field is consistent with a Kolmogorov spectral slope of −5/3 and in tension with the −2 slope predicted by the critical balance conjecture (Goldreich & Sridhar 1995); in addition, the perpendicular transport in the Goldreich & Sridhar (1995) was found to be inefficient in the perpendicular diffusion of fast particles (Fraschetti 2016a,b). The spectral slope found by PSP confirms previous findings at 1 AU from the Wind spacecraft, restricted to fast solar wind (Telloni et al. 2019), and a number earlier analyses (e.g., Jokipii & Coleman 1968).
The total magnetic field is decomposed as where the large-scale component B 0 (x) is the 3D magnetic field generated by the 3D-MHD wind simulations (see Section 2); the random component δB = δB(x, y, z) has a zero mean ( δB(x) = 0). As for the turbulent environment of TRAPPIST-1 (Fraschetti et al. 2019), the fluctuation δB(x, y, z) is calculated as the sum of plane waves with random orientation, polarization, and phase following the prescription in Giacalone & Jokipii (1999); Fraschetti & Giacalone (2012), with an inertial range k min < k < k max , with k max /k min = 10 2 , where k max is the magnitude of the wavenumber corresponding to a turbulence dissipation scale.
The advantage of the test-particle approach used here is that particle trajectories enable tracking of the pitchangle scattering, of the perpendicular diffusion, and also of the transport across field-lines both in the quiescent and CME-disrupted winds; such effects are known to contribute significantly to particle transport in the heliosphere (e.g., Dröge et al. 2010;Fraschetti & Jokipii 2011;Gómez-Herrero et al. 2015) but are often neglected for analytic tractability. Moreover, a 1D-spatial transport equation approach cannot be applied to a wind disrupted by the CME passage as the radial scaling of the diffusion coefficient, κ, typically inferred from the radial scaling of the large scale magnetic field, is reshuffled in the post-CME wind: the expected strong angular dependence of κ cannot be included in a semi-analytic model.
A second parameter of the stellar wind magnetic turbulence is the correlation length L c , i.e., the outer scale of turbulence injection. Due to the lack of observational constraints on L c , and the likely observational inaccessibility to L c in the near future, in Fraschetti et al. (2019) we used for TRAPPIST-1 a range of values of L c , each one kept uniform throughout the simulation box, and found no significant difference in the spatial distribution of EPs at the distances corresponding to the semi-major axes of the planets in that system. Likewise, for AU Mic we adopt here the uniform value L c = 10 −5 AU throughout the simulation box. Such a value warrants that the resonant scattering condition holds with good approximation during the EP propagation throughout the wind. The resonance condition reads kr g (x)/2π = r g (x)/L c < 1 for each wave-number k within the inertial range; here, r g (x) = p ⊥ c/eB 0 (x) is the gyroradius of a proton with p ⊥ momentum perpendicular to the unperturbed and space-dependent magnetic field B 0 (x), e is the proton electric charge and c is the speed of light in vacuum. As for the case of TRAPPIST-1, the combined effect of a high surface stellar magnetic field strength and its decrease with radius make L c = 10 −5 AU a reasonable value within the assumed circular orbits of AU Mic b and c, for the particle energies considered.
The power of the magnetic fluctuation δB(x) relative to B 0 (x) is defined as (1) Given the current lack of any observational constraint of the magnetic turbulence around AU Mic, it seems reasonable to assume a uniform σ 2 , following Fraschetti et al. (2018). Here, σ 2 is assumed to be independent of space throughout the simulation box. The solar wind measurements between 0.3 and 4 AU yield for the turbulence amplitude δB a power-law dependence on heliocentric distance with a comparable slope (−2.2) at a variety of helio-latitudes (Horbury & Tsurutani 2001). In the steady-state reconstructed 3D magnetic fields used here (see Sect.2), the spherical average of the unperturbed field B 0 (x) Ω drops with radius R as ∼ R −2.2 (see also the case of TRAPPIST-1 in Fraschetti et al. 2019). The high anisotropy of the post-CME stellar wind due to the CME eruption causes deviations from the monotonic scaling of B 0 (x) Ω , but it is conceivable that the level of small-scale turbulence is not significantly altered (see Sect.4.4 and Kilpua et al. (2021)). The turbulence within the young and active M dwarf magnetosphere is likely to be much stronger than in the solar wind (σ 2 0.1, Burlaga & Turner 1976), hence the broader σ 2 -range 0.01 − 1.0 is spanned here. The interpretation of our simulations makes use of the scattering mean free path, λ , given by quasi-linear theory (Jokipii 1966), that reads (Giacalone & Jokipii 1999;Fraschetti et al. 2018) The choices of uniform L c and σ 2 imply that λ depends on spatial coordinates only via r g (x), i.e., B 0 (x). Fitz Axen et al. (2021) investigated the transport of < 1 GeV protons within protostellar cores by implementing scattering off magnetic turbulence via a Monte Carlo algorithm that neglects perpendicular transport; thus, cross-field diffusion and consequent longitudinal spread cannot be incorporated. It is noteworthy that assigning L c and σ 2 for a given particle energy (namely λ ) does not uniquely describe the particle transport due to the increase of the perpendicular diffusion as σ 2 increases (Giacalone & Jokipii 1999;Fraschetti & Giacalone 2012;Dröge et al. 2016). Therefore the EP trajectory needs to be calculated step by step in the given magnetic turbulence without assuming that EPs follow the magnetic field lines.

RESULTS
Here, we discuss the results from injecting 0.1 and 1 GeV kinetic energy protons at R s = 2R = 0.0069 au and at R s = 5R = 0.017 au for σ 2 = 0.01, 1.0. EPs are propagated in our simulations throughout the astrosphere until either they collapse back to the star or hit (for the first time along their trajectory) the spherical surface at R = R b and R = R c . A small fraction of EPs hit at first the R b -sphere at a latitude different from the geometrical cross-section of the planetary orbit, then backscatter and in their subsequent star-ward propagation hit the R b -(R c -) sphere again at the latitude of the planet orbital plane. Such a fraction is < 1% at most, so is neglected here and EP trajectories are followed within the region R < R c .

Quiescent stellar wind
In Fig. 1 open/closed magnetic field lines are seeded through the orbits of planets b and c. The unperturbed magnetic field lines that approximately track the motion of EPs have a predominantly dipolar structure (Alvarado-Gómez et al. 2022).
For the quiescent stellar wind numerically reconstructed from the ZDI map from Klein et al. (2021b), Figure 2 shows the 2D histogram in spherical coordinates at two distinct radii (R = R b and R c ) of the first-crossing points for 1 GeV kinetic energy protons in the case of strong turbulence, i.e., σ 2 = 1, injected at R s = 5R . The corresponding stellar wind magnetic field is shown in Fig.1.
Depleted regions (in blue) reached by no EPs are found at both distances R b and R c . Such regions broaden progressively as the distance from the host star increases. The azimuthal oscillation of the depleted regions maps the slow speed wind and the stellar current sheet, both shown in 2D spherical projections of the flow speed and the magnetic field strength in Fig. 3; similar correspondence between the EP 2D histogram and the current sheet was found in the case σ 2 = 1 for the HZ of the TRAPPIST-1 system (Fraschetti et al. 2019). As for TRAPPIST-1, the depleted regions result from the perpendicular transport, enhanced in the case of σ 2 = 1, of EPs at the boundary between open and closed field lines that favours a net migration of EPs across the large scale magnetic field from open to closed lines, due to the larger scattering mean free path in the weaker B-field of that region (see Eq.2): after transferring onto a closed line, EPs precipitate in a nearly scatter-free regime along the closed lines toward the star surface. Due to the larger B-field (smaller mean free path at fixed σ 2 ) in the open lines region, the number of EPs migrating via perpendicular transport in the opposite direction, from closed to open, is smaller, as some can backscatter and return to a closed line. Along the open lines the EPs proceed in an outward trajectory toward the planet. The effect of a weaker turbulence (σ 2 = 0.01) is discussed below in this section.
As found in the case of TRAPPIST-1 (Fraschetti et al. 2019), for AU-Mic the EP-depleted regions are explained by the combined effect enhanced perpendicular diffusion and B-dependence of λ . The planet AU-Mic b is further out from the host star than TRAPPIST-1e, the closest HZ planet therein, i.e., 0.066 au compared with 0.029 au, but closer in units of stellar radius, 19 R compared with 52 R ; however, such differences do not lead to significant differences in the 2D histogram. Figure 4 shows that, if EPs are injected further in (R s = 2 R ), the larger fraction of closed-to-open magnetic field lines in the inner wind increases the likelihood for EPs to be captured by the closed lines, hence a smaller EP flux at the planet. The upper panels show the decrease with radius of EPs due to the trapping by closed lines.
The combination of open field lines and strong turbulence focusses EPs into caps far out, allowing to reach the planetary ecliptic. These caps track the high B-field projected region (see Fig.3). The high B-field shrinks the particle gyroscale, keeping them confined within a limited angle defining the caps. The magnetic fieldrotation axis inclination angle of 19 • (Klein et al. 2021b) causes the caps to be tilted toward the equatorial plane. Likewise, a magnetic field/rotation axis inclination angle of ∼ 40 • in TRAPPIST-1 focusses EP caps toward the equatorial plane (Fraschetti et al. 2019)]. In case of alignment of the B-field and rotation axes, the projection of the current sheet would appear closer to an horizontal stripe in Fig.3 and the EP caps would be expected to be closer to the polar region.
The caps cross the equatorial plane, i.e., planetary orbit (Plavchan et al. 2020), implying a modulation in the bombardment of the planet and in the consequent atmospheric ionization rate. Likewise, a modulation of the EP flux at the HZ planet TRAPPIST-1e was found to be up to ∼ 4 − 5 orders of magnitude greater than experienced by Earth (Fraschetti et al. 2019), and its    implications for the EP penetration depth were outlined in (Fraschetti et al. 2021).
Also of interest is the timescale of particle modulation due to orbital motion and stellar rotation. The rotation period of AU Mic is shorter than the orbital period of the planets (8.46 days for AU Mic b Klein et al. 2021b). The stellar rotation relative to the orbital motion will then sweep the EP caps over the planet with an effective period of approximately 11 days. The change in EP flux from EP-depleted to EP-enhanced regions occurs over an azimuthal angle range of a few 10s of degrees, such that the EP flux variation timescale would be of the order of a day. This is relevant for the recovery timescale of a planetary atmosphere to EP ionization events, and whether or not the atmosphere would be in a perturbed equilibrium state or subject to strong secular variation (Herbst et al. 2019;Chen et al. 2021).
The azimuthal variation of the EP flux impinging on the planet along its orbit around the star also strongly depends on the strength of the magnetic turbulence, as can be seen by comparing the case of strong (Fig.2) and weak turbulence (Fig.5). In the former case, the planet crosses two discernible caps, whereas in the latter the planet orbits into an azimuthally homogeneous , but much more sparse, distribution of EPs with far lower EP flux, evident also from the uniform color distribution in the latter. In the weak turbulence case, the distribution is closer to homogeneity as a result of the homogeneous injection on the R s -sphere and the boundary between closed and open field lines does not act as an attractor for EPs toward the stellar surface as in the case σ 2 = 1; a comparably homogeneous distribution was found in the case of the TRAPPIST-1 system (Fraschetti et al. 2019). The difference between the distributions in Figs. 5 and 2 emphasizes the role of the diffusion in the direction perpendicular to the large-scale unperturbed field that is not accounted for in a purely Monte-Carlo approach to scattering off the magnetic fluctuations.

Energy dependence of particle propagation
The spatial distribution of EPs is fairly independent of particle energy, as shown by comparing the 2D histograms for 0.1 GeV EPs (see Fig. 6) with 1 GeV (see Fig.2) at two distinct radii, for the same injection radius R s = 5R . As a consequence an energy-dependence of the diffusion coefficient does not alter the EPs 2D histograms. A comparison of the EP energy spectrum at various astrospheric distances requires spanning a wide range of EP energies, in order to include the effect of the perpendicular transport that solar in-situ measurements suggest contribute significantly to circumsolar events (Fraschetti & Jokipii 2011;Gómez-Herrero et al. 2015).
This work in hand focuses on the spatial distribution of EPs throughout the astrosphere to investigate the effect of the magnetic connection source-planet on the EP propagation. Transport might steepen the momentum spectrum of EPs at high energy, as shown with a pure scattering model with no perpendicular diffusion (Li & Lee 2015); however, spectral modifications are not investigated herein because the source of EPs is not localized to an individual shock with specified parameters, i.e., fixed spectral shape.

EP injection by flares in the stellar corona
Although the structure of a flaring loop within a highresistivity stellar corona cannot be produced by our ideal MHD simulations, we can mimick the flare-produced EPs by releasing them within 1 R from the stellar surface. Figure 7 depicts for the quiescent SW the 2D histogram at four distinct locations (R = 0.25 R b , R = 0.5 R b , R b and R c ) of EPs injected within the lower corona (R s = 1.2R ), that represents flare-emitted EPs. The same EPs pattern as at larger injection radius R s is found, including depleted regions, as wide as ∆φ ∼ 100 • . In the coronal region the magnetic field suppresses λ by about a factor of 10, favouring EPs precipitation to the star from the open/closed lines boundary, as discussed above. The EP flux to the planet is as intense as for greater R s , although is limited over two smaller azimuthal intervals, i.e., ∆φ ∼ 40 • or ∼ 1 (or 2) days along planet b (c) orbits.

Post-CME stellar wind
In this section we present for the first time the propagation of EPs within the wind of a highly magnetized and active star 90 minutes after the eruption of a very energetic CME (see Fig.8). The CME is initiated by coupling the Alfvén Wave Solar Model (AWSoM, van der Holst et al. 2014) and the Titov & Démoulin (1999) flux-rope eruption model, jointly used, for instance, to study fast CME-driven shocks with associated solar EP events (Jin et al. 2013). The initial conditions for the stellar wind are the same as considered for the quiescent state in Sect.4.1 from (Klein et al. 2021b). In addition, we have used as initial condition the ZDI maps derived in Kochukhov & Reiners (2020), and show here only the result for the post-CME case. In particular, EPs are released on a spherical surface at radii R s = 2 and 5 R , mimicking travelling shocks, at 90 minutes past the CME onset, as the CME front has crossed the entire simulation box.
EPs are assumed to diffuse into a turbulence with the same spectral index as the quiescent wind, as justified below. The CME is likely to stir up the parameters of  the quiescent stellar wind turbulence more significantly the higher the CME kinetic energy.
In the case of heliospheric CMEs, the power spectrum of the magnetic turbulence in the CME sheath (region between the shock front and the front of the CME driving it, crossed by a spacecraft typically in several hours) has been measured in-situ at 1 au by the Wind spacecraft and analyzed by Kilpua et al. (2021); a steepening of about 5 − 10% of the inertial range power law index for an interval of 2 hours was revealed, with a considerable data spread. However, the same statistical analysis has not been carried out for post-CME front turbulence, needed for a lag of a few hours in the present analysis.
In the case of EPs released at R s = 5 R ∼ 2.5 × 10 11 cm the wind advection time to a certain radius R is ∆t(R) = (R − R s )/Ū , whereŪ is an average wind speed. Since the wind speed is highly variable in this region between ∼ 1, 000 and ∼ 5, 000 km/s, at the planet AU Mic b the time-lapse along the stellar wind ∆t(R b ) is between 25 and 125 minutes and at the box boundary R box = 120 R = 6.0×10 12 cm the ∆t(R box ) is between 3 and 16 hours; thus, the parcel of wind plasma where EPs are released 90 minutes after the CME onset is likely to arrive to the box boundary between 3 and 16 hours after the CME front. As mentioned above, the heliospheric turbulence past the CME front has not been accurately investigated, so the turbulence seen by the EPs at the injection is not constrained by solar measurements. It is reasonable to assume that the pre-CME turbulence conditions (Kolmogorov isotropy) are restored in the parcel of gas, and at the time, of EP release. An increase of the wind total magnetic field magnitude, with a wider Figure 8. Left Snapshot in a 175 R field of view of the bi-lobed plasma density isosurface (n(x, t)/nSS(x) = 10.0, where n(x, t) is the density at location x and time t and nSS(x) is the steady-state density at that location) of the propagating CME within the wind reconstructed from ZDI maps in Klein et al. (2021b) at a time 90 minutes past the CME initiation. The sphere at R b /2 is color coded by the B-field strength. The two cyan circles represent the orbits of planet b (solid) and c (dotted). The magenta lines indicate selected closed magnetic field lines. Open magnetic field lines are in green if one end is attached to the star surface and in black if none of the ends is tracked, respectively. Right Same as the Left panel with a 60 R field of view.
spread, in heliospheric CMEs has also been measured by (Kilpua et al. 2021); however, this effect is already accounted for in our 3D-MHD simulations. Figure 9 (all panels in the left column) shows for weak turbulence (σ 2 = 0.01) a significantly different pattern from the nearly homogeneous distribution of the quiescent case: Figure 5 shows depleted regions partially corresponding with the current sheet (hereafter CS) silhouette and not significantly broadening between R b and R c . As in the quiescent case of TRAPPIST-1 (Fraschetti et al. 2019), the near-homogeneity of the 2D histogram at R = R b (for σ 2 = 0.01) reflects the homogeneity of the injection of EPs on the injection sphere at R s . In the post-CME case, the EPs distribution in the southern hemisphere mirrors the region of minimal B-strength (blue in Fig. 10) only in a narrow depleted and tilted segment (110 • < φ < 210 • , 50 • < θ < 70 • ) in the midand bottom panels of Fig.9, in contrast with the quiescent case: the CME disrupts the large scale structure of the B-field by distorting and pushing the CS away from its original location (compare the CS in Fig.1 with Fig.11, top panels) and compressing the field strength by a factor up to a few tens from its quiescent value. Comparison of Fig. 10 and Fig.3 shows such a CME-driven compression by a factor > 10 throughout the astrosphere. The magnetic compression reduces throughout the astrosphere the EPs λ (by a factor ∼ 100 1/3 ∼ 4.6, see Eq. 2).
Upon comparison of the right column, mid-panel, of Fig. 9 for the post-CME case (at R = R b ) with the case of quiescent wind in Fig. 2, EPs are channelled in the latter case into two polar caps connected by an axis tilted by ∼ 10 − 20 • from the ecliptic plane. The asymmetry of the caps in the post-CME case, both largely in the southern hemisphere, results from the change in the large-scale B-field caused by the CME eruption. Fig. 11 maps the 3D spatial location of the EPs hitting (spherical) points at the R b -sphere. Comparison of the EP locations with the B-field strength in the two top panels shows that the region with relatively large and CMEcompressed magnetic field on the sphere is EPs-depleted (cfr. also with the EP-depleted regions in Fig.9) and corresponds to the launched high-density CME-lobes (bottom panel in Fig.11). On the back-side of the lobes the EPs fill the region with lower magnetic field. We conclude that the rising CME inflates closed field lines in the direction of its motion (the CME cannot break field lines in our non-resistive MHD simulations) so that closed lines extend out to larger radii and expand sideway; as seen in the previous section, those closed lines cause the back-precipitation of EPs to the star surface. On the back side (no CME lobes, no magnetic field compression) EPs follow the open lines and escape toward the planets. Figure 9 shows also that the maximum EP flux is greater in the post-CME scenario than in the quiescent wind (cfr. Fig. 2). Top row: For a stellar wind 90 minutes past the eruption of a 10 36 erg kinetic energy CME, spherical-coordinates 2D EP histogram at radius R = R b /2 (top row), R = R b (mid row), R = Rc (bottom row) of the hitting points for 1 GeV kinetic energy protons, for σ 2 = 0.01 (left column) and σ 2 = 1 (right column), injected at Rs = 5R ; here Lc = 10 −5 AU. The quiescent stellar wind is constructed from the ZDI radial field map from Klein et al. (2021b). The x (y) axis indicates the azimuthal (polar) coordinate on that sphere. The log-scale colorbar counts logarithmically the number of EPs. Bottom row: Same as top row for R = R b .

DISCUSSION
In Fig.12 the scattering mean free paths in units of the stellar radius (λ /R ) as a function of R/R are compared for the innermost planets in AU Mic and for the innermost HZ TRAPPIST-1 planet. The smaller (by a factor ∼ 7) star radius but the greater (by a factor ∼ 3) surface average magnetic field (hence the smaller λ by a factor ∼ 3 1/3 ) for TRAPPIST-1 combine so that at R = R s the mean free paths λ /R have comparable values (∼ 0.02), thereby explaining the comparable angular size of the depleted regions, i.e., vanishing EP flux, in the two systems.
The choice of a 90 minute post-CME snapshot does not require severe constraints on the flaring rate or more generally on the time lag between two transient events, i.e., very energetic CMEs or coronal flares. Consistently with the 90 minutes chosen, TESS data indicate for AU Mic a flaring rate of ∼ 1 flare every 3.8 hours (Gilbert et al. 2021) for a flux at AU Mic b of 1.8 × 10 4 erg cm −2 s −1 sr −1 sr and 4.4 × 10 5 erg cm −2 s −1 sr −1 sr, respectively; this value is up to 30% the solar flux on Earth (i.e., 1, 373 W m −2 ).
A comparison of the spatial distribution of the EPs at the distances R b /2, R b and R c for σ 2 = 0.01 in the quiescent and in the post-CME case (Figs. 9 and 5), shows that, if a CME escapes, the closed magnetic field lines are inflated and hence trap efficiently EPs (in the region θ > 50 • , φ < 100 • and φ 280 • ); moreover, the reshuffling of the large-scale magnetic field caused by the CME opens regions of the southern hemisphere to EPs, de- flecting them toward hot spots (in red at 30 • < θ < 60 • , 0 • < φ < 40 • and 50 • < θ < 70 • , φ > 280 • ) from the equatorial plane. This trapping effect by the expanding CME is shown for the strong turbulence case by Fig. 11. The angular location of the EPs hot spots (red regions in Fig. 9) depends on the spatial initialization of the CME, whose choice herein refers to Au Mic observations (Wisniewski et al. 2019): a different CME initialization might drive a more intense or weaker EP flux toward the planets. This result might seem in contrast with the expectation that the CMEs open and stretch out magnetic field lines providing additional routes for EPs to reach the planets; however, a resistive non-ideal MHD simulation would be necessary to overcome such a limitation. We have carried out multiple runs with distinct single realizations of the magnetic turbulence with no significant deviation from the conclusion above. Figure 13 (left panel), shows that the EP number at radius R relative to N inj at each distance is enhanced by the prior passage of a CME for σ 2 = 0.01, and lowered for σ 2 = 1, for the magnetic field reconstruction in Klein et al. (2021a) map. This inversion can be explained as follows. The CME inflates the closed field lines out to a large distance from the star and re-shuffles the closed field lines (see Fig.11, top row) in a pattern dependent on the location of the CME initialization region with respect to the current sheet. If σ 2 = 0.01, EPs follow the field lines with little scattering (see Fig.12), reach larger distances along the closed lines, i.e., travel a longer time, and are therefore more likely migrate to open lines and escape to the planets; thus, the EP flux is greater than the quiescent case. In the case σ 2 = 1, the number of EPs arriving to planets does not increase with respect to the case σ 2 = 0.01 as much as in the quiescent case: the migration from closed to open field lines due to the perpendicular diffusion is suppressed as the post-CME chaotic structure of the field lines dominates over the transport (see Sect. 4.1). Figure 13 (right panel) shows the relative EP number obtained with identical spatial CME initialization to the left panel and the magnetic field map in Kochukhov & Reiners (2020).
In addition, a different turbulence realization (with the same statistical properties) from the left panel was used (Fraschetti & Giacalone 2012) to show that the N (R)/N inj increase with σ 2 is not dependent on the particular details of the turbulence (simulations using the turbulence ensemble average are not shown as EP distribution is nearly homogeneous with a smaller EP flux onto the planets, thus not relevant to the present study). For the Kochukhov & Reiners (2020) map, the slope of the EP number vs σ 2 is comparable to the slope for the Klein et al. (2021a) field because the chaotic large scale structure of the post-CME field dominates over transport. This comparison shows that different conditions can lead to very different EP num- Figure 11. Top Left In the 90-minutes post-CME case based on the Klein et al. (2021b) magnetogram as in Fig. 8, EP hitting points on the R b -sphere seen from the side of the two CME expanding lobes. The strength of the unperturbed wind magnetic field B0 color-codes the R b -sphere. The two magenta circles represent the orbits of planet b (solid) and c (dotted). The purple lines indicate selected closed magnetic field lines. Open magnetic field lines are in green or black if one end is attached to the star surface or if none of the ends is tracked, respectively. Top Right Same as Top Left from the back-side. Bottom Left Snapshot in a 180 R field of view of the bi-lobed plasma density isosurface (n(x, t)/nSS(x) = 10.0) of the propagating CME within the wind reconstructed from ZDI maps in Klein et al. (2021b) at 90 minutes past the CME initiation. The spherical dots at R b /2, R b and Rc mark the EP hitting points and are color-coded by the B-field strength. The two cyan circles represent the orbits of planet b (solid) and c (dotted). Open magnetic field lines are in black. Q, R s = 2, R b /4 Q, R s = 2, R b /2 Q, R s = 2, R b Q, R s = 2, R c Q, R s = 5, R b /2 Q, R s = 5, R b Q, R s = 5, R c CME, R s = 5, R b /2 CME, R s = 5, R b CME, R s = 5, R c 0.01 0.1 1 0.1 1 σ 2 CME, R s = 2, R b /4 CME, R s = 2, R b /2 CME, R s = 2, R b CME, R s = 5, R b Figure 13. Left: Fraction of 1 GeV EPs reaching a given sphere at radius R to the total number of injected EPs, Ninj as a function of σ 2 in the stellar wind reconstructed from the ZDI maps in Klein et al. (2021b). The quiescent case (labelled as "Q") is compared with the 90-minutes post-CME case ("CME"). The red-filled circles refer to the flare case (Rs = 1.2 R , Fig. 7) corresponding, from lowest to highest N (R)/Ninj to R = 0.15R b , R b /4, R b /2, R b . Right: Same as left panel for the stellar wind reconstructed from the ZDI maps in Kochukhov & Reiners (2020), case 1 from (Alvarado-Gómez et al. 2022), in the 90-minutes post-CME case. ber at a given radius, and likely to a different planet bombardment, after the passage of the CME. As pointed out above, the perhaps unlikely CME escape from such a magnetically confining star, as well as the technical limitations in confirming CMEs from active stars, led to extrapolations of the coronal flare/CME relation from the solar system to M dwarfs. However, the tension between low mass-loss rate associated with M dwarfs (Wood et al. 2021) and the wind flux required to support very energetic CME (Drake et al. 2013) seems to indicate that flares should be more common than CMEs, hence a large fraction of EPs impinging onto planets might be released very close to the stellar surface by coronal flares rather than from CMEs. The lack of radio bursts resembling the solar type II bursts (Villadsen & Hallinan 2019) supports such a conclusion. In addition, in CME simulations shocks are generated further out in the corona, where densities are considerably smaller than in the solar region of type II burst formation and frequencies below detection threshold (Alvarado-Gómez et al. 2020). This effect is partially compensated by the relatively small flux of EPs emitted in the lower corona and reaching the planets (Fig.13, left panel, red-filled circles for σ 2 = 1).
The time-variability of the flux of stellar EPs and its effect on the planetary atmosphere (Fraschetti et al. 2019;Herbst et al. 2019), as well as on the ionization of proto-planetary disks around young stars (Fraschetti et al. 2018;Rodgers-Lee et al. 2017;Padovani et al. 2018), have been under increasing scrutiny in the past few years. The evolution of planetary atmospheres can be also affected by Galactic Cosmic Rays (GCRs), that are likely unmodulated by the stellar wind at energies > 1 − 10 GeV. Several works have considered the effect of stellar wind on the energy spectrum of GCRs impinging onto the planet, e.g., Archean Earth (Cohen et al. 2012) or exoplanets (Herbst et al. 2020;Mesquita et al. 2022). The ∼ 0.3 GeV EP propagation and modulation has been long known to be dominated by drifts in the solar system (Jokipii et al. 1977): a calculation of stellar modulation in such an energy range needs to account for drifts (Mesquita et al. 2022). However, for the solar system the flux of protons at 0.1 − 0.7 GeV during GLE events exceeds typically by 1 or 2 orders of magnitude the GCRs flux, at 1 AU. For active stars, likely producing many more energetic events than the Sun, the flux of stellar EPs at distance < 1 AU (planets b and c in AU Mic and HZ planets in TRAPPIST-1) is likely to be much higher than the local GCRs flux, even including the adiabatic losses due to the radially expanding wind (Youngblood et al. 2017). The higher energy range of GCRs (> 1 − 10 GeV) than stellar EPs might partially compensate the much lower flux of GCRs in the effect on the atmosphere evolution at the inner planets. The impact the global planet atmosphere (Segura et al. 2010;Airapetian et al. 2016) has to be investigated in further detail. Figure 14 shows that the EP flux along the planetary orbit undergoes orders of magnitude fluctuations, as was also shown in the TRAPPIST-1 case (Fraschetti et al. 2019). Likewise, we calculate here the flux of EPs impinging onto the planet by using the estimate of > 10 MeV proton flux inferred for GJ 876 by (Youngblood et al. 2017). The flux on the ecliptic plane is determined within a ring of semi-latitude aperture 5 • (despite the near-complanarity of the planet), to determine the flux of EPs in the planetary environment. The cases of quiescent wind and 90-minutes post CME are compared in Fig. 14 for 1 GeV protons at the AU-Mic b orbit. In the case of the post-CME, only ∼ 5% (∼ 1%) of the total injected EPs hit the ring enclosing the planet orbit for strong (weak) turbulence as the large part of the EPs travel toward other latitudes. The low EP flux in the plane of the planetary orbits results from the particular geometry of the fluxtube setup. This is due to initialization of the CME at low latitude, suggested by observations (Wisniewski et al. 2019), perhaps contrary to expectation: the expansion of the CME inflates closed field lines over a vast angular region that includes the equatorial plane, preventing escape of EPs toward the planets. It is conceivable that with a smaller misalignemnt between the B-field and stellar rotation axis, CME lobes travel poleward and subsequently emitted EPs might more easily be magnetically connected to the planets along open field lines.

CONCLUSION
We have carried out numerical simulations of the propagation of ∼GeV protons out to the two innermost planets in the reconstructed astrosphere of the dM1e star AU Microscopii, for the first time both in the quiescent and CME-disrupted state. Energetic particles are injected at a variety of distances from the star on spherical surfaces with an isotropic velocity distribution and diffuse in the turbulent stellar magnetic field.
The post-CME wind is likely to be the most common stellar wind configuration of very active stars encountered by propagating EPs due to the very high flaring rate; however, large stellar magnetic fields hamper CME escape and observational constraints on the rate of escaped CME are currently lacking. We determine the spherical pattern of EPs reaching the distances of planets b and c; the projection of the current sheet at the planetary distance maps the back-precipitation of EPs to the star and is enhanced by perpendicular diffusion in the strong turbulence regime.
The CME eruption re-shuffles the dipolar structure of the large scale magnetic field and dominates over the magnetic turbulence in controlling the EP flux at least 90 minutes after its eruption; as a result, the bombardment of planets by the EPs released after the CME passage can be suppressed or enhanced by the CME. A stronger turbulence leads instead in all cases to a larger EP flux at the planets. We emphasize that, even for very energetic and wide-front CMEs such as the one examined here, the EP flux along the planetary orbits depends on the region of the CME initialization, similar to the case of solar CMEs.
The effect of EPs released by CME-driven shocks localized to small spatial regions has not been considered here but merits future investigation.
We thank the referee for useful and constructive comments. FF was supported, in part, by NASA through Chandra Theory Award Number T M 0 − 21001X, T M 6 − 17001A issued by the Chandra X-ray Observatory Center, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of NASA under contract NAS8-03060; FF is also partially supported by NASA under Grants NNX16AC11G and 80NSSC18K1213 and by NSF under grant 1850774. JJD was supported by NASA contract NAS8-03060 to the Chandra X-ray Center and thanks the Director, Pat Slane, for continuing advice and support. OC is supported by NASA XRP grant 80NSSC20K0840. Resources supporting this work were provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center.