Simulating the Space Weather in the AU Mic System: Stellar Winds and Extreme Coronal Mass Ejections

The dissipation of Alfvén waves in the upper chromosphere constitutes the main mechanism behind the coronal heating and the acceleration of the stellar wind within AWSoM. These contributions appear as additional source terms in the nonideal magnetohydrodynamic (MHD) momentum and energy relations, which are solved alongside the magnetic field induction and mass conservation equations. Electron heat conduction and radiative losses are also incorporated in the simulation.

As described by van der Holst et al. (2014), AWSoM-specific proportionality constants, together with the base plasma density (n₀ = 2 × 10¹⁰ cm⁻³) and temperature (T₀ =5 × 10⁴ K), are required as boundary conditions. For the models discussed here, we assume the same set of AWSoM parameters ¹⁰ as in the M-dwarf simulations presented in Alvarado-Gómez et al. (2020b), adjusting the basic stellar properties to AU Mic (M_★ = 0.5 M_⊙, R_★ = 0.75 R_⊙, P_rot = 4.85 days; Kiraga 2012; Plavchan et al. 2020).

In order to calculate the required propagation and reflection of Alfvén waves, AWSoM also uses the distribution of the radial magnetic field (B_R) on the stellar surface as input. For simulations in the stellar regime, this information can be retrieved from spectropolarimetric observations and Zeeman–Doppler Imaging (ZDI) magnetic field reconstructions (Donati & Brown 1997, Piskunov & Kochukhov 2002).

In the case of AU Mic, two recent observational studies have provided information on the properties of its large-scale magnetic field. Kochukhov & Reiners (2020) obtained the first ZDI map of this star, revealing a nonaxisymmetric field with a (signed) average strength of ∼88 G and values up to ∼184 G in local regions (Figure 1, Case 1). However, the authors caution that due to the particular line-of-sight orientation of AU Mic (close to 90° of inclination), ZDI reconstructions employing circular polarization alone would not be sensitive to components antisymmetric with respect to the stellar equator (e.g., a dipole aligned with the stellar rotation axis). To circumvent this limitation, Kochukhov & Reiners (2020) analyzed the Zeeman signatures from AU Mic in linear polarization, which indicates the presence of an additional axisymmetric large-scale dipole field of ∼2.0 kG (Figure 1, Case 2). Employing a different set of circular polarization observations, Klein et al. (2021) reported a second ZDI map of AU Mic. This map is characterized by a 450 G large-scale dipole field inclined by 19° with respect to the rotation axis, displaying an average longitudinal magnetic field strength of 475 G (Figure 1, Case 3). The effects from the equator-on orientation of AU Mic were not discussed by Klein et al. (2021) in their ZDI reconstruction. The radial magnetic field maps, incorporated at the inner boundary of our simulations, are presented in the different panels of Figure 1.

Zoom In Zoom Out Reset image size

The grid in our simulations is spherical, with the inner boundary representing the stellar surface (∼1 R_★), and the rotation axis of the star aligned with the z Cartesian direction. The spatial resolution is higher at the surface (ΔR = 0.025 R_★, ΔΦ = 1.4° ¹¹ ), with a radial stretching factor $\propto \mathrm{ln}(R)$. Given the varying maximum field strength of the input ZDI maps (see Figure 1), a different domain size is employed for each run: 80 R_★ (Case 1), 200 R_★ (Case 2), and 120 R_★ (Case 3). This is done to guarantee a closed Alfvén surface (AS) of the stellar wind, necessary for a correct description of the escaping wind solution within each three-dimensional domain. From an initial finite-difference potential field extrapolation (Tóth et al. 2011), the simulations evolve until a steady-state solution is reached.

The SWMF implementation of the Titov & Démoulin (TD; 1999) flux rope eruption model is used to perform our CME simulations. This model is based on a prescribed magnetic tension imbalance between the ambient magnetic field configuration (provided by the AWSoM steady-state solution) and a twisted magnetic loop inserted at the inner boundary of the domain (see, e.g., Tóth et al. 2007). The initialization of the TD flux rope requires seven parameters describing its location (2), orientation, size (2), loaded mass (M^FR), and magnetic free energy available in the eruption (${E}_{{\rm{B}}}^{\mathrm{FR}}$).

Table 1 contains the model parameters employed in our TD flux rope CME simulation on AU Mic. Both tilt angle and longitude were chosen in such a way that the inserted flux rope would lie in the vicinity of a polarity inversion line on all three surface-field configurations. For the anchoring latitude, we used one of the values of the starspot locations on AU Mic derived from the photometric light-curve modeling performed by Wisniewski et al. (2019). This selection was made so that the emerging CMEs would have a higher probability of impacting the planets lying in the equatorial plane of the system. The remaining parameters were adjusted so that the inserted TD flux rope would be able to power a CME with mass and kinetic energy consistent with the best candidate event observed in this star so far (M^CME ∼ 10²⁰ g, ${E}_{{\rm{K}}}^{\mathrm{CME}}\,\sim {10}^{36}$ erg; Katsova et al. 1999).

Table 1. Initialization Parameters of the TD Flux Rope CME Simulation on AU Mic

Parameter	Value	Unit
Tilt angle a	215.0	deg
Longitude	330.0	deg
Latitude	9.6	deg
Radius (RFR)	20.0	Mm
Length (LFR)	99.5	Mm
Mass (MFR)	1019	g
Magnetic energy (${E}_{{\rm{B}}}^{\mathrm{FR}}$)	4.2 × 1037	erg

Note.

^a Measured with respect to the stellar equator in the counterclockwise direction.

Download table as:  ASCIITypeset image

After the TD flux rope insertion, the evolution of the resulting CME on each surface magnetic field configuration is followed by a time-dependent simulation covering 90 minutes (real-time), with full-domain snapshots extracted at a cadence of 1 minute. Note that while the spatial resolution in our simulations is not able to capture the internal structure of the flux rope (six resolution elements between the flux rope ends at the surface), it is sufficient to describe its global evolution and the subsequent CME propagation within the domain.

The three panels of Figure 2 contain the steady-state stellar wind solutions obtained for AU Mic, employing observed large-scale magnetic field distributions of this star as boundary conditions (Section 2). For each case we include the orbits of planets b and c (Section 1), as well as the AS of the stellar wind. This structure is defined by the spatial locations with an Alfvénic Mach number M_A = 1, where the wind velocity matches the Alfvén speed V_A (i.e., ${M}_{{\rm{A}}}\equiv U/{V}_{{\rm{A}}}=U\sqrt{4\pi \rho }/B=1$, where ρ, U, and B correspond to the local stellar wind density, velocity, and magnetic field strength, respectively). The AS establishes the boundary for the escaping stellar wind (M_A > 1) and the magnetically coupled outflows in the corona (M_A < 1) that could still fall back to the star. The AS also separates the sub- and super-Alfvénic regimes of the stellar wind, which will influence dramatically the structure of any magnetosphere/ionosphere around an orbiting planet in the system (see Cohen et al. 2014 and references therein).

Zoom In Zoom Out Reset image size

In line with other stellar wind studies (e.g., Cohen & Drake 2014; Vidotto et al. 2014; Réville et al. 2015), stronger surface magnetic fields lead to faster stellar winds and larger AS structures. The terminal radial wind speeds ¹² within the respective simulation domains were roughly 1150 km s⁻¹ (at 80 R_★, Case 1), 2200 km s⁻¹ (at 200 R_★, Case 2), and 2000 km s⁻¹ (at 150 R_★, Case 3). Similarly, we obtain average AS radii of 27.8 R_★ (Case 1), 106.4 R_★ (Case 2), and 73.1 R_★ (Case 3). In terms of the stellar wind mass-loss rates (${\dot{M}}_{\unicode{x02605}}$), our models yield values between ∼5 and 10 ${\dot{M}}_{\odot }$. ¹³

Unfortunately, there are no observational constraints on the stellar wind properties of AU Mic. Applying the observed coronal/disk properties in this system to the results of theoretical, numerical, and semiempirical studies from the literature (Augereau & Beust 2006; Strubbe & Chiang 2006; Cranmer et al. 2013; Schüppler et al. 2015) leads to mass-loss rate estimates ${\dot{M}}_{\unicode{x02605}}\sim 10\mbox{--}300\,{\dot{M}}_{\odot }$. Very recently, Kavanagh et al. (2021) simulated the stellar wind of AU Mic using the ZDI map of Klein et al. (2021) as a boundary condition (Case 3). They report solutions with low (${\dot{M}}_{\unicode{x02605}}=27\,{\dot{M}}_{\odot }$) and high (${\dot{M}}_{\unicode{x02605}}=590\,{\dot{M}}_{\odot }$) mass-loss rates. ¹⁴ It is good to note that our relatively low ${\dot{M}}_{\unicode{x02605}}$ values appear close to estimates from stellar wind astrospheric detections for other M-dwarf stars (Wood et al. 2021). Similarly, our resulting range of mass-loss rate per unit surface area on AU Mic (${\dot{M}}_{\unicode{x02605}}/{A}_{\unicode{x02605}}\,={\dot{M}}_{\unicode{x02605}}/4\pi {R}_{\unicode{x02605}}^{2}\simeq 9-18\,{\dot{M}}_{\odot }/{A}_{\odot }$) is consistent ¹⁵ with the Lyα astrospheric constraint available for EV Lac (M3.5V, P_rot = 4.38 days, F_X = 1.56 × 10⁷ erg cm⁻² s⁻¹, ${\dot{M}}_{\unicode{x02605}}/{A}_{\unicode{x02605}}\,\simeq 9.8\,{\dot{M}}_{\odot }/{A}_{\odot };$ Wood et al. 2005), which displays very similar rotation and coronal activity to AU Mic (P_rot = 4.85 days, F_X = 1.60 × 10⁷ erg cm⁻² s⁻¹).

Both transiting planets in the AU Mic system lie very close to the equatorial plane (i.e., i ≃ 90°; Martioli et al. 2021). Figure 2 shows the projection of the AS (white solid line) as well as its average size on this plane ${\left\langle R\right\rangle }_{\mathrm{AS}}^{\mathrm{eq}}$. In all three cases, we find that for the majority (if not the totality) of their orbit, both planets remain inside the sub-Alfvénic regime of the stellar wind. As described by Cohen et al. (2014), a planetary global magnetosphere within such a sub-Alfvénic stellar wind will tend to be open (similar to the Alfvén wing configuration observed in the Galilean moons; see Neubauer 1998). Under those conditions, the magnetic field of the star and the planet could reconnect directly, impeding the generation of a bow-shock structure (see Strugarek 2021). This could modify drastically the amount of stellar wind energy dissipated into the planetary ionosphere/atmosphere (Cohen et al. 2018), as well as the observable patterns of evaporative outflows from the exoplanets in the system (Harbach et al. 2021).

Note also that the 2 kG dipole aligned with the stellar rotation axis in Case 2 generates a smaller ${\left\langle R\right\rangle }_{\mathrm{AS}}^{\mathrm{eq}}$ compared to the 450 G inclined dipole from Case 3. This follows from a general result of magnetically driven stellar winds models in which, for simple large-scale magnetic field geometries (i.e., dipole, quadrupole), the largest AS lobes are found right above the regions with the strongest magnetic field, while the smallest portions of the AS will emerge in the vicinity of the polarity inversion lines of the surface field (Garraffo et al. 2015; Vidotto et al. 2015; Alvarado-Gómez et al. 2016). In this way, a stellar dipole field aligned with the z-axis (as in the Case 2 simulation) will have a polarity inversion line coincident with the equatorial plane so that the resulting stellar wind AS will be the smallest at latitude 0°. On the other hand, the ∼20° inclination of the large-scale dipole field in Case 3 is sufficient to make the large lobes of the AS cross the equatorial plane, making it much bigger at that particular latitude despite the weaker surface field. This example clearly illustrates that the magnetic field strength alone is insufficient to characterize the stellar wind environment of a planet-hosting star. It is also worth bearing in mind that the secular change of the surface magnetic field due to the growth and decay of active regions, and perhaps also to a magnetic cycle, can change the location of the current sheet and the shape of the AS.

The visualizations of Figure 2 also display the stellar wind dynamic pressure (P_dyn = ρ U²), normalized to typical conditions experienced by Earth ¹⁶ (${P}_{\mathrm{dyn}}^{\oplus }$). Within the three considered cases, we find values of P_dyn along the orbit of AU Mic b ranging between ∼1150 and 13,300 ${P}_{\mathrm{dyn}}^{\oplus }$, going down to ∼450–6800 ${P}_{\mathrm{dyn}}^{\oplus }$ for AU Mic c. Similar conditions have been predicted for the habitable-zone planet Proxima b (Garraffo et al. 2016), as well as for the exoplanet candidate Proxima d (Alvarado-Gómez et al. 2020a). Still, we find that the planetary environment around AU Mic is much harsher than for Proxima Centauri. The reason for this is that the relatively strong large-scale surface magnetism of AU Mic leads to magnetic pressure values (P_mag = B²/8π ∝ 1/R⁶) greater than the associated P_dyn at the exoplanet orbits (which is not the case for Proxima b). At the distance of AU Mic b, we find that, on average, P_mag appears larger than P_dyn by factors of 1.4 (Case 1), 11.1 (Case 2), and 5.2 (Case 3). Likewise, along the orbit of AU Mic c, Cases 2 and 3 show average P_mag values greater than P_dyn by factors of 2.3 and 2.0, respectively. Only in the simulation with the weakest surface magnetic field (Case 1) does the average P_dyn dominate over P_mag at the distance of AU Mic c (by a factor of ∼2.7). The quiescent environment obtained for the AU Mic planets then appears more comparable to the extreme conditions expected for the planets of the TRAPPIST-1 system (Garraffo et al. 2017).

The steady-state stellar wind solutions discussed in the previous section serve as the initial state for our time-dependent CME simulations. The latter are driven by a single erupting TD flux rope, initialized with the parameters provided in Section 2.2. As expected, while the properties of the inserted flux rope are identical, the varying magnetic/stellar wind conditions among the three cases lead to quite different CMEs.

We follow the procedure introduced by Alvarado-Gómez et al. (2019a, 2020b) to identify and characterize our simulated CME events. For each solution, we compute the 3D density contrast n(t)/n^SS, with n(t) and n^SS as the instantaneous and pre-CME (steady-state) particle density distributions, respectively. The escaping CME front is traced by the collection of points on the time-evolving isosurface n(t)/n^SS = 3.0 with radial speeds greater than the local escape velocity. ¹⁷ As different sections of the CME face dissimilar coronal and stellar wind conditions, the front expansion is highly inhomogeneous. To estimate the instantaneous maximum radial speed of the CME as a whole (${U}_{{\rm{R}},{\rm{T}}}^{\mathrm{CME}}$), we take the mean position among the 1% outermost grid points forming the escaping front and determine the relative variation of this average value between consecutive frames (1 minute cadence). The escaping CME front also serves as the limit of integration for a volume integral of the plasma density, allowing the determination of the total mass (${M}_{{\rm{T}}}^{\mathrm{CME}}$) and kinetic energy (${E}_{{\rm{K}},{\rm{T}}}^{\mathrm{CME}}$) of the eruption at each time step.

This analysis revealed that during their evolution, the emerging CME structures enclosed two separate regions where the density contrast increased by more than one order of magnitude (i.e., n(t)/n^SS > 10). Figure 3 shows the CME fragmentation in Case 3, but it was also observed in the other runs as well. As can be seen in the visualization, despite having a launching latitude of only ∼10° (see Table 1), one of the CME fragments follows a high-latitude trajectory (close to the north pole of the star) whereas the second one propagates closer to the equator. Similar behavior was obtained on the strongest CME events analyzed in Alvarado-Gómez et al. (2018), where the large-scale stellar field only provided weak magnetic confinement. The Case 2 simulation, in which the large-scale dipole field is the strongest of all cases (see Figure 1), shows that the two CME fragments are slowly guided in opposite poleward directions. Interestingly, only in the simulation for Case 1, where the surface field lacks a large-scale dipole field (see Section 2.1), is the propagation of both CME fragments centered along the equatorial plane where the planets reside. These results are better illustrated in Figures 4 and 5, where the different snapshots correspond to the arrival time of the event at the orbit of AU Mic b (Figure 4) and AU Mic c (Figure 5).

Zoom In Zoom Out Reset image size

To quantify parameters such as the speed (${U}_{{\rm{R}},1-2}^{\mathrm{CME}}$), mass (${M}_{1-2}^{\mathrm{CME}}$), and kinetic energy (${E}_{{\rm{K}},1-2}^{\mathrm{CME}}$) of these CME fragments, we followed the methodology applied to the event as a whole, taking in this case an n(t)/n^SS = 10 isosurface as the boundary of reference. Table 2 contains time averages of the CME properties for the entire event, as well as values associated with the individual CME fragments on each run. Notice that in all cases the global values of $\left\langle {M}_{{\rm{T}}}^{\mathrm{CME}}\right\rangle$ are close to two orders of magnitude lower than the CME candidate event parameters reported on AU Mic by Katsova et al. (1999) and roughly five times lower than the amount of mass originally inserted in the TD flux rope (see Table 1). This indicates that an important fraction of the flux rope eruption does not escape and remains magnetically confined in the low corona (see Alvarado-Gómez et al. 2019a). Despite this, the escaping eruptions display $\left\langle {E}_{{\rm{K}},{\rm{T}}}^{\mathrm{CME}}\right\rangle$ values consistent with the observed CME candidate event in this star.

Table 2. Properties Obtained in the Simulated CME Events on AU Mic

Case	# of CME	$\left\langle {M}_{{\rm{T}}}^{\mathrm{CME}}\right\rangle$	$\left\langle {U}_{{\rm{R}},{\rm{T}}}^{\mathrm{CME}}\right\rangle$	$\left\langle {E}_{{\rm{K}},{\rm{T}}}^{\mathrm{CME}}\right\rangle$	$\left\langle {M}_{1}^{\mathrm{CME}}\right\rangle$	$\left\langle {U}_{{\rm{R}},1}^{\mathrm{CME}}\right\rangle$	$\left\langle {E}_{{\rm{K}},1}^{\mathrm{CME}}\right\rangle$	$\left\langle {M}_{2}^{\mathrm{CME}}\right\rangle$	$\left\langle {U}_{{\rm{R}},2}^{\mathrm{CME}}\right\rangle$	$\left\langle {E}_{{\rm{K}},2}^{\mathrm{CME}}\right\rangle$
	Fragments a	(1018 g)	(km s−1)	(1035 erg)	(1018 g)	(km s−1)	(1035 erg)	(1018 g)	(km s−1)	(1035 erg)
1	2	2.15	10146	9.45	1.46	9924	6.24	0.56	4112	0.41
2	2	1.94	5102	3.51	0.95	4078	0.73	0.27	4968	0.65
3	2	2.51	9653	12.8	2.19	8278	8.53	0.09	9631	0.39

Notes. Parameters for the events as a whole (T) as well as for the individual CME fragments (1, 2) are included (see text for details). Values listed in columns 3–11 correspond to temporal averages of instantaneous maxima during the entire event evolution.

^a Identified by spatially disconnected escaping perturbations with n(t)/n^SS ≥ 10.

Download table as:  ASCIITypeset image

In further agreement with the CME magnetic suppression results of Alvarado-Gómez et al. (2018), the mean total mass of the eruptions is comparable between Cases 1 and 2, while the presence of a strong surface large-scale magnetic field in the latter case leads to an important reduction in the CME radial speed (∼50%) and kinetic energy (∼65%). Interestingly, the Case 3 simulation (driven by a different geometry, which includes an inclined 450 G large-scale dipole field; see Section 2.1) only shows a minor reduction in $\left\langle {U}_{{\rm{R}},{\rm{T}}}^{\mathrm{CME}}\right\rangle$ compared to Case 1 (by ∼5%). While further investigation is required, this suggests that the large-scale dipole field inclination will strongly influence the dynamical properties of an emerging CME.

In terms of the CME fragmentation, our calculations indicate that at least one fragment propagates with a similar speed to the global CME, and that combined they carry a significant fraction of the total escaping CME mass (between 60% and 90%). However, the mass distribution among them is not uniform, with one fragment being considerably more massive than the other (see Table 2). For Cases 1 and 3, this mass discrepancy translates to a difference of more than one order of magnitude in the average kinetic energy of each CME fragment. The situation is more balanced in Case 2, as the lighter CME fragment compensates by traveling at a higher speed.

As can be seen from Figures 4 and 5, the impact of the CME on the exoplanets in the system will be mediated by the fragmentation of the eruption. Depending on the large-scale magnetic field geometry, both exoplanets can be strongly affected by one or two CME fragments (Cases 1 and 3) or can receive a relatively weak perturbation compared with the global energetics of the entire event (Case 2). In the following section, we evaluate how the space weather environment at the exoplanet orbits is affected by the arrival and passing of these CMEs.

Figure 6 shows the conditions along the orbit of AU Mic b during the evolution of the exoplanet-effective CMEs emerging from each analyzed case. We include the behavior of local dynamic pressure (P_dyn, left) as well as its ratio to the magnetic pressure (P_dyn/P_mag, right). An analogous diagram along the orbit of AU Mic c is presented in Figure 7. The steady-state solutions (Figure 2, Section 3.1) provide the conditions at t = 0 minutes in these time-phase representations.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Apart from differences in the CME speeds, the visualizations of Figures 6 and 7 reveal various noteworthy features. The CME arrival in all cases is characterized by an extremely sharp increase in P_dyn with respect to the pre-CME state —roughly of two to four orders of magnitude for planets b and c, respectively. While the CME conditions are harsher in absolute terms for the planet closer to the star, the milder steady-state stellar wind environment surrounding the outermost planet leads to a larger P_dyn jump at the CME arrival.

In the denser regions of the expanding CME front (forming the cores of the CME fragments reaching the planets; ¹⁸ see Figures 4 and 5), P_dyn rises to values up to 10⁸ times the dynamic pressure experienced by Earth on nominal solar wind conditions. Even the arrival of a Carrington-type CME event at Earth's orbit would induce a P_dyn transient five orders of magnitude smaller. ¹⁹ While such P_dyn values also exceed by far the total pressure (P_dyn + P_mag) predicted for planets orbiting low-mass stars at much closer distances compared to the AU Mic planets (e.g., Vidotto et al. 2015; Garraffo et al. 2017), these extreme CME-driven conditions persist only for a relatively short time (tens of minutes for the events considered here). Therefore, their global impact on a given atmospheric and magnetospheric structure of an exoplanet will mainly depend on the effective escaping CME rate of the star (mediated by the stellar large-scale magnetic field).

The bulk of the eruption strongly affects a large fraction of the orbit of both exoplanets—between 60% and 100% depending on the case. In those regions the dynamic pressure ranges from ${10}^{6}\,{P}_{\mathrm{dyn}}^{\oplus }$ down to a few ${10}^{4}\,{P}_{\mathrm{dyn}}^{\oplus }$. As illustrated in the right panels of Figures 6 and 7, these values are sufficient to overcome the general dominance of P_mag in this system. Furthermore, the CME-driven enhancement in the local plasma density and velocity will shift the conditions along the affected orbital phases from the sub-Alfvénic to the super-Alfvénic regime (see Section 3.1). However, this situation is expected to be temporary, as the sub-Alfvénic conditions along the orbit return once the stellar wind has relaxed into a new steady state (not captured in our simulations).

One interesting and counterintuitive consequence of this CME-driven transition is that the eruption passage would shut down any possible planet-induced stellar radio emission via electron cyclotron maser instability (ECMI; see Saur et al. 2013; Kavanagh et al. 2021). This type of emission is expected to occur in star–planet interactions (SPIs) within the sub-Alfvénic regime of the stellar wind, where Alfvén waves traveling toward the star would carry enough energy to induce radio emission via ECMI in the stellar corona (Turnpenney et al. 2018). As the sub-Alfvénic regime is temporarily lifted by the CME passage, the dimming of this radio signal (assuming that a pre-CME baseline could be detected and isolated from other possible stellar sources) may potentially serve to find candidate CME events on the host star—analogous to coronal dimming events observed in the Sun and other stars (see Mason et al. 2016; Veronig et al. 2021). Furthermore, if such a CME event would be associated with a flare, the radio-dimming time lag should be able to yield an estimate of the CME speed (provided that the orbit of the exoplanet is known). While their characterization is beyond the scope of this paper, these SPI-ECMI radio-dimming features could constitute an alternative pathway to study the eruptive behavior of exoplanet hosts. This is not only relevant as the information on stellar CMEs in general is very scarce (see Argiroffi et al. 2019; Moschou et al. 2019) but also because the radio properties of cool stars at low frequencies are starting to be unveiled using ground-based observations (Vedantham et al. 2020; Callingham et al. 2021).

Finally, as noted in Section 1, observations have shown considerable flare activity in AU Mic, with estimated rates in the range of 5.54–6.35 flares/day for events with bolometric energies E^FL < 10³² erg, down to one event every 10 days for large flares with E^FL ∼ 10³⁴ erg (Martioli et al. 2021; Gilbert et al. 2022). As such, on a given rotation AU Mic displays on the order of 25–30 flaring events with E^FL < 10³² erg and about 5 flares that reach 10³³ erg in bolometric energy. Even though CMEs associated with this flare activity will suffer strong magnetic suppression, a quiescent state in AU Mic (described in Section 2.1 and also by other authors in previous studies) might not represent the most common conditions around this flare star and its exoplanets. Such conditions are expected to be characterized by a rapidly changing corona and stellar wind environment, with considerable variability on a subrotational period timescale.

While the flare-CME candidate inspiring our simulations will not be as frequent as those milder events, ²⁰ its evolution reveals how much the space weather of the exoplanets and the entire inner astrosphere can be affected by these energetic transients. The resulting CME fragmentation leads to a wide range of latitudes, including the polar regions, being disturbed by the expanding structures (Figures 4 and 5). Changes in the geometry of the AS as well as the complete disruption of the current sheet (see Figures 6 and 7) are just two examples of the CME influence on the circumstellar environment. As the typical timescale involved in all of these CME-related affectations is comparable to the exoplanets' transit times (3.56 hr for AU Mic b and 4.42 hr AU Mic c, Gilbert et al. 2022), any observational and/or modeling efforts performed on this object should include time-dependent effects in their analysis. We will pursue this research direction in two complementary numerical studies on AU Mic, evaluating the influence of this extreme CME on atmospheric outflow patterns from the innermost planet (O. Cohen et al. 2022, in preparation), as well as investigating the pre- and post-CME energetic particle environments in this system (F. Fraschetti et al. 2022, in preparation).