K2 Ultracool Dwarfs Survey. III. White Light Flares Are Ubiquitous in M6-L0 Dwarfs

Ultracool dwarfs (hereafter UCDs) are stellar or substellar objects with effective temperatures of no more than 2700 K (Kirkpatrick et al. 1999; Martín et al. 1999). Understanding the nature of magnetic dynamos in UCDs has been very challenging. UCDs show deviations from age, rotation, and activity relations which are seen in early-M and mid-M dwarfs. The rotation rate plays a significant role in shaping the magnetic dynamo of early-M and mid-M dwarfs. When the stars are young, their rapid rotation rates empower strong magnetic dynamos. As they evolve, magnetic braking slows the rotation, which in turn decreases the magnetic activity (Gershberg 2005; Telleschi et al. 2005; Donati & Landstreet 2009). UCDs (which include both young brown dwarfs and old low-mass stars), despite being rapid rotators, display a poor relation between rotation and magnetic activity. Some of the usual indicators of activity (X-ray and Hα emission) weaken significantly among UCDs (Gizis et al. 2000; Berger et al. 2010; Schmidt et al. 2015). This could be due to a number of factors, including cool atmospheres which have reduced amounts of ionized gas (Mohanty et al. 2002), or atmospheres which are undergoing centrifugal coronal stripping (Jardine & Unruh 1999; James et al. 2000; Berger et al. 2008). On the other hand, observations of radio emission show that strong magnetic fields exist in UCDs (Route & Wolszczan 2012, 2016; Williams & Berger 2015). One interpretation of these data is that turbulent dynamos in UCDs may be producing both large- and small-scale magnetic fields (Reiners & Christensen 2010; Yadav et al. 2015). Alternatively, X-rays and Hα emission may be powered by fast magnetic reconnection, whereas radio emission may be generated by electrons which emerge from slow magnetic reconnection (Mullan 2010).

Stellar flares are transient events which are caused by sudden releases of magnetic energy in the upper atmosphere of the star. During this process, energy which was previously stored in magnetic form is converted (by reconnection) in part to kinetic energy of electrons and ions, in part to bulk flow of ejected matter (coronal mass ejections: CMEs; Benz & Güdel 2010; Mullan 2010), and in part to thermal energy. White light flares (hereafter WLFs) are assumed to be produced when non-thermal electrons accelerated after reconnection hit a cold, thick target in the lower chromosphere or upper photosphere. The precipitated electrons cause the formation of hot "chromospheric condensations" which emit a white light continuum: the continuum has a wavelength dependence in visible photons which approximates that of a blackbody with a temperature of order 10⁴ K (Kowalski et al. 2015 and references therein). During a WLF, a faint star can become significantly brighter in optical light, by as much as several magnitudes. Two huge WLFs on UCDs will be discussed later in this paper. Estimates of surface areas of flares show that the WLFs in UCDs cover larger fractional areas of the surface than do flares on bright stars such as the Sun (Kowalski et al. 2010; Walkowicz et al. 2011). Pineda et al. (2017) suggest that flares on brown dwarfs could result from planet-like auroral emissions produced by large-scale magnetospheric currents.

The discovery of the TRAPPIST-1 planetary system (Gillon et al. 2016, 2017; Luger et al. 2017) around an M8.0 dwarf suggests that more such planetary systems around other M stars or possibly brown dwarfs may be discovered. Thus, there will be increased chances of finding Earth-like planets in the habitable zones of such stars. To know the habitable conditions of such planets, it will be important to know flare rates on the host stars: ultraviolet (UV) photons from huge flares may have significant impact on the chemical evolution of the atmospheres of the planets (Segura et al. 2010; Grenfell et al. 2014; Armstrong et al. 2016; Owen & Mohanty 2016). In this regard, studies of the WLF rates of UCDs will contribute to the conclusions.

The Kepler mission (Koch et al. 2010) was originally aimed at finding more Earth-sized planets. But it has also turned out to be useful for studying stellar properties, including WLF rates, astroseismology, etc. WLF rates of several early-M and mid-M dwarfs were estimated using Kepler data (Martín et al. 2013; Ramsay et al. 2013; Hawley et al. 2014; Davenport 2016). The occurence of WLFs on L dwarfs and young brown dwarfs (Gizis et al. 2013, 2017a, 2017b; Schmidt et al. 2016) is proof that WLFs are common in some UCDs. In this paper, we continue our monitoring of various UCDs using Kepler K2 (Howell et al. 2014) data, with a goal of studying the WLF rates of various UCDs. Our targets, which include mainly late-M dwarfs and early L dwarfs, were monitored by Kepler K2 mission during Campaigns 3, 4, 5, 6, 10, and 12. Some targets in our sample have spectral type M6. We consider objects with spectral types ≥M6 as UCDs in this paper. The results of work presented here are an important step that in future works will be helpful to understand different flare properties (e.g., flare energy, duration, rate, etc) in UCDs and how these properties depend on spectral type, age, mass, etc. In cases where rotation periods and ages of targets are known, our study may shed light on the rotation–age–activity relationships in UCDs. We also include in this paper our own analysis of TRAPPIST-1 flares which were previously discussed by Vida et al. (2017) and Davenport (2017).

In this paper, we discuss flare properties in the context of the flare frequency distribution (FFD). For each star, the FFD is assumed to be fitted (over a range of energies) by a power law (Gershberg 1972; Lacy et al. 1976):

$$\begin{eqnarray}&&\mathrm{log}\tilde{\nu }={\alpha }_{o}-\beta \mathrm{log}E\end{eqnarray} \tag{ 1 }$$

where $\tilde{\nu }$ is cumulative (or integrated) flare frequency, i.e., the number of flares with energies ≥E which were detected per unit observation time. The constant α_o represents the intercept at zero energy, and the constant β represents the slope of the FFD. In stars which are found to have spectral index β > 1, the weakest flares contribute most to the total energy emitted by flares. In the stars with β < 1, the strongest flares contribute most to the total energy emitted by flares. Here, the phrase "total energy emitted by flares" refers to energy of all flares which were detected during a given observation time. The FFD can also be expressed in terms of differential form as

$$\begin{eqnarray}&&{dN}\propto {E}^{-\alpha }{dE},\end{eqnarray} \tag{ 2 }$$

where dN is the number of flares having energies in the range E and E + dE. The indices in Equations (1) and (2) are related as α = β + 1. Many studies have been undertaken to compute the FFD in the Sun and in early-M and mid-M dwarfs. Kurochka (1987) reported the value of the spectral index β to be ∼0.80 for the energy distribution of 15,000 solar flares observed during 1978 and 1979. Hilton (2011) calculated β = 0.73 ± 0.1 in the U -band energy range 10^27.94 ≤ E_U ≤ 10^30.60 erg for four M6–M8 dwarfs. Likewise, using a simple linear fit, Gizis et al. (2017a) reported β = 0.59 ± 0.09 in the energy range 10³¹ erg to 2 × 10³² erg for a field L1 dwarf WISEP J190648.47+401106.8 (hereafter W1906+40) and β = 0.66 ± 0.04 in the energy range 4 × 10³¹ erg to 1.1 × 10³³ for a 24 Myr brown dwarf 2MASS J03350208+2342356 (hereafter 2M0335+2342). Gizis et al. (2017a) used a maximum likelihood estimation (MLE) to obtain α = 1.6 ± 0.2 and 1.8 ± 0.2 for W1906+40 and 2M0335+2342 respectively. Assuming that all the flares with energy in the range from E_min to E_max follow a FFD with a uniform power law, the total energy of all flares during the observation time T can be computed using the spectral index β. This total energy is expressed as (Gershberg & Shakhovskaia 1983)

$$\begin{eqnarray}&&\varepsilon =T\times {10}^{\alpha }\beta ({E}_{\max }^{1-\beta }-{E}_{\min }^{1-\beta })/(1-\beta ).\end{eqnarray} \tag{ 3 }$$

In this paper, we report in Section 2 on Kepler K2 photometry of our 10 UCD targets, and we use the photometric data to estimate the energies of each flare. In Section 3, we discuss artificial flare injection and recovery. In Section 4, we present estimates of flare rates. In Section 5, we concentrate on the detailed properties of two superflares in our sample. A discussion of our results is presented in Section 6.

2.1. Targets

There are 10 UCDs in our sample. In Table 1 we list, for each target, the full name, the Kepler ID (EPIC), the Kepler magnitude, the 2MASS J magnitude, the K2 campaign number in which the target was observed, the tangential velocity, the optical spectral type, and the distance to the star. To calculate the tangential velocity for each target, we used the relation V_tan = 4.74 dμ, where d (in parsecs) is the distance of each target as given in Table 1 and μ is the proper motion. For all but one star, we used proper motions from Gagné et al. (2015): the exception is 2M0326+1919, for which we used Schneider et al. (2016). Our sample contains two M6, three M7, three M8, and two L0 dwarfs. Most of the targets are old low-mass stars and some may be brown dwarfs. The distances of two targets 2M2228-1325 and TRAPPIST-1 are taken from the literature. For the remaining stars, distances are estimated using either the M_J/ST (ST = spectral type) relationship from Dupuy & Liu (2012) or the i–z/M_i relationship from S. J. Schmidt et al. 2018 (in preparation). The i–z/M_i relation is an updated version of that from Schmidt et al. (2010). It is based on a linear fit to the color–magnitude diagram of 64 M5-L8 dwarfs and has typical uncertainties of ∼12%. The ninth column in Table 1 gives the references we used to identify spectral types of each target. The tenth column indicates whether the distances were obtained from the literature or were estimated using photometry: the reader is referred to the notes to the table for explanations of "mjst" and "miiz" in Column 10.

Table 1.  Target Properties

Name	EPIC	$\tilde{{K}_{p}}$	J	Cam. #	Vtan (kms−1)c	Spt.	Distance (pc)	Ref.	Remarks
2MASS J22285440-1325178	206050032	14.66	10.77	3	59	M6.5	11.26 ± 0.62	1	a
(LHS 523, GJ 4281, LP 760-3)
2MASS J22021125-1109461	206135809	16.72	12.36	3	25	M6.5	22.57 ± 4.16	2	mjst
2MASS J08352366+1029318	211332457	17.55	13.14	5	23	M7	32.26 ± 5.95	3	mjst
2MASS J22145070-1319590	206053352	17.74	13.46	3	55	M7.5	33.93 ± 6.26	4	mjst
2MASS J13322442-0441126	212826600	16.93	12.37	6	5.0	M7.5	20.54 ± 3.79	2	mjst
2MASS J23062928-0502285	200164267	15.91	11.40	12	61	M8	12.10 ± 0.40	5	b
(TRAPPIST-1)
2MASS J12215066-0843197	228754562	17.93	13.52	10	29	M8	32.11 ± 5.93	4	mjst
2MASS J03264453+1919309	210764183	18.08	13.12	4	61	M8.5	24.68 ± 4.55	6	mjst
2MASS J12212770+0257198	201658777	18.40	13.17	10	14	L0	19.66 ± 2.67	6	miiz
2MASS J12321827-0951502	228730045	18.85	13.73	10	28	L0	26.41 ± 4.87	6	mjst

Notes: mjst = Distance estimated using a combination of 2MASS J and spectral type based on the M_J/ST relationship from Dupuy & Liu (2012).

 miiz = Distance estimated using a combination i–z color based on the i–z/M_i relationship from S. J. Schmidt et al. (2018, in preparation).

^aDistance taken from Henry et al. (2004). ^bDistance taken from Gillon et al. (2016). ^cV_tan was calculated using the relation V_tan = 4.74dμ.

References: (1) Giampapa & Liebert (1986); (2) Cruz et al. (2003); (3) Burgasser et al. (2002); (4) Faherty et al. (2009); (5) Gillon et al. (2016); (6) Reid et al. (2008)

Download table as:  ASCIITypeset image

2.2. K2 Photometry

2.3. Flare Detection

Flare detection in the light curve of targets was a multi-step process. The initial step was to remove any periodic features in the light curve, which might be due to systematic or astrophysical variability. These features add complexity to the light curve, and alter the morphology and duration of flares. We began by smoothing the original light curve of each target using the Python package "pandas.rolling_median" (McKinney 2010) to remove the long-term trends mainly caused by systematic errors in the light curve (Davenport 2016). We used a window of w = 3 day data points (Handberg & Lund 2014). We fitted this smoothed light curve with a third-order polynomial (as suggested by Davenport 2016) and subtracted from the original light curve. We then followed a similar method as described in Osten et al. (2012) to identify the flares in the smoothed light curve. We calculated the relative flux F_rel,i for each data point in the smoothed light curve, defined as

$$\begin{eqnarray}&&{F}_{\mathrm{rel},i}=\displaystyle \frac{{F}_{i}-{F}_{\mathrm{mean}}}{{F}_{\mathrm{mean}}},\end{eqnarray} \tag{ 4 }$$

where F_i is the flux in ith epoch and F_mean is the mean flux of the entire light curve of each target. This relative flux was used to identify the flare candidates. We used a Lomb–Scargle periodogram to examine any other periodic features, which we expect to be mainly due to astrophysical variability, e.g., due to the presence of starspots. If any periodic feature was detected, we fitted the smoothed light curve with a sinusoidal function using the dominant period, and subtracted from the smoothed light curve. In this way, we prepared the detrended light curve for our targets. We then calculated a statistic ϕ_ij for each consecutive observation epoch (i, j) as

$$\begin{eqnarray}&&{\phi }_{{ij}}=\left(\displaystyle \frac{{F}_{\mathrm{rel},i}}{{\sigma }_{i}}\right)\times \left(\displaystyle \frac{{F}_{\mathrm{rel},j}}{{\sigma }_{j}}\right),j=i+1\end{eqnarray} \tag{ 5 }$$

where σ_i is the error in the flux associated with the ith epoch. This statistic, defined in Welch & Stetson (1993) and Stetson (1996), was used to study variable stars using automated searches. It was later used by Kowalski et al. (2009) and Osten et al. (2012) for flare searches in different stars. In order to identify the possible flare candidates in the light curve, we used the false discovery rate (FDR) analysis described in Miller et al. (2001). This method uses a critical threshold value of the ϕ_ij statistic which is different for each target. To calculate this critical value of ϕ_ij, we first discarded all those epoch pairs for which ϕ_ij > 0 but F_rel,i,j < 0. We then divided the remaining ϕ_ij distribution in two distributions: the null distribution, for which ϕ_ij < 0, and the possible flare candidate distribution, for which ϕ_ij > 0. The absolute value of the null distribution was fitted by a Gaussian function. The parameters of this Gaussian function were then used to calculate the p-values of each ϕ_ij in the flare candidate distribution. We then followed each step as described in Appendix B of Miller et al. (2001) to calculate the critical p-value and hence the critical ϕ_ij. Epochs with ϕ_ij greater than this critical ϕ_ij value were considered to be better flare candidates in the light curve. The value of variable α, which was used in the Miller et al. (2001) FDR analysis, was chosen to be 0.05 for the Kepler data (based on a private communication from R. Osten). This value of α signifies that no more than 5% of the epochs with ϕ_ij greater than the critical ϕ_ij are false positives (mostly due to noise in the data). We used additional criteria, namely that the detrended flux should exceed the photospheric level by 2.5σ. This decreased the number of flare candidates in our data set to a few hundred. The final flares were chosen by inspecting the data by eye. In this way, even for the weakest flares we ensured that there was at least a pair of epochs for which F_rel,i,j > 0: by this means, we were able to exclude any flares which had only a single measurement of flux brightening. For strong flares, there were multiple consecutive epoches with F_rel,i > 0. A more detailed explanation regarding this method of flare detection can be found in Osten et al. (2012).

2.4. Calibration of Equivalent Duration and Calculation of Flare Energy

To calculate the flare energies, we first estimated the equivalent duration (hereafter ED) of each flare. This depends on the filter used but is independent of the distance to the flaring object and so it is widely used for determining flare energies. The ED of a flare is expressed as

$$\begin{eqnarray}&&{ED}=\int [({F}_{f}-{F}_{c})/{F}_{c}]{dt}\end{eqnarray} \tag{ 6 }$$

where F_f is the flare flux and F_c is the continuum flux (i.e., when the star is in its quiescent state). It has units of time and gives the area under the flare light curve. It is the equivalent time during which the star (in its quiescent state) would have emitted the same amount of energy as the flare actually emitted (Gershberg 1972). We follow the method described in Gizis et al. (2017a, 2017b) to calibrate the ED in terms of energy. Kepler measures photometry in the wavelength range 430–900 nm. In the case of UCDs which have lower effective temperatures, a significant part of the flux is contributed by the longer-wavelength part of this range but the WLF radiation contributes flux throughout the whole range of wavelengths in the Kepler band. This means that a given number of WLF counts measured in the Kepler band will have higher mean energy than the same number of photospheric counts from the UCD (Gizis et al. 2013). For each target, we estimate the photospheric spectrum using the matching late-M or L dwarf template spectrum (Bochanski et al. 2007; Schmidt et al. 2014) normalized to match the Pan-STARRS i-band photometry (Tonry et al. 2012; Chambers et al. 2016; Magnier et al. 2016).⁵ We compute the photospheric specific flux of a 10,000 K blackbody which is normalized to have the same count rate through the Kepler filter as the photosphere of each target. (Using an 8000 K blackbody gives values only 2% lower; this reduction is much less than other sources of uncertainty.) We multiply this photospheric specific flux by the full width half maximum (FWHM) of the Kepler band pass (4000 Å), 4π d² (d = distance of target), and the Kepler ED to obtain the Kepler flare energy E_Kp. Since a 10,000 K blackbody is more energetic for the same count rate, we apply a correction factor of 1.3 to get the final estimate of E_Kp. The photospheric specific flux of a 10,000 K blackbody corresponding to each target is given in Table 5. Figure 1 shows the optical and near-infrared (near-IR) spectral energy distribution of TRAPPIST-1 and a 10,000 K flare with the same count rate through the Kepler filter.⁶ Computation of total flare energies integrated over the UV, visible, and IR is useful to compare the results with those obtained using other surverys. For the 10,000 K blackbody flare model we adopt, this bolometric energy (UV/optical/IR) is 3.1 times the E_Kp we report. Likewise, for an 8000 K blackbody the factor is 2.5. Gizis et al. (2013) argued that this range is simlar to that seen in an M dwarf flare by Hawley & Pettersen (1991). In Table 3, we list the time at which peak flare emission occurred, the ED, and the Kepler energy for all of the flares which we identified in our targets.

Zoom In Zoom Out Reset image size

Table 3.  Flare Properties

EPIC	Peak Flare Time	ED	log EKp
	(BJD—2454833)	s	erg
206050032	2145.4316	1.1e+00	29.5
206050032	2145.5278	6.3e-01	29.3
206050032	2146.7530	2.4e+02	31.9
206050032	2147.2441	3.1e+01	31.0
206050032	2147.7570	1.5e+00	29.7

Note. Only a portion of this table is presented here to show its form and content.

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

In order to obtain an estimate of the minimum flare energy that could be detected by our algorithm, we generated artificial flares of randomly chosen amplitude and duration using the Davenport (2014; hereafter D14) model. This was done via a slight modification of a similar module used in the software package known as "appaloosa" (Davenport et al. 2016). Then we injected the artificial flares at random times to the detrended light curve with 1σ noise. To prepare this detrended light curve, we followed the detrending process described in Section 2 and masked all other fluxes greater than 1σ level from the median flux. For simplicity, all injected artificial flares were single-peak "classic flares." Care was taken to avoid any overlapping of the injected flares. We injected 10 artificial flares at once and used our algorithm to detect them. We kept track of times at which the artificial flares were injected and their EDs. We repeated this process 1000 times, so a total of 10,000 artificial flares were generated; however, due to our restriction to non-overlapping events, some fraction of the 10,000 could not be injected. We then calculated the flare energies of the injected artificial flares and compared these values with the energies we recovered by means of our algorithm. We found that weak flares having energies less than a certain energy were not detected by our algorithm. The light curves of different targets had different noise levels, so the minimum detectable energy of weak flares as found by our algorithm was different for each target. To estimate this minimum energy, we repeated the above process separately for each target. A list of minimum energies of artificial flares injected and later detected by our algorithm is given in Table 4. For flare analysis, we discarded all flares (if any) having energies less than the minimum energy of the artificial flares detected by our algorithm. Figure 2 shows a detrended light curve with 1σ noise (upper panel), and with artificial flares injected at random times (lower panel) for one of our targets: 2M2228-1325.

Zoom In Zoom Out Reset image size

Table 5 lists various properties of flares on sample targets. In this table, "specific flux" is the flux per cm² per second per Angstrom for a 10,000 K blackbody which has the same number of counts as the target in Kepler bandpass. The column "Periodic Feature" gives information if any periodic feature is seen in the light curve of each target after removing the long-term trends using a rolling median method. If any periodicity is seen, the dominant period is listed. The periodicity might be due to the presence of starspots, and most probably gives an indication of the rotation period of the target. N is the number of flares observed on a given target during the entire interval T of K2 observations of that target (for T values, see Table 6). Likewise, E_Kp,min is the minimum Kepler energy of all the identified flares on a given target and E_Kp,max is the equivalent maximum.

Table 6.  Power-law Fits to FFDs of Targets

Name	T	β	α = β + 1	αo	$\hat{\alpha }$	Emin	Emax
	(days)					(1030 erg)	(1030 erg)
2M2228-1325	66.8	1.03 ± 0.05	2.03	−1.14 ± 0.05	1.90 ± 0.18	3.6	32
2M2202-1109	67.5	0.65 ± 0.02	1.65	−1.50 ± 0.02	1.49 ± 0.07	0.64	100
2M0835+1029	73.3	0.65 ± 0.05	1.65	−2.11 ± 0.03	1.43 ± 0.14	1.1	15
2M2214-1319	66.9	0.70 ± 0.02	1.70	−1.45 ± 0.02	1.54 ± 0.11	2.2	160
2M1332-0441	76.4	0.57 ± 0.03	1.57	−1.81 ± 0.02	1.49 ± 0.09	0.69	38
TRAPPIST-1	70.6	0.61 ± 0.02	1.61	−1.94 ± 0.01	1.47 ± 0.08	0.21	18
2M1221-0843	53.4	0.75 ± 0.02	1.75	−1.37 ± 0.02	1.66 ± 0.12	1.4	49
2M0326+1919	69.2	0.45 ± 0.04	1.45	−1.98 ± 0.03	1.32 ± 0.08	0.67	44
2M1221+0257	53.5	1.04 ± 0.08	2.04	−2.34 ± 0.02	1.85 ± 0.28	0.48	2.3
2M1232-0951	53.5	0.34 ± 0.04	1.34	−2.16 ± 0.03	1.37 ± 0.12	0.82	130
2M0335+2342	69.0	0.64 ± 0.03	1.64	−1.23 ± 0.05	1.77 ± 0.19	15	320
W1906+40	83.0	0.52 ± 0.05	1.52	−2.12 ± 0.03	1.37 ± 0.09	0.56	19

Note: The intercept α_o reported in this table gives the flare rate for a flare energy of 10³⁰ erg, not the zero energy. See the text for further description.

Download table as:  ASCIITypeset image

Table 6 lists the values of fitted parameters for the FFD of each of our targets. In this table, T is the total observation time for a given target, and β and α_o are the fitted values of parameters in Equation (1). Likewise, E_min and E_max are the minimum and maximum Kepler energies used for fitting. The total observation time T (in seconds) is computed by counting the total number of good (Quality = 0) data points and multiplying by 58.85 s which is the correct exposure time equivalent to one short-cadence. We used the maximum-likelihood method described in Hogg et al. (2010) and implemented in the routine known as "emcee" (Foreman-Mackey et al. 2013) to fit a straight line to our data (in log scale) and hence obtain the optimal values of parameters β and α_o. Here we report the intercept α_o corresponding to energy 10³⁰ erg, not the zero energy. It will be helpful in comparing the flare rates of our targets with previously reported flare rates of other targets observed by Kepler and K2, most of which have flare energies greater than 10³⁰ erg (Ramsay et al. 2013; Hawley et al. 2014; Davenport 2016; Gizis et al. 2017a). The routine emcee uses the standard Metropolis–Hastings Markov chain Monte Carlo procedure for marginalization and uncertainty estimation. We neglected the highest observed energy for fitting the line to reduce any bias in the analysis. Since 2M2228-1325 has a heavy-tailed distribution, we need to select a minimum value of energy to be considered for fitting. This minimum value was chosen on the basis of fitting for a broken power law discussed below. For targets which did not have a heavy tailed distribution, all flares, even those with the lowest energies, were considered for fitting. We also used the analytic solution to get the estimate of parameter α which is derived in Clauset et al. (2009) and references therein. Here we denote this estimate as $\hat{\alpha }$ which is the MLE of the true parameter α and is expressed as

$$\begin{eqnarray}&&\hat{\alpha }=1+n{\left[\displaystyle \sum _{i=1}^{n}\mathrm{ln}\displaystyle \frac{{E}_{i}}{{E}_{\min }}\right]}^{-1},\end{eqnarray} \tag{ 7 }$$

with error

$$\begin{eqnarray}&&\sigma =\displaystyle \frac{\sqrt{n+1}(\hat{\alpha }-1)}{n}.\end{eqnarray} \tag{ 8 }$$

Here, n is the number of flares and E_i, i = 1...n are the observed values of energies E such that E_i ≥ E_min. The MLE solution $\hat{\alpha }$ is an unbiased estimator of α in the asymptotic limit of large sample size, n $\to \,\infty$. In addition, a more reliable estimate for parameter α can be obtained for sample size n ≳ 50 (Clauset et al. 2009). For small sample sizes, Arnold (2015) suggests the value of $\hat{\alpha }$ obtained using Equation (7) can be multiplied by a factor of (n − 2)/n to make the result unbiased. The unbiased values of $\hat{\alpha }$ for each target are listed in Table 6. The same energy range used for fitting the power-law FFD is used for estimating $\hat{\alpha }$. As the number of flares on most of our targets is unfortunately rather small for a more accurate estimation of parameter α using the analytic solution, we will use the results obtained using the emcee routine for discussion and comparison of our results with previous works. Table 6 also lists the slopes of FFDs for 2M0335+2342 and W1906+40, obtained using the methods described above. It should be noted that they are slightly different from the values reported in previous papers because somewhat different energy intervals were chosen for fitting.

Figure 3 shows the FFD of each target in the sample. The FFD for each star is presented separately in each panel, and the panels are arranged in order of increasingly late spectral type. The FFDs of other targets in the sample are also plotted in the background of each panel to enable the reader to compare the FFD of a given target with others in our sample. Figure 4 compares the fitted FFDs of all the targets in the sample. It should be noted that the fitted line of 2M1221+0257 covers a very small range of energy, less than one order in magnitude. Likewise, Figure 5 shows a comparison of the FFD for a young brown dwarf 2M0335+2342 with the FFD for an L0 dwarf 2M1232-0951.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

4.1. Possibility of a Broken Power Law or a Power Law with Exponential Cutoff in the FFDs of Some UCDs

The presence of long tails at high energies in the FFDs suggests that a single power law might not be the optimal fit to the FFDs of all UCDs. In the case of two of our targets (2M2228-1325 and 2M0326+1919), we tried to fit two different models: (i) a broken power law, and (ii) a power law combined with an exponential cutoff. The results are shown in Figures 6 and 7. The broken power-law model can be expressed as

$$\begin{eqnarray}&&f(E)=\left\{\begin{array}{c}A{(E/{E}_{\mathrm{break}})}^{-{\alpha }_{1}}\,:E\lt {E}_{\mathrm{break}}\\ A{(E/{E}_{\mathrm{break}})}^{-{\alpha }_{2}}\,:E\gt {E}_{\mathrm{break}}\end{array}\right.,\end{eqnarray} \tag{ 9 }$$

while the power law with exponential cutoff can be expressed as

$$\begin{eqnarray}&&f{(E)=A(E/{E}_{0})}^{-\alpha }\exp (-E/{E}_{\mathrm{cutoff}}).\end{eqnarray} \tag{ 10 }$$

We used the astropy-affiliated packages "modeling.powerlaws.BrokenPowerLaw1D" and "modeling.powerlaws.ExponentialCutoffPowerLaw1D" to estimate the model parameters. We also used the astropy-affiliated package "modeling.fitting.LevMarLSQFitter" to fit the data. The last package uses the Levenberg–Marquardt algorithm and a least-squares statistic for fitting. We used a likelihood ratio test to decide which model better fits the observed data, by considering the broken power-law model as an alternative, and the power law with exponential cutoff as the null model. In the case of 2M2228-1325, we found that the broken power law provides a better fit to the data with α₁ = 0.20, α₂ = 0.83, E_break = 3.25 × 10³⁰ erg and A = 29.86. The p-value of the likelihood-ratio test statistic is 0.005 (i.e., we reject our null model for a significance level of 0.05). We used the energies greater than E_break for determining the slope of the FFD. But in the case of 2M0326+1919, we found that the power law with exponential cutoff provides a better fit to the data. In this case, the fitted parameters are A = 15.15, E_cutoff = 3.83 × 10³¹ erg, E₀ = 1.40 × 10³⁰ erg and α = 0.23. The p-value of the test statistic is 0.60 (i.e., we accept our null model for a significance level of 0.05). The possible break in the power law could be due to a number of factors, including the sensitivity of the instrument in observing weak flares, or saturation, or an upper limit on the energy which any flare on a given target is able to release (Gershberg 2005). As far as we have been able to determine, the FFDs which follow broken power-law are seldom discussed in the literature.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

One of the advantages of short-cadence Kepler K2 data is that they help to study the timescales associated with the rapid rise, rapid decay, and gradual decay phases of superflares (flares with energies ≥10³³ erg). Here we present the photometry of two superflares observed on our targets.

5.1. Photometry of Superflare Observed on L0 Dwarf 2MASS J12321827-0951502

We observed a superflare on L0 dwarf 2M1232-0951 at Kepler time 2811.766232, when the flare flux rose to 51,529 counts s⁻¹ ($\tilde{{K}_{p}}$ =13.44, ${\rm{\Delta }}\tilde{{K}_{p}}$ = −5.41): this count rate is ∼144 times larger than the photospheric level (358 counts s⁻¹). To determine the relevant timescales of the flare, we used the D14 model to fit the flare light curve. This yields values for the timescale t_1/2 associated with the rise of flux from and return to half maximum flux. The model also yields values for two decay timescales, one rapid (close to the flare maximum), and the other gradual (later in the flare). The D14 model uses a flare template that was based on the flare properties of the M4 dwarf GJ 1243 using Kepler short-cadence data. In this model, the rapid rise phase is best fitted using a fourth-order polynomial, and the decay phase is best fitted using sum of two exponentials which can be expressed as:

$$\begin{eqnarray}&&{\rm{\Delta }}F=A({\alpha }_{i}{e}^{-{\gamma }_{i}{\rm{\Delta }}t/{t}_{1/2}}+{\alpha }_{g}{e}^{-{\gamma }_{g}{\rm{\Delta }}t/{t}_{1/2}}),\end{eqnarray} \tag{ 11 }$$

where ΔF is flare-only flux, A is flare amplitude and Δt = t−t_f (t_f is the peak flare time). In the D14 template, the values of the different parameters of Equation (11) are α_i = 0.6890(±0.0008), α_g = 0.3030(±0.0009), γ_i = 1.600(±0.003), and γ_g = 0.2783(±0.0007). Schmidt et al. (2016) used this model to estimate the value of t_1/2 to be in the range 3 (best fit) to 6.2 (minimal fit) min for a superflare observed on the L0 dwarf ASASSN-16ae. Likewise, Gizis et al. (2017b) estimated t_1/2 = 6.9 min for the L1 dwarf W1906+40 using the same model. Slightly different values of the parameters were used by Gizis et al. (2017b) to fit a superflare on the L1 dwarf SDSSp J005406.55-003101.8 to estimate t_1/2 = 7.8 min. Here, the estimation of t_1/2 for 2M1232-0951 was done using Kepler short-cadence data, so it is more accurate than those reported for ASASSN-16ae and SDSSp J005406.55-003101.8. To fit the flare of interest to us here (i.e., the superflare on 2M1232-0951), we started by examining whether the D14 model could fit the observed data. However, we found that it did not provide an especially good fit for the late decay phase. To get an initial estimation of the decay phase timescale, we fitted the late decay phase separately by a single exponential curve and used the parameters obtained in this way to fit the entire curve. Figure 8 shows the observed flux and fitted model.⁷ At first glance in Figure 8(a) there seems to be a good agreement between the observed data and the fitted model. But the discrepancy can be clearly seen in the log–log version of same plot as shown in Figure 8(b). The values of the fitted parameters for this flare are: α_i = 0.9691, α_g = 0.0310, γ_i = 0.4551, and γ_g = 0.01543. Likewise, the other fitted parameters are A = 51005.76, t_1/2 = 1.054 min, and time of flare = 2811.7661 days. The observation shows that the flux decreases from its maximum value of 51,529 counts s⁻¹ to about one-half of its maximum value, which is 27,370 counts s⁻¹ in an interval of about 1 min. Since the best cadence time for gathering K2 data is also about 1 min, we estimate that a more accurate value of t_1/2 could be around 2 min for this superflare. The ED of this flare is 11.4 hr and its total bolometric (UV/optical/IR) energy is 3.6 × 10³³ erg. The total flare duration is 4.2 hr.

Zoom In Zoom Out Reset image size

Using the above fitted parameters and solving for the value of Δt = t−t_f at which impulsive decay switches to gradual decay, we get Δt = 7.83 t_1/2 after the peak flare, which occurs at Δt = 0. Now, using this time reference, we find that the rise phase contains 22.15%, the impulsive decay phase 38.70%, and the gradual decay phase 39.15% of the total energy.

5.2. Photometry of Superflare Observed on M7 Dwarf 2MASS J08352366+1029318

We also observed a superflare on M7 dwarf 2M0835+1029 at Kepler time 2379.882883, when the flare flux was 53,213 counts s⁻¹ ($\tilde{{K}_{p}}$ = 13.26, ${\rm{\Delta }}\tilde{{K}_{p}}$ = −4.29): this peak flux is ∼60 times larger than the photospheric level (891 counts s⁻¹). For this flare also, we found that the original parameters of the D14 model do not fit the data particularly well. We considered that we needed to make slight changes to the D14 values. We used a fitting procedure similar to that used for the 2M1232-0951 superflare discussed in Section 5.1. We found the following values for the fitted parameters of the M7 superflare: α_i = 0.8182, α_g = 0.1818, γ_i = 0.6204, and γ_g = 0.08467. Likewise, we found the other fitted parameters to be A = 58480, t_1/2 = 1.73 min, and time of flare = 2379.8825 days. Both the observed flux and fitted model for this superflare are shown in Figure 9(a). This figure and its corresponding log–log version (Figure 9(b)) show that there is better agreement between the observed flare curve and fitted model than in case of the superflare observed on the L0 dwarf 2M1232-0951. The ED of this flare is 7.7 hr and its total bolometric (UV/optical/IR) energy is 8.9 × 10³³ erg. The total flare duration is 4.5 hr.

Zoom In Zoom Out Reset image size

In this superflare, we found that that the value of Δt = t−t_f at which the transition between impulsive decay and gradual decay takes place is Δt = 2.81 t_1/2 after the peak flare (Δt = 0). Now, using this time reference, we find that the rise phase contains 19.57%, the impulsive decay phase contains 25.86%, and the gradual decay phase 54.58% of the total flare energy.

We identified a total of 283 white light flares in our sample of 10 UCDs. The flares we detected have Kepler energies in the range log E_Kp ∼ (29–33.5) erg, and the flares follow power-law distributions with slopes −α lying in the range from −1.3 to −2.0. The values of these slopes were determined by using a maximum-likelihood method of fitting a straight line, as implemented in the emcee software. We also estimated the values of the slopes using an analytical method. The presence of noise in the data makes it difficult to detect the weakest flares in UCDs. For this reason, the minimum detectable energy is >10²⁹ erg in our targets. The late-M dwarfs in our sample have FFDs comparable to those of active mid-M and late-M dwarfs studied by Hilton (2011). Compared to that work, the flares observed by Kepler have higher energy and are less frequent. Yet, we find similar slopes (−α). The slopes of late-M dwarfs lie between −1.5 and −1.8 except for the M6.5 dwarf 2M2228-1325 for which −α = −2.0. We should note that the slopes also depend on the range of energy chosen for fitting. We also analyzed flares on TRAPPIST-1 using the official K2 pipeline reduced data. We identified 39 good flares on it with bolometric (UV/optical/IR) flare energies in the range from 6.5 × 10²⁹ to 7.2 × 10³² erg. We find that its FFD has a slope of −α = −1.6. Previously, Vida et al. (2017) published the FFD of TRAPPIST-1 using raw data and estimated a similar value of the slope but in a slightly different energy range: 1.3 × 10³⁰ − 1.2 × 10³³ erg. In comparison to the other M8 dwarfs in our sample, TRAPPIST-1 has a steeper slope than the M8.5 dwarf 2M03264+1919 (−α = −1.5) and shallower slope than another M8 dwarf, 2M1221-0843 (−α = −1.8), in slightly different energy ranges.

In Figure 10, we compare the average flare rates of UCDs with a range of spectral types. The average flare rates are computed by mixing both active and less active targets of similar spectral types. For the L0+L1 catagory, we used the average flare rate of 2M1232-0951 (L0) and W1906+40 (L1). For comparison with previous works, we compare the flare rates of the M4 dwarf GJ 1243 and M5 dwarf GJ 1245 AB taken from Hawley et al. (2014) and Lurie et al. (2015), who also report the slopes of FFDs using Kepler energies of flares. We can see that the flare rates of GJ 1243 are higher than those of late-M dwarfs and early-L dwarfs by about two orders of magnitude for flares with energy 10³⁰ erg. In addition, Figure 10 shows that cool stars tend to have shallower slopes. The β values of averaged FFDs in Figure 10 are 0.84, 0.64, 0.60, and 0.43 for M6, M7, M8, and L0+L1 spectral types respectively.

Zoom In Zoom Out Reset image size

The L dwarfs have the lowest flare rates for weak flares. The FFD of the L0 dwarf 2M1232-0951 has a very shallow slope (−α = −1.3) compared to other UCDs in sample. The fitted line covers three orders in magnitude. A shallower slope signifies that the occurrence rate of bigger flares is higher in this target than others in our sample. The FFD of another L dwarf, W1906+40, as obtained by Gizis et al. (2013) also suggests that L dwarfs have shallower slopes. Figure 5 compares the FFD of 2M1232-0951 with that of the young BD 2M0335+2342. The convergence of two fitted lines at the observed high flare energies suggests that both targets have a comparable flare rate for such energies. For flare energy log E_Kp ∼ 29.0 erg, 2M1232-0951 has a lower flare rate with the difference being more than one order of magnitude. The convergence of lines is also clearly seen in Figure 10 in case of spectral types ≥M5, which again indicates comparable flare rates for observed higher flare energies. A more fascinating fact about 2M1232-0951 is that it has a flare with the second highest energy among all the targets in our sample. Why does it have better efficiency in converting magnetic energy to higher energy flares than those with low energies? We have to understand the relationship between the magnetic field, volume, and rate at which the magnetic field lines are stretched to store energy. The stretching of magnetic field lines may also depend on the speed of convective flows.

The steeper slope (−α = −2.0) which we have obtained in case of another L0 dwarf (2M1221+0257) is valid for a very narrow energy range (log E_Kp ∼ 29.7–30.4 erg). There is one high-energy flare on this target but its energy was not included for fitting purpose. It will be inappropriate to compare the flare rates of two L dwarfs based on the slopes obtained here as the energy ranges considered for obtaining the slopes are very different. As of now, white light flares are observed on five L dwarfs. Two such flaring L dwarfs are discussed in this paper. The remaining three are discussed in Gizis et al. (2013), Schmidt et al. (2016), and Gizis et al. (2017b). Figure 11 shows the bolometric (UV/optical/IR) flare energy distribution of the biggest flares observed on all five L dwarfs. ASASSN-16ae has the largest bolometric flare energy, which is >6.2 × 10³⁴ erg (Schmidt et al. 2016).

Zoom In Zoom Out Reset image size

The fraction of L0 and L1 dwarfs having chromospheric Hα emission is ∼90% and ∼67%, respectively, with a decline in H_α activity in comparison to earlier spectral types (Schmidt et al. 2015). This may be due to the lower effective temperatures and hence less ionization, reducing the effectiveness of the interaction between the magnetic field and gas (Mohanty et al. 2002). In addition, the L0−L1 dwarfs do not have clearly developed rotation−activity connections despite being rapid rotators (Reiners & Basri 2008). We do not have proper information about the rotation periods, or ages, or activity levels for the two L dwarfs 2M1232-0951 and 2M1221+0257. But what we do know is that 2M1232-0951 has a longer timescale variability (Koen 2013), and 2M1221+0257 has a variable Hα emission with equivalent width 25.65 Å and logL_Hα/L_bol = −4.18 (Reiners & Basri 2008). This limited information is not enough to interpret the observed results for L dwarfs. One possible physical process which might contribute to shallower FFD slopes in L dwarfs is discussed by Mullan & Paudel (2018). The model is based on production of flares by instability in coronal magnetic loops in which the footpoints of magnetic flux ropes are subject to random walk due to convective flows. Shallower slopes of the FFDs are found to arise in the presence of reduced electrical conductivity in the coolest stars.

We observed superflares on two targets: 2M1232-0951 and 2M0835+1029. Those flares have total bolometric (UV/optical/IR) energies 3.6 × 10³³ erg and 8.9 × 10³³ erg. They have very short FWHM timescales of ∼2 min. In the case of 2M1232-0951, the superflare brightened by a factor of ∼144 relative to the quiescent photospheric level. Likewise, the superflare observed on 2M0835+1029 brightened by a factor of ∼60 relative to the quiescent photospheric level. Another superflare with bolometric (UV/optical/IR) energy 2.6 × 10³⁴ erg was observed on a 130 Myr old brown dwarf CFHT-PL-17 (Gizis et al. 2017a). We have also detected the most powerful white light flare on another young brown dwarf (Paudel et al. 2018). The presence of superflares on so many cool targets (the flares discussed here, plus those discussed in Schmidt et al. 2016 and Gizis et al. 2017b) suggests that this is a universal phenomenon which is commonly observed on solar-type stars (Maehara et al. 2012; Notsu et al. 2013; Shibayama et al. 2013). It also suggests that very low-mass stars and brown dwarfs are magnetically active despite having cool atmospheres. The underlying nature of the magnetic dynamo is still unknown in such objects.

If we make use of the tangential velocity estimates as age indicators, we find that our targets have different ages. TRAPPIST-1 and 2M0326+1919 seem to be the oldest targets and 2M1332-0441 seems to be the youngest. In support of our claim of a great age for TRAPPIST-1, we may cite Burgasser & Mamajek (2017) who have reported an age of 7.6 ± 2.2 Gyr. Despite having the same tangential velocity and spectral type, the FFD of TRAPPIST-1 has a slightly steeper slope than that of 2M0326+1919; moreover, TRAPPIST-1 has a higher occurrence rate of low-energy flares than the same M8.5 dwarf. The FFD of 2M2228-1325 has a steeper slope than that of 2M2202-1109, although 2M2228-1325 has a larger tangential velocity and the same spectral type. Likewise, 2M1332-0441 has a shallower slope than other targets of similar spectral type, despite having the lowest tangential velocity, and is therefore presumably the youngest. In Figure 12, we compare the flare rates of targets having similar ages with different spectral types. The targets are categorized as very young, young, or old according to the tangential velocity estimates we have. The average tangential velocities and average slopes of the FFDs of the three categories of targets are listed in Table 7. 2M1221+0257 is not included for this purpose because its observed flare energies lie in a very narrow range. While there is a slight difference in the slopes of the FFDs of targets of various ages in Table 7, we cannot conclude anything definitive because of the small sample size.

Zoom In Zoom Out Reset image size

There is diversity both in age and spectral type of our targets. The diversity in the slopes of FFDs and flare rates observed in our targets suggests that the flare rate of a given target may depend on many factors, rather than just age and spectral type (effective temperature). Some of these factors may be rotation rate, magnetic field topology, the number of spots, etc. Due to the small sample size and smaller number of flares detected on faint targets, we cannot conclude any relation between the spectral type or age and the slope −α of the FFD. We need to study the FFDs for a larger sample of UCDs over a wide range of energies to see if any such relations exist. It will be interesting to compare the fl are rates of some targets which will again be observed by K2 in future campaigns. If we can use new data for the same targets, we may be able to confirm if the flare rates of UCDs remain constant, change slightly, or change drastically within short intervals of time. The results of Lee et al. (2016) suggest that the occurrence rate (R) of big solar flares and front-side halo CMEs are higher during the descending phase of the solar cycle than in other phases. They report a strong anti-correlation between R and annual average latitude of sunspot groups. Is this just a coincidence or due to some underlying physical phenomenon? We may also consider the same scenario for our targets. Route (2016) has reported some evidence of magnetic cycles in UCDs. The flare rate on UCDs might depend on the phase of their magnetic cycle. More information about rotational velocities of UCDs will be helpful to understand such correlations (if they exist).

The FFDs of two of our targets seem to show deviations from a single power-law dependence. Using the likelihood ratio test, we found that the FFD of one target (2M2228-1325) seems to follow a broken power-law distribution while that of another target (2M0326+1919) seems to follow a power-law distribution with exponential cutoff. Unfortunately, since the number of flares observed on both targets is small, we cannot conclude if such deviations of FFDs from regular power laws are due to instrumental sensitivity or due to saturation at large energies. Gershberg (2005, p. 227) mentions that the curvatures seen in FFDs of some targets were absent when they were observed again, and the number of observed flares was increased. Curved FFDs can be seen in EQ Peg, UV Cet, and AD Leo in Figure 38 of Gershberg (2005, p. 224). Other examples of departures of FFDs from single power laws are provided by GJ 1243 and GJ 1245 AB (Hawley et al. 2014), and by KIC 11551430 (Davenport 2016).

The flare rates of UCDs reported in this paper should be helpful in predicting the number of flares on targets of similar types, which could be observed by future photometric surveys. They are also very important for gyrochronology and studying planets in the habitable zones of stars like TRAPPIST-1. The biggest flares might be capable of damaging atmospheric chemistry and other habitable conditions of the planets. However, a detailed discussion about the possible impact of flares of TRAPPIST-1 on its planetary atmospheres is beyond the scope of this paper. One can find the necessary information in the recent papers of Roettenbacher & Kane (2017), Tilley et al. (2017), etc.

We thank Rachel Osten, Eric D. Feigelson and James R. A. Davenport for their valuable suggestions on some methods and statistical tools used for flare analysis and FFD fitting. We also thank the anonymous referee for providing valuable suggestions to improve this paper. The material is based upon work supported by NASA under award Nos. NNX15AV64G, NNX16AE55G and NNX16AJ22G. P.K.G.W. and E.B. acknowledge support for this work from the National Science Foundation through Grant AST-1614770. Some/all of the data presented in this paper were obtained from the Mikulski Archive for Space Telescopes (MAST). STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. Support for MAST for non-HST data is provided by the NASA Office of Space Science via grant NNX09AF08G and by other grants and contracts. This paper includes data collected by the Kepler mission. Funding for the Kepler mission is provided by the NASA Science Mission directorate. The Pan-STARRS1 Surveys (PS1) and the PS1 public science archive have been made possible through contributions by the Institute for Astronomy, the University of Hawaii, the Pan-STARRS Project Office, the Max Planck Society and its participating institutes, the Max Planck Institute for Astronomy, Heidelberg and the Max Planck Institute for Extraterrestrial Physics, Garching, The Johns Hopkins University, Durham University, the University of Edinburgh, Queen's University Belfast, the Harvard-Smithsonian Center for Astrophysics, the Las Cumbres Observatory Global Telescope Network Incorporated, the National Central University of Taiwan, the Space Telescope Science Institute, the National Aeronautics and Space Administration under Grant No. NNX08AR22G issued through the Planetary Science Division of the NASA Science Mission Directorate, the National Science Foundation Grant No. AST-1238877, the University of Maryland, Eotvos Lorand University (ELTE), the Los Alamos National Laboratory, and the Gordon and Betty Moore Foundation.

Facility: Kepler - The Kepler Mission.

Software: astropy (Astropy Collaboration et al. 2013), emcee (Foreman-Mackey et al. 2013).