Stellar energetic particles in the magnetically turbulent habitable zones of TRAPPIST-1-like planetary systems

TL;DR

This study uses simulations to analyze how energetic particles propagate around TRAPPIST-1-like systems, revealing high proton fluxes on habitable planets and potential natural containment mechanisms.

Contribution

It provides the first detailed simulation-based analysis of energetic particle propagation in TRAPPIST-1-like systems, highlighting the effects of turbulence and magnetic focusing.

Findings

Few percent of particles escape near the star surface.

Escaping particles are focused onto two caps and the equatorial plane.

Proton flux on TRAPPIST-1e could be up to 10^6 times Earth's flux.

Abstract

Planets in close proximity to their parent star, such as those in the habitable zones around M dwarfs, could be subject to particularly high doses of particle radiation. We have carried out test-particle simulations of ~GeV protons to investigate the propagation of energetic particles accelerated by flares or travelling shock waves within the stellar wind and magnetic field of a TRAPPIST-1-like system. Turbulence was simulated with small-scale magnetostatic perturbations with an isotropic power spectrum. We find that only a few percent of particles injected within half a stellar radius from the stellar surface escape, and that the escaping fraction increases strongly with increasing injection radius. Escaping particles are increasingly deflected and focused by the ambient spiralling magnetic field as the superimposed turbulence amplitude is increased. In our TRAPPIST-1-like simulations, regardless of the angular region of injection, particles are strongly focused onto two caps within the fast wind regions and centered on the equatorial planetary orbital plane. Based on a scaling relation between far-UV emission and energetic protons for solar flares applied to M dwarfs, the innermost putative habitable planet, TRAPPIST-1e, is bombarded by a proton flux up to 6 orders of magnitude larger than experienced by the present-day Earth. We note two mechanisms that could strongly limit EP fluxes from active stars: EPs from flares are contained by the stellar magnetic field; and potential CMEs that might generate EPs at larger distances also fail to escape.

Figures

Equations

$$\mathrm{k}{r}_g(\mathbf{x})/2\pi =r_g(\mathbf{x})/L_c<1,$$

$${\sigma }^2=(\delta B(\mathbf{x})/B_0(\mathbf{x}){)}^2.$$

$${\lambda }_{\parallel }(\mathbf{x})\simeq 4.8(r_g(\mathbf{x})/L_c{)}^{1/3}{L}_c/{\sigma }^2.$$

$${\mathbf{B}(\mathbf{x})=\mathbf{B}}_0(\mathbf{x})+\delta \mathbf{B}(\mathbf{x}),$$

$$F_{\mathrm{inj}}(R_s)={\left(\frac{r_{876}^{HZ}}{R_s}\right)}^2{F}_{876}^{max}\simeq 10^3\times F_{876}^{max}\simeq 10^6\frac{\mathrm{protons}}{{\mathrm{cm}}^2\mathrm{s}\mathrm{sterad}}.$$

$$F(R_p)=\frac{N_{R_p}^{\prime }}{N_{\mathrm{inj}}}\frac{F_{\mathrm{inj}}(R_s)}{A}$$

Full text

…

††journal: ApJ

1 Introduction

The definition of planet habitability has been based in the last decades on the orbital distance (or habitable zone, hereafter HZ, Kasting et al., 1993) at which the steady stellar irradiation allows for a temperature consistent with the presence of liquid water on the planetary surface. However, charged energetic particles (hereafter EPs) produced by stellar flares or shock waves driven by Coronal Mass Ejections (hereafter CMEs) and travelling into the interplanetary medium may significantly impact the conditions for life to exist in planets beyond the solar system (Segura et al., 2010; Ribas et al., 2016; Lingam & Loeb, 2018).

In the case of the solar wind, in-situ measurements of EP irradiation are used to assess shielding requirements for astronauts at $1$ AU (Mewaldt, 2006; Mewaldt et al., 2007). Multi-spacecraft observations of solar eruptive events during the solar maximum of cycle $23$ ($2002-2006$) show that between $0.4$ and $20\%$ of the kinetic energy of CMEs in the energy range $10^{31}-10^{32}$ erg (in the solar wind frame) is expended in accelerating solar EPs (Mewaldt et al., 2008; Emslie et al., 2012).

Stellar EPs are in some cases expected to cause depletion of planetary ozone layers (Segura et al., 2010; Tilley et al., 2017). Such depletion allows penetration of UV radiation with consequent degradation of proteins (Kerwin & Remmele, 2007) but also, in contrast, catalysis of pre-biotic molecules (Airapetian et al., 2016; Lingam et al., 2018). Loyd et al. (2018) note that ozone depletion by photolysis alone was expected to be significant only for very major flares expected to occur monthly or yearly, but note that effects of very commonly occurring weaker flares in their study could be enhanced by EPs. Such multiple lines of evidence suggest that EPs are a component of the star/planet interaction worthy of detailed investigation in relation to habitability.

Propagation of EPs from the injection location to a planet is mediated by the large-scale and the turbulent components of the stellar magnetic field. Studies of the effect of EPs on the ionization of protoplanetary disks (Turner & Drake, 2009) or on the synthesis of short-lived nuclides in the early solar system (see, e.g., Dauphas & Chaussidon, 2011) assumed that EPs propagate rectilinearly, unimpeded by the magnetic field structure. However, both the components of the magnetic field have been shown to lead to an efficient confinement of EPs close to young active stars (see, e.g., Fraschetti et al., 2018).

M dwarfs, the most abundant and long-lived stars in the Milky Way, are currently among the primary targets in exoplanet searches. This is largely due to their small radius that increases the likelihood of detecting orbiting Earth-sized planets with transit techniques, or due to their low mass compared with other spectral types that increases a planet-induced radial velocity Doppler shift in the stellar spectrum.

Youngblood et al. (2017) have recently used the MUSCLES (Measurements of the Ultraviolet Spectral Characteristics of Low-mass Exoplanetary Systems) Treasury Survey (France et al., 2016) to determine that large flares on M dwarfs, i.e., with a soft $X$-ray (hereafter SXR) peak flux $\geq 10^{-3}$ W m*-2* at $1$ AU or class X10.0 in the GOES (Geostationary Operational Environmental Satellite) classification, lead to a $>10$ MeV proton flux on planets in the HZ up to $\sim 4$ orders of magnitude higher than the present-day Earth.

Likewise, the assumption of a solar-like correlation for T Tauri stars between peak emission of large flares ($X$-ray luminosity $>10^{30}$ erg s*-1*) and energetic proton enhancements (Feigelson et al., 2002; Turner & Drake, 2009) leads to suggest an enrichment by $\sim 4$ orders of magnitude over the present-day proton density at $1$ AU. These fluxes imply that the ionization of protoplanetary disks can locally exceed ionization due to stellar $X$-rays as a result of EPs being channeled and concentrated by magnetic turbulence (Fraschetti et al., 2018).

Such cases show that the EPs emitted by stars more active than the Sun can play a crucial role in the evolution of the circumstellar medium, or inner “astrosphere” (here within $\sim 100$ stellar radii), and potentially in the habitability of exoplanets. However, while active stars might generate copious EPs, it is necessary to understand how they propagate within the stellar and interplanetary magnetic field in order to assess their potential impact.

The seven Earth-sized transiting exoplanets recently discovered in the TRAPPIST-1 system (Gillon et al., 2017) are surprisingly packed within a distance of $0.062$ AU from the host star (Delrez et al., 2018). Three planets (TRAPPIST-1e, f, g) have been found to orbit the HZ, that spans the range $\sim 0.029-0.047$ AU (Delrez et al., 2018), raising the question whether the enhanced EP flux at such a close distance affects the atmosphere and planetary habitability.

In this work we determine the flux of EPs impinging onto the HZ planets in the TRAPPIST-1 system by using a realistic and turbulent magnetized wind model of an M dwarf star proxy for the yet poorly-constrained wind of TRAPPIST-1. We adopt the extended magnetic field structure computed using a three-dimensional magnetohydrodynamic (MHD) model previously calibrated to the solar wind and recently applied to study the coronal structure, winds, and inner astrospheres of Sun-like stars (Alvarado-Gómez et al., 2016a, b) and M-dwarfs (Garraffo et al., 2016, 2017), together with the propagation of EPs in stellar turbulence (Fraschetti et al., 2018). We directly solve for the propagation of individual EPs in the turbulent inner astrosphere of an M dwarf wind. The turbulence is calculated via the prescription defined in Giacalone & Jokipii (1999); Fraschetti & Giacalone (2012).

In section 2, general properties of the MHD model simulations are outlined. Section 3 describes the assumptions adopted regarding EP propagation and the magnetic turbulence. Section 4 presents the numerical model. Section 5 contains the main results and 6 quantifies the flux impinging on the HZ planets in the TRAPPIST-1 system. Discussion and conclusion are in Sections 7 and 8, respectively.

2 TRAPPIST-1 MAGNETOSPHERIC MODEL

TRAPPIST-1 is a low-mass M dwarf ($0.089M_{\odot}$) with a $3.3$ day rotation period and a radius $R_{*}\sim 0.114~{}R_{\odot}$ according to the latest observations (Luger et al., 2017). It was confirmed to host seven planets orbiting in a co-planar system (within $\sim 30$ arcmin) viewed nearly edge-on (Gillon et al., 2017). All the planets reside close to the host star, with semi-major axes from 0.01 AU to 0.062 AU (Mercury orbits at 0.39 AU), with orbital periods from 1.5 days to 20 days.

As a background medium for studying the propagation of EPs within the TRAPPIST-1 system, we adopt the wind and magnetosphere model computed by Garraffo et al. (2017) using the 3D MHD code Block Adaptive Tree Solarwind Roe Upwind Scheme (BATS-R-US, Powell et al., 1999; Tóth et al., 2012), in the version that incorporates the Alfvén Wave Solar Model (AWSoM) (van der Holst et al., 2014). A data-driven global MHD method is used that was initially developed to reconstruct the solar atmosphere and the solar wind. BATS-R-US employs a radial field magnetogram as a boundary condition for the stellar photospheric magnetic field. In the case of application to the Sun, this is a solar magnetogram but stellar magnetograms obtained using the Zeeman-Doppler Imaging technique (Donati & Brown, 1997) can also be used.

Zeeman-Doppler Imaging is presently limited to luminous, fairly rapidly rotating stars. TRAPPIST-1, despite its relatively fast spin, is optically faint ($M_{v}=18.8$, Gillon et al., 2017) and out of reach of current Zeeman-Doppler Imaging capabilities. Unfortunately both the distribution of the magnetic field on its surface and the direction of the rotation axis are unknown due to the extreme faintness of the star; moreover, both are subject to change in time with time scales of years to greater, due to the periodic change of magnetic polarity and to the axis precession, respectively. Its average magnetic field, however, has been estimated to be $\sim 600$ G using Zeeman broadening (Reiners & Basri, 2010). There is growing agreement that the geometry of the magnetic field depends on the rotation period and spectral type of the star (Vidotto et al., 2014; Garraffo et al., 2015; Réville et al., 2015; Finley & Matt, 2018). Garraffo et al. (2017) therefore used as a proxy for TRAPPIST-1 the magnetogram observed for GJ 3622 (Morin et al., 2010), an M4 dwarf with a rotation period of $1.5$ days. The field on its surface reaches a maximum of $1.4$ kG, yielding an average field of $\sim$600 G, consistent with the TRAPPIST-1 observations. The magnetic structure is not expected to change significantly between stars with periods of 1 to 3 days. We note that our approach is different to that of Dong et al. (2018), who estimated the ion escape rate in the seven planets using a wind model based on a solar magnetogram under solar minimum conditions, rescaled to a magnetic field strength more like typical M-dwarf values (Morin et al., 2010).

The GJ 3622 magnetic field is vaguely dipolar with a notable misalignment between the rotation axis and the magnetic field amounting to a few tens of degrees ($\sim 40^{\circ}$ - $50^{\circ}$). The wind and magnetosphere model are illustrated in Figure 1.

3 Stellar energetic particles in the TRAPPIST-1 environment

3.1 General assumptions on EPs: origin and propagation

Our general goal here is to explore the effect of small-scale magnetic turbulence on the propagation of EPs through the magnetosphere of the host star TRAPPIST-1, and as far as the outermost planet located at a distance of $\sim 0.062$ AU. In particular, we focus on a comparison of the EP flux generated at the star itself with that which propagates out to planets 1b, 1e and 1h.

Two processes are assumed to produce the non-thermal particles (Fraschetti et al., 2018): 1) shock waves driven by CMEs, travelling in the interplanetary medium and therein accelerating and releasing EPs; 2) flares occurring within the stellar corona and releasing EPs within a small distance from the stellar surface ($\sim 0.5R_{*}$). Both such processes are assumed to produce the $\sim$ GeV kinetic energy protons studied here. This assumption can be justified by a solar analogy: former GOES measurements correlating solar proton enhancements at $1$ AU with SXR flares do not unequivocally pinpoint the flares as the only sources of particle acceleration as CME-driven shocks are consistent with such a correlation as well (Belov et al., 2007).

In our simulations only the location of injection of EPs (at a distance $R_{s}$ from the star), rather than the acceleration mechanism, is assigned. As for the abundance of accelerated particles in the circumstellar medium at a given distance from the host star, 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). This scaling provides a time-averaged EP enrichment for time scales comparable with a statistically typical flare duration (Vida et al., 2017) .

We calculate the propagation of the EPs using a test particle approach within a realistic representation of the interplanetary medium that includes magnetic field fluctuations. The large-scale structure used here for the TRAPPIST-1 magnetic field (see Fig. 1) has an approximately dipolar structure with no significant field lines wrapping around the star as might be expected for T Tauri stars and some fast rotators (see, e.g., Gregory et al., 2009; Cohen et al., 2010; Fraschetti et al., 2018). Nevertheless, it is still uncertain whether the average $\sim$ kiloGauss magnetic field of TRAPPIST-1 allows for CME escape and the outward driving of EPs accelerated at shocks (Drake et al., 2016; Osten & Wolk, 2015). Under the assumption that EPs can be steadily supplied by flares and CMEs, the dominant magnetic effects we are concerned with for EP propagation in TRAPPIST-1 are expected to be scattering and perpendicular diffusion in the turbulent stellar field.

The MHD wind solution and the magnetic turbulence are stationary on the time-scale of EP propagation to a good approximation. The EPs travel at speed $\simeq c$, whereas the stellar rotation speed close to the surface is $\sim 2$ km s*-1*, and the Alfvén wave speed in the circumstellar medium is $\sim 10^{4}$ km s*-1* ($\sim 10^{3}$ km s*-1*) at a distance $\sim 10R_{*}$ ($110R_{*}$, semi-major axis of the outermost planet) from the host star. This holds for M dwarfs in general. The visible light periodograms of M dwarfs (with radii in the range $0.08-0.6R_{\odot}$)—presumably dominated by rotational modulation signatures—typically peak at a few days over a range of periods $\sim 1$–100 days (Hawley et al., 2014), with a corresponding surface rotation speed over a range $0.04-30$ km/s (Barnes et al., 2014; Jeffers et al., 2018). Only the earliest M dwarfs ($0.6R_{\odot}$) with rotation periods $\leq 3$ days have surface rotation speeds $>10$ km s*-1*. Dynamical timescales are therefore much longer than the EP travel time in our simulations (typically $<1$ hour).

3.2 Turbulent stellar magnetic field

In analogy with the measurements of interplanetary magnetic turbulence (e.g., Jokipii & Coleman, 1968), and of interstellar density turbulence (Armstrong et al., 1995), we prescribe a magnetic turbulence power spectrum having the shape of a power-law (Kolmogorov) in the 3D turbulent wavenumber ${\rm k}$ (see Fig. 2). Scale-dependent anisotropic turbulence (á la the Goldreich & Sridhar 1995, e.g. model) explaining the origin of the solar wind MHD-scale turbulence anisotropy (e.g., Horbury et al., 2008), has unsettled theoretical transport properties (Laitinen et al., 2013; Fraschetti, 2016a, b) and would require a more cumbersome numerical code.

The test-particle simulations presented here track naturally the pitch-angle scattering and cross-field motion of EPs caused by the small-scale turbulence. An alternative approach to EP transport involves Monte Carlo simulations that reproduce the pitch-angle scattering and neglect perpendicular transport (see, e.g., Ellison et al., 1981). The nearly radial spread-out of the open magnetic field lines of the astrosphere used here leads to an observable consequence (see Sect. 5) of the turbulent transport across field lines (Fraschetti & Jokipii, 2011; Strauss et al., 2017). In contrast, in the case of the T Tauri star studied in Fraschetti et al. (2018) the wrapping of magnetic field lines around the star prevented an assessment of the effect of the transport across field lines.

Due to the lack of observational estimates of the correlation length, or injection scale, $L_{c}$, of the magnetic turbulence within the circumstellar medium (see Fig. 2), we adopt the uniform value $L_{c}=10^{-5}$ AU throughout the simulation box. A simulation set carried out with a smaller uniform $L_{c}=10^{-6}$ AU shows that the statistical properties of EPs are not significantly affected by the choice of $L_{c}$, provided that the resonance condition is satisfied. In this regard, $L_{c}=10^{-5}$ AU is a reasonable value for the quite small range in radial distance of the planets in the TRAPPIST-1 system, within $0.062$ AU. The chosen value of $L_{c}$ ensures resonance with turbulent inertial scales at each EP energy considered (see Fig. 2) during their entire propagation. Such a condition reads

[TABLE]

for each wave-number ${\rm k}$ within the inertial range; here, $r_{g}({\mathbf{x}})=p_{\perp}c/eB_{0}({\mathbf{x}})$ is the gyroradius of a proton with momentum $p_{\perp}$ perpendicular to the unperturbed and space-dependent magnetic field $B_{0}({\mathbf{x}})$ of TRAPPIST-1, $e$ the proton electric charge and $c$ the speed of light in vacuum.

The power of the magnetic fluctuation $\delta B(\mathbf{x})$ relative to $B_{0}(\mathbf{x})$ is defined as

[TABLE]

Here, $\sigma^{2}$ is assumed to be independent of space throughout the simulation box as well. The spherical average of the unperturbed field $\langle B_{0}({\bf x})\rangle_{\Omega}$ produced by the 3D-MHD simulations (see Sect. 2) drops with radius $R$ from $2R_{*}$ as $\sim R^{-2.2}$. On the other hand, the solar wind measurements yield for the turbulence amplitude $\delta B$ between $0.3$ and $4$ AU a power-law dependence on heliocentric distance with a very similar index ($\simeq−2.2$) at a variety of helio-latitudes (Horbury & Tsurutani, 2001). Thus, in the lack of any current measurement of the magnetic turbulence around TRAPPIST-1, it seems reasonable to assume a uniform $\sigma^{2}$, following Fraschetti et al. (2018).

The turbulence might be generated by the stirring of the plasma at the outer scale $L_{c}$, followed by a cascade, or by plasma instabilities at kinetic scales generated, e.g., by streaming of EPs along the field; we neglect the latter here as we are restricted to the test-particle limit. The turbulence within the violently active M dwarf magnetosphere is likely to be much stronger than that in the solar wind ($\sigma^{2}$ not greater than $0.1$, Burlaga & Turner, 1976). Thus, we considered values of $\sigma^{2}$ spanning the range $0.01-1.0$. The interpretation of our simulations makes use of the scattering mean free path, $\lambda_{\parallel}$, given by quasi-linear theory (Jokipii, 1966), that reads (Giacalone & Jokipii, 1999; Fraschetti et al., 2018)

[TABLE]

The choices of uniform $L_{c}$ and $\sigma^{2}$ imply that $\lambda_{\parallel}$ depends on spatial coordinates only via $r_{g}(\mathbf{x})$ (i.e., $B_{0}(\mathbf{x})$).

4 Numerical Method

In our numerical experiments, we have directly integrated the trajectories of $\sim 10^{4}$ energetic protons propagating in a turbulent magnetic field that can be decomposed as

[TABLE]

where the large-scale component, ${\bf B}_{0}({\bf x})$, is the 3D magnetic field generated by the 3D-MHD simulations as calculated in Garraffo et al. (2017) and described in Section 2; the random component ${\bf\delta B}={\bf\delta B}(x,y,z)$ has a zero mean ($\langle\delta{\bf B(x)}\rangle=0$). Here ${\bf\delta 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). We use an inertial range $k_{\rm min}<k<k_{\rm max}$, with $k_{\rm max}/k_{\rm min}=10^{2}$, where $k_{\rm min}=2\pi/L_{c}$ and $k_{\rm max}$ is the magnitude of the wavenumber corresponding to some turbulence dissipation scale. In Fraschetti et al. (2018) we verified that an inertial range extended by one decade to smaller scales does not substantially change the resulting distribution of a large number of EPs hitting a protoplanetary disk, despite being computationally much more expensive; we assume that a larger inertial range is not relevant for the M dwarf circumstellar turbulence either.

The turbulence power spectrum within the inertial range (Fig. 2) is assumed to be a three-dimensional Kolmogorov power-law (index $-11/3$). At scales larger than $k_{\rm min}^{-1}$ ($k_{0}<k<k_{\rm min}$), the power spectrum is taken as constant (see, e.g., Jokipii & Coleman (1968) for the solar wind case).

In our simulations, the EPs are injected uniformly on spherical surfaces at a variety of radii, $R_{s}$, with a velocity distribution isotropic in pitch-angle. The number of particles is then rescaled by using the enhancement in EP flux inferred at a given distance from the star in Youngblood et al. (2017). After propagation through the inner astrosphere, the EP angular location is recorded on spherical surfaces at distances $R_{p}$. We verified that the particle energy is conserved to a relative accuracy of  $10^{-3}-10^{-4}$.

5 Results

5.1 Particle Trajectories

Figures 3 to 5 show the trajectories of selected individual EPs injected at $R_{s}=10R_{\star}=0.0056$ AU with $\sigma^{2}=1.0$. All EPs are allowed two possible fates in our simulations: hitting (in blue) the spherical surface at $R_{p}=R_{b,e,h}$, where $R_{b,e,h}$ equals the semi-major axes of the planets TRAPPIST-1b, e, h (respectively $20.4R_{*}=0.0115$ AU, $51.8R_{*}=0.02916$ AU and $110.R_{*}=0.0617$ AU (Delrez et al., 2018)), or collapsing (in red) back to the star.

5.2 Weak Turbulence Case

Figures  6, 7 and 8 show the spherical coordinates of the hitting points for 1 GeV (Fig. 6) and 10 GeV (Figs. 7 and 8) kinetic energy protons, injected at $R_{s}=10R_{\star}$ (Figs. 6, 7) and $R_{s}=5R_{\star}$ (Fig. 8) recorded at the spheres $R_{p}=R_{b,e}$. The total number of injected EPs ($N_{\rm inj}$) is the same in all cases. Different rows correspond to different values of $\sigma^{2}$, increasing from top to bottom; different columns correspond to a different planet, 1b (left) and 1e (right). The colorbar is scaled to the maximum number of EPs per pixel, and varies strongly between panels; thus, the same color in different panels does not indicate the same absolute number of EPs. The plane $\theta^{\prime}=\theta+90^{\circ}=90^{\circ}$ perpendicular to the stellar rotation axis, where $-90^{\circ}<\theta<90^{\circ}$ is the latitude, marks the plane of the planetary orbits within $30$ arcmin (Delrez et al., 2018).

In Figs. 6 and 7, for weak turbulence ($\sigma^{2}=0.01$, upper row), the distribution of hitting points spreads fairly uniformly over the $R_{p}$-sphere. Such a distribution mirrors the uniform distribution of the injection points of EPs and results from the EPs propagation outward close to the scatter-free limit, i.e., uniform and static electric and magnetic field, along the open field lines intercepted on the sphere at $R_{s}$ (greater $\lambda_{\parallel}$ for small $\sigma^{2}$, from Eq. 3).

The perpendicular diffusion coefficient $\kappa_{\perp}$ grows, regardless of the model, as $\kappa_{\perp}\sim\sigma^{2}$ (Giacalone & Jokipii, 1999; Fraschetti & Jokipii, 2011; Strauss et al., 2017) leading to a negligible decorrelation of EPs, for small $\sigma^{2}$ from the direction of the average magnetic field. Thus, the resulting distribution of hitting points at $R_{p}$ is close to the injection distribution at $R_{s}$ and the trajectories nearly map the unperturbed magnetic field ${\mathbf{B}}_{0}$. However, we note that the ratio of the number of EPs at $R_{p}$-sphere ($N_{R_{p}}$) to $N_{\rm inj}$ is limited to $20-25\%$ (see also Fig. 9, left panel), as a large fraction ($75-80\%$) collapse back to the star. The latter EPs are released on closed field lines that are prevalent at $R_{s}=10R_{*}$ (see Fig. 1), and propagate along those closed field lines back to the star, due to the large $\lambda_{\parallel}$ (see Eq. 3) and negligible perpendicular diffusion.

We also note in Fig. 9, left panel, that for each value of $\sigma^{2}$ the ratio $N_{R_{p}}/N_{\rm inj}$ decreases for greater $R_{p}$, i.e., decreasing from 1b (red) to 1h (blue). This occurs because some EPs that propagate past an inner $R_{p}$-sphere undergo pitch-angle diffusion that leads them to move backward and to collapse to the star without reaching the outer $R_{p}$-sphere. In addition, Fig. 9, left panel, shows a smaller difference for each value of $\sigma^{2}$ between the blue and green curves as compared with green and red ones: this change results from the transition of the large-scale $B_{0}$-field structure from closed/open to prevalently open field lines between the 1b (red) and 1e (green), whereas between 1e and 1h (blue) all field lines are open (cfr. Fig. 1), so that no significant difference is expected between the green and the blue curves. We note that the likelihood of backward trajectories decreases further out due to the increase of the mean free path: $\lambda_{\parallel}$ increases outward as $r_{g}^{1/3}\propto B_{0}^{-1/3}$ (see Eq. 3) for $B_{0}$ decreasing outward in a uniform $\sigma^{2}$, so most EPs channelled onto an open line that reach 1e will also reach 1h.

We have run an additional set of simulations with $R_{s}=1.5R_{*}$, i.e., at a distance of $0.5R_{*}$ from the stellar surface, for particles with $E=0.3$ GeV. For these simulations, negligible turbulence was adopted ($\sigma^{2}=10^{-8}$) since within the chosen turbulent inertial range the EPs would not scatter resonantly as $r_{g}$ is suppressed by the strong $B_{0}$ field close to the surface. We find that the ratio $N_{R_{p}}/N_{\rm inj}$ is in the range $3.0-3.7\%$ for $R_{p}=R_{b}$ or $R_{h}$.

5.3 Effect of Stronger Turbulence

The histogram on the $R_{p}$-sphere changes dramatically in the presence of stronger turbulence ($\sigma^{2}=0.1$, $1.0$, middle and lower row in Figs. 6, 7 and in Fig. 8): EP hitting points on the $R_{p}$-sphere are confined to equatorial caps. We find a depleted region, in white, that is barely discernible at $R_{p}=R_{b}$ but conspicuous at $R_{p}=R_{e}$, and that azimuthally oscillates in the middle and bottom rows in Figs. 6, 7 and 8. This arises from the inclination of the magnetic axis to the rotation axis, and traces the azimuthal variation of the slow wind (see the spherical map of the wind speed, upper row in Fig. 10).

Inspection of the structure of the average magnetic field (see Fig. 1) confirms that closed (open) field lines populate mainly the slow (fast) wind region. Moreover, a comparison of the middle row of Fig. 7 with Fig. 8 shows that injection further out ($R_{s}=10R_{*}$ rather than $5R_{*}$) reduces the chances of intercepting a closed field line due to the opening of field lines in the slow wind region as one proceeds outward. Consequently, the depleted white regions narrow down as the injection radius is increased from $R_{s}=5$ to $10$.

The broadening of the depleted regions as $\sigma^{2}$ increases, shown in the bottom rows of Figs. 6 and 7, can be explained as follows. A greater amplitude of magnetic fluctuation, i.e., greater $\sigma^{2}$, leads to a reduced $\lambda_{\parallel}$ (see Eq. 3) and to an enhanced perpendicular diffusion: EPs more frequently decorrelate via cross-field transport. Near the boundary between open and closed field lines, a fraction of particles diffusing from open onto closed field will collapse back to the star, depleting the region corresponding to the current-sheet. There is then a net migration from open to closed field due to this loss of particles at the stellar surface.

The diffusive motion in the opposite direction, i.e., from a closed field lines near the boundary to an open line, and subsequent escape is less likely due to smaller $B_{0}$ of the closed line regions (see Fig. 10, lower row), i.e. larger $\lambda_{\parallel}$, that might lead EPs rapidly to the stellar surface. Indeed, EPs can travel a short distance before falling to the star as the path length of the closed field lines is only a few times $\lambda_{\parallel}$ (from Eq. 3, a 10 GeV proton at $R_{s}=10R_{*}$, with $r_{g}/L_{c}\sim 0.1$, for $\sigma^{2}=0.1$ has $\lambda_{\parallel}\simeq 3.3\times 10^{9}$ cm $\simeq 0.5R_{*}$ that increases outward as shown in Sec. 5.2). We note that for the case of weak turbulence ($\sigma^{2}=0.01$, Figs. 6 and 7, upper row) the depleted regions seen at higher $\sigma^{2}$ are not visible on the $R_{p}$-sphere as on the spheres at $R_{p}=R_{b},R_{e}$ the points intercepted by open field lines are approximately uniform and closed lines do not reach such distances.

As for the escaping EPs, once they are channelled into the fast wind region, the large $B_{0}$ (see Fig. 10, lower row) keeps them confined and focussed toward the caps, where $B_{0}$ is larger and hence $r_{g}$ smaller.

Particularly relevant to the influence of EPs on planets in our simulated magnetic field configuration is the approximate symmetry of the caps (see Sect. 6) with respect to the equatorial plane ($\theta^{\prime}=90^{\circ}$); such a pattern results within the fast wind region from the approximately symmetric and greater $B_{0}$ (lower row in Fig. 10) that reduces $r_{g}$ thus favouring the confinement and focussing EPs within the caps.

In the case of a Sun-like $B_{0}$-field, i.e., approximate alignment of ${\bf B}_{0}$ with the rotation axis, with $\sigma^{2}\simeq 1$ (within the solar system typically $\sigma^{2}<0.1$), EPs would be directed preferentially into the polar regions, leaving planets relatively unaffected. The solar wind latitudinal dependence of EPs in large events is, however, poorly constrained due to the limited number of events with high-latitude in-situ measurements (see Sect. 7).

Surprisingly we find that EPs are focussed toward the equatorial plane even when injected at high latitude, i.e., close to the pole. Such an effect is shown in Fig. 11 where EPs are injected, with isotropic velocity distribution, in the latitudinal ring in the upper hemisphere close to the geographic north pole with $\theta^{\prime}=160-170^{\circ}$. In this case, EPs are focused on the $R_{p}$-sphere within $40^{\circ}$ from the equatorial plane mostly in the upper hemisphere, except for a few points in the lower hemisphere ($180^{\circ}<\phi<230^{\circ}$) due to an additional diffusion in the azimuthal direction.

We note that, despite the reduced filling factor of the EP caps for greater values of $\sigma^{2}$ shown in Figs. 6 and 7, that would seem to suggest a smaller $N_{R_{p}}$, the ratio $N_{R_{p}}/N_{\rm inj}$ actually increases for greater $\sigma^{2}$ (see Fig. 9). This effect results again from (1) a more efficient perpendicular diffusion at the boundary between open and closed field lines and from (2) the increase of $\lambda_{\parallel}$ with distance from the star ($\lambda_{\parallel}\propto B_{0}(r)^{-1/3}$). For most EPs injected on open field lines near the boundary, the former enhances the frequency of decorrelation from a given field line, as discussed above, and the latter favours EPs moving outward with an increasing $\lambda_{\parallel}$ rather than back to the star. Such combined effects ultimately prevent most particles from collapsing to the star and allow them to propagate outward toward the equatorial caps.

At larger EP energy, the escape of EPs injected at the open/closed field line boundary is favoured, as suggested by Fig. 9, right panel: $10$ GeV protons arrive more copiously on the $R_{p}$-spheres than $1$ GeV ones. This is a result of a larger perpendicular transport coefficient at larger energy, regardless of the particular model.

Finally, the features in the bottom rows of Fig. 7 protruding out of the caps toward greater $\phi$, and also present to a lesser extent in Fig. 6, map the stripe at constant latitude of maximal wind flow visible in red in Fig. 10, lower panels. On the other hand, the EP caps are shifted to smaller $\phi$ as a result of the stellar rotation.

6 Energetic Particle Flux within the TRAPPIST-1 System

The total output of EPs from M dwarf stars cannot be measured directly at present. A possible approach to estimate the EP abundance relies on the solar correlations between the observed properties of coronal flares and in-situ spacecraft measurements of EP fluxes at $1$ AU. GOES observations of $800$ SXR solar flares ($1.5-12.4$ keV) at the Sun and measurements of the associated $>10$ MeV energetic protons events have shown an approximately linear correlation of the far-UV emission line flux to the proton flux (Belov et al., 2007).

Youngblood et al. (2017) found two correlations: (1) between SXR peak flux and the flux of $>10$ MeV protons from GOES data only; (2) between SDO/EVE He II 304 $\AA$ emission line fluence during the entire durations of flares and $>10$ MeV GOES protons fluence. By using a sample of stellar flares observed by the Hubble Space Telescope (HST) and Chandra/ACIS, Youngblood et al. (2017) finally inferred the proton enhancement for other stars. The He II 304 $\AA$ ($41$ eV) flare fluence was related to the HST far-UV ($7.3-13.6$ eV) fluence with the M dwarf synthetic spectrum created in Fontenla et al. (2016). The solar flaring rates for M- and X-class (corresponding to a SXR peak flare flux of $10^{-5}$ and $10^{-4}$ W/m2 at $1$ AU in the $[1-8]\,\textup{\AA}$ band in the GOES classification, respectively) are estimated to be $0.02$ hr*-1* and $2.3\times 10^{-5}$ hr*-1*, respectively, based on flare observations in the period 1976-2000 (Veronig et al., 2002). Therefore, the estimated rates for M- and X-class flares on the M4 dwarf GJ 876 are $\sim 0.4$ hr*-1* (Youngblood et al., 2017), $20$ and $1.7\times 10^{5}$ times more frequent than the Sun for M- and X-class, respectively. The rescaling to the average HZ radius $r_{876}^{HZ}\sim 0.18$ AU (Youngblood et al., 2017, via the empirical scaling in Kopparapu et al. (2014)), leads to an increase of the flux by a factor $30$ for the HZ of GJ 876 (a flaring rate $600$ and $5\times 10^{5}$ times higher for M- and X-class, respectively); it should also be noted that, due to the closer HZ, M-class flares are scaled up to X10. Therefore, Youngblood et al. (2017) estimate that large GJ 876 flares (SXR peak flux $\geq 10^{-3}$ W m*-2*) lead to a $>10$ MeV proton flux ($F_{876}^{max}$) on HZ planets up to $10^{3}$ protons cm*-2* s*-1* sr*-1*, and enhanced up to $\sim 4$ orders of magnitude higher than for the present-day Earth by both the higher flaring rate and closer distance.

Since the Youngblood et al. (2017) scaling applies to EPs of any energy $>10$ MeV, it should be noted that here we implicitly assume a uniform EP energy spectrum, although different spectral shapes, e.g., power-law or log-parabola, normalized to $>10$ MeV could be used.

The TRAPPIST-1 HZ is dramatically closer to the host star ($R_{e}=0.029$ AU) than the GJ 876 HZ, leading to a much higher EP flux. Rescaling the flux from $r_{876}^{HZ}=0.18$ AU to the injection radius in our simulations, $R_{s}=10R_{*}=0.0056$ AU, we find an EP flux enhancement

[TABLE]

The relation above holds for very intense flares.

By using the maximal EP flux in Eq. 5, we can determine the flux $F(R_{p})$ of EPs impinging on the planet 1e along its 6 day orbital motion around the star. The EP flux impinging on a ring of the $R_{p}$-sphere with semi-aperture $\Delta\theta^{\prime}=5^{\circ}$ centered on the equatorial plane is given by

[TABLE]

where $N^{\prime}_{R_{p}}$ is the number of EPs hitting the ring and we have used $A=\int_{85^{\circ}}^{95^{\circ}}{\rm sin}\,\theta^{\prime}d\theta^{\prime}=0.17$.

The flux of $10$ GeV EPs with $\sigma^{2}=1$, $R_{s}=10R_{*}$ along the orbit of planet $1$e is shown in Fig. 12. The maximal flux, $\sim 1.2\times 10^{5}\,\frac{{\rm protons}}{{\rm cm}^{2}\,{\rm s}\,{\rm sterad}}$, exceeds by roughly $6$ orders of magnitude the EP abundance at the present-day Earth. However, such an estimate is subject to several caveats, which we discuss in the following section.

7 Discussion

The results described in Sect. 5 show that the magnetic fluctuations not only affect the small-scale particle motion but change drastically the behaviour of EPs over the entire inner astrosphere.

7.1 The spatial distribution of propagating EPs

The EP-depleted angular regions on the $R_{p}$-sphere track the slow wind populated by closed field lines that lead to EPs being trapped and lost due to their trajectories leading back to the stellar surface. For relatively large values of $\sigma^{2}$, particles are lost due to enhanced perpendicular diffusion into the closed field region (see Fig.s 6 and 7). The opening of the closed field lines further out results in the narrowing of the depleted regions for larger particle injection radii $R_{s}=10R_{*}$ as compared to $R_{s}=5R_{*}$ (see Fig. 8).

The stronger unperturbed magnetic field in the fast wind region on the equatorial plane (see Fig. 10, lower row) favours EP focussing. The EP caps are centered in the region of fast wind speed at $\sim 800-1,000$ ($\sim 950-1,100$) km$/$s at the planet 1b (1e).

A key characteristic of the GJ 3622 proxy magnetogram we adopted for TRAPPIST-1 is its resemblance to a tilted dipole. This gives rise to the focus of EPs at low latitudes, and into the planetary orbital plane. The location of the spherical caps of EPs hitting the $R_{p}$-sphere has potentially important consequences for the energetic particle flux experienced by the planets in our TRAPPIST-1-like system (the TRAPPIST-1 planets themselves are all in coplanar orbits to within $30$ arcmin). It should be noted that the locations of the EP caps are subject to shift both along the orbital plane due to differences in the stellar rotation and planetary orbital periods, and in latitude due to the evolution and probable cyclic behavior of the stellar surface magnetic field. Both times scales associated with these processes are much greater than the EP propagation time scale. We investigate the EP flux variation planets could experience below.

We also point out that the EP focussing onto planets seen in our simulations is not expected to occur in a stellar wind driven by a dipolar magnetic field closely aligned with the stellar rotation axis (such as the solar wind), where the wind is fast at high latitudes (see Fig. 10, lower row). Moreover, $\sigma^{2}$ might attain values greater than $0.1$ only in transients, such as CME-driven shocks, or corotating interaction regions. In-situ solar wind measurements following large solar flares ($>10^{30}$ erg) do not strongly constrain the latitudinal dependence in EP intensity: for instance, in the Bastille day event (Zhang et al., 2003) Ulysses high heliolatitude EP intensity, in the fast wind, was measured at $3.2$ AU distance from the Sun whereas lower latitude intensity, in the slow wind, was measured at a different distance ($1$ AU).

The spatial distribution of EPs centered on the equatorial plane might raise the question of a possible relation with the spatial distribution of CMEs in active M dwarfs found in numerical simulations by Kay et al. (2016): regardless of the latitude of injection, CMEs are deflected further out ($\sim 60R_{*}$) along the near-equatorial current sheet, where the $B$-field is minimum and therefore CME expansion encounters the minimal magnetic confinement as the ratio between the CME ram pressure to the stellar magnetic pressure is highest. In our simulations, EPs are unleashed from the bulk motion of CME-driven shocks at the initial time, so their motion is independent of the subsequent CME trajectory. We expect that in a stellar wind with a highly tilted magnetic-to-rotation axis, such as the one in Fig. 1, particles emitted at $R>5R_{*}$ by CMEs along the current sheet (blue-purple stripe in Fig. 10) will be transported toward the fast wind region for $\sigma^{2}>0.1$ (cfr. Fig. 11). However, for $\sigma^{2}<0.1$ we expect that the fewer escaping EPs will concentrate along the current-sheet stripe.

7.2 On the absolute EP flux and trapping of EPs and CMEs

Since EPs can be trapped by close field line regions, they can also be liberated from these regions when the closed field is perturbed or broken open. Such a disruption to the stellar magnetic $B_{0}$-structure can result from a CME-driven shock (not accounted for in our static solution MHD simulations), increasing the chances for EPs to fill the depleted regions on the $R_{p}$-sphere.

On the other hand, EPs accelerated and injected directly by coronal flares at $R_{s}<2R_{*}$, rather than by the travelling shock scenario considered in Figs. 6, 7, 8, are efficiently trapped by the very intense stellar magnetic field and by the closed field lines. Figure 9, left panel, shows that doubling $R_{s}$ approximately doubles $N_{R_{p}}$. The low $N_{R_{p}}/N_{\rm inj}$ ($3.0-3.7\%$) for $R_{s}=1.5R_{*}$ described in Sect. 5, might be considered a lower limit if disturbances of the $B_{0}$ topology by flares or CMEs can enable a larger $N_{R_{p}}/N_{\rm inj}$.

These results indicate that a fairly simple dipole-like magnetic field structure on a magnetically active star prevents coronal flares from contributing significantly to the steady abundance of EPs further out. Thus, at face value in the undisrupted magnetic topology used here, CME-driven shocks might be expected to be the dominant supplier of EPs within the interplanetary medium of a very active star.

In this context, the underlying assumption that CMEs can successfully escape the strong magnetic confinement of the stellar magnetic field to drive shock waves that accelerate EPs is uncertain and needs further investigation. Drake et al. (2016) presented a preliminary simulation of what would have been a large CME on the Sun induced on the surface of the very active K dwarf AB Dor, and found the event to be entirely contained by the strong overlying magnetic field. Indication that a $75$ G dipolar field prevents the escape from the stellar corona of CMEs with kinetic energy $<10^{32}$ erg has also been found by Alvarado-Gómez et al. (2018) based on a number of detailed numerical CME simulations.

There are thus two potentially powerful mechanisms that could strongly limit EP fluxes from active stars: EPs from flares are contained; and CMEs that might generate EPs at larger distances also fail to escape.

The morphology of $N_{R_{p}}/N_{\rm inj}$ in Figs. 6, 7 are, to a good approximation, independent of the EPs energy. In addition, the Youngblood et al. (2017) correlation is determined for $>10$ MeV protons, with an unspecified EP energy-dependence. Regardless of the specific shape, we expect EPs energy spectrum to decrease at larger energy; thus, the EP flux $\sim 10^{5}\,\frac{{\rm protons}}{{\rm cm}^{2}\,{\rm s}\,{\rm sterad}}$ impinging on 1e (see Sect. 6 and Fig. 12) will be lower at $\gg 10$ MeV. We will investigate this effect in a forthcoming work.

We emphasize that our estimated number of injected EPs (Sect 6) is based on strong flares in SXR observed from GJ 876 and classified as large, i.e., time-integrated SXR flux larger than $10^{29}-10^{30}$ erg, due to the small distance to the star. The extrapolation of the correlation between SXR and EP fluence to such large events is uncertain due to the scatter of the observations and to the fact that no solar events beyond a certain energy have been observed ($>$ X10, Hudson, 2007; Drake et al., 2016). However, Kepler-2 constraints (Vida et al., 2017) on TRAPPIST-1 white light flares lead to an estimated total flare energy (in the optical) between $10^{31}$ and $10^{33}$ erg, similar to other very active M dwarfs (Hawley et al., 2014) and beyond the total estimated energy of the Carrington event ($10^{32}$ erg, Carrington, 1859) that is among the most energetic geomagnetic storms ever recorded on Earth. Thus, we argue that the dramatic EP enhancement in the HZ of M dwarfs like TRAPPIST-1 or GJ 876, as compared to present-day Earth, might be not uncommon. Such EP fluxes could have a significant impact on exoplanet atmospheric ionization.

We do not consider the spatial distribution of the EP hitting points on the planetary surface or through the planetary atmosphere, since they depend strongly on the propagation through the planetary magnetosphere and atmosphere: the magnetospheric properties of the TRAPPIST-1 HZ planets—or any other exoplanets—are at present unknown. The effect of EPs on the atmospheric evolution also depends on the atmospheric mass and chemical composition, which are also unknown for TRAPPIST-1. Lyman $\alpha$ detection of variability during transits (observed for planets 1b and 1c, but not 1e, Bourrier et al. 2017), could be useful for further atmospheric characterization, although more detailed constraints will likely have to await observations by next generation facilities.

By using preliminary 3D-MHD simulations here, we instead consider simply the geometrical flux impinging onto a latitudinal ring, centered on the equatorial plane. We have integrated fluxes over a $5^{\circ}$ semi-aperture, which is much broader than the dispersion of the planetary orbits, in order to obtain sufficient signal from our test particle results (see Fig. 12).

8 Conclusions

We have carried out numerical test-particle simulations to calculate for the first time the propagation of stellar energetic particles through a realistic and turbulent magnetic field of an M dwarf star and its wind. Our simulations have been tailored to a proxy for TRAPPIST-1, and we have investigated the flux of energetic particles throughout the habitable zone of the TRAPPIST-1 system to the outermost planet. Particle acceleration by flares close to the stellar surface and further out by CME-driven shocks is mimicked here by injecting particles at various distances from the star over the full sphere and with an isotropic velocity distribution. We highlight three important aspects of the results.

Particles injected close to the stellar surface, regardless of their energy, are trapped within the strong stellar magnetic field. In our simulations, only a 3–4% of particles injected within half a stellar radius from the surface escape. The escaping fraction increases strongly with increasing injection radius: Particles accelerated further from the stellar surface have a much greater chance of escaping the closed stellar magnetic field.

Particles are increasingly focussed and directed toward the equator and toward open field fast wind regions with increasing turbulence amplitude. This results from asymmetric perpendicular diffusion from stronger to weaker field regions. In our TRAPPIST-1 proxy, strong turbulence produces two concentrated polar streams $180^{\circ}$ apart of energetic particles in the fast wind region focussed on the planetary orbital plane, regardless of the angular location of the injection. Based on the scaling relation between far-UV emission and energetic protons for solar flares by Youngblood et al. (2017), we estimate that the innermost putative habitable planet, TRAPPIST-1e, is bombarded by a proton flux up to 6 orders of magnitude larger than experienced by the present-day Earth. Such a bombardement of planets in this study is found to result largely from the misalignment of the $B$-field/rotation axis assumed for the star-proxy. Since the exact magnetic morphology and alignment of the magnetic field is currently unknown for TRAPPIST-1, and for M dwarfs in general, our results indicate that determination of these quantities for exoplanet hosts would be of considerable value for understanding their radiation environments.

The trapping of EPs produced close to the stellar surface suggests that particles directly accelerated in flares do not generally escape, and that the ambient energetic particle environment of planets is dominated by particles accelerated in CME shocks. However, recent findings that CMEs can be strongly suppressed by strong stellar magnetic fields (Drake et al., 2016; Alvarado-Gómez et al., 2018) point to a consequent large uncertainty in our understanding of the EP fluxes that exoplanets around active stars sustain.

We thank the referee whose constructive comments helped to improve significantly the clarity of the presentation. We are very grateful to Drs G. Ballester, K. France, M. Lingam and A. Loeb for comments. The work of FF was supported, in part, by Scholarly Studies Award 40488100HH00181 at the Harvard-Smithsonian Center for Astrophysics, by NASA under Grants NNX15AJ71G and 80NSSC18K1213, by NASA through Chandra Award Number $TM6-17001B$ 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 and by NSF under grant 1850774. JJD was funded by NASA contract NAS8-03060 to the CXC and thanks the Director, Belinda Wilkes, and the CXC science team for continuing advice and support. JDAG was supported by Chandra grants AR4-15000X and GO5-16021X. Resources supporting this work were partially provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center.