KEPLER FLARES. II. THE TEMPORAL MORPHOLOGY OF WHITE-LIGHT FLARES ON GJ 1243

The first step in our flare selection from the short-cadence light curve was an unsupervised detection of candidate flare events. The starspots produce ∼3% flux variations, which are smooth sinusoidal features that evolve slowly (timescales ≳100 days). A great many flares exist in the light curve with amplitudes smaller than the starspot, and as such a simple flux threshold selection would only be useful in detecting the largest flares.

Instead, we subtracted the starspot features using a custom smoothing function. The smoothing function used a variable span smoothing method, inspired by the Supersmoother algorithm (Friedman 1984), and with a three-pass iterative approach. The light curve was first smoothed with a large boxcar kernel of 75 data points. All data falling more than 1σ away from the boxcar smoothed light curve were removed, and the entire light curve was then fit using a cubic spline. This process was repeated twice more, using progressively smaller boxcar kernels. The resulting smoothed model light curve was only minimally affected by the presence of large flares, and effectively traced the starspot. The starspot curve was then subtracted from the original light curve. We emphasize that while this iterative "sigma-clipping" method is not always the favored approach, it provided a rapid and effective means to identify flares, whose properties were then measured on the original data.

From the flattened (starspot-free) data, we selected all single epochs with positive flux excursions greater than 2.5 times the standard deviation (σ) of the light curve. Examples of these epochs are highlighted in Figure 2 as green diamonds. Flare candidates were required to have at least two consecutive epochs that passed this threshold cut. To avoid inadvertently splitting flares with complex morphologies (e.g., multi-peaked events) in to separate events, neighboring candidate events were merged in to a single flare event if they were only separated by one epoch below the 2.5σ threshold. We initially found that the recovered flare candidates systematically had a shortened duration compared to those defined by visual inspection of the first month of data examined. This was due to over-smoothing during the starspot removal. A final empirical correction to the end time of the candidate flares accounted for this truncation, which extended the candidate flare duration by ∼57%. The effect of this correction is evident in the flares shown in Figure 2, with the event duration extending beyond the epochs highlighted with green diamonds. The completeness of recovering flares with this automatic algorithm is presented in detail in Section 5.

Zoom In Zoom Out Reset image size

A crude flare-type classification was also automatically assigned. By default, all flare candidates were considered "classical" (one peak) unless they passed some conservative criteria. Candidate events with durations shorter than 4 minutes were labeled as "maybe" flares, indicating uncertainty in their identification. Events with a duration of at least 20 minutes, and having a secondary peak with amplitude greater than 30% the maximum peak (either before or after the maximum peak) were labeled as "complex." This simple type classification was not robust enough for scientific analysis, frequently missing true complex events, but was helpful in the human validation stage.

The sample of automatically detected flare candidates was then validated by manual inspection. To facilitate this validation, we developed a tool called Flares By EYE (hereafter FBEYE), an IDL suite of programs to view light curves and identify and classify flares. FBEYE also contains the auto-finding routines described above. Users interacted with FBEYE using an interactive graphical window, shown in Figure 2.

Briefly, the work flow for processing data with FBEYE is as follows: users first load the light curve of interest, typically a single month of short-cadence data. FBEYE runs the smoothing and auto-finding algorithms described in Section 2.1, and identifies candidate flares. Then, starting at the beginning of the light curve, the user will scan through the entire month of data, inspecting each flare. The user can choose to look at the smoothed, starspot-subtracted light curve, or the original data (shown in Figure 2). For each time window, the user will validate all the flares present. Flares that have been incorrectly selected are deleted. Events with incorrect start or end times are modified as needed. Flares that were missed entirely by the auto-finder are then defined. Flares with the wrong type classification, e.g., classical versus complex, have the correct classification assigned. The user then steps forward by half-window increments and repeats the validation processes. Every operation with FBEYE is automatically saved, updating the master flare list for the month.

Typical users would validate a single month of GJ 1243 data in 1–2 hr, while the auto-finder takes only a few seconds to run on a typical workstation. This highlights the great expense in human labor needed to validate these events in a data set as rich as Kepler, and for a star as active as GJ 1243. Each of the 11 months of short-cadence GJ 1243 data was inspected by at least five separate users, with a maximum of eight users, producing over 37,000 individual (though not necessarily unique) flare identifications.

Once the data were manually validated, we selected the final flares from a composite of all the user identifications. Figure 3 shows a representative day of data from GJ 1243, with epochs manually selected as containing a flare highlighted, and agreement between users indicated with colors. Many flares show a color gradient, particularly in the decay phase, indicating the regions of strongest disagreement between users. We selected our final flare start and end times to include epochs with at least two users indicating the presence of a flare. As with the auto-finding procedure, flare events were defined as continuous sequences of epochs. Events that were separated by a single non-flaring epoch were merged. The final flare classifications were simplified in to "classical" and "complex." Events that had been tagged as "maybe," either by hand or automatically, where called "classical," while the few "unusual" events were called "complex." The mode of the users' flare classifications was used for each event. The maximum number of users who agreed on each event, as well as the number of users who inspected the corresponding month's data, were also saved in the flare database.

Zoom In Zoom Out Reset image size

The final sample of flares from our FBEYE analysis of GJ 1243 contained 6107 unique events. Of these, 5162 were classified as "classical," and 945 as "complex" (15.5%). The general properties of these flares, such as the correlation between flare duration and energy from ∼1000 flares (two months of data) are presented in Paper 1. The sample of >6000 flares presented here is the largest catalog of stellar flares for a single object (excluding the Sun) that we are aware of, and will be a benchmark for future studies of stellar activity across the main sequence. In Table 1, we present the total number of flares and number of users who examined each month. While the month-to-month number of flares varies considerably, no overall trend in the flare rate was observed.

In Paper 1, we found the quiescent luminosity for GJ 1243 to be log L = 30.66 erg s⁻¹ in the Kepler bandpass using ground-based spectrophotometry. By integrating the fractional flux under each flare light curve, we computed the equivalent duration (hereafter ED; Gershberg 1972; Hunt-Walker et al. 2012), which has units of seconds. Multiplying the ED by the quiescent luminosity for GJ 1243, we found that the largest flares in our sample had energies of ∼10³³ erg, while the smallest were ∼10²⁷ erg.

The median flare template can be described to a very high precision with low-order analytic functions, which in turn makes fitting real flares straight forward. We have divided the flare template into two regimes: the rise and decay phases, or equivalently in Figure 4, t ⩽ 0 and t > 0, respectively. The resulting continuous functions can also be compared with theoretical predictions for flare heating and cooling curves.

The observed rise time for most flares in our sample was between one and five minutes (see Paper 1), yielding very little information about the morphology of the rise phase with the Kepler one minute cadence for any particular flare. However, with our high-fidelity template, details about the rise phase do emerge as shown in Figure 5(a). We fit the rise phase with a fourth-order polynomial⁷ spanning the time region −1 < t_1/2 ⩽ 0, chosen to reproduce the initial slow rise and nearly linear behavior near the peak. The fit was forced to go through the point (0, 1), in other words to equal a relative flux of 1 at maximum light. The best-fit solution took the form

$$\begin{eqnarray} F_{\rm {rise}} &=& 1 + 1.941(\pm 0.008) t_{1/2} -0.175(\pm 0.032) t_{1/2}^2 \nonumber \\ &&-2.246(\pm 0.039) t_{1/2}^3 - 1.125(\pm 0.016) t_{1/2}^4, \end{eqnarray} \tag{ 1 }$$

where values in parentheses indicate uncertainties on each fit coefficient. This analytic model is shown in Figure 5(a), and we emphasize it is valid only for the time region −1 < t_1/2 ⩽ 0. At times before t_1/2 = −1, the analytic model was forced to equal zero relative flux.

Zoom In Zoom Out Reset image size

The initially gradual rise, followed by a very rapid climb toward peak flux is similar to the morphology seen in ground-based white-light photometry for many of the "impulsive flares" observed in Kowalski et al. (2013). However, we were not able to resolve any significant "roll over" or flattening of the median flare shape as it approached maximum light, as has been observed with much higher time cadence observations (e.g., Kowalski et al. 2011). This is due to both the coarse sampling of our light curve compared to the duration of this peak phase, as well as our method of stacking flares based on their observed peak time (the mid-point of the exposure of maximum flux). We believe the latter reason is the primary cause for the "sharpness" of our empirical template at maximum light. Our flare model could be used in future studies to provide a robust peak time estimate for flares, and thus investigate the detailed morphology at the light curve maximum with higher time resolution data.

Previous efforts to fit the decay phases (often termed "cooling curves") for flares have used a variety of parameterizations. A generic exponential decay with time is frequently assumed in searches for flares in large catalogs (Walkowicz et al. 2011; Parke Loyd & France 2014). A departure from a single exponential (or "teapot") cooling curve toward a more detailed two-phase model for large M dwarf flares was explored as early as Andrews (1965), who parameterized the two phases as a steep linear decline, followed by a gradual inverse square shape. The linear regime was then postulated to be due to bremsstrahlung radiation of fully (or near fully) ionized hydrogen in a small region, while the gradual phase indicated radiative recombination. Hilton (2011) parameterized the decay phase observationally with an initial linear decline followed by an exponential profile. Recently, Kowalski et al. (2013) determined that the impulsive cooling phase in the white light was dominated by the decay of a 10⁴ K blackbody, which cools rapidly to ∼8000 K. A red continuum component, dubbed the "Conundruum," is often present during this impulsive cooling phase, and begins to dominate during the gradual cooling phase. Balmer continuum radiation is also observed throughout the flare cooling, but at a lower level.

The decay portion of our median flare template exhibited two clear exponential regimes, seen in Figure 5(b): an initially rapid decay, followed by a longer timescale gradual phase. The flare template is displayed using the natural log of the flux in Figure 5(b), so that straight lines correspond to exponential functions of the form $F(t) = a \, e^{\,(b\, t_{1/2})}$. A single exponential curve is clearly inappropriate in representing the flare decay. The shape was also not well fit using a single power law.

We fit the decay profile using two separate parameterizations. Our initial model used least-squares minimization to fit the median flare template with two exponential curves. The fit was computed in two well-separated time regions: 0 < t_1/2 < 0.5 and 3 < t_1/2 < 6. The resulting exponential fits were

$$\begin{eqnarray} F_1 = & 0.948\, e^{-0.965\, t_{1/2}}, \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray} F_2 = & 0.322\, e^{-0.290\, t_{1/2}}, \end{eqnarray} \tag{ 3 }$$

respectively. The two curves intersected at t_1/2 = 1.60, and relative flux = 0.20 (ln flux =−1.595), where the flare is presumed to instantaneously switch between decay profiles. This corresponds to the time of transition between the impulsive and gradual decay phases, and occurs at nearly the same relative flux value found by Hilton (2011). Such an abrupt threshold between the impulsive and gradual phases would represent a state change in the flare cooling, or rapid change between the dominant emission components.

The energy budget of the template flare shape can then be divided into three regimes: the impulsive rise, impulsive decay, and gradual decay, with the transition between decay phases set using the intersection of the two exponential curves defined above. Integrating the flare template within each regime, we found that the rise phase contains 19.9%, the impulsive decay phase 41.1%, and the gradual decay phase 38.9% of the total energy. By combining the impulsive rise and decay phase values, the fraction of total energy emitted during the impulsive phase observed in the Kepler band is 61%. This is close to the V- and R-band continuum energy fractions found in a very large flare on AD Leonis (Table 6 of Hawley & Pettersen 1991).

However, the template flare decay profile in Figure 5(b) indicates a smoother transition between decay phases. We thus fit the entire decay profile of the flare template with a continuous function, using the sum of two exponential curves:

$$\begin{eqnarray} F_{\rm {decay}} &=& 0.6890(\pm 0.0008) \, e^{-1.600 (\pm 0.003)\, t_{1/2}} \nonumber\\ && +0.3030(\pm 0.0009)\, e^{-0.2783 (\pm 0.0007)\, t_{1/2}}. \end{eqnarray} \tag{ 4 }$$

This parameterization would represent two physically distinct regions, each with its own exponential cooling profile, radiating throughout the entire flare decay. The initial decay would then be dominated by a brighter (presumably hotter) region that cools more quickly, and the gradual decay to a cooler region that cools more slowly, similar to the physical description of the emission components from the spectroscopic analysis of Kowalski et al. (2013). The transition time would correspond to the time when the two regions were equal in luminosity.

Likely neither of these simple parameterizations is entirely correct for describing the evolution of white-light emission region(s) in a flare. Thus, we have provided both sets of equations for use in comparing to future flare atmosphere models. For the remainder of our paper we adopt the latter parameterization, given in Equation (4).

The analytic (classical) flare template we have developed can be used as a model or basis function to decompose complex events. In this ansatz, complex events are described as the superposition of several classical flares. Linearly adding a series of our models then reproduces the observed complex morphology. For a given complex event, the task is to determine both the specific properties of each constituent flare, namely (t_peak, amplitude, t_1/2), as well as the number of classical flares needed to describe the event.

We fit model flares to every event in the flare sample with a (by-eye determined) duration of 10 minutes or greater. This sub-sample included 3737 individual flare events. For each event, we fit 1 to 10 individual flare models. This was done iteratively, and the final "best-fit" number of flares was determined using the Bayesian information criterion (BIC). Two example flare events fit by this procedure, which we describe in detail below, are shown in Figure 6.

Zoom In Zoom Out Reset image size

Every flare event was fit in fractional flux units relative to the global median flux (Figure 1), and using the detrended starspot-removed light curve. The flare region to fit for each event included 3 minutes before the user-defined t_start and 15 minutes after t_stop, to ensure that the flare fully returned to the quiescent flux level within the fitting window. As was done in constructing the median flare model, we subtracted any residual starspot signal using a linear fit to small windows of time before and after each flare region.

For the first pass (n = 1 model), we fit each event with a single flare template. We again utilized the IDL nonlinear least-squares fitting suite MPFIT (Markwardt 2009) to minimize the empirical model to the flare event. We seeded this minimization using the peak flux amplitude and time, and 15% of the full duration for t_1/2. The resulting n = 1 best fit was subtracted from the observed event, and the maximum peak in the residual flux was assumed to be an additional flare. To seed the n = 2 model minimization, we used the best-fit n = 1 model values, as well as the amplitude and time of the largest positive flux residual, and again 15% of the full event duration for t_1/2. The n = 2 model was then minimized, allowing both the primary and secondary flare to be fully solved to the observed data. We repeated the process of fitting and subtracting for n = 1 through n = 10 models on each flare event, solving for every flare component in each iteration.

Boundary conditions on the fit parameters were also imposed. Each component flare amplitude was required to be larger than twice the average photometric uncertainty within the flare event. We also required t_peak for each component to occur within the boundaries of the flare time window, and t_1/2 for each component to be larger than 1 minute and smaller than 50% of the observed flare duration. No priors on the relationships between flare amplitude, t_1/2, and t_peak were included in our fitting.

We then choose the "best" model to represent each event using the BICs. This statistic attempts to determine the improvement of the fit (decreased χ²) while penalizing the increasing number of free parameters used in subsequent models. We computed the BIC for all n = 1 through n = 10 models, which took the form

$$\begin{equation*} {\rm BIC}_n = \chi ^2 + k_n \ln (M)\,, \end{equation*}$$

where M was the number of observed data points that fell within the flare event window, and k_n the number of degrees of freedom in the nth model. The "best-fit" solution was then selected as the nth model with the smallest BIC parameter, where we additionally required the BIC to have decreased by at least 10% from the previous (n − 1) model. Finally, complex flares were defined to be any event best fit with a n > 1 model. The choice of a 10% BIC improvement threshold was determined by manual inspection of repeated fits to complex flares, and ensures we are not over-fitting these events.

In Figure 6, we show two examples of this fitting procedure. The top panel demonstrates a classical flare, very well fit by the n = 1 solution shown in red. The bottom panel contains an example of a very complex event, with four clearly separated peaks in the light curve. Our best model decomposed this event into seven distinct components, shown individually as colored lines. While the total morphology of this event is well fit by our procedure, a few discrepancies are apparent. For example, the broadest component flare (purple) does not account for the four or so low-amplitude secondary flares in the gradual decay phase. Additionally, the t_peak for this component aligns with a much smaller amplitude component (yellow), which may not be physically realistic. This broad component flare over-estimates the flux near the end of the flare event, further indicating that this component should have a shorter timescale, and additional component flares should be fit in the decay phase. This example illustrates both the utility of our template in decomposing these events and classifying them as complex, as well as the intrinsic degeneracies in such a procedure.

We have thus far described the light curve structure of complex flares using superpositions of our empirical flare template. This methodology implies the varied substructure seen in complex events is the result of additional flares, which have the same morphological properties as classical events. From this approach two physical interpretations for complex flare events are possible: either (1) the multiple component flares seen in complex events are physically associated within a single or nearby active regions on the star, or (2) complex events occur due to random superpositions of unassociated flares from separate active regions on the star. The former interpretation may include phenomena such as homologous flares, induced or sympathetic flares, and large tangled or multi-loop structures, all of which are seen on the Sun. Large complex flare events on M dwarfs, such as the one in Kowalski et al. (2010), have previously been interpreted in this manner (Anfinogentov et al. 2013). Given the high level of magnetic activity and flaring on GJ 1243, the latter interpretation that some complex flares are chance superpositions likely accounts for at least some of the observed complex events. Larger energy, longer duration flares clearly have a higher likelihood of overlapping other non-associated flare events.

Our user-validated flare sample included a flare classification of complex versus classical, as described in Section 2.2. However, deciding if an event is complex requires that the two (or more) underlying flares be sufficiently separated in time, at least by a few minutes, such that the peaks are distinct. Additionally, the amplitudes of secondary flares within complex events must be high enough to clearly be distinguished above the morphology of the primary flare. These two considerations mean that the complex flare rate determined by eye is a lower limit. In Figure 7 (top), we show the distribution of flare durations for all events in our by-eye validated sample. The fraction of events classified by eye as complex as a function of their duration is presented in Figure 7 (bottom, solid gray line). As naively predicted above, the rate of complex flares does increases with the duration.

Zoom In Zoom Out Reset image size

As described in Section 4.1, we fit all 3737 flare events with durations of 10 minutes or longer with n = 1 through n = 10 template fits. This produced a sample of 1141 flares (30.5%) that were best fit by an n > 1 model, as in Figure 6 (bottom), which we then classified as complex events. The fraction of model-fit complex events as a function of duration is also shown in Figure 7 (bottom, black points). The uncertainties were computed from the binomial error within each duration bin. The rate of complex flares from the model was higher than in our by-eye sample, as expected. Flares with durations greater than ∼50 minutes (log duration ∼1.7) show more than 80% chance of being complex events from the model, but only 60% from the users.

To compare with the observed fraction of complex flares as a function of duration, we generated two Monte Carlo models for creating complex flares in our data. In the first model, we generated random flares with durations 10 minutes or greater, drawn from a single power law with slope −1.99 ± 0.02, which was determined by fitting the observed flare duration distribution for events longer than 20 minutes, shown in Figure 7 (top, dashed blue line). A single power-law fit for flare energy distributions is commonly used (see Paper 1). In the second model, we used a broken power law, with a slope of −0.25 ± 0.03 for flares with duration between 10 and 20 minutes, and −1.99 ± 0.02 for flares greater than 20 minutes to more closely fit the data, as shown in Figure 7 (top, dashed and dotted blue lines). This two component, broken power-law model is similar to the flare energy distribution model of Kashyap et al. (2002), though we emphasize we do not extend our simulation to unobservable "microflares."

The procedure for both the single and broken power-law Monte Carlo models was the same. We began each model with a series of 50 trials, simulating a fixed number of flare events with durations ranging between 10 and 200 minutes. For these 50 trials, we adjusted the number of simulated flares between 3000 and 7500 events, with the first trial having 3000 flares and each subsequent trial increasing the number of flare events by 90. In each trial the specified number of flares was drawn from the respective duration distribution, and then placed at random start times throughout a blank light curve with the same time sampling as our data set. Any flare events that overlapped in time were combined, and the resulting event was classified as complex. We saved the resulting total number of events, both complex and classical, for each trial.

We used these 50 trials to determine the number of simulated flares required to reproduce the observed total number of flare events. A second-order polynomial was fit to the number of resulting versus number of simulated flares. For these models we required the resulting number of flares to match observations for the number of events (both classical and complex) with durations of 20 minutes or greater, which included 1750 flares from our by-eye sample. For the single power-law model, generating this number of ⩾20 minute events required an input of 6587 simulated flares, while the broken power-law model required 4147 flares.

Using the respective number of input flares needed to generate the observed number of ⩾20 minute events, we repeated this procedure 1000 times for both models. The fraction of complex events as a function of the resulting flare durations were recorded for each of the 1000 trials. The averaged complex flare fractions for both the single and broken power-law models are shown in Figure 7 (bottom, blue and red dashed lines, respectively).

Our data appear to rule out a single power law for the underlying duration distribution. While the single power-law model more closely reproduced the template-fit complex flare fraction curve, and was tuned to match the number of flares with durations 20 minutes or greater, it required far more flares than were seen in our data. The total number of short duration events (less than 20 minutes) resulting from this model exceeded our observation by over 1000. According to Paper 1, flares with durations between 10 and 20 minutes would be expected to have energies of log E ≈ 30.5–31 erg, and are easily detectable in the Kepler data for GJ 1243.

Instead, a broken power law is favored by our observations. Using the slopes fit from the observed duration distribution in Figure 7, the broken power-law model was able to reproduce the user-selected complex flare fractions. This model, however, somewhat under-produced the 10–20 minute duration flares as compared to our sample. However, we have made very simple assumptions about the duration distribution. The true, underlying power-law slopes for the distribution of individual flare durations are almost certainly different than what we observe, as many of the flares overlap and thus skew the observed duration distribution from its true shape.

The preferred (broken power-law) model did not fully reproduce the model-fit complex fraction. We believe this indicates that some of the complex structure our model-fitting scheme recovered must come from actual sympathetic flaring, for example, from the same or other nearby active regions. We note that for a given complex event, from the Kepler observations alone we cannot differentiate whether the complex structure was due to random superpositions or sympathetic flaring. The overall rate of sympathetic flaring is likely represented by the excess complex flare fraction above the broken power-law model. Observations with higher time cadence, and/or bluer wavelength coverage, may be able to detect additional lower amplitude complex flare structure, particularly for shorter duration flares. Further study on the complex flare fraction for other active stars would also be useful to constrain these model results. Studies of less active stars, where random superpositions are less frequent, will also help constrain the rate of sympathetic flaring.