INFLUENCE OF STELLAR FLARES ON THE CHEMICAL COMPOSITION OF EXOPLANETS AND SPECTRA

Olivia Venot; Marco Rocchetto; Shaun Carl; Aysha Roshni Hashim; Leen Decin

doi:10.3847/0004-637X/830/2/77

1. INTRODUCTION

M stars are the most common stars in our Galaxy, and thus are likely to harbor most of the planetary systems. A particularity of M stars is their high activity. This most active class of stars undergoes great stellar variability, such as star spots, granulation, and flares. The latter are violent and unpredictable outbursts of the photosphere that are caused by sudden changes in magnetic flux profiles. During such events, the emitted radiative energy can increase by several orders of magnitude over the complete wavelength range, including prominent H_α emissions.

A few studies on the effect of flares on exoplanets have been conducted to date. Segura et al. (2010) focused on the composition of Earth-like planets. They simulated an Earth-like planet located in the habitable zone of the active M dwarf AD Leo and tracked the evolution of the abundances of water and ozone in the planetary atmosphere during a flare event. They found that the enhanced UV radiation emitted during a flare would not affect the habitability of the planet. They did not address the effect of the flare on the planetary spectra. Complementary to this latter paper, Tofflemire et al. (2012) studied the consequences of M dwarf stellar flares on the observation of exoplanets. They found that, if a flare with an energy of 10³¹–10³² erg—i.e., a standard flare occurring in M-type dwarfs—occurs during a transit, it would have no impact on the infrared (IR) detection or characterization of exoplanets above a level of ∼5–11 mmag. The most energetic flares (E > 10³⁴ erg) would be easily detected in IR observations, as they would produce an increase in magnitude significantly greater than the photometric noise obtained with most observational instruments. Here, again, the authors did not consider the effect of flares on the atmospheric composition of exoplanets. Recently, Rugheimer et al. (2015) studied the atmospheres and spectra of Earth-like planets orbiting different kinds of M dwarfs (both active and inactive) leading to a very useful spectral database. They did not, however, vary the stellar flux to mimic stellar activity during their simulations.

Far from astrobiological consideration, in this paper, we simulate two exoplanets (sub-Neptune/super-Earth) orbiting the star AD Leo, and study the consequences of a stellar flare on both the chemical composition of the atmospheres and the synthetic spectra. We model two cases: a warm and a hot planet, with effective temperatures of 412 K and 1303 K, respectively. In the context of future space missions (James Webb Space Telescope (JWST), PLAnetary Transits and Oscillations of stars (PLATO), etc.), the aim is to determine whether or not planetary spectra acquired during a stellar flare can be differentiated from "quiescent" spectra, and whether or not false conclusions can be drawn, i.e., we aim to study the extent to which a stellar flare can influence the chemical composition of the atmosphere and the resulting synthetic spectra. In Section 2, we describe the one-dimensional (1D) kinetic model and the radiative transfer model, along with the adopted thermal profiles and stellar fluxes. We present our results in Section 3, and discuss them in Section 4. Our conclusions are given in Section 5.

2. MODELS

2.1. 1D Neutral Chemical Model

We used our thermo-photochemical model, developed especially for warm atmospheres (Venot et al. 2012). This code has been used to study exoplanets (Venot et al. 2012, 2013, 2014, 2015; Agúndez et al. 2014), as well as the deep atmosphere of giant planets in our solar system (Cavalié et al. 2014; Mousis et al. 2014). We use the numerical solver DLSODES.⁴ It can be found in ODEPACK, a collection of powerful numerical solvers for ordinary differential equation systems. To resolve stiff problems (i.e., those containing processes having a great disparity in rates) as systems of time-dependent chemical rate equations, DLSODES uses the backward differentiation formula (BDF) method. This multiple-step method was implemented for the first time by C. W. Gear (Gear 1971), and is known to yield very accurate solutions in such cases. The superiority of LSODES over other ODE solvers, in resolving astrochemical kinetics systems, has been shown by Nejad (2005). Also, the advantage of using DLSODES instead of the Cranck–Nicholson method (a common method used in photochemical models) is presented in Dobrijevic et al. (2010). More details about this solver can be found in the work of Hindmarsh (1983) and Radhakrishnan & Hindmarsh (1993).

In our study, the time step for the solver starts at 10⁻¹⁰ s and increases continuously as chemistry evolves.

One of the strengths of our model is its robust chemical scheme, which has been validated against experiments over a large range of pressures (0.01–100 bar) and temperatures (300–2500 K). This reaction network has been constructed thanks to a close collaboration with specialists in combustion (see details in Venot et al. 2012 and references therein). This C₀–C₂ scheme describes the kinetics of species with up to two carbon atoms. The network contains 105 neutral species involved in 960 reversible reactions. To model the effect of UV irradiation, we added to this network a set of 55 photodissociations. References to the associated cross-sections and quantum yields used can be found in Venot et al. (2012) and Dobrijevic et al. (2014). For the cross-sections of CO₂ (Venot et al. 2013) and NH₃ (O. Venot et al. 2016, in preparation), we used our recent measurements at 500 K.

We recently developed an extended chemical scheme (C₀–C₆) that is able to reproduce the kinetics of species with up to six carbon atoms (Venot et al. 2015). This chemical network is very useful for studying atmospheres rich in carbon atoms. As we have shown in Venot et al. (2015), the C₀–C₂ scheme is sufficient for atmospheres with solar C/O ratios. It yields the same results and is significantly more efficient to compute. Here, we have not not explored C/O ratios different than solar, so our C₀–C₂ scheme is sufficient. In order to parameterize the vertical mixing, we choose a constant eddy diffusion coefficient K_zz = 10₈ cm² s⁻¹, which corresponds to the lower value employed by Fortney et al. (2013) (from which we take our pressure–temperature profiles, see Section 2.3) and which is a commonly used value in studies of exoplanetary atmospheres (e.g., Lewis et al. 2010; Line et al. 2011; Moses et al. 2011; Venot et al. 2013, 2015). The planets have a radius R_p of 0.238 R_Jup (2.61 R_Earth) and a mass M of 0.02 M_Jup (6.35 M_Earth).

2.2. Stellar Fluxes

AD Leo (also called Gl 388) is a M3.5V star (Shkolnik et al. 2009) located at 4.9 pc from the Sun (Rojas-Ayala et al. 2012). Its luminosity is 0.024 L_☉ (Pettersen & Coleman 1981) and its effective temperature is 3390 ± 19 K (Rojas-Ayala et al. 2012). It is one of the most active flare stars known. The great flare affecting AD Leo in 1985 was observed by Hawley & Pettersen (1991) in the visible and UV spectral region. Unfortunately, the time distribution of this acquisition was not well-distributed (especially in the UV spectral region). From these original data, Segura et al. (2010) constructed a time-dependent sequence of spectra that can be used in photochemical models. A complete description of the method used can be found in their paper.

To model the flare event, we used the data between 100 nm and 444 nm. For wavelengths less than 100 nm, we use the mean of the Sun spectra at maximum and minimum activity (Gueymard 2004), multiplied by a factor of 100 in order to reproduce the high X-ray luminosity of AD Leo (L_X ≃ 4 × 10²⁸ erg s⁻¹ (Favata et al. 2000) versus L_X☉ ≃ 10²⁶–10²⁷ erg s⁻¹ for the Sun (Judge et al. 2003)). For wavelengths larger than 444 nm, we used the quiescent spectrum of AD Leo from Segura et al. (2005). To calculate the steady states before and after the flare event, we used the same spectrum—except that the quiescent spectrum of AD Leo (Segura et al. 2005) was used down to 100 nm. We note that, for AD Leo, log[L_X/L_bol] = −3.3, which is close to the saturation level (∼−3) that very active stars reach (e.g., Vilhu & Walter 1987; Stauffer et al. 1994). The quiescent and flares stellar fluxes are shown in Figure 1. They are separated into three phases: the first impulsive phase from 0 to 800 s, during which the stellar flux increases rapidly by one order of magnitude; the second impulsive phase, from 912 to 1497 s, during which the stellar flux, after reaching the maximum at 912 s, decreases by less than one order of magnitude over approximately the same timescale as the first phase; and finally, the gradual phase during which the stellar flux continues to decrease, but more slowly.

**Figure 1.** Stellar spectra of AD Leo at a distance of 1 AU from the star. In the [1–99] nm and [445–900] nm ranges, all spectra are identical. Spectra at different time steps during the flare are represented by different colors, as labeled on the plots, and are used as input for the photochemical model at the corresponding time step. In each figure, the quiescent spectra is represented in black. The three phases of the flare are represented: the first impulsive (top), the second impulsive (middle), and the gradual (bottom).
Download figure:
Standard image High-resolution image

2.3. Thermal Profile

To model our hypothetical planet's atmosphere, we use the thermal profiles from Fortney et al. (2013). Their fully non-gray atmosphere code has been adapted to exoplanetary atmospheres (e.g., Fortney et al. 2005, 2008; Morley et al. 2013). The code considers both the incident radiation from the star and the thermal radiation from the planet's atmosphere, which are resolved using the algorithm developed by Toon et al. (1989). Opacities are tabulated using the correlated-k approximation. Using this code, they computed several temperature profiles for GJ1214b-like planets, corresponding to different intensities of irradiation of the planet. All atmospheric models consider an enrichment of 50 times solar metallicity, an internal temperature of 60 K (and not 75 K as it is misprinted in the legend of their Figure 3, J. Fortney 2016, private communication), and a Bond albedo of ∼0.05. Here, the albedo is very low due to the absence of clouds in the model. However, clouds, which have been detected in GJ1214b (Kreidberg et al. 2014), can significantly increase albedo (up to 0.4–0.6 according to Charnay et al. 2015). More details about the data used to calculate these profiles can be found in Fortney et al. (2013), and references therein.

We consider that using thermal profiles of GJ1214b as analogs to model our hypothetical "AD Leo b" planets is appropriate, given that AD Leo and GJ 1214 are both M stars, with similar effective temperatures (respectively, 3 390 ± 19 K and 3 245 ± 31 K, Rojas-Ayala et al. 2012). We used the profiles corresponding to the effective temperatures of 412 and 1303 K (see Figure 2). These profiles originally end at 10⁻³ mbar, with a temperature increasing at low pressures, which is not a physical effect but rather an artifact due to boundary conditions of the model. In order to model completely the photochemical processes, we extrapolated the profiles to lower pressures, using an isothermal profile from 10⁻² to 10⁻⁴ mbar to avoid this artifact. The distances between the star and the planets necessary to be consistent with these thermal profiles are 0.069 AU and 0.0069 AU, for 412 K and 1303 K, respectively. The hottest profile could raise doubt about the possible existence of a sub-Neptune/super-Earth so close to its star. Planets observed at a distance less than 0.01 AU are usually smaller than our super-Earth models (i.e., R ∼ 0.07 R_J, Charpinet et al. 2011; Muirhead et al. 2012). However, Sanchis-Ojeda et al. 2015 reported very recently the discovery of an exoplanet with similar properties to our model (R = 0.22 R_J at a distance of 0.008 AU) and they suggested that this planet was a disintegrating rocky planet with a cometary head and tail. This kind of object can be seen as a final step in the evolution of a planet being destroyed by tidal forces (S. Raymond 2016, private communication).

**Figure 2.** Thermal profiles of the hypothetical AD Leo b planets.
Download figure:
Standard image High-resolution image

2.4. Spectra Model

Synthetic transmission spectra in the wavelength range 0.8–20 μm for the different atmospheric compositions were obtained using the forward model described in Waldmann et al. (2015), based on the Tau code by Hollis et al. (2013) and modified to allow variable pressure-dependent temperature profiles. The 1D radiative transfer model is based on a code that calculates the optical path through the planetary atmosphere, resulting in a transmission spectrum of transit depth as a function of wavelength. The wavelength and temperature-dependent absorption cross-sections for the absorbing molecules were computed using the line lists from ExoMol (Tennyson & Yurchenko 2012), HITRAN (Rothman et al. 2009, 2013), and HITEMP (Rothman et al. 2010). We only considered the following molecules: H₂O, CO₂, OH, HCN, NH₃, CH₄, and CO. Although these molecules represent only a fraction of the 105 molecules considered in the thermo-photochemical model, they are the most abundant; therefore, they dominate the spectral features. We also included collision-induced absorption cross-sections for H₂–H₂ and H₂–He (Borysow et al. 2001; Borysow 2002). The atmosphere is assumed to be cloud-free.

3. RESULTS

3.1. Chemical Composition

Although many species are impacted by the increase of UV irradiation due to the flare, we describe in the following the evolution of the mixing ratios of only six of them: hydrogen (H), the amidogen radical (NH₂), ammonia (NH₃), carbon dioxide (CO₂), nitric oxide (NO), and the hydroxyl radical (OH). These species are among the most abundant, and undergo deeper and more significant changes than the other chemical species in our model. Moreover, three of them are considered in our radiative transfer code to calculate synthetic spectra. Their abundances during the different phases of the flare are shown on Figures 3, 4, 6, and 7. For both thermal profiles, most species see their mixing ratios globally increase during the flare event, with a maximum abundance around 912 s (corresponding to the peak of the stellar flux). Exceptions are NH₃, CO₂, and NO, which can have a lower abundance during the stellar flare than at steady state. The other species considered in this calculation, but not plotted (H₂O, HCN, and CH₄), experienced changes only in the very upper atmosphere (P < 10⁻¹ mbar and P < 4 × 10⁻⁴ mbar for the 1303 K and 412 K profiles, respectively). Their abundances globally decrease by ∼1–5 orders of magnitude. The abundances after steady state, during the "return to quiescence" phase are shown in Figures 5 and 8. Figure 9 shows a comparison between the initial and final steady states.

**Figure 3.** Evolution of H, NH₂, and NH₃ mixing ratios during the flare event with the thermal profile corresponding to T_eq = 412 K. The legend for all figures is in the upper panels.
Download figure:
Standard image High-resolution image

**Figure 4.** Evolution of CO₂, NO, and OH mixing ratios during the flare event with the thermal profile corresponding to T_eq = 412 K. The legend for all figures is in the upper panels.
Download figure:
Standard image High-resolution image

**Figure 5.** Evolution of H, NH₂, NH₃, CO₂, NO, and OH mixing ratios after the flare event with thermal profile corresponding to T_eq = 412 K. The legend for all figures is in the upper panels.
Download figure:
Standard image High-resolution image

**Figure 6.** Evolution of H, NH₂, and NH₃ mixing ratios during the flare event with thermal profile corresponding to T_eq = 1303 K. The legend for all figures is in the upper panels.
Download figure:
Standard image High-resolution image

**Figure 7.** Evolution of CO₂, NO, and OH mixing ratios during the flare event with thermal profile corresponding to T_eq = 1303 K. The legend for all figures is in the upper panels.
Download figure:
Standard image High-resolution image

In several instances, the most important reactions controlling changes in species concentrations were identified. This was achieved using a module of our chemical pathway analysis software suite, which essentially keeps a tally of the cumulative reaction flux for each reaction in the scheme, such that the total and net reaction flux for any species over any time period can be identified. The net reaction fluxes were checked for equality against net changes in concentration.

3.1.1. Thermal Profile T_eq = 412 K

Hydrogen: The abundance of this species changes between 3 × 10⁻⁴ and 7 × 10² mbar. At higher or lower pressures, the mixing ratio of hydrogen (y_H) remains identical to that at the steady state. Globally, during the impulsive phase (0–912 s), the amount of H increases with time (by a factor 25 at 1 mbar). We note that, between 3 × 10⁻³ and 3 × 10⁻² mbar, the evolution of the mixing ratio does not evolve linearly with time—one can notice some overlap and crossing between the different temporal profiles. This is caused by the peak of destruction between these two pressures, which moves downward to higher pressures (lower altitudes). The reactions responsible for the increase of y_H depend on the pressure level. For pressures lower than 10⁻² mbar, H₂O + hν $\to$ H + OH and OH + H₂ $\to$ H + H₂O are the main production pathways. Around 10⁻² mbar, the photodissociation NH₃ + hν $\to$ NH₂ + H is the principal source of hydrogen, with a contribution of the two aforementioned reactions. Then, around 10 mbar, the main reactions producing hydrogen are NH₂ + H₂ $\to$ NH₃ + H and NH₃ + hν $\to$ NH₂ + H, which form a limited catalytic cycle. Between 912 and 2586 s, we see a decrease of the abundance of H, toward the initial steady state value. This is due to the decrease in efficiency of production routes that can no longer compete with the destruction pathways: N₂H₂ + H (+M) $\to$ N₂H₃ (+M) and N₂H₃ + H $\to$ N₂H₂ + H₂ for P ∼ 10⁻² mbar, and CH₃ + H (+M) $\to$ CH₄ (+M) and CH₄ + H $\to$ CH₃ + H₂ for P > 0.1 mbar.

Amidogen: This species experiences abundance variations from the top of the atmosphere, down to 6 × 10² mbar. The mixing ratio of NH₂ increases during the impulsive phase and decreases after 912 s. Its time evolution is similar to that of hydrogen because this species is mainly produced by one of the reactions responsible for the increase of H; namely, NH₃ + hν $\to$ NH₂ + H. We also note the important contribution of the reaction N₂H₃ + H $\to$ NH₂ + NH₂ in the production of amidogen. The reactions that destroy NH₂ are NH₂ + H (+M) $\to$ NH₃ (+M), NH₂ + NH₂ (+M) $\to$ N₂H₄ (+M), and NH₂ + H₂ $\to$ NH₃ + H. Their relative efficiencies depend on the pressure level and phase of the flare (their contributions become more important after 912 s). We note that, interestingly, at very low pressure (∼4 × 10⁻⁴ mbar), the decrease of NH₂ after 912 s is due to the reaction NH₂ + OH $\to$ NH₂OH that becomes dominant because of the increase of y_OH.

Ammonia: Unlike the other species considered here, ammonia globally decreases with time during the first 1038 s, down to 3 × 10⁻² mbar. At 5 × 10⁻³ mbar, ${y}_{{\mathrm{NH}}_{3}}$ decreases by a factor of 3. At lower pressures (P ∼ 1.5 × 10⁻⁴ mbar), the abundance of ammonia decreases by a factor of ∼20. The loss of NH₃ is due mainly to photodissociation; NH₃ + hν $\to$ NH₂ + H. At some levels, we also note the contribution of other reactions: NH₃ + OH $\to$ NH₂ + H₂O at very low pressure (∼4 × 10⁻⁴ mbar), and NH₃ (+M) $\to$ NH₂ + H (+M) around 3 × 10⁻⁴ mbar. From 1130 s, the evolution of NH₃ is more moderate. Around 4 × 10⁻⁴ mbar, and for 8 × 10⁻⁴ < P < 9 × 10⁻³ mbar, ${y}_{{\mathrm{NH}}_{3}}$ increases together with time. This increase (which is most visible around 3 × 10⁻³ mbar) is due to the reaction NH₂ + H (+M) $\to$ NH₃ (+M). Apart from these pressures, the change in mixing ratio of ammonia is relatively small: we only observe a slight decrease between 3 × 10⁻⁴ and 2 × 10⁻⁴ mbar.

Carbon dioxide: CO₂ is affected by the stellar flare only in the very upper atmosphere, down to 6 × 10⁻⁴ mbar. In the first part of the flare, we observe two different evolutions for the abundance of this species. For pressures lower than 3 × 10⁻⁴ mbar, it decreases with time (a factor ∼15 at ∼1.5 × 10⁻⁴ mbar). Then, for pressures between 3 × 10⁻⁴ and 6 × 10⁻⁴ mbar, CO₂ sees its abundance slightly increasing (by a factor ∼1.3). During the second part of the flare, the abundance of CO₂ hardly changes. We observe a very small decrease of ${y}_{{\mathrm{CO}}_{2}}$ for P < 3.5 × 10⁻⁴ mbar. The loss of CO₂ is largely due to the photodissociation of CO₂. The branching to O(¹D) is dominant throughout, but the relative yield of O(³P) increases at lower pressures (∼4 × 10⁻⁴ mbar). The production of carbon dioxide is due to CO + OH $\to$ CO₂ + H. We will see in Section 3.4 that the variation of abundance of this species is quite small, but sufficient to affect the synthetic spectra.

Nitric oxide: The abundance of this species changes from the top of the atmosphere down to 0.5 mbar. We observe different regimes depending on the pressure. During the first 1038 s, for pressures lower than 2.5 × 10⁻⁴ mbar, y_NO decreases slightly with time, from 0 to 1038 s. The decrease is mainly due to the photodissociation NO + hν $\to$ N(⁴S) + O(³P), which is somewhat tempered by the formation reactions N(⁴S) + OH $\to$ NO + H. We also observe a decrease of the abundance of NO for 4 × 10⁻⁴ < P < 3 × 10⁻³ mbar (up to a factor of 3), and for 2 × 10⁻² < P < 0.5 mbar (also up to a factor of 3), due to the reaction NH₂ + NO $\to$ N₂ + H₂O. Between these two regimes, y_NO increases with time because it is produced by HNO + H $\to$ NO + H₂. In the second part of the flare, the abundance of NO varies less. There is a small decrease of its abundance between 3 × 10⁻³ and 0.5 mbar, and also around 2 × 10⁻⁴ mbar.

Hydroxyl radical: The abundance of OH varies from the top of the atmosphere down to ∼500 mbar. From 0 to 912 s, we globally observe an increase of the abundance of OH (by a factor of 100 at 10⁻³ mbar), except for pressures lower than 3 × 10⁻⁴ mbar. From 912 to 2586 s, we observe a decrease of the abundance of OH toward the initial steady state. This species evolves similarly to hydrogen because the abundance of OH is strongly connected to H through the dominant forming reaction H₂O + hν $\to$ H + OH. The reaction that mainly destroys OH is OH + H₂ $\to$ H + H₂O.

3.1.2. Thermal profile T_eq = 1303 K

Hydrogen: The amount of hydrogen varies from 30 mbar up to 10⁻³ mbar; y_H increases significantly up to 1038 s. The maximum enhancement occurs at 3 mbar (a factor of ∼4.5). The reactions responsible for the increase of hydrogen are different depending on the pressure level. We found that, in the very upper atmosphere (down to 10⁻¹ mbar), the increase of y_H is due to the photodissociation of water (H₂O + hν $\to$ H + OH) and to the reaction OH + H₂ $\to$ H + H₂O. At higher pressures (∼1 mbar), the production of hydrogen is due mainly to NH₂ + H₂ $\to$ NH₃ + H and NH₃ + hν $\to$ NH₂ + H, with a small contribution of CO + OH $\to$ CO₂ + H. Around 3 mbar (where we observe the maximum change in H abundance), the reactions that produce the most hydrogen are OH + H₂ $\to$ H + H₂O, OH + H₂ $\to$ OH + H, and CN + H₂ $\to$ HCN + H.

From 1130 to 2586 s, the variations of the abundances are small. We observe a transition between two regimes around 0.5 mbar. For smaller pressures, y_H keeps increasing slightly (by a factor less than ∼1.25 between 1130 s and 2586 s), but for higher pressures, the amount of hydrogen decreases with time (by a factor of ∼1.5 at 10 mbar).

Amidogen: As for H, the abundance of NH₂ is affected down to 30 mbar. We observe three regimes of variations. For pressures greater than ∼0.2 mbar, the abundance increases with time (up to 912 s). Between 1 × 10⁻³ and 2 × 10⁻³ mbar, there is a production peak of NH₂. Here and above, after a sudden increase of ${y}_{{\mathrm{NH}}_{2}}$ (by a factor of ∼3 at 1 × 10⁻² mbar) from 0 to 100 s, the amount of amidogen decreases with time. Finally, between 2 × 10⁻³ and 0.2 mbar, the temporal evolution of the abundance of NH₂ is more variable—but remains greater than the initial steady state. From 1038 s, ${y}_{{\mathrm{NH}}_{2}}$ diminishes with time at all pressures. We analyzed the chemical pathways involved in the production and destruction of NH₂, and found that, over the entire duration of the flare, the production of amidogen was due to NH₃ + hν $\to$ NH₂ + H, so the production rate of NH₂ follows closely that of ammonia. The dominant removal process of NH₂ is by reaction with H₂: NH₂ + H₂ $\to$ NH₃ + H. However, we note that, around 2 mbar, the main reactions involved in the formation/destruction of amidogen are different than at lower pressures: NH₂ + H $\to$ NH + H₂ for the destruction; HNCO + H $\to$ NH₂ + CO and OH + HCN $\to$ NH₂ + CO for the production.

Ammonia: This species sees its abundance affected by the flare from the top of the atmosphere down to ∼0.2 mbar. As observed for NH₂, ${y}_{{\mathrm{NH}}_{3}}$ undergoes an initial increase (by a factor of 15 at 2 × 10⁻³ mbar), due to the fact that the production of ammonia via the reaction NH₂ + H₂ $\to$ NH₃ + H dominates over the process of destruction through photolysis. At longer times, the amount of ammonia decreases, ending up lower than the initial steady-state at pressures greater than 10⁻² mbar for t > 1864 s. At this stage, NH₃ removal (dominated by the photodissociation NH₃ + hν $\to$ NH₂ + H) is dominant over ammonia production.

Carbon dioxide: The abundance of CO₂ changes down to 3 × 10⁻³ mbar. During the flare, one can observe a decrease of the mixing ratio by four orders of magnitude at 9 × 10⁻⁴ mbar, between the initial steady state and 1038 s. In the warm atmosphere case, this decrease is due to the photodissociation of CO₂: CO₂ + hν $\to$ CO + O(¹D). In the second part of the flare, there is almost no change in the abundance of CO₂, due to the balance between the production and destruction processes.

Nitric oxide: NO is affected by the stellar flare down to 30 mbar. Its evolution depends on the pressure level. For the first part of the flare (0–1038 s), from 30 to 0.2 mbar, y_NO increases as a function of time until the inflection point of the abundance profile (located at ∼0.2 mbar). At these lower altitudes, where NO increases, NO is produced mainly by thermal dissociation of HNO (HNO + H $\to$ NO + H₂) and reaction with nitrogen atom (N(⁴S) + OH $\to$ NO + H). The maximum enhancement of the mixing ratio (a factor of 10) is found at 5 mbar. Above this point, the abundance profile of NO decreases with time. At 1.2 × 10⁻³ mbar (pressure corresponding to the production peak of NO), the abundance of NO decreases by a factor of 4000 between 0 and 1038 s. We also notice that this peak is shifted downward with time. At these higher altitudes, the NO decrease is due mainly to photodissociation (NO + hν $\to$ N(⁴S) + O(³P)) and, at the highest altitudes, also to reactions with carbon atoms (C + NO $\to$ CO + N(⁴S) and C + NO $\to$ CN + O(³P)), produced by photodissociation of CO. From 1130 s, we observe a transition in the abundance at 4 mbar. Below this pressure, y_NO decreases with time, and for higher pressures, y_NO increases with time.

Hydroxyl radical: This species sees its abundances change down to ∼20 mbar. From 0 to 1038 s, y_OH increases with time from 20 mbar up to 2 × 10⁻³ mbar, the pressure at which there is a production peak. We observe an increase by a factor of four at 3 mbar between the abundances at 0 and at 1038 s.

At these lower altitudes (∼2 mbar) the reversible reactions O + H₂ ↔ OH + H and H + H₂O ↔ OH + H₂ are finely balanced. Interestingly, production of OH via photolysis of water (H₂O + hν $\to$ H + OH) tips the balance in favor of overall OH production, even at these high pressures. At lower pressures, in addition to the photodissociation of water, OH produced by the reaction O(³P) + H₂ $\to$ OH + H also affects the balance. Finally, at high altitudes, OH decreases due to the dominance of OH photodissociation and its reaction with carbon and oxygen atoms; namely, OH + C $\to$ CO + H and OH + O(³P) $\to$ O₂ + H. From 1130 s, we observe a decrease of the abundance with time for pressures between 20 and 1 mbar. For pressures between 1 and 2 × 10⁻³ mbar, y_OH increases slightly. Between 2 × 10⁻³ and 10⁻³ mbar, the abundance of OH slightly decreases with time. Finally, for pressures lower than 10⁻³ mbar, y_OH increases with time.

3.2. Long-term Effect of the Flare

After 2586 s, the stellar flux returns to quiescence. We study how long it takes for all species to return to steady-state. For the two atmosphere cases considered here, the final steady-states are reached after 10¹² s (between 10⁹ and 10¹² s, the variations are very minor). However, we see in Figures 5, 8, and 9 that some species exhibit differences between their initial and post-flare steady-states in the upper atmospheres (and down to 2 mbar for NO, in the case of the warm planet). These variations are more important for the T_eq = 1303 K atmosphere. In the hot atmosphere, the major species under consideration in this study—except for H, which hardly changes concentration—are less abundant after the flare event than before. This is not the case for the warm atmosphere, where we observe abundances changes in either direction, depending on species and pressure level. We suggest that these new steady-states are due to the permanent modification of the photolysis rates and the coupling between the different layers. The intense flare effectively bleaches the upper atmosphere, thus exposing lower levels also to more intense radiation. The resulting induced chemistry partly maintains this increased atmospheric transparency even when photon fluxes return to background levels after the flare. We ran numerous tests (i.e., cut off the photolysis processes after the flare event, apply the flare event from thermochemical equilibrium, change the initial conditions) that confirm our hypothesis. In a previous model, in which the mixing ratio of dihydrogen was 2000 times higher than what is shown here (i.e., 10% of y_H), the initial and post-flare steady-states were identical. This shows that the differences expected between the two steady-states depend also on the chemical composition of the upper atmosphere.

**Figure 8.** Evolution of H, NH₂, NH₃, CO₂, NO, and OH mixing ratios after the flare event with the thermal profile corresponding to T_eq = 1303 K. The legend for all figures is in the upper panels.
Download figure:
Standard image High-resolution image

**Figure 9.** Atmospheric composition at the initial steady-state (dashed line) and final steady-state (solid line), with the thermal profile corresponding to T_eq = 412 K (left) and T_eq = 1303 K (right). They correspond to the compositions obtained 10¹² s after the flare event.
Download figure:
Standard image High-resolution image

3.3. Recurrence of Flares

A flare is not a unique event in the lifetime of an active star. Pettersen et al. (1984), who studied the flare activity of AD Leo during observations of this star between 1971 and 1984, saw no obvious periodicity in the occurrence of the flares, but did determine that the duration between two flares ranged from 1 to 137 minutes. Such a random distribution of flare events has also been observed in other active stars (e.g., Lacy et al. 1976).

To account for the recurrence of flares, we ran two tests:

1. We re-ran our time-dependent model from the post-flare steady-state, up to a new steady-state. This case is not realistic because it implies that flares occur every 10¹² s, but can be considered as an extreme case for low-activity stars. Here again, we observed variations between the initial and final steady-states, but to a lesser extent than for our nominal case (Figure 10). We note that hydrogen returns to the initial steady-state (i.e., the previous post-flare steady-state) after this second flare event.

**Figure 10.** Atmospheric composition at the initial steady-state (dotted line), at the first final steady-state (after one flare—dashed line), and at the second final steady-state (after a second flare—solid line) with the thermal profile corresponding to T_eq = 1303 K. They correspond to the compositions obtained 10¹² s after each flare event.
Download figure:
Standard image High-resolution image

2. We ran our time-dependent model from the initial steady-state, imposing a period between flares of 2 × 10⁴ s (i.e., ∼5 hr), which is a more realistic timescale and close to the maximum duration between flares observed by Pettersen et al. (1984) on AD Leo. We chose this upper value, which is more likely to be representative of the occurrence of flares observed by Tofflemire et al. (2012) and Hunt-Walker et al. (2012), as well as the one in lower spectral type (M0-M3) stars, due to their lower magnetic activity (Kowalski et al. 2009; Davenport et al. 2012). We simulate a repetition of 47 flare cycles (representing ∼10⁶ s or ∼11 years). We found that the abundance of species continuously evolves at each flare event, but it seems that the maximum (and minimum) abundances reached at each cycle tend toward a limiting value (Figure 11) The number of cycles necessary to regulate the variations of abundances depends on both species and pressure level. Hydrogen reaches the same maximum abundance of 7 × 10⁻² at 0.022 mbar after only six cycles, whereas methyl radical (CH₃) needs about 20 cycles to reach its maximum abundance of 9 × 10⁻¹⁴ at 0.295 mbar. For ammonia, more than 47 cycles are necessary to reach a limiting value.

**Figure 11.** Temporal evolution of the molar fractions of H (top left), CH₃ (top right), and NH₃ (bottom) at different pressure levels, with the thermal profile corresponding to T_eq = 1303, assuming that the host star underwent 47 flares.
Download figure:
Standard image High-resolution image

3.4. Synthetic Spectra

We computed synthetic transmission spectra during the different phases of the flare and compared them to that of the initial steady state. Figure 12 shows the instantaneous transmission spectra for the initial steady state for the planets with T_eq = 412 K and T_eq = 1303 K. The planetary radius is defined at 1 bar. The spectra are binned at a constant resolution of R = 300. Figure 13 shows the differences between the initial steady state and the transmission spectra obtained during the different phases of the flare. It can be seen that, although small changes occur in the warm atmosphere, significant changes occur in the hotter planet's atmosphere. The strongest changes are seen in the CO₂ features at 4.6 and 14 μm in the hotter planet case, with total amplitude variations of, respectively, 220 and 200 ppm in the impulsive phase, and up to 500 ppm in the return to quiescence phase. Other changes with amplitude smaller than about 50 ppm are seen in other parts of the spectra, especially in the return to quiescence phase. Weaker variations of only 40 ppm are seen in the 412 K planet, especially in the CO/CO₂ features at 4.6 and 14 μm. Other very small changes of about 20 ppm are seen in other parts of the spectra.

**Figure 12.** Synthetic transmission spectra for the initial steady state for T_eq = 412 K (left) and T_eq = 1303 K (right) cases. The resolution is 100 across the entire spectrum.
Download figure:
Standard image High-resolution image

**Figure 13.** Difference in relative absorption between the initial steady state and the instantaneous transmission spectra, obtained during the different phases of the flare, for the T_eq = 412 K (top) and T_eq = 1303 K (bottom) cases. The left plots refer to the impulsive and gradual phases, whereas the right plots depict the return to quiescence phase. The color legends of the right plots are represented with a logarithm scale. The resolution is 100 across the entire spectrum.
Download figure:
Standard image High-resolution image

Two effects contribute to the observed changes. First, changes in the molecular abundances vary the strength of the absorption features specific to the different molecules. Second, these changes modify the mean molecular weight (μ) for the different atmospheric layers, amplifying or reducing the amplitude of the absorption features throughout the entire spectrum. This is because the scale height, i.e., the altitude at which the atmospheric pressure decreases by 1/e, is inversely proportional to the mean molecular weight (H = kT/μg; where: k is the Boltzmann constant, T is the atmospheric temperature, g is the gravity acceleration, and μ the mean molecular weight). Relatively small variations in mean molecular weight are seen in the warm planet case as the flare evolves. For the hotter planet case, we found that μ decreases (and hence H increases) at high altitudes (i.e., low pressures) during the flare, indicating that the density of the upper atmosphere decreases. However, we found that the decrease in abundances in the atmosphere outweighs the effect of an increased scale height on the transmission spectra.

We note that, although the mean molecular weight for each atmospheric layer is calculated taking into account the entire set of 105 molecules contained in the chemical model, we only compute the spectral opacities for a set of seven molecules that are the most abundant. The changes observed in the transmission spectra are, therefore, limited to these molecules. Other absorbers that might contribute to additional observable features will, however, show much weaker—if not unobservable—features.

3.4.1. Effects on Current and Future Observations

Transmission spectra of exoplanetary atmospheres are obtained by observing simultaneously the transit of the planet at different wavelengths. By tracing the transit depth as a function of wavelength, it is possible to reconstruct the transmission spectrum of the planet. Different molecules absorb differently in different spectral regions, making the planet appear smaller or larger at different wavelengths. A varying atmospheric composition during the transit event will introduce a time- and wavelength-dependent shift to the transit light curve.

Figure 14 shows the expected change in transit depth as a function of wavelength and time, with respect to the transit depths measured during the initial steady state. It must be noted that the timescale of the flare is on the same order of magnitude as the typical transit duration (∼10⁴ s). Therefore, if a transit is observed during the flare, the transit depth would change during the transit—shifts of up to 500 ppm are expected for the hotter planet. Clearly, these changes will be wavelength-dependent, as shown by Figure 12. Changes of only ≈40 ppm would be seen in the warm planet.

The bulk of recent observations of exoplanetary atmospheres have been obtained from space using the Hubble Space Telescope Wide Field Camera 3 (HST WFC3), covering the 1.1–1.7 μm range. Observations of hot Jupiters with WFC3 claim an uncertainty in transit depth per wavelength bin of 40–100 ppm (Deming et al. 2013; Fraine et al. 2014; Knutson et al. 2014; Tsiaras et al. 2015), for resolutions of about 50. For smaller super-Earths, such as GJ1214b, an uncertainty of 60 ppm was achieved (Kreidberg et al. 2014); however, this was only accomplished by combining observations of 12 transits, whereas for 55 Cancri e, a relative uncertainty of 22 ppm was achieved for a single transit observation (Tsiaras et al. 2016).

Other space-based instruments that have been used to obtain exoplanetary spectra include IRAC, IRS, and MIPS onboard the Spitzer Space Telescope (e.g., Richardson et al. 2006; Tinetti et al. 2007; Grillmair et al. 2008), as well as NICMOS, STIS, and COS onboard HST (e.g., Sing et al. 2008; Swain et al. 2008; Linsky et al. 2010). Several ground-based facilities have also been used (e.g., Bean et al. 2010; Snellen et al. 2010; Swain et al. 2010; Mandell et al. 2011; Sedaghati et al. 2015). Future missions, such as JWST, will achieve a similar or greater level of precision (Barstow et al. 2015).

The predicted changes in the spectra during a flare are, therefore, comparable to the sensitivity of current and future instruments, especially for the 1303 K planet case, where relative changes of up to 500 ppm are seen. We note, however, that the current instrumentation covers spectral ranges that are either outside the spectral features with the strongest changes (HST/WFC and STIS), or have only photometric channels with relatively low sensitivity (Spitzer/IRAC). Therefore, only with future instruments, such as JWST, which will cover a spectral range that includes the CO/CO₂ features that show the strongest variations, will it be possible to observe changes of spectral features of planetary atmospheres due to stellar activity in the form of flares.

3.4.2. Recurrence of Flares

Figure 15 shows the shifts as a function of wavelength and time that would occur over a period of ten days, assuming that the flare recurs every five hours. It can be seen that, compared to a single flare event, the changes in the transmission spectrum are even more dramatic. We see variations of up to 1200 ppm after one day, which would be detectable with future instruments onboard JWST. This shows even more clearly how strong stellar activity influences the atmospheric chemistry and the resulting spectra.

4. DISCUSSIONS

In this study, we focused on the effect of stellar flares on the chemical composition of the planet's atmosphere. One major simplification in this paper is that we did not change the thermal profile during the flare event, nor during the recovery phase. To be rigorous, the temperature of the atmosphere should be modified by both the increase of the irradiation (that could possibly also take place in the IR wavelength range) and the change of the atmospheric composition. However, if the change in temperature is about the same order as the one found in Segura et al. (2010) (i.e., less than 8 K), the chemical composition and thus the findings of this paper will not be significantly changed. However, because of the shorter radiative timescale of warmer planets (Showman et al. 2008), the change of temperature might be higher. A more complete study will be necessary to draw firm conclusions—for that purpose, we are currently developing a 3D self-consistent model including dynamics and chemistry. Thus, a more detailed study will be produced in a forthcoming paper.

Another aspect of stellar flares that is not discussed here is the effect of energetic particles. Indeed, a flare event can be accompanied by an ejection of energetic particles (Segura et al. 2010). In Earth-like planets, a high flux of energetic particles indirectly destroys ozone, because it favors the formation of nitrogen oxides, but in warm hydrogen-rich atmospheres, like the ones studied in this study, ozone is not present. Most studies dealing with the thermospheres of exoplanets (e.g., Yelle 2004; García Muñoz 2007; Koskinen et al. 2007) suggest that the charged particles will affect the atmospheric part above the one we studied here (P ≲ 10⁻³/10⁻⁴ mbar). However, the work of Koskinen et al. (2010) suggests that the ionization peak in the atmosphere of planets orbiting young host stars with a strong X-ray emission could be located down to 1 mbar. Thus, it would be interesting to quantify the effect of the possible increase of energetic particles expected during a stellar flare using one of the rare complex ion-neutral chemical scheme developed for warm exoplanet atmospheres (Lavvas et al. 2014; Rimmer & Helling 2015).

Concerning the spectra, we assumed a cloud-free atmosphere, in accordance to the thermal profiles and the atmospheric model. However, if clouds were considered, their effect would be to flatten the spectra, especially at small wavelengths. However, a cloud deck would not hide the strongest features at long wavelengths (such as the ones caused by CO₂). Given that the spectral changes are seen at these strong features, the actual differences would stay relatively constant.

Finally, although investigation of this effect is beyond the scope of this paper, we also note that the spectrum of a star varies as the flare evolves, introducing another systematic and time-dependent change on the transit light curves. However, stronger variations will only be seen in the bluest bands (<450 nm), compared to the red and in infrared, because emission in the blue is due to material heated by magnetic reconnection up to tens of thousands of degrees (Davenport et al. 2012).

5. CONCLUSIONS

We have studied the effect of a stellar flare on the chemical composition of sub-Neptune/super-Earth-like planets orbiting active M stars. As a case study, we chose AD Leo and used data recorded during a stellar flare event this star underwent in 1985. Two hypothetical planets have been modeled with different thermal profiles (T_eq = 412 K and T_eq = 1303 K). Using the 1D photochemical model of Venot et al. (2012), we followed the evolution of the atmospheric compositions during the flare. In both cases, these compositions change with time. Prominent species with relative abundances that are the most affected are H, OH, CO₂, NH₂, NO, and NH₃. We found that a single flare event could irreversibly alter the chemical composition of warm/hot planetary atmospheres, with final steady-state (reached 10¹² s after the end of the flare, i.e., ∼30000 years) being significantly different than the initial steady-state in both cases. Even if the frequency at which flares occur is random, it has been observed on AD Leo that the typical time between two flares is about a few hours (Pettersen et al. 1984; Hunt-Walker et al. 2012), which is far less than the time needed to reach the post-flare steady-state. To discuss this fact, we simulated a series of flares occurring about every five hours and found that the chemical abundances of species oscillate at each flare, around a mean value that evolves with time toward a limiting value. The number of flare events required to reach a limiting value depends not only on the species but also on pressure (and thus, also local temperature). This result suggests that planets around very active stars (with frequent flares) could never be at a steady-state, but rather are constantly and permanently altered by flare events.

Using the radiative transfer model of Waldmann et al. (2015), we computed synthetic transmission spectra for the two planet cases during the different phases of the flare (including the return to quiescence phase) and found variations in both cases. We see relative variations of up to 500 ppm in the CO₂ bands at respectively 4.6 and 14 μm in the hotter planet case, and changes of up to 40 ppm in the CO₂/CO at 4.2 μm feature in the the warm planet case. In particular, we note that the transmission spectra of the initial and final steady states are different, especially for the closer-in planet. Stronger variations of up to 1200 ppm are seen in the hotter planet case, assuming that the flare comes back every five hours over a period of ten days. These variations, although difficult to detect with current instrumentation, could easily be detected using JWST.

We thank Sarah Rugheimer for the quality of her referee work. Our paper has been greatly improved due to her useful comments. We also thank Antígona Segura and Lucianne Walkowicz for providing the stellar fluxes of AD Leo. We thank Ludmila Carone and Ingo Waldmann for their insightful comments; Vivien Parmentier and Kevin Heng for useful discussions about thermal profiles and radiative transfer; and Jonathan Fortney for providing the thermal profiles files. OV acknowledges support from the KU Leuven IDO project IDO/10/2013 and from the FWO Postdoctoral Fellowship program. L.D. acknowledges support from the ERC consolidator grant 646758 AEROSOL and the FWO Research Project grant G024112N. This work was supported by STFC (ST/K502406/1) and the ERC projects ExoLights (617119) and ExoMol (267219).

INFLUENCE OF STELLAR FLARES ON THE CHEMICAL COMPOSITION OF EXOPLANETS AND SPECTRA

Article metrics

Permissions

Author e-mails

Author affiliations

Dates

ABSTRACT

1. INTRODUCTION