A ΔR ∼ 9.5 mag Superflare of an Ultracool Star Detected by the SVOM/GWAC System

As one of the main ground facilities of the SVOM ⁹ mission (Wei et al. 2016), the GWAC system located at Xinglong observatory of NAOC is an optical transient survey that images the sky in the optical down to R ∼ 16.0 mag at a cadence of 15 s. It aims to detect various short-duration astronomical events, including the electromagnetic counterparts of gamma-ray bursts (Wei et al. 2016) and gravitational waves (Turpin et al. 2020), and stellar flares. The main characteristic and the survey strategy of GWAC is presented below. More detailed information on GWAC is provided by Wang et al. (2020).

The effective aperture size of each GWAC JFoV camera is 18 cm. The f-ratio is f/1.2. Each camera is equipped with a 4096 × 4096 E2V back-illuminated CCD chip. The wavelength range is from 0.5 to 0.85 μm. The field of view for each camera is 150 deg² and the pixel scale is 117. For GWAC, each mount carries four JFoV cameras (called a unit in the GWAC system). The total field of view for each unit is ∼600 deg². Currently, four units have been installed at Xinglong Observatory, Chinese Academy of Sciences, China. More units will be installed before the launch of the SVOM mission in 2022, aiming to cover about 5000 deg² simultaneously. During the survey, each unit is assigned to a given grid, which is predefined for the whole sky according to the field of view of each unit. Areas of sky at a Galactic latitude of b < 20° as well as the grids near the Moon are set with a lower priority since the detection efficiency of any transient observation in these areas will be reduced by the higher star density or higher background noise.

A dedicated rapid follow-up system has been developed for each candidate by using two Guangxi–NAOC 60 cm optical telescopes (F60A and F60B) deployed beside GWAC with a typical delay time of one minute (Xu et al. 2020). More deep imaging and spectroscopy can be carried out through Target of Opportunity observations by the 2.16 m telescope (Fan et al. 2016) at Xinglong Observatory and by the 2.4 m telescope at Gaomeigu Observatory, China. The high cadence, moderate detection limit, self-automatic trigger capability, and its dedicated rapid follow-up telescopes enable the GWAC system to detect a great number of stellar flares and to capture events similar to the superflare ASAS-SN-16ae (ΔV < 11 mag, Schmidt et al. 2016) with greater temporal resolution.

On 2018 December 29 UT10:42:51, an alert was generated by the GWAC online pipelines for a very bright optical transient (GWAC 181229A) during a survey for one predefined field from 10:03:07.8 to 14:55:21.0 UT on the same night.

The detection magnitude was 13.5 mag in the R band measured by the real-time pipelines. The coordinates of the new source measured from the GWAC images are R.A. = 01:33:33.08, decl. = 00:32:23.02 (J2000). The corresponding astrometric precision is about 20 typically (1σ). This source was not detected in the reference image, which was obtained by stacking 10 images taken at around 10:04:21 UT, i.e., about 38 minutes before the trigger time. The finding charts of the detection image and the reference image observed by GWAC are shown in Figure 1. The candidate shows a stellar profile, indicating that it likely did not originate from a hot pixel, fast moving objects, or ghosts in the GWAC system. No apparent motion was obtained by the pipeline for the transient. No known minor planet or comet brighter than V = 20.0 mag was found in the 150 region around the transient. ¹⁰ All this information indicates that the transient is a real astronomical event with a high level of confidence.

Zoom In Zoom Out Reset image size

The online data processing showed that a fading of the transient by 0.9 mag can be seen in all single exposures within 2.5 minutes of the first detection by GWAC. The detection limit of all these single exposures was R ∼ 15.0 mag at a significance level of 3σ.

We re-performed an offline pipeline with standard aperture photometry at the location of the transient and for several nearby bright reference stars by using the IRAF APPHOT package, including corrections for bias, dark, and flat-field in a standard manner. After differential photometry, the finally calibrated brightness of the transient was obtained by using the Sloan Digital Sky Survey (SDSS) catalogs through the transformation of Lupton (2005). ¹¹

When the flare was triggered by the GWAC real-time pipeline, it was immediately followed-up by F60A ¹² in the standard Johnson–Cousins R band via a dedicated real-time automatic transient validation system (RAVS, Xu et al. 2020), which is developed to confirm candidates triggered by GWAC and to obtain an adaptive light-curve sampling for an identified target. With RAVS, the exposure time can be dynamically adjusted automatically based on the evolution of brightness of an object. For the case of GWAC 181229A, the range of exposure time is from 30 to 150 s. The follow-up observations by F60A started 2 minutes after the trigger, and stopped at the time when the object was fainter than the detection limit of ∼19.0 mag, which corresponds to a total duration of about 120 minutes.

The raw images were reduced by following the standard routine in the IRAF ¹³ package, including bias and flat-field corrections. The correction for dark current was not made since the impact on the photometry can be negligible if the CCD is cooled down to −70°C. After aperture photometry, absolute photometric calibration was performed with several nearby comparison stars with the transformation of Lupton (2005) from the SDSS Data Release 14 catalog to the Johnson–Cousins system. ¹⁴

Figure 2 compares the SDSS image centered on the target to the images obtained by F60A, in which there is a faint red counterpart within a distance of 0697 between the locations measured by F60A and reported by the SDSS Stripe 82 catalog (SDSS J013333.08+003223.7, Annis et al. 2014). Its brightness is r = 24.05 ± 0.15 mag (Annis et al. 2014), which is taken as the quiescent brightness for further analysis.

Zoom In Zoom Out Reset image size

One long-slit spectrum was obtained by the NAOC 2.16 m telescope (Fan et al. 2016) by using the Beijing Faint Object Spectrograph and Camera (BFOSC) ¹⁵ via a Target of Opportunity request. The observation start time for the spectrum was at 11:21:51.0 UT, 39 minutes after the trigger time. The exposure time was 30 minutes. The coverage of the exposure time during the flare is shown by the yellow shaded area in Figure 3. With a slit width of 18 oriented in the south–north direction, the spectral resolution is ∼10 Å when grating G4 was used, which results in a wavelength coverage of 3850–8000 Å. The wavelength calibration was carried out with iron–argon comparison lamps. Standard procedures were adopted to reduce the two-dimensional spectra by using the IRAF package, including bias subtraction and flat-field correction. The extracted one-dimensional spectrum was then calibrated in wavelength and in flux by the corresponding comparison lamp and standard calibration stars.

Zoom In Zoom Out Reset image size

In order to further investigate the nature of this object, it is crucial to analyze its properties in the quiescent state. We retrieved photometry from the SDSS (York et al. 2000), Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010), Pan-STARRS DR1 catalog (PS1, Chambers et al. 2016), and other catalogs based on a coordinate cross-match through the VizieR Service. ¹⁶ Each catalog returns only one source, named SDSS J0133, within our search radius of 2''. Some of the queried parameters are shown in Table 1.

At the beginning, based on the color–magnitude transformations given by Lupton (2005), ¹⁷ we estimate a quiescent brightness in the R band of 23.03 mag, which results in a flare magnitude as large as ΔR = 9.5 mag. The derived quiescent flux is F_R,q = 1.4 × 10⁻¹⁸ erg cm⁻² s⁻¹ Å⁻¹ by converting the quiescent magnitude above with the zero flux and the transformation for the R band (Bessell et al. 1998). Ahmed & Warren (2019) reported that the quiescent counterpart is of spectral type M9. Due to the faintness of this source, there is no parallax or other report about its distance, including in the Gaia DR2 catalog (Gaia Collaboration et al. 2018). With the corresponding SDSS i- and z-band magnitudes, based on the relation between color i − z and absolute magnitude provided by Bochanski et al. (2010, 2012), an absolute magnitude of M_r = 17.7 mag is derived for the quiescent counterpart. Consequently, a distance of d ∼ 155.8 pc can be calculated with the estimation of the absolute magnitude and the apparent magnitude above. The effect of reddening could be neglected for the above colors and the derived spectral type, since the extinction in the Galactic plane along the line of sight is not significant with E(B − V) = 0.021. ¹⁸ This distance is roughly consistent with the value of 144.6 pc reported by Ahmed & Warren (2019). The mean value of the distance of 150 pc will be used hereafter for further analysis.

However, it is noted that a spectral type of M7 would be obtained if the estimation were based on the i − z value provided by the PS1 catalog. The difference in the derived spectral type is possibly caused by the difference between Pan-STARRS and SDSS filters. The alternative possibility is that SDSS J0133 is active with a low amplitude at the PS1 survey time. Another clue to activity is the blue WISE infrared color of ∼−0.15 with W1 (15.366 ± 0.049) and W2 (15.517 ± 0.152) (Cutri et al. 2013), which is slightly bluer than the expectation (W1 − W2 ∼ 0.2) from the empirical relationships for UCDs reported in Schmidt et al. (2015).

According to the relation between metallicity and color of late-type stars (Equation (3) in West et al. 2011), the metallicity-dependent parameter ζ is estimated to be 0.859, which is slightly larger than the criterion for the classification of a subdwarf (ζ < 0.825, Lépine et al. 2007).

Figure 3 shows the optical light curve of GWAC 181229A, in which the data taken by GWAC and by F60A are shown by blue and red points, respectively. The horizontal red line marks the brightness level of the quiescent counterpart. The inset shows the GWAC data around the peak time. Before the first detection, the long-term monitors give an upper limit of 15.3 mag in the R band. In the late phase, there are some fluctuations at low confidence since the signal-to-noise ratio decreases with time. The vertical error bars are measurement-by-measurement estimates of the photon statistical error including instrumental characteristics. The horizontal error bars correspond to an exposure of 10 s duration.

With a cadence of 15 s, the first detection of GWAC 181229A shows that the brightness of the object was 13.9 mag in the R band, and the second one reaches the peak with a brightness of 13.5 mag. The brightness then falls to less than half the maximum in only two images, taking 30 s. The total duration of the flare from the onset to the quiescent flux level is estimated to be about 14,465 s by assuming that the brightness fades with a constant slope determined by fitting the late data as shown in Figure 3.

In order to have a more precise description of the morphology of the flare of GWAC 181229A, we fit the light curve for the decay phase by following the procedure of Davenport et al. (2014, hereafter D14), who tried to build a template from the single-peak flares detected in active flare star GJ 1243. Their procedure is as follows. For each flare, the flux and time after the onset are normalized to the quiescent level and the full time width at half the maximum flux (t_1/2), respectively. The key parameter t_1/2 can be obtained by (1) fitting the light curve as a free parameter, or (2) estimating in advance if the sampling of the light curve around the peak is dense enough. The decaying light curve is described by a sum of two exponential curves as in Equation (4) in D14, standing for the two components: the impulsive decay phase and the gradual decay phase.

For the case of GWAC 181229, the uncertainty in peak time is less than 7.5 s due to the GWAC's short cadence of 15 s. By assuming that the peak magnitude we detected is the real peak brightness of the flare, the amplitude of ΔR ∼ 9.5 mag corresponds to the relative flux of F_amp = 6500, which will be fixed during the analysis in our work. We here model the rising and the decaying phases separately as follows.

4.3.1. Rising Phase

In the template of D14, the rising phase is fitted with a fourth-order polynomial. However, for the case of GWAC 181229A, most of the observation data before the peak time are upper limits, except for one real detection. The behavior could not be well constrained with a fourth-order polynomial such as the template of D14. Here we have only to describe the rising phase of the flare briefly by assuming that this part follows a straight line for a few detections.

$$\begin{eqnarray}&&{F}_{\mathrm{decay}}/{F}_{\mathrm{amp}}={k}_{0}+{k}_{1}t\end{eqnarray} \tag{ 1 }$$

where F_decay is the relative flux and F_amp the peak relative flux that is fixed to be 6500. The values of k₀ and k₁ are calculated to be 0.69 and 0.02, respectively. The uncertainties of two parameters cannot be well estimated since there is only one positive detection before the peak. The uncertainties of these values are about 10% if only the precision of photometry measurements is taken into account. With this model, the onset time for the flare is about 35 s before the first detection, or 50 s before the peak time.

4.3.2. Decaying Phase

After modeling of the rising phase, we started by examining whether the D14 model can fit the observed data in the decaying phase. In D14, a sum of two exponential laws as shown in the Equation (2) was adopted to describe the light curve.

$$\begin{eqnarray}&&{F}_{\mathrm{decay}}/{F}_{\mathrm{amp}}={k}_{1}{e}^{-{\alpha }_{1}t/{t}_{1/2}}+{k}_{2}{e}^{-{\alpha }_{2}t/{t}_{1/2}}\end{eqnarray} \tag{ 2 }$$

where k₁ = 0.6890 ± 0.0008, k₂ = 0.3030 ± 0.0009, α₁ = 1.600 ± 0.003, and α₂ = 0.2783 ± 0.0007 as given in D14 are fixed in the subsequent modeling. By setting the peak flux (F_amp) and the timescale t_1/2 as free parameters, the best fitting returns F_amp = 3059 ± 63.6 and t_1/2 = 517.4 ± 12.0 s. The reduced χ²/dof = 3.63 with a degree of freedom of 54. The large χ² indicates that the template of D14 does not provide a good fit to the data, especially near the peak time as shown in the left panel in Figure 4. In fact, by checking the light curve by eye, the real t_1/2 should be around 30 s due to the sharp curve around the peak.

Zoom In Zoom Out Reset image size

To improve the fitting, we set the parameters in Equation (2) to be free except for F_amp = 6500, t_1/2 = 1. The modeled values are tabulated in Table 3, and the reduced χ²/dof = 2.65 with a degree of freedom of 52. The fitting results are shown in the right panel in Figure 4. In the upper panel of the figure, the total fitting result is displayed by the red line, and the two components by the blue and green lines, respectively. The time at which the two components have equivalent contributions is 793 s after the peak time. The lower panel shows the residual data that are obtained by a subtraction of the total fitting result from the observation data. The data near the peak time are still poorly reproduced, indicating that they might be from a new, steeper component that is not included in Equation (2).

In order to reproduce the light curve around the peak, we then model the light curve in the decaying phase by adding an exponential component:

$$\begin{eqnarray}&&{F}_{\mathrm{decay}}/{F}_{\mathrm{amp}}={k}_{1}{e}^{-{\alpha }_{1}t/{t}_{1/2}}+{k}_{2}{e}^{-{\alpha }_{2}t/{t}_{1/2}}+{k}_{3}{e}^{-{\alpha }_{3}t/{t}_{1/2}}.\end{eqnarray} \tag{ 3 }$$

A much better fitting with a reduced χ²/dof = 1.15 with a degree of freedom of 50 can be determined from Figure 5. The modeled parameters are again listed in Table 2. This good fitting suggests that there are three components in the decay phase. After the peak time, there is a very sharp decay component. At around 75 s, the light curve transfers to the second gradual component. After about 1500 s, the third, shallow decay is dominant until the end of the flare.

Zoom In Zoom Out Reset image size

Table 2. Parameters of the Modeled Decaying Light Curve of GWAC 181229A

k1	k2	k3	α1	α2	α3
Two-component model
0.444 ± 0.002	0.145 ± 0.007	...	0.005 ± 0.001	0.0005 ± 0.0001	...
Three-component model
0.373 ± 0.016	0.128 ± 0.008	2.248 ± 1.061	0.106 ± 0.008	0.014 ± 0.001	2.946 ± 0.895

Note. α₃ is for the first impulsive decay phase. α₁ and α₂ stand for the gradual phase and shallow phase, respectively.

Download table as:  ASCIITypeset image

A Bayesian information criterion (BIC) is used to test whether the three-component model used in the fitting is required or results from overfitting the data. The BIC values are 522.46, 660.03, and 602.56 for the three-component model, D14 model, and two-component model, respectively. These BIC values are also summarized in Table 3. This result confirms that the three-component model is more reasonable for the data.

Table 3. BIC for Three Models

Model	BIC
D14 model	660.03
Two-component model	602.56
Three-component model	522.46

Download table as:  ASCIITypeset image

Although some complex light curves have been observed (e.g., Kowalski et al. 2010), previous works found that the morphology of flare light curves is typically divided into two phases: an impulsive phase and a gradual decay phase (e.g., Moffett 1974; Moffett & Bopp 1976; Hawley & Pettersen 1991; D14). However, for the case of GWAC 181229A, three phases are needed to describe the high-cadence light curves well. The initial decay lasts 20 s after the first detection (5 s after the peak time), and is likely dominated by a brighter, hotter region that cools very quickly; then there is a gradual decay phase from about 20 to 350 s, which corresponds to a cool region in which the radiation cools slowly. Finally, the event moves into the shallower decay phase, lasting from about 350 s to the quiescent state.

4.3.3. Ratio of Decay Indices

Figure 6 shows the spectrum taken by the 2.16 m telescope. A series of strong emission lines such as Hα, He i λ5876, Hβ, Hγ, and Hδ are marked on the spectrum. The fluxes measured by a direct integration are presented in Table 4. After excluding the regions with strong emission lines, we modeled the underlying continuum by a blackbody in the wavelength range 4000–8000 Å, which returns a temperature of T_bb = 5340 ± 40 K.

Zoom In Zoom Out Reset image size

These emission lines are commonly detected during a dMe flare (e.g., Kowalski et al. 2013) and thought to be associated with chromospheric temperatures. By summarizing the flux of these strong emission lines as in Table 4, the total energy in the emission lines of 4.8 × 10⁻¹⁴ erg s⁻¹ cm⁻² in our range of observation wavelengths could be derived. The total emission of 5.13 × 10⁻¹³ erg s⁻¹ cm⁻² for the continuum within the wavelength range from 4000 to 8000 Å can also be measured. The ratio of the energy in the emission lines to that in the underlying continuum is about ∼9.3% for GWAC 181229A, which is higher than the percentage (∼4%) in the impulsive phase (Hawley & Pettersen 1991) and is significantly smaller than the values (17%–50%) in the gradual decay phase reported in the literature (e.g., Hawley & Pettersen 1991; Hawley et al. 2007).

Previous works in the literature show that the temperature in the gradual phase is lower than the values obtained at the peak time (e.g., Fuhrmeister et al. 2008; Schmitt et al. 2008). Our measured temperature of ∼5340 K in the shallow decay phase is similar to the reported temperature of 5500–7000 K in the decay phase of a flare event presented by Mochnacki & Zirin (1980), but is slightly higher than the values reported in the decay phase (Fuhrmeister et al. 2008; Schmitt et al. 2008), where a blackbody temperature of 3200–5600 K was given after measuring the continuum shape in their red higher-cadence spectra.

The equivalent duration (ED) of a flare is defined to be the time needed to emit all the flare energy at a quiescent flux level (e.g., Kowalski et al. 2013). By integrating the model of the light curve over the range of the light curve from the start to the end of the flare, the ED is estimated to be ∼2.584601 × 10⁶ s, or 29.9125 days for GWAC 181229A. Following the method of Kowalski et al. (2013), the total energy E_R in the R band can be calculated with the equation E_R = 4π r² × F_R,q × ED, where the quiescent flux F_R,q = 1.4 × 10⁻¹⁸ erg cm⁻² s⁻¹ Å⁻¹ and the distance is r = 150 pc; the energy E_R is found to be 1.54 × 10³⁴ erg. ¹⁹

To estimate the bolometric energy, one needs to know the effective temperature. In this work, our spectrum during the decay phase gives a temperature of 5430 ± 40 K by fitting a blackbody spectrum. On the other hand, the temperature at the peak time for a dMe flare could be as high as T_eff = 10⁴ K (e.g., Kowalski et al. 2013). More evidences indicates that the temperature will evolve during the flare from the peak time to the gradual decay phase (e.g., Hawley & Pettersen 1991; Hawley & Fisher 1992). Here, for simplicity, the bolometric energy will be estimated based on two effective temperatures: one is T_eff = 10⁴ K and the other is T_eff = 5340 K. By integrating the spectrum of a blackbody shape with these effective temperatures over the wavelength range from 1 nm to 3000 nm, and calibrating the energy with the R-band flux, bolometric energies E_bol of 9.25 × 10³⁴ erg and 5.56 × 10³⁴ erg for T_eff = 10⁴ K and T_eff = 5340 ± 40 K could be obtained, respectively. With the same method, the U-band energy of the flare is E_U ∼ 1.5 × 10³⁴ erg and E_U ∼ 3.6 × 10³³ erg for the two temperatures, respectively. Such a large amount of energy makes this flare comparable to the flare events SDSS J0221 (E_U = (3.2–5.5) × 10³⁴ erg) reported by Schmidt et al. (2016) and CZ Cnc reported by Schaefer (1990), and one of the largest energy events from UCDs.

The flare emission at optical and UV wavelengths is believed to be contributed by two major components. The dominant one is a hot blackbody emission (continuum emission) with a template of about T ∼ 10,000 K (e.g., Hawley & Fisher 1992), which is considered to be produced at the bottom of the stellar atmosphere near the footpoints of the magnetic field loops. The second component is the atomic emission lines (e.g., Fuhrmeister et al. 2010) and hydrogen Balmer continuum (Kunkel 1970). The proportion of the two contributors changes with the evolution of the flare. Near the peak time, the continuum emission could contribute more than 90% (Hawley & Pettersen 1991) of the total energy of the flare. In the gradual phase, the fraction of the continuum can drop to 69% (Hawley & Pettersen 1991) or even to 0% (Hawley et al. 2003).

The filling factor X_fill is the fraction of the area of the projected visible stellar disk that emits flare continuum emission, which allows us to understand what type of heating distribution is responsible for the observed light curve (Kowalski et al. 2013). Following the method of Hawley et al. (2003), X_fill in the impulsive and gradual phase can be deduced from

$$\begin{eqnarray}&&{F}_{\lambda }={X}_{\mathrm{fill}}\displaystyle \frac{{R}^{2}}{{d}^{2}}\pi {B}_{\lambda }(T)\end{eqnarray} \tag{ 4 }$$

where R is the stellar radius, d the distance, and T the characteristic temperature of the blackbody emission. F_λ is the flare flux observed at Earth at wavelength λ, which can be measured from the optical spectrum within a range of wavelength free of emission lines.

For the case of GWAC 181229A, only one spectrum was obtained at about 54 minutes after the event (mid-time of the exposure as presented in Figure 6). The continuum flux level is measured to be 1.8 × 10⁻¹⁶ erg cm⁻² s⁻¹ Å⁻¹ within the wavelength range 6800–7200 Å. There are no apparent emission lines within this wavelength range. Adopting R = 0.1 R_⊙ for the typical radius of an M9 brown dwarf (Baraffe et al. 2015), d = 150 pc, and a blackbody temperature of T_bb = 5340 K yields X_fill ∼ 19.3% for the decay phase, by assuming that all the emission measured within the wavelength range is blackbody emission.

Although no spectrum was obtained near the peak time, the temperature and the corresponding filling factor X_fill can be estimated as follows. Assuming that 95% observed peak emission is contributed by continuum emission, a critical temperature T_c = 10,000 K of blackbody emission is deduced, which corresponds to a filling factor of 100% of the surface of the object, indicating that the temperature of the blackbody emission near the peak time is much higher than T_c . Further calculations are made with T = 16,000 K, T = 20,000 K, T = 30,000 K, and T = 35,000 K to estimate X_fill, which results in X_fill of 36%, 24%, 13%, and 10%, respectively. We noted that Kowalski et al. (2013) reported that the temperature of the blackbody body is from T = 9800 to 14,100 K for the peak of the flares of the mid-M dwarf. If this is true for the later-M dwarf in GWAC 181229A, the value of X_fill is at the level of ∼30% at the peak time.

The maximum magnetic field strength B_z ^max associated with the superflare observed on GWAC 181229A could be estimated with the scaling relation in Aulanier et al. (2013) and Paudel et al. (2018) by assuming that the flare on GWAC 181229A is similar to solar flares:

$$\begin{eqnarray}&&{E}_{\mathrm{bol}}=0.5\times {10}^{32}{\left(\displaystyle \frac{{B}_{z}^{\max }}{1000{\rm{G}}}\right)}^{2}{\left(\displaystyle \frac{{L}^{\mathrm{bipole}}}{50\mathrm{Mm}}\right)}^{3}\,\mathrm{erg}\end{eqnarray} \tag{ 5 }$$

where E_bol is the bolometric flare energy and L^bipole is the linear separation between bipoles. Since the filling factor X_fill is at the level of 30% in the early phase, we could take L^bipole as π R as the maximum distance between a pair of magnetic poles on the surface of GWAC 181229A. With these parameters, a strong magnetic field of 3.6–4.7 kG is deduced. This magnetic strength is at the level of the saturated value of 3–4 kG (Reiners et al. 2009), and slightly smaller than the reported values of 7.0 kG for WX Ursae Majoris (Shulyak et al. 2017) and 5 kG for the M8.5 brown dwarf LSR J1835+3259 (Berdyugina et al. 2017).