Modeling a Carrington-scale Stellar Superflare and Coronal Mass Ejection from $\kappa^{1}Cet$

TL;DR

This paper models a massive stellar superflare and associated coronal mass ejection from $$ Cet, revealing how magnetic energy buildup can produce events comparable to historic solar flares, with significant implications for exoplanet habitability.

Contribution

It introduces a 3D magnetohydrodynamic model of a Carrington-scale superflare on a young solar-like star, linking magnetic shear to extreme energy release and space-weather effects.

Findings

Superflare energy comparable to the 1859 Carrington event.

Synthetic emission and shock signatures suggest extreme space-weather impacts.

The 1986 Robinson-Bopp superflare may have been as extreme as the Carrington event for $$ Cet.

Abstract

Observations from the Kepler mission have revealed frequent superflares on young and active solar-like stars. Superflares result from the large-scale restructuring of stellar magnetic fields, and are associated with the eruption of coronal material (a coronal mass ejection, or CME) and energy release that can be orders of magnitude greater than those observed in the largest solar flares. These catastrophic events, if frequent, can significantly impact the potential habitability of terrestrial exoplanets through atmospheric erosion or intense radiation exposure at the surface. We present results from numerical modeling designed to understand how an eruptive superflare from a young solar-type star, $\kappa^{1}Cet$, could occur and would impact its astrospheric environment. Our data-inspired, three-dimensional magnetohydrodynamic modeling shows that global-scale shear concentrated near the radial-field polarity inversion line can energize the closed-field stellar corona sufficiently to power a global, eruptive superflare that releases approximately the same energy as the extreme 1859 Carrington event from the Sun. We examine proxy measures of synthetic emission during the flare and estimate the observational signatures of our CME-driven shock, both of which could have extreme space-weather impacts on the habitability of any Earth-like exoplanets. We also speculate that the observed 1986 Robinson-Bopp superflare from $\kappa^{1}Cet$ was perhaps as extreme for that star as the Carrington flare was for the Sun.

Equations

$$\frac{\partial \rho }{\partial t}+\nabla \cdot \left(\rho \mathbf{V}\right)=0,$$

$$\frac{\partial }{\partial t}\left(\rho \mathbf{V}\right)+\nabla \cdot \left(\rho \mathbf{V}\mathbf{V}\right)+\nabla P=\frac{1}{4\pi }\left(\nabla \times \mathbf{B}\right)\times \mathbf{B}-\rho \mathbf{g},$$

$$\frac{\partial \mathbf{B}}{\partial t}=\nabla \times \left(\mathbf{V}\times \mathbf{B}\right),$$

$$\frac{V_{\mathrm{sw}}^2}{c_0^2}-\ln \left(\frac{V_{\mathrm{sw}}^2}{c_0^2}\right)=-3+4\ln \left(\frac{r}{r_c}\right)+4\frac{r_c}{r},$$

$${\mathbf{V}}_{\mathrm{shear}}(\theta ,\varphi ,t)=V_0{f}_T{f}_B\left[\hat{\mathbf{r}}\times {\nabla }_{\mathbf{t}}{B}_r(\theta ,\varphi )\right],$$

$${\nabla }_{\mathbf{t}}=\frac{1}{R_{\odot }}\frac{\partial }{\partial \theta }\hat{\mathbf{\theta }}+\frac{1}{R_{\odot }\sin \theta }\frac{\partial }{\partial \varphi }\hat{\mathbf{\varphi }},$$

$$f_{T}(t)=\left\{\begin{array}[]{ll}\frac{1}{2}-\frac{1}{2}\cos{\left[\pi\frac{(t-100)}{5}\right]}&{\rm for}\;100\leq t<105,\\ 1&{\rm for}\;105\leq t<140,\\ \frac{1}{2}-\frac{1}{2}\cos{\left[\pi\frac{(145-t)}{5}\right]}&{\rm for}\;140\leq t<145,\end{array}\right.$$

$$f_{B}(B_{r})=\left\{\begin{array}[]{ll}\sin{\left[2\pi\frac{(B_{r}+8)}{7}\right]}&{\rm for}\;-8.0\leq B_{r}\leq-1.0,\\ \sin{\left[2\pi\frac{(B_{r}-1)}{7}\right]}&{\rm for}\;+1.0\leq B_{r}\leq+8.0.\end{array}\right.$$

$${\Phi }_{\mathrm{rxn}}=\sum _{j,k}|B_r(R_{\odot },{\theta }_j,{\varphi }_k)|d{A}_{jk},$$

$${\epsilon }_J(\theta ,\varphi )=\sum _{i}d{\epsilon }_J({\mathbf{r}}_{ijk})=\sum _i⟨J^2{⟩}_{ijk}$$

$$⟨J^2{⟩}_{ijk}=\frac{1}{L}{\int }_0^{L}d{\ell }^{\prime }|J(\mathbf{r}({\ell }^{\prime })){|}^2$$

$${\epsilon }_N(\theta ,\varphi )=\sum _{i}d{\epsilon }_N({\mathbf{r}}_{ijk})=\sum _i{n}^2(r_i,{\theta }_j,{\varphi }_k).$$

$$I_{\{J,N\}}=\sum _{j,k}{\epsilon }_{\{J,N\}}({\theta }_j,{\varphi }_k)d{A}_{j,k},$$

$$f_{pe}=8.98\times 10^3{n}_e^{1/2}\left[\mathrm{Hz}\right].$$

Full text

magnetohydrodynamics (MHD) – stars: magnetic field – stars: solar-type – Sun: flares – Sun: coronal mass ejections (CMEs) – solar-terrestrial relations

1 Introduction

There is growing appreciation that planetary atmospheric chemistry, and even the retention of an atmosphere in many cases, depends critically upon the high-energy radiation and particle environments around the host star (Segura et al., 2005; Cohen et al., 2014; Domagal-Goldman et al., 2014; Rugheimer et al., 2015). This has led to a number of increasingly sophisticated modeling efforts to characterize the space environment of exoplanets and stellar-wind/exoplanet interactions (Cohen et al., 2011; Vidotto et al., 2011, 2012; Garraffo et al., 2016, 2017; Garcia-Sage et al., 2017). The cumulative effects of both steady-state and extreme space weather from active stars, including intense X-ray and extreme ultraviolet (EUV) radiation, large fluxes of highly energetic particles, and frequent exoplanet interactions with stellar coronal mass ejections (CMEs), will have a significant impact on the exoplanets’ atmospheric evolution, and ultimately on their habitability (Lammer et al., 2007, 2009; Drake et al., 2013; Osten & Wolk, 2015; Kay et al., 2016; Airapetian et al., 2019b). Additionally, stellar activity from young solar-type stars, such as $\kappa^{1}Cet$, is now being considered in investigations of the evolution of our own solar system, where enhanced levels of extreme space weather may have profoundly affected the chemistry and climate of the early Earth (Airapetian et al., 2016b, 2017; Dong et al., 2017; Garcia-Sage et al., 2017).

Solar flares—the explosive release of energy in the solar atmosphere across a wide range of electromagnetic wavelengths—occur due to the rapid release of free energy stored in the strong sheared and/or twisted magnetic fields typically associated with sunspots and active regions (Forbes, 2000; Fletcher et al., 2011; Kazachenko et al., 2012). Solar flares and CMEs are widely, although not universally, accepted as being driven by magnetic reconnection (Klimchuk, 2001; Lynch et al., 2008; Karpen et al., 2012). The long-standing CSHKP model (Carmichael, 1964; Sturrock, 1966; Hirayama, 1974; Kopp & Pneuman, 1976) for eruptive solar flares explains many of their generic observed properties (e.g. Janvier et al., 2015; Lynch et al., 2016; Welsch, 2018, and references therein).

The total energy released during solar flares typically ranges over $E\sim 10^{29-32}$ erg (Emslie et al., 2012), whereas stellar flares can extend to a much higher range of energies $E\sim 10^{32-36}$ erg (Shibata et al., 2013; Maehara et al., 2015; Notsu et al., 2019). Kepler observations of superflaring solar-type stars indicate that large starspots typically are associated with their flares, the frequency and maximum energy of the flares depend critically upon the age of the star, and younger stars exhibit greater maximum flare energies and higher flare frequencies (Maehara et al., 2012; Notsu et al., 2013, 2019; Shibayama et al., 2013).

Schaefer et al. (2000) presented some of the earliest superflare observations associated with solar-type stars, including $\kappa^{1}Cet$, a G5 young solar analog aged $\sim 0.7$ Gyr. They estimated the $\kappa^{1}Cet$ superflare energy as $E\sim 2\times 10^{34}$ erg, based on observations by Robinson & Bopp (1987) of He I emission. $\kappa^{1}Cet$ is reported as having magnetic cycles and magnetic field strengths in the kiloGauss range (Saar & Baliunas, 1992) and showing evidence of starspots that rotate differentially across the stellar disk (Rucinski et al., 2004; Walker et al., 2007). The star’s quasi-steady wind and radiation environments are objects of ongoing study for their effects on planetary habitability (e.g. Ribas et al., 2010; do Nascimento et al., 2016).

In this paper, we present results from a three-dimensional numerical magnetohydrodynamics (MHD) simulation of a global stellar superflare and its CME, based on observations by Rosén et al. (2016) of the magnetic field of $\kappa^{1}Cet$. Although stellar active-region fields associated with large starspots are thought to contain the magnetic free energy necessary to power superflares and their associated eruptions, stellar magnetogram observations tend to leave such strong-field regions unresolved. Our approach is to model self-consistently the gradual accumulation of free magnetic energy via the introduction of large-scale stresses to the observed global stellar field. In this way, we maximize both the amount of free energy available for release during the eruptive flare and the spatial scale of the stellar CME, without making any assumptions about the existence or strength of unresolved starspots on $\kappa^{1}Cet$. Therefore, our global-scale eruption represents the most extreme stellar space-weather event possible within the observed constraints on the surface magnetic flux distribution, and sets a baseline for comparison with observations of superflares on $\kappa^{1}Cet$ and similar stars.

Our paper is organized as follows. In section 2 we present the MHD numerical model. Section 3 describes the pre-eruption phase of the simulation, including the magnetic field configuration ($\S$3.1), the stellar-wind outflow ($\S$3.2), and the self-consistent, slowly driven energy-accumulation phase ($\S$3.3). In section 4 we present the eruption results and their analysis: the global-eruption overview and energy evolution ($\S$4.1), the stellar-flare reconnected flux ($\S$4.2), the synthetic X-ray and EUV emission proxies ($\S$4.3), and the properties of the CME-driven shock ($\S$4.4). Section 5 concludes with a discussion of our results in the context of future modeling and observations of stellar space weather.

2 Adaptively Refined MHD Solver (ARMS)

The ARMS code, developed by (DeVore & Antiochos, 2008) and collaborators, calculates solutions to the 3D nonlinear, time-dependent MHD equations that describe the evolution and transport of density, momentum, energy, and magnetic flux throughout the system. The numerical scheme used is a finite-volume, multidimensional flux-corrected transport algorithm (DeVore, 1991). ARMS is fully integrated with the adaptive-mesh toolkit PARAMESH (MacNeice et al., 2000) to handle dynamic, solution-adaptive grid refinement and enable efficient multiprocessor parallelization.

For our simulation, we use ARMS to solve the following ideal MHD equations in spherical coordinates

[TABLE]

[TABLE]

[TABLE]

where all the variables retain their usual meaning, solar gravity is $\boldsymbol{g}=g_{\odot}(r/R_{\odot})^{-2}\boldsymbol{\hat{r}}$ with $g_{\odot}=2.75\times 10^{4}$ cm s*-2*, and we use the ideal gas law $P=2(\rho/m_{p})k_{B}T$. Given the isothermal approximation used in the construction of our background solar wind, we do not solve an internal energy or temperature equation. The plasma temperature remains uniform throughout the domain for the duration of the simulation.

Additionally, while there is no explicit magnetic resistivity in the equations of ideal MHD, necessary, and stabilizing numerical diffusion terms introduce an effective resistivity on very small spatial scales, i.e., the size of the grid. In this way, magnetic reconnection can occur when current sheet features and the associated gradients of field reversals have been distorted and compressed to the local grid scale.

The spherical computational domain uses logarithmic grid spacing in $r$ and uniform grid spacing in $\theta,\phi$. The domain extends from $r\in\left[1R_{\odot},30R_{\odot}\right]$, $\theta\in\left[11.25^{\circ},168.75^{\circ}\right]$ ($\pm 78.75^{\circ}$ in latitude), and $\phi\in\left[-180^{\circ},+180^{\circ}\right]$ (longitude). The initial grid consists of $6\times 6\times 12$ blocks with 83 grid cells per block. There are two additional levels of static grid refinement. The highest refinement region (level 3) is $r\in\left[1R_{\odot},5.485R_{\odot}\right]$ for all $\theta,\phi$, and the level 2 refinement extends from $r\in\left[5.485R_{\odot},9.650R_{\odot}\right]$. The maximum resolution is therefore $192\times 192\times 384$ with the level 3 grid cells having an angular width of $0.820^{\circ}\times 0.938^{\circ}$ in $\theta,\phi$ and a radial extent of $\Delta r=0.01787R_{\odot}$ at the lower boundary. After the solar wind relaxation phase we turn on the adaptive-mesh refinement criteria. The maximum refinement remains at level 3 but as the eruption progresses and regions of high electric current density evolve, the computational grid refines and de-refines to track the evolution of the strong currents. The refinement criteria are described in Karpen et al. (2012).

3 Pre-Eruption Stellar Corona

3.1 Magnetic Field Configuration

For decades, the fossil magnetic fields of massive, early-type stars were analyzed assuming simple dipole or dipole-plus-quadrupole magnetic field geometries. For late-type active stars, the development of Zeeman-Doppler imaging (ZDI; Donati & Brown, 1997; Piskunov & Kochukhov, 2002) has made it possible to invert time series of high-resolution spectroscopic and spectropolarimetric Stokes observations (circularly and linearly polarized stellar light) into surface maps of parameters such as brightness, temperature, elemental abundance, and complex magnetic field geometry (Lüftinger et al., 2010a, b). ZDI has further matured and can now be used to determine the strength, distribution, polarity, and polarity reversals of surface magnetic fields (Lüftinger et al., 2015; Kochukhov, 2016). The availability of stellar magnetic field maps has significantly advanced our capacity for sophisticated numerical modeling of stellar coronae, winds, and star-planet interactions (Cohen et al., 2011; Vidotto et al., 2011; do Nascimento et al., 2016; Garraffo et al., 2016).

In our numerical simulation, we use the $\kappa^{1}Cet$ surface field distribution from Rosén et al. (2016) during the epoch 2012.9 (mid-to-late August) derived from the ZDI analysis of data from PolarBase (Petit et al., 2014). Figure 1 shows this $B_{r}(\theta,\phi)$ distribution at the lower radial boundary. Positive radial field is denoted as red, negative as blue, and the black line indicates the $B_{r}=0$ polarity inversion line (PIL). We generate the potential field source surface (PFSS; Wang & Sheeley, 1992) solution which is used to initialize the magnetic field throughout the computational domain.

3.2 Stellar Wind Outflow

The stellar properties of $\kappa^{1}Cet$ are discussed by do Nascimento et al. (2016) in a study where they present results from a data-driven MHD simulation utilizing a polytropic ($\gamma=1.1$) stellar-wind model (Vidotto et al., 2012) to show a mass-loss rate about 50 times greater that that of the Sun. While more advanced stellar-wind treatments are under development (e.g. Airapetian et al., 2019a), our focus is on the storage and release of magnetic energy in the stellar corona during the eruptive superflare. Since our background stellar wind needs only to create the distinct open- and closed-flux systems characteristic of solar and stellar coronae, here we use the isothermal Parker (1958) wind solution for a uniform temperature, $T_{0}=2\times 10^{6}$ K.

The solar wind is initialized in ARMS by first solving the one-dimensional Parker (1958) equation for a spherically symmetric isothermal corona,

[TABLE]

where the base number density, pressure, and temperature are $n_{0}=\rho_{0}/m_{p}=9.05\times 10^{8}$ cm*-3*, $P_{0}=0.5$ dyn cm*-2*, and $T_{0}=2.0\times 10^{6}$ K, respectively. Here, $c_{0}=(2k_{B}T_{0}/m_{p})^{1/2}=181.7$ km s*-1* is the thermal velocity at $T_{0}$ and the location of the critical point is $r_{c}=GM_{\odot}/2c_{0}^{2}=2.87R_{\odot}$. These parameters yield a solar wind speed at the outer boundary of $V_{\rm sw}(30R_{\odot})\approx 550$ km s*-1*.

At time $t=0$ s we impose this Parker $V_{\rm sw}(r)$ profile and use it to set the initial mass density profile $\rho(r)$ from the steady mass-flux condition ($\rho V_{\rm sw}r^{2}=$ constant) throughout the computational domain. We then let the system relax to time $t_{\rm rel}=3.6\times 10^{5}$ s (100 hr). The initial discontinuities in the PFSS magnetic field solution at the source surface ($r=2.5R_{\odot}$) propagate outwards and eventually through the outer boundary. During the relaxation phase, the highest layers of the closed streamer belt flux are carried outward by the stellar wind flow, setting up the condition for the transverse pressure from the open fields to push in behind the expanding streamer belt structure forming the elongated current sheet. Eventually, the numerical diffusion allows magnetic reconnection between the elongated streamer belt field lines and gives the system the opportunity to adjust the amount of open flux relative to the new pressure balance associated with the background stellar wind outflow. The inner-boundary mass source allows material to accumulate in the closed-field regions and sets up steady radial outflow along open field lines. Figure 2 shows a snapshot of the radial velocity (top panel) and the logarithmic number density (bottom panel) at the end of the relaxation phase ($t_{\rm rel}=100$ hr). The Alfvén surface is shown as a white contour and representative magnetic field lines are shown in each plot.

Figure 3 plots the global energy evolution during the relaxation phase: internal $E_{\rm int}$ (dash-dotted line), kinetic $E_{\rm K}$ (dashed line), gravitational $E_{\rm grv}$ (dotted line), magnetic $E_{\rm M}$ (solid thin line), and total $E_{\rm total}$ (solid thick line). The global energy curves reflect the initial PFSS magnetic configuration and the isothermal Parker wind outflow equilibrating, as described above. By the end of the solar wind relaxation phase, the global energies are $E_{\rm int}=3.24\times 10^{33}$ erg, $E_{K}=4.07\times 10^{33}$ erg, $E_{\rm grv}=8.26\times 10^{33}$ erg, $E_{M}=1.32\times 10^{34}$ erg, and $E_{\rm total}=2.88\times 10^{34}$ erg. We can estimate a globally-averaged plasma beta as $\langle\beta\rangle=E_{\rm int}/E_{M}=0.245$. We note, however, that the plasma $\beta$ is on the order of $10^{-3}$ in the strong-field regions of our stellar corona. Similarly, a globally averaged Alfvén speed can be calculated as $\langle{V_{A}}\rangle=\sqrt{2E_{M}/M_{\rm tot}}$ where $M_{\rm tot}$ is the total mass of the system. Once the mass density profiles along open-field flux tubes reach quasi-equilibrium at the end of the relaxation phase, $M_{\rm tot}(t_{\rm rel})=6.548\times 10^{18}$ g, yielding $\langle V_{A}\rangle=635$ km s*-1*. While the stellar wind outflow is both supersonic and super-Alfvénic for $r\gtrsim 5-10R_{\odot}$ (Figure 2), the globally averaged Alfvén speed can be compared to the CME and CME-driven shock velocities showing that our entire global CME eruption and its evolution proceeds at or exceeds $\langle V_{A}\rangle$.

3.3 Accumulation of Magnetic Free Energy

Idealized surface flows are imposed on the lower radial boundary of our simulation to accumulate magnetic free energy in our stellar corona system. The boundary flows are constructed to follow the contours of $B_{r}(\theta,\phi)$ exactly so that the stellar radial flux distribution remains constant throughout the simulation as in DeVore & Antiochos (2008). This ensures that the potential magnetic energy of our system remains constant during the simulation and therefore any increase in the magnetic energy represents free magnetic energy that will be available for the stellar eruption. While these flows are obviously simpler than the complex photospheric motions observed on our Sun (Li et al., 2004), the cumulative effect of both the largest-scale photospheric motions of differential rotation and meridional flows (van Ballegooijen et al., 1998; Sheeley, 2005; Yeates, 2014) acting on emerged active-region fields and the smallest-scale granulation, diffusion, and chaotic rotational motions of helicity condensation (Antiochos, 2013; Knizhnik et al., 2017) is to form elongated, global-scale PILs with highly concentrated, non-potential, sheared field structures.

The boundary shearing flow profiles are as given as

[TABLE]

where $\boldsymbol{\nabla_{t}}$ is the tangential gradient operator,

[TABLE]

acting on $B_{r}$ at the $r=R_{\odot}$ boundary. The temporal dependence is given by

[TABLE]

with $t$ in units of simulation hours and represents a smooth ramp-up period followed by a uniform driving period and a smooth ramp-down back to zero for $t\geq 145$ hr. The function $f_{B}(B_{r})$ defines the spatial extent over the surface by smoothly enforcing the range of radial field magnitudes over which to calculate the flow profiles to $|B_{r}(R_{\odot},\theta,\phi)|\in[1,8]$ G as

[TABLE]

The coefficient $V_{0}=\pm 4\times 10^{15}$ cm2 s*-1* G*-1* yields a maximum magnitude $\boldsymbol{V}_{\rm shear}\approx 20$ km s*-1*. This maximum flow speed ensures that the evolution during the energization phase is quasi-static, i.e. much less than both the Alfvén speed in the vicinity of the global PIL, $V_{\rm shear}/V_{A}\lesssim 5$%, and the sound speed, $V_{\rm shear}/c_{0}\approx 10$%.

Figure 4(a) plots the magnitude of the surface velocity components, $(V_{\theta}^{2}+V_{\phi}^{2})^{1/2}$, on the $r=R_{\odot}$ lower boundary at $t=110$ hr during the uniform shearing phase. The $\boldsymbol{V}_{\rm shear}$ distribution traces the global PIL underneath the stellar streamer belt. The small velocity magnitudes (dark blue) in the polar regions are the non-radial components of the steady-state stellar-wind outflow. Figure 4(b) plots representative magnetic field lines at $t=135$ hr, late in the shearing phase. The large-scale, sheared-arcade field structure above the global PIL is a common feature of extended filament channels on the Sun (Mackay et al., 2010; Pevtsov et al., 2012). Our field lines develop a weak twist from the structure of the boundary flows and form the characteristic dips found in many prominence field models and observations (DeVore & Antiochos, 2000; Parenti, 2014). The distribution of mass density is shown in the plane of the sky to highlight the closed-flux-streamer belt region of the stellar corona.

4 Simulation Results

4.1 Carrington-scale Eruptive Stellar Flare and CME

The eruption process in our simulation follows the CSHKP scenario described by Lynch et al. (2016) and references therein. The energized field slowly rises due to the force imbalance set up by the increased magnetic pressure of the sheared-flux core. A radial current sheet forms underneath the expanding sheared flux. Eventually, fast magnetic reconnection sets in at the current sheet facilitating the rapid release of stored magnetic energy, ejecting a coherent twisted flux rope, and rebuilding the closed-flux system as the post-eruption flare arcade. Due to the global scale of our eruption, the post-eruption flare arcade becomes the closed-flux streamer belt encircling the entire star. The global eruption and evolution of the super-CME are shown in Figure 5 and its animation. The 3D structure of the erupting stellar CME flux rope can also be visualized from different vantage points in the ecliptic plane corresponding to the viewpoint of white-light coronagraph observations. Figure 6 shows a snapshot of our stellar CME flux rope at $t=150$ hr from three different meridional perspectives in the ecliptic plane: (a) $\phi=+120^{\circ}$, (b) $\phi=0^{\circ}$, and (c) $\phi=-120^{\circ}$. The plane-of-the-sky contour plots show the electric current density magnitude on a logarithmic scale to illustrate the CME cross sections. Representative field lines are chosen to illustrate the transition of the magnetic structure of the pre-eruption sheared flux into the erupting CME flux rope and its propagation through the outer stellar corona. Figure 6 is included as an animation.

The typical three-part structure of CMEs in coronagraph observations—a bright leading circular front, a dark circular cavity, and a bright central or trailing core region—is one of the best proxy measures of the magnetic structure of flux-rope CMEs (Vourlidas et al., 2013). Figure 6 shows this characteristic magnetic structure in the plane-of-the-sky cross-sections of our global CME eruption. The figure also shows the magnetic structure of the flux rope propagating towards the observer in the ecliptic plane—the configuration of most halo-CME eruptions that impact the Earth and cause significant geomagnetic responses (Zhang et al., 2007). Understanding this connection between the pre-eruption magnetic configuration of the CME source region and the CME’s structure and evolution during the eruption and propagation through the heliosphere is of critical importance to terrestrial space-weather forecasting (Palmerio et al., 2018) and will play an increasingly important role in characterizing exoplanetary space weather (Airapetian et al., 2016b; Cohen et al., 2018).

Figure 7 plots the evolution of the global magnetic energy ($E_{M}$, black) and kinetic energy ($E_{K}$, red). Here the energization phase is indicated with the light gray line. $E_{M}$ increases as the applied shearing flows drive the accumulation of free magnetic energy. The light-blue shaded region, defined as $t_{\rm pre}\leq t\leq t_{K}$ where $t_{\rm pre}=143.58$ hr is time of maximum $E_{M}$ and $t_{K}=150.83$ hr is the time of maximum $E_{K}$, indicates the impulsive phase of the eruption where the energy conversion is most rapid. Defining $\Delta E_{M}(t)\equiv E_{M}(t_{\rm pre})-E_{M}(t)$ and $\Delta E_{K}(t)\equiv E_{K}(t)-E_{K}(t_{\rm pre})$, the total magnetic energy released by the end of the simulation, $t_{f}=163.17$ hr, is $\Delta E_{M}(t_{f})=7.13\times 10^{33}$ erg, the maximum increase in kinetic energy is $\Delta E_{K}(t_{K})=2.84\times 10^{33}$ erg, and the magnetic-to-kinetic energy conversion ratio is $\Delta E_{K}(t_{K})/\Delta E_{M}(t_{K})=68.5\%$ during the impulsive phase and $\Delta E_{K}(t_{K})/\Delta E_{M}(t_{f})=39.9\%$ over the entire eruption process. The global scale of our eruption means that the stellar flare current sheet is reasonably well-resolved; consequently, our energy ratio is very similar to the $\approx$ 30% found by Karpen et al. (2012) for a global eruption with adaptive mesh refinement. We note that the CME associated with the famous Carrington event of 1859 was estimated to have a kinetic energy $\Delta E_{K}\sim 2\times 10^{33}$ erg (Cliver & Dietrich, 2013).

4.2 Stellar Flare Reconnection Flux

The flare ribbons and reconnection flux in the simulation data are calculated using a version of the Kazachenko et al. (2017) methodology used to characterize two-ribbon flares in SDO data. However, instead of using SDO/AIA 1600Å emission thresholding to determine whether a pixel is “flaring,” here we use the change in field line length $\Delta L\equiv L(t)-L(t-\Delta t)$ as a proxy for the rapid geometric re-configuration of the field-line connectivity between simulation data output intervals of $\Delta t=20$ min. We create a $192\times 384$ uniform grid in $(\theta,\phi)$ at the $r=R_{\odot}$ lower boundary and trace magnetic field lines from these footpoints for each simulation output time. If a long field line becomes significantly shorter by $\Delta L\leq-3R_{\odot}$ over the $\Delta t$ interval between consecutive output files, then we consider that pixel to have undergone reconnection. The reconnection pixels are accumulated in time to create the cumulative ribbon area map.

Figure 8(a) shows the time evolution of the area on the stellar surface swept out by the two-ribbon flare. The color scale indicates the first time the magnetic flux bundle at a given pixel has reconnected through the flare current sheet and become part of the post-eruption flare arcade. Large two-ribbon flares are often characterized by the “zipper effect,” where the ribbons form and rapidly grow parallel to the source region PIL and then move more slowly away from the PIL in the perpendicular direction (Moore et al., 2001; Qiu, 2009; Linton et al., 2009; Aulanier et al., 2012; Priest & Longcope, 2017). This effect is clearly seen in the Figure 8 animation.

From the time series of ribbon area pixel mask shown in Figure 8(a), we can calculate the stellar reconnection flux as

[TABLE]

where the pixel area is the usual $dA_{jk}=R_{\odot}^{2}\sin{\theta_{j}}\Delta\theta\Delta\phi$. The reconnection rate is then calculated as the central-differenced time derivative.

Figure 9(a) shows the unsigned reconnection flux $\Phi_{\rm rxn}$ (black squares) and the reconnection rate $\dot{\Phi}_{\rm rxn}$ (blue diamonds) versus time, along with the global kinetic energy ($E_{K}$) shown in gray. The total unsigned reconnection flux at the end of the simulation ($t_{f}=162.67$ hr) is $\Phi_{\rm rxn}(t_{f})=2.26\times 10^{23}$ Mx, and the maximum reconnection rate is $\dot{\Phi}_{\rm rxn}=8.0\times 10^{18}$ Mx s*-1*. The impulsive phase of the eruptive flare, defined previously, is highlighted by the light-blue shading. The stellar CME is driven by the eruptive flare reconnection, as seen by the substantial $\dot{\Phi}_{\rm rxn}$ increase before the global kinetic energy increase.

The Kazachenko et al. (2017) power-law relationship between GOES X-ray flux and reconnection flux was determined from 3137 solar flares of classes $\geq$C1.0 analyzed by in the SDO observations between 2010–2016. It is given by $I_{\rm Xray}=\alpha\left(\Phi_{\rm rxn}/10^{21}\right)^{1.454}$, where $\alpha=2.19\times 10^{-6}$ W m*-2* Mx*-1.454*. Figure 10 plots each of the Kazachenko et al. (2017) $(I_{\rm Xray},\Phi_{\rm rxn})$ data points (gray diamonds) as well as the power-law fit (solid red line) for the range of their observations. The power law is continued (dotted red line) beyond the solar observations, and we have shown the location of our simulation reconnection flux on it as the red square. The calculated X-ray flux of $I_{\rm Xray}=5.8\times 10^{-3}$ W m*-2* corresponds to a X58 class flare comparable to the estimate for the 1859 Carrington event of an X45($\pm$5) class flare (Cliver & Dietrich, 2013).

4.3 Synthetic Flare X-Ray and EUV Emission

To compare with general morphology and evolution of solar flare observations, we construct two synthetic emission proxies representing hot (10 MK) and ambient (1 MK) emission that we will refer to as synthetic X-ray and EUV emission, respectively.

The first emission proxy uses the method developed by Cheung & DeRosa (2012). The synthetic differential emissivity proxy $d\varepsilon_{J}$ at a three-dimensional point in space $\boldsymbol{r}_{ijk}$ is taken to be proportional to the average squared current density over the magnetic field line $\langle J^{2}\rangle$, normalized by the total field-line length, $L$. The total synthetic emission intensity $\varepsilon_{J}$ is then the line-of-sight integral over this differential emissivity proxy. For simplicity, here we calculate the (spherical) differential emissivity cube as a $100\times 192\times 384$ grid in $(r,\theta,\phi)$ and take the line-of-sight to be in the radial direction for $r\in\left[1R_{\odot},2R_{\odot}\right]$. The synthetic emission is thus

[TABLE]

where

[TABLE]

for a magnetic field line that passes through $\boldsymbol{r}_{ijk}$ and is parameterized by a differential arc length $d\ell^{\prime}$. Figure 8(b) plots the logarithm of line-of-sight-integrated emission $\varepsilon_{J}$ at $t=162$ hr corresponding to the spatial extent of the stellar flare arcade, which has become the entire closed-flux region of the stellar corona after the eruption. We note that $\varepsilon_{J}$ is qualitatively most similar to the SDO/AIA 131Å emission where the filter bandpass captures lines of Fe VIII associated with transition region temperature of $\log{T}\sim 5.6$ but also contains a significant flaring-corona contribution from Fe XXI at $\log{T}\sim 7.0$ (Lemen et al., 2012). Therefore, we refer to our synthetic emission distribution and light curves based on $\varepsilon_{J}$ as synthetic X-ray emission, representing hot (10 MK) plasma temperatures.

Our second synthetic emission proxy, $d\varepsilon_{N}$ is taken as proportional to $n_{e}^{2}$ which we equate with the single-fluid plasma number density $n^{2}$ at point $\boldsymbol{r}_{ijk}$. This yields

[TABLE]

Figure 8(c) plots the logarithim of the line-of-sight-integrated emission $\varepsilon_{N}$. Since $\varepsilon_{N}$ is qualitatively similar to SDO/AIA 171Å emission, which is primarily an Fe IX contribution from upper transition-region/quiet or ambient coronal temperatures of $\log{T}\sim 5.8$, our light curves based on $\varepsilon_{N}$ represent ambient (1 MK) emission and which we refer to as synthetic EUV emission.

The Figure 8 animation also shows the spatiotemporal evolution of our synthetic flare emission proxies. We note that by $t=162$ hr, the global flare arcade encompasses the entire pre-eruption closed-flux region of the stellar corona.

To compare with general morphology and evolution of stellar flare light curve observations, we integrate the Figure 8(b) and 8(c) spatial emission distributions. The surface-integrated light curves derived from our synthetic X-ray and EUV proxies are calculated as

[TABLE]

and normalized by the pre-eruption values at $t=130$ hr. The total surface-averaged intensity light curves ($I_{J}$, $I_{N}$) are plotted in Figures 9(b) and 9(c) as black squares, respectively. The maximum (radial) line-of-sight-integrated emission curves in Figure 9(b) and 9(c) are taken as $\max\left[\varepsilon_{J}dA\right]$ and $\max\left[\varepsilon_{N}dA\right]$, respectively, and are plotted as solid lines (X-ray, red; EUV, green).

Our light curves all show a qualitative transition during the impulsive phase of the global eruptive flare, but each curve has a slightly different character. The X-ray mean intensity shows a $\approx$ 30% pre-eruption dimming and a $\approx$ 200% post-eruption brightening—reminiscent of the $\Phi_{\rm rxn}$; the maximum X-ray emissivity shows a highly variable pre-eruption enhancement transitioning to a less-variable post-eruption dimming of $\approx$ 40%. In contract, the EUV mean intensity shows a pre-eruption dimming of $\approx$ 15% that transitions during the flare to a modest enhancement of $\approx$ 18% afterwards; and the maximum EUV emissivity shows the sharpest rise just after the onset of the flare, an $\approx$ 80% enhancement, and the clearest post-eruption decay. We note that, while solar and stellar X-ray and EUV observations measure the combined temperature and density evolution during flares, the synthetic emission from our isothermal simulation is determined solely by the density component.

Stellar flare observations in X-ray, UV, optical, and radio wavelengths are generally consistent with the CSHKP understanding of eruptive solar flares (Hawley et al., 1995; Güdel et al., 2002; Osten & Wolk, 2015). Harra et al. (2016) discussed the solar SDO/EVE disk-integrated light curves of Fe ion spectral lines for signatures that could be applied to stellar flare observations. They found the lower ionization states of Fe show prompt dimming curves, whereas the higher ionization states show a rapid increase followed by a slow decay. The stellar observations of pre-flare dimming have competing physical interpretations (Leitzinger et al., 2014; Osten & Wolk, 2017). The origin of the pre-flare dimming in our mean intensity light curves can be understood from the spatiotemporal evolution of the emission distributions shown in the Figure 8 animation. As the pre-eruption closed-flux regions expand and the outer layers open into the stellar wind, the emission from the streamer belt decreases gradually. This is followed by the rapid CME-related dimming due to the eruption/opening-up of most of the remaining closed flux once the flare reconnection begins in earnest. The eruption-related dimming is coincident with the initial formation of the flare arcade that has significantly enhanced emission compared to the pre-flare configuration. As the post-eruption flare arcade grows, its enhanced emission becomes the dominant feature of the global mean intensity light curves.

4.4 Stellar CME and CME-driven Shock Parameters

CMEs are responsible for some of the most geoeffective space-weather impacts at Earth and other solar system bodies. The combination of southward-directed $B_{z}$ in the sheath and ejecta flux rope with the increased dynamic pressure in dense sheath regions driven by fast events can cause significant geomagnetic responses (Zhang et al., 2007; Li et al., 2018). Fast CMEs drive shocks, and these coronal and interplanetary shocks are often sites of substantial energetic particle acceleration (Lario et al., 2016; Luhmann et al., 2018). Our simulation results can be used to begin an investigation of stellar CMEs and CME-driven shocks, as a precursor to estimating exoplanetary space-weather.

Figures 11(a) and 11(b) show the distributions of velocity magnitude and number density, respectively, in the meridional plane at $\phi_{\rm cut}=+90^{\circ}$ during the global eruption. The black lines indicate the radial sampling we use to quantify the CME-driven shock properties. Figure 11(c) shows the shock density compression ratio (black) and velocity magnitude (red) for the $\phi_{\rm cut}=+90^{\circ}$ meridional plane. The density compression ratio is calculated as $d=n(r,t)/n(r,130)$ so that the upstream, unperturbed values $\approx 1$. The shock location is determined along the radial sampling trajectory by calculating the radial gradient of the compression ratio $\partial d(r,t)/\partial r$ and choosing the maximum (negative) value at the largest radial distance—shown as the vertical blue dotted line. The peak number density of the stellar CME sheath region is found as the maximum number density within a spatial window just downstream of the shock location (shown as the solid blue line).

CME-driven shock strengths can be estimated from white-light coronagraph observations. Ontiveros & Vourlidas (2009) examined a number of shock fronts driven by fast ($\geq$1500 km s*-1*) CMEs and showed their compression ratios ranged from 1.5–3 in the LASCO C3 field of view. Kwon & Vourlidas (2018) used multi-viewpoint STEREO observations of two fast halo CMEs to show how the compression ratios were strongest at the nose of the shock and weaker at the flanks. Manchester et al. (2008) presented a comparison of the CME-driven shock in their numerical simulations of the $\approx$ 2000 km s*-1* halo CME of 28 Oct 2003 with the coronagraph observations, showing good agreement with the observed compression ratio of $\approx$ 5 at $r\sim 15R_{\odot}$. In Figure 11(c), our compression ratio reaches $\approx$ 4 by $25R_{\odot}$, indicative of a strong shock.

Figure 12 shows the spatiotemporal evolution of physical properties at the location of the peak compression ratio of the CME-driven shock for the six radial cuts through our simulation at longitudes $\phi_{\rm cut}\in\{\pm 30^{\circ},\pm 90^{\circ},\pm 150^{\circ}\}$. The latitude of each radial cut was chosen to approximate the center of the erupting CME flux rope cross-section in each meridional plane.

Figure 12(a) shows the shock velocities, 12(b) the shock compression ratios, and 12(c) the shock densities and plasma frequencies. In (b) at some locations, the compression ratio significantly exceeds 4, which is the upper bound for a plasma with adiabatic index $\gamma=5/3$. This can occur in our isothermal model, for which $\gamma=1$ and the compression ratio is theoretically unlimited, as well as in the more general polytropic models, for which $1<\gamma\ll 5/3$ and the compression ratio has a large but finite upper bound. Irrespective of the energy model that is used, large compression ratios reliably indicate the occurrence of strong shocks with Mach numbers $M\gg 1$—the key requirement for shock acceleration of energetic particles.

The shock strength (compression ratio), the magnetic field strength, and the shock orientation angle with respect to the magnetic field are the main inputs into models for the shock-acceleration of solar or stellar energetic particles (SEPs; e.g. Zank et al., 2007). A greater shock strength and $B$ magnitude result in a higher energy roll-over of the spectrum (a “harder” spectrum) that increases the fluxes of the highest energy particles, e.g., protons $\gtrsim 50$ MeV. The ambient particles that are swept up by the shock and then accelerated are referred to as the seed particle population. The upstream density and turbulence determine the number of seed particles and energizes their distribution, respectively. The downstream turbulence/waves generated in quasi-parallel shocks and low cross-field diffusion at quasi-perpendicular shocks both contribute to greater SEP energization.

The higher density and velocities in the ambient solar and stellar winds of the young Sun and analogs like $\kappa^{1}Cet$, can drive shocks in the stream interaction regions, providing greater and more highly energized seed particle populations (Airapetian & Usmanov, 2016) that will enhance SEP production by CME-driven shocks. Stellar energetic particle fluxes may significantly impact exoplanetary atmospheres. For example, Airapetian et al. (2016a) modeled the atmospheric chemistry of the early Earth under large SEP event conditions, showing that precipitating energetic protons enhanced chemical reactions that convert molecular nitrogen, carbon dioxide, and methane into nitrogen oxide and hydrogen cyanide.

Solar flares and CMEs often generate radio emission signatures. For example, CME-driven shocks can produce Type II radio bursts (Gopalswamy et al., 2005; Schmidt & Cairns, 2016) during their propagation through the corona and interplanetary space and flare-accelerated electron beams often produce a variety of Type III bursts (Reid & Ratcliffe, 2014; McCauley et al., 2018).

To estimate the radio emission from electrons at our CME-driven shocks, we use the standard formulation (Huba, 2013) where the electron plasma frequency $f_{pe}$ depends on $n_{e}$ $\left[{\rm cm}^{-3}\right]$ as

[TABLE]

Figure 12(c) shows the number density and corresponding plasma frequency and we obtain a range of characteristic frequency drifts of $\Delta f/\Delta t\approx$ 0.6–2.0 MHz hr*-1* with initial plasma frequency emission in the 5–20 MHz range. Given the apparent difficulty in observing radio signatures of stellar CMEs thus far (e.g. Crosley & Osten, 2018; Mullan & Paudel, 2019), the use of modeling to constrain the parameter space of stellar wind and CME densities and velocities that should result in observable signatures will be increasingly important. We note that this is a zeroth-order estimate of radio emission frequencies associated with our MHD fluid properties. In general, in order to estimate synthetic radio fluxes that can be directly compared to observations, a number of additional factors are required including source geometry considerations, signal propagation effects, and the specific details of various non-thermal emission mechanisms, as discussed by Moschou et al. (2018) and references therein.

5 Discussion

Our simulation self-consistently models both gradual accumulation and rapid release of free magnetic energy during the eruption of a global-scale stellar CME from the measured background magnetic field of $\kappa^{1}Cet$. The total magnetic flux and energy contained in the equatorial streamer belt of a star provides a baseline estimate for the anticipated strengths of its eruptive superflares/superCMEs. The energization of the background field of $\kappa^{1}Cet$ and the ensuing impulsive phase of stellar flare reconnection in our simulation created and accelerated a global CME eruption with kinetic energy $\Delta E_{K}\approx 3\times 10^{33}$ erg and duration $\Delta t_{K}\approx$ 10 hr. This energy is similar to the estimated energies of both the extreme 1859 Carrington event from the Sun (Cliver & Dietrich, 2013) and a 2016 naked-eye superflare from the nearby M-dwarf star Proxima Centauri (Howard et al., 2018). Our results show that even the comparatively weak background magnetic fields of such stars, relative to the estimated strengths of the fields in starspots, can store sufficient energy to power detectable superflares. We found that the free energy that could be stored in our pre-eruptive configuration for $\kappa^{1}Cet$ reached $E_{M}(t_{\rm pre})/E_{M}(t_{\rm rel})-1\approx 63\%$. This is very close to the free energy ($66\%$) required to open to infinity all of the magnetic field lines of an elementary dipole configuration (Mikić & Linker, 1994). The amount of energy that was converted to kinetic energy of our eruption is $\Delta E_{K}(t_{K})/E_{M}(t_{\rm rel})\approx 22\%$, or approximately one-third of the stored free energy.

The energetic superflare on $\kappa^{1}Cet$ observed by Robinson & Bopp (1987) was estimated by Schaefer et al. (2000) to have released about $2\times 10^{34}$ erg, an order of magnitude more energy than our simulated event from the recently measured background magnetic field of the star (Rosén et al., 2016). It is interesting to note that the energy of the observed 1986 $\kappa^{1}Cet$ superflare is roughly twice the energy contained within its background field, $2\times 10^{34}$ erg vs. $1\times 10^{34}$ erg (cf. Figure 3). Correspondingly, the estimated energy of the historic Carrington 1859 solar flare (Cliver & Dietrich, 2013) and the energy of the Sun’s present-day background field (Yeates et al., 2018) are each smaller by about one order of magnitude, $2\times 10^{33}$ erg vs. $1\times 10^{33}$ erg. The similarity in the flare-to-background energy ratio suggests that the Robinson-Bopp superflare on $\kappa^{1}Cet$ may have been as extreme an event for that star as the Carrington flare was for the Sun. Verifying this correspondence would require obtaining and analyzing a long-duration database of flare observations of $\kappa^{1}Cet$ and similar stars. Encouragingly, Notsu et al. (2019) and collaborators are investigating precisely this subset of superflaring solar-type stars in the Kepler data.

There are two primary areas where our idealized modeling could be improved upon in future numerical simulations. First, the magnetic flux distribution in the ZDI stellar magnetogram captures the global structure of the background magnetic field but does not resolve the field strengths or spatial scales associated with starspots or active regions. Second, the isothermal model for the stellar atmosphere does not account for the full thermodynamic evolution of the plasma, which affects the density and temperature distribution of the background wind and also means that the flare-related energy deposition into bulk plasma heating during the eruption is not captured in the simulation. We discuss each of these issues, in turn, below.

Our energization process imparts shear and twist to the coronal field structure in exactly the place the Sun requires (i.e., localized above the polarity inversion line). This enables the system both to gradually accumulate free energy and to transition rapidly to an unstable, runaway eruption that removes this localized stressed field from the closed-field corona via magnetic reconnection. Thus, the energization and eruption-triggering processes in our simulation are completely generic and should be universally applicable, including to more realistic solar and stellar magnetic field configurations.

More realistic stellar active region fields will be stronger (100s–1000s G) and significantly more localized spatially, although potentially over areas larger than observed on the Sun. For example, Rucinski et al. (2004) used the Microvariability and Oscillations of Stars (MOST) observations of $\kappa^{1}Cet$ to derive the sizes of two large starspots which covered 1.4% and 3.6% of the stellar surface, respectively. These areas are roughly 10 times greater that the largest sunspot/AR areas ever observed on the Sun (Hoge, 1947). The magnetic energy estimates used in stellar flare analyses are typically of the form $E_{M}\sim(1/(8\pi))B_{AR}^{2}A_{AR}^{3/2}$ (e.g. Shibata et al., 2013; Maehara et al., 2015; Notsu et al., 2019). If we were to scale the magnetic energy of our system, $E_{M}(t_{rel})\sim 1.3\times 10^{34}$ erg, to a starspot area estimate of $A_{AR}\sim f(\pi R_{\star}^{2})$ where $f=0.002-0.02$, then we obtain a range of areas $A_{AR}=3\times 10^{19-20}$ cm2 and a resulting range of starspot/AR field strengths, $B_{AR}\sim 250-1400$ G. These field strengths are comfortably within the range of both solar observations (e.g. Figure 10 of Kazachenko et al., 2017) and the Saar & Baliunas (1992) estimate of $350-500$ G for the $\kappa^{1}Cet$ surface-averaged field magnitudes.

We emphasize that our approach represents an attempt to model the most extreme stellar space-weather event possible within the observational constraints imposed by the surface magnetic-flux distribution from the ZDI reconstruction. Despite the global spatial scale and the stellar magnetogram’s unresolved starspot/AR flux distribution, the magnetic energy stored and released in our simulation’s eruptive superflare is compatible with the order-of-magnitude estimate one obtains from typical starspot/AR areas and magnetic field strengths. On the other hand, if much stronger unresolved starspots/AR flux distributions were assumed, the resulting superflare and super-CME energies could be much higher. For example, recent analyses by Okamoto & Sakurai (2018) of Hinode spectropolarimetry data of NOAA AR 11967 in 2014 Feb 1–6 showed a peak magnetic field strength of 6.2 kG. This solar magnetic field magnitude over the range of starspot/AR areas above would result in much larger magnetic energy estimates of $E_{M}\sim 2.5\times 10^{35}-8\times 10^{36}$ erg.

Our use of an isothermal stellar atmosphere neglects the internal energy equation in the MHD system and consequently, the evolution of our gas pressure is determined solely by the variation in mass density. In more complex thermodynamic MHD models, the internal energy equation typically includes field-aligned heat conduction, radiative cooling, magnetic and viscous dissipation, and a source term for the local contribution from coronal heating. One such parameterized coronal heating model represents the dissipation of energy associated with Alfvénic turbulence, which contributes both to local plasma heating and, through the wave pressure gradient, to the acceleration of the solar wind (van der Holst et al., 2014; Lionello et al., 2014; Oran et al., 2017). While a polytropic stellar wind outflow naturally creates a significant mass density gradient between the open-flux regions (coronal holes) and the closed-field corona (the helmet streamer belt), a more realistic treatment of the internal energy equation is expected to increase this density contrast for both the quiet-Sun and active regions (e.g. Török et al., 2018). The isothermal model for our stellar atmosphere has two main consequences for our simulation results: first, on the structure and intensities associated with our synthetic EUV and soft X-ray emission profiles; second, on the interaction between the CME and the background wind, including the properties of the CME-driven shock.

As discussed above, approximately one-third of the free magnetic energy stored prior to the eruption is converted to kinetic energy of the CME. About one-sixth of the stored energy remained in the closed coronal field after the eruption; that is, the magnetic field did not relax completely to the initial, minimum-energy state of the system. The remaining energy released by the magnetic field during the eruption, $\Delta E_{H}/E_{M}(t_{\rm rel})\approx 32\%$ of the initial magnetic energy, was not captured by our simple isothermal model. This remnant released energy $\Delta E_{H}$ would appear as thermal energy due to magnetic and viscous dissipation, heating the stellar plasma to high flare temperatures and being radiated away into space. With a more realistic temperature structure in the dynamic formation and evolution of the flare current sheet and the post-eruption arcade (e.g. Reeves et al., 2010; Lynch et al., 2011), one could improve the synthetic emission calculations of $\S$4.3 by calculating the density- and temperature-dependent emission intensities in various spectral lines and convolving these with instrument response functions (e.g. Lionello et al., 2009; Reeves et al., 2010; Shen et al., 2013; Oran et al., 2017; Jin et al., 2017).

Relaxing the isothermal assumption would result in a modified stellar wind profile and, thus, quantitative differences in the plasma properties of the CME–stellar wind interaction region. Qualitatively, the overall picture will remain the same as in Figure 11, i.e., the energetic eruption driven by the flare reconnection will create and accelerate a rapidly expanding, highly-magnetized (low $\beta$) flux rope ejecta which will generate a shock and a dense compression region at the CME’s leading edge. A modified upstream stellar wind profile will obviously impact the absolute number densities and an improved thermodynamic treatment will allow for additional localized, compressional heating. The shock is expected to remain Alfvénic, however the detailed shock parameters such as the compression ratio, field strength, and the shock normal orientation will have different quantitative values. Even with the simple isothermal model used here, the compression ratios generated in our simulation are broadly consistent with those determined from the white-light coronagraph observations of large solar events.

Therefore, limitations notwithstanding, our simulation results represent an important step toward understanding observations of superflares from active solar-type stars and characterizing CMEs from magnetically active stars across the K–M dwarf-star spectrum. We analyzed the properties of CME-induced shocks, including the frequency and duration of associated Type II events that may be detected in future low-frequency (10 MHz or lower) radio observations of magnetically active stars. The derived properties of CMEs and their associated shocks can provide inputs to models of stellar SEP energization and transport via the diffusive shock acceleration mechanism (Li et al., 2013; Hu et al., 2018). Knowledge of the SEP energy spectrum and particle fluence is critical for evaluating biogenic conditions on terrestrial-type exoplanets around active stars, as well as for the early Earth and Mars. Airapetian et al. (2016b) have shown that large SEP fluxes can increase the production rates of nitrous oxide, a powerful greenhouse gas, and hydrogen cyanide, a feedstock molecule for prebiotic synthesis. Characterizing the SEP environment will help to determine the boundaries of the planetary “biogenic” zone (Airapetian et al., 2019b) and will be important in specifying the efficiency of ozone destruction and surface dosages of ionizing radiation that are damaging to life on planetary surfaces. Future modeling efforts that describe the conditions of CME initiation in more realistic stellar magnetic field distributions can help to understand observations of G, K, and M dwarf stars using the currently implemented international multi-observatory program (TESS, Hubble, XMM-Newton, Apache Point Observatory) and to prepare for next-generation JWST observations.

BJL, MDK, and WPA acknowledge support from NASA NNX17AI28G, NSF AGS-1622495, and the Coronal Global Evolutionary Model (CGEM) project NSF AGS-1321474. VSA acknowledges support from the NASA Exobiology 80NSSC17K0463 and the TESS 1 Cycle 1 program. CRD acknowledges support from NASA H-SR and LWS TR&T programs. TL acknowledges support via the Austrian Space Application Programme (ASAP) of the Austrian Research Promotion Agency (FFG) within ASAP11 and the FWF NFN project S11601-N16. OK acknowledges support by the Knut and Alice Wallenberg Foundation (project grant “The New Milky Way”), the Swedish Research Council (project 621-2014-5720), and the Swedish National Space Board (projects 185/14, 137/17). The computational resources for this work were provided to BJL by the NASA High-End Computing Program through the NASA Center for Climate Simulation at Goddard Space Flight Center.