Broadening and Redward Asymmetry of Hα Line Profiles Observed by LAMOST during a Stellar Flare on an M-type Star

The continuum levels of the original Hα spectra vary with time, e.g., in the fourth and fifth exposures the continuum levels are about twice as large as in other exposures. This variation is too big to be explained by a white-light flare. Instead, it is more likely caused a variation in seeing. Due to the lack of absolute flux calibration, we thus assumed that the continuum level remains unchanged before, during, and after the flare. However, we also point out that the Hα flare might be associated with a white-light flare, especially considering that the broadened Hα profiles observed in the flare (see Figure 2) indicate a nonthermal heating process (see details in Section 3.3 and also Namekata et al. 2020), but we cannot confirm whether a white-light flare exists from the available data. Before further analysis, we need to normalize the original spectra to the continuum intensities. To obtain the continuum intensities, we used a summation of a Voigt function, a Gauss function (only for the last four spectra), and a constant background to fit each Hα profile. The constant was taken as the continuum level. The Gauss function was introduced here to account for the additional bulk flow component, which can be inferred from the redward asymmetry. The Voigt function is a convolution of a Lorentz function and a Gauss function. A Voigt function fitting has often been applied to analyze the spectral broadening of hydrogen Balmer lines (i.e., Hα, H_β ) (e.g., Luque et al. 2003; Tremblay & Bergeron 2009; Falcon et al. 2015) and helium lines (e.g., Kanodia et al. 2022; Schad et al. 2021) in astrophysical and laboratory plasmas. Here, we used it to fit the Hα line profile because the chromospheric line broadening is likely dominated by a combination of pressure broadening and Doppler broadening during flares, which can be described as a Voigt function (Zhu et al. 2019). The total fitting function can be written as follows (Ginsburg & Mirocha 2011):

$$\begin{eqnarray}&&I(\lambda )=\displaystyle \frac{{\rm{Re}}\left[{\rm{wofz}}\left(\tfrac{\lambda -{\lambda }_{{\rm{cen}}}+i\gamma }{\sqrt{2}\sigma }\right)\right]}{\sqrt{2\pi }\sigma }{I}_{0}+{{ae}}^{-{\left(\tfrac{\lambda -m}{s}\right)}^{2}}+C,\end{eqnarray} \tag{ 2 }$$

where λ and λ_cen are the wavelength and centroid wavelength of the Hα line, respectively. Re represents the real part of a complex number. And wofz is the Faddeeva function, which satisfies ${\rm{wofz}}(z)={e}^{-{z}^{2}}{\rm{erfc}}(-{iz})$.

There are seven free parameters: I₀ characterizes the total flux of the Voigt function; γ is the half width at half maximum (HWHM) of the Lorentzian component. A larger γ means that the Voigt function is more similar to a Lorentz profile, indicating a larger effect of pressure broadening. The parameter σ is the standard deviation of the Gaussian component in the Voigt function. The height, centroid, and width of the Gauss function are represented by a, m, and s, respectively. The continuum level is represented by C. By dividing the original spectrum taken during each exposure by C, we obtained the normalized spectrum I_n (λ). All normalized spectra are shown in Figure 2.

We used the same function to fit the normalized spectra. After normalization, the fitting function becomes

$$\begin{eqnarray}&&{I}_{n}(\lambda )=\displaystyle \frac{{\rm{Re}}\left[{\rm{wofz}}\left(\tfrac{\lambda -{\lambda }_{{\rm{cen}}}+i\gamma }{\sqrt{2}\sigma }\right)\right]}{\sqrt{2\pi }\sigma }{I}_{n,0}+{a}_{n}{e}^{-{\left(\tfrac{\lambda -m}{s}\right)}^{2}}+1.\end{eqnarray} \tag{ 3 }$$

The subscript n means normalized intensity. The normalized spectrum, fitting function, and the residual between fitting and observation for each exposure are shown in Figure 3. Note that we removed the continuum just for the sake of simplicity. The variations of I₀, widths of the Gaussian and Lorentzian components of the Voigt–Gaussian fitting function, as well as the fitting parameters of the Gaussian function are shown in Figure 4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

From Figures 3 and 4, we noticed the following characteristics:

• 1.  
  The total flux of the Voigt function, I₀, increases quickly between 70 and 117 minutes, indicating the occurrence of a flare.
• 2.  
  The HWHM of the Gaussian component in the Voigt–Gauss fitting function is close to 0 for three spectra taken at 70, 93, and 117 minutes. At these occasions, the wings of the Hα line profiles are clearly characterized as a Lorentz function. Interestingly, this change of spectral profile starts as soon as the rapid increase in I₀. When the rapid increase in I₀ stops, the HMHW of the Gaussian component in the Voigt–Gauss fitting function begins to increase.
• 3.  
  a_n keeps on increasing after 93 minutes, meaning that the redward asymmetry of the Hα line profile becomes more and more significant.
• 4.  
  m slowly decreases after 93 minutes, indicating that the moving plasma responsible for the redward asymmetry of Hα possibly decelerates after 93 minutes (see Figure 4). From 93 minutes to 163 minutes, the average deceleration rate was estimated as 7.6 m s⁻².

The time variation of the integrated Hα flux can reflect the change in radiation power through the Hα line. We first integrated each normalized spectrum from 6557 to 6572 Å, then removed the continuum. The remaining integral is the equivalent width of the Hα line. Using the fitting parameters of the Gauss function, we can also obtain the equivalent width of the redshifted enhancement responsible for the line asymmetry in the last four exposures. The time variations of the equivalent widths of the Hα line and the redshifted enhancement are shown in Figure 5. We can see that the equivalent width of the Hα line increases quickly between 70 and 117 minutes, and decreases slowly after 117 minutes, which validates the conclusion from the analysis of the total flux of the Voigt function. At the same time, the equivalent width of the redshifted enhancement simultaneously grows by a factor of ∼6 from 93 to 163 minutes.

Zoom In Zoom Out Reset image size

A technique called red–blue (RB) asymmetry analysis has been often used to analyze the asymmetries of solar spectral lines both under the quiet-Sun condition (e.g., De Pontieu et al. 2009; Tian et al. 2011) and during solar eruptions (Tian et al. 2012). The technique can reveal where the asymmetry is most obvious. Here we used this method to analyze the asymmetry of the Hα line profiles. The RB asymmetry for an offset velocity can be expressed as

$$\begin{eqnarray}&&\mathrm{RB}({u}_{c})={\int }_{{u}_{c}-\delta u/2}^{{u}_{c}+\delta u/2}I(u){du}-{\int }_{-{u}_{c}-\delta u/2}^{-{u}_{c}+\delta u/2}I(u){du},\end{eqnarray} \tag{ 4 }$$

where u_c , δ u, and I(u) respectively represent the velocity from the line center, the velocity range over which the RB asymmetry is determined, and the ratio between spectral intensity and peak intensity. Here the Hα line center and peak intensity are λ_cen and I₀, respectively, as calculated in Section 2.1. We found that a smaller δ u leads to too much fluctuation, and a bigger δ u ignores too much data in the spectra. By trial and error, we found that 0.5 Å is just appropriate for our calculation. So we let δ u be the velocity corresponding to a Doppler shift of 0.5 Å, which is around 22.8 km s⁻¹. We calculated RB as a function of u_c , where u_c = 22.8 × (0.5 + n) km s⁻¹, n = 0, 1, 2, ⋯ .

The results of the RB asymmetry analysis for the spectra taken during the last four exposures are shown in Figure 6. We can see that an obvious redward asymmetry, which is indicated by the positive peak of an RB asymmetry profile, is clearly present after 93 minutes, when the Hα line intensity shows an impulsive increase. The redward asymmetry becomes more evident toward the end of the observation. The peaks of the RB asymmetry profiles indicate a redshifted velocity of 100–200 km s⁻¹, far below the escape velocity of the star (550 km s⁻¹).

Zoom In Zoom Out Reset image size

We tried to estimate the radiation power of the flare. We first estimated the luminosity of the continuum (L_c ) under the Hα line. Since the white-light continuum of a star generally can be approximated by the Planck function, L_c can be calculated through the following equation:

$$\begin{eqnarray}&&\begin{array}{c}{L}_{c}=4{\pi }^{2}{R}^{2}{\int }_{6557\mathring{{\rm{A}}}}^{6572\mathring{{\rm{A}}}}\displaystyle \frac{1}{{e}^{{hc}/\lambda {kT}}-1}\displaystyle \frac{2{{hc}}^{2}}{{\lambda }^{5}}d\lambda \\ =\,7.58\times {10}^{28}\,{{\rm{e}}{\rm{r}}{\rm{g}}\,{\rm{s}}}^{-1},\end{array}\end{eqnarray} \tag{ 5 }$$

where T represents the effective temperature of the star, which is 3404.72 K. L_c represents the luminosity of the continuum between 6557 and 6572 Å. R is the radius of the star, which is 0.411 R_☉.

Considering that the equivalent width of the normalized continuum in the wavelength range 6557–6572 Å is 15 Å and the equivalent width of the Hα line has a peak of 28.36 Å, we can calculate the peak luminosity (L_f ) of the flare through the following relation:

$$\begin{eqnarray}&&\displaystyle \frac{{L}_{f}}{{L}_{c}}=\displaystyle \frac{28.36}{15}\end{eqnarray} \tag{ 6 }$$

and L_f was estimated to be 1.43 × 10²⁹ erg s⁻¹.

We can also estimate the total energy of the flare from the integrated flux calculated in Section 2.2. This observation did not cover the whole process of the flare. We need to remove the Hα contribution from non-flare state of the star. We used the average of first three spectra to estimate it, and the result in the form of equivalent width is 5.43 Å. During the observation, from 70 minutes to 163 minutes, the flare energy radiated through Hα is 4.45 × 10³² erg. If we assume that the power of the flare and the equivalent width of Hα in the last three spectra decay at the same speed as the last three spectra, we estimate that the flare lasts for 418 minutes, and the total flare energy radiated through Hα was estimated to be 2.42 × 10³³ erg.

We have tried to estimate the mass of the moving plasma. We decided to use the method that was first introduced by Houdebine et al. (1990) and recently used by Koller et al. (2021) to estimate the mass. We used the following formula to obtain a lower limit on the mass (Houdebine et al. 1990):

$$\begin{eqnarray}&&M\geqslant \displaystyle \frac{4\pi {d}^{2}{F}_{m}({N}_{{\rm{tot}}}/{N}_{j}){m}_{{\rm{H}}}{\eta }_{{\rm{OD}}}}{h{\upsilon }_{j-i}{A}_{j-i}}=\displaystyle \frac{{L}_{m}({N}_{{\rm{tot}}}/{N}_{j}){m}_{{\rm{H}}}{\eta }_{{\rm{OD}}}}{h{\upsilon }_{j-i}{A}_{j-i}},\end{eqnarray} \tag{ 7 }$$

where d is the distance from the star and F_m is the integrated excess flux of the moving plasma. Since we did not have absolute flux calibration, we used L_m , the luminosity of the moving plasma, to estimate 4π d² F_m . N_tot/N_j is the ratio of the number density of hydrogen atoms and the number density of hydrogen atoms at excitation level j. m_H is the mass of a hydrogen atom. η_OD is the opacity damping factor. h is the Planck constant. υ_j−i and A_j−i are the frequency and Einstein coefficient of the spectral line, respectively.

As for other parameters in the formula above, we used the same settings as Koller et al. (2021) used. This flare event under study was detected from the Hα line (transition from j = 3 to i = 2). However, they could not find an appropriate calculation to estimate N_tot/N₃, so they decided to use the relevant data for the Hγ line (j = 5, i = 2) for their calculation. Houdebine et al. (1990) performed nonlocal thermal equilibrium (non-LTE) modeling for AD Leo and estimated N_tot/N₅ = 2 × 10⁹. They set η_OD to 2, which indicates that half of the radiation escapes the moving feature (Leitzinger et al. 2014). For Hγ, A_j−i is 2.53 × 10⁶ (Wiese & Fuhr 2009).

We used the same method mentioned in Section 3.1 to estimate the radiation power of the moving plasma through the red-wing enhancement of the Hα line L_m,α . In Section 2.3, we have calculated the equivalent width of the redshifted enhancement, and the equivalent width reaches a maximum in the last exposure. The equivalent width is 3.29 Å, which means ${L}_{m,\alpha }/{L}_{c}=3.29/15$. In this way we found L_m,α =1.66 × 10²⁸ erg s⁻¹. To obtain the radiation power of the Hγ line L_m,γ , we need to divide L_m,α by a factor of ∼3 (Butler et al. 1988), just as Koller et al. (2021) did. Finally, we obtained a lower limit for the mass of the moving plasma, which is 3.2 × 10¹⁵ kg.

From 70 to 117 minutes, the Hα line appears to reveal a Lorentz profile with a symmetric line broadening. Similar line broadening phenomena have been reported in observations of solar flares (e.g., Kowalski et al. 2017; Zhu et al. 2019) and also in studies of stellar flares (Fuhrmeister et al. 2018; Namekata et al. 2020), and have been interpreted as nonthermal broadening or Stark (pressure) broadening. There are two classical theories for the explanation of Stark broadening: the impact theory (Weisskopf 1932) and the statistical theory. The former results in a Lorentz profile and the latter leads to a Holtsmark profile (Holtsmark 1919). Since the observed Hα spectra can be well fitted by a Voigt–Gauss function that includes a Lorentz profile, here we take the former to explain our current observation. According to the impact theory, each time the radiating atom collides with a charged particle, the radiation is interrupted. After a Fourier transform, this interrupted radiation forms a Lorentz profile in the frequency domain. Moreover, the RADYN numerical simulation of M-dwarf flares from Namekata et al. (2020) recently suggested that, due to the non-LTE formation properties and highly variable opacity of the Hα line, the broadened Hα profile observed during flares is not only determined by the Stark broadening but also sometimes relates to "opacity broadening." Indeed, through changing the energy spectra from soft- into hard-energy input cases, their simulation results revealed that once hard high-energy electron beams appear during flares and deposit their energy in the deep chromosphere, strong self-absorptions and a Stark effect can simultaneously cause a wider line broadening. Therefore, this radiative-transfer-induced broadening process may also be important for the broadened flaring Hα line, because stellar and solar flares are often associated with intense nonthermal radiation/heating.

Clear Lorentz profiles were observed at 70, 93, and 117 minutes, which correspond to the impulsive increase in the Hα line intensity. On the basis of the theory mentioned above, pressure broadening implies that collisions between electrons and hydrogen atoms become more frequent, that is to say, electron density increases significantly. It is well known that in the impulsive phase of solar flares, energetic electrons are accelerated in the corona, then propagate downward along flare loops and finally collide with hydrogen atoms in the (deep) chromosphere. Meanwhile, as revealed in Namekata et al. (2020), if a strong nonthermal heating occurs in this stellar flare, the formation height of the Hα line can further extend into the lower chromosphere and then a larger opacity (or strong self-absorption) would also contribute to a wider line broadening.

After 117 minutes, the Hα line intensity stops increasing, indicating the end of the impulsive phase. From Figure 4 we can see that the HWHM of the Gaussian component in the Voigt–Gauss fitting function becomes nonzero again, meaning that the spectral line profiles start to deviate from a Lorentz function and the electron density begins to decrease, but the spectral line profiles are still wider than their pre-flare state. This likely indicates that nonthermal heating, i.e., energetic electron beams, is still continuing but with a softer/weaker energy injection. Indeed, Namekata et al. (2020) recently reported a good correspondence between the white-light emission and the variation of Hα line width in stellar flares of the M-type star AD Leo, suggesting that Hα line broadening has a fairly good relationship with nonthermal heating. In their observed flares, a similar decrease in the Hα line width from the impulsive to the decay phase was also found. Such a scenario is in accord with the physical explanation of the decaying phase of the flare in the standard solar flare scenario (e.g., Benz and Güdel 2010; Shibata & Magara 2011).

On the Sun, numerous spectroscopic observations of flares have revealed that chromospheric lines often exhibit line asymmetries during flares. Redward asymmetries indicative of velocities of several tens of kilometers per second have been frequently reported during the impulsive phase of the flare. This is commonly thought to be induced by chromospheric condensation (e.g., Fisher et al. 1985), the downward movement of the chromospheric plasma caused by the high pressure in the chromosphere. Downflows of dense plasma draining back to the solar surface along post-flare loops also lead to redward asymmetries of chromospheric spectral lines (e.g., Tian et al. 2015). Blue asymmetries or blueshifts of spectral lines have also been detected in the early phase of solar flares, which indicate upflows with speeds of a few to a few hundred kilometers per second in different spectral lines. This has often been attributed to the chromospheric evaporation, namely warm/hot plasma upflows driven by intense chromospheric heating. In the meantime, CMEs may also produce absorption line asymmetries if they contain erupting filament materials. For example, Den & Kornienko (1993) reported a strong mass ejection during a solar flare through obvious blue-wing absorption in Hα, indicating a radial velocity of over 600 km s⁻¹. Similarly, Ding et al. (2003) investigated a series of filament-induced absorptions in the Hα blue wing, revealing a good temporal correspondence between the acceleration process of a filament and the impulsive phase of a two-ribbon solar flare.

Doppler shifts or asymmetries of spectral lines have also been frequently detected during flares on other stars (e.g., Houdebine et al. 1990; Gunn et al. 1994; Crespo-Chacón et al. 2006; Vida et al. 2016; Fuhrmeister et al. 2018; Cao et al. 2019). Using X-ray observations, Argiroffi et al. (2019) recently reported a long-lasting blueshift of ∼90 km s⁻¹ of cool plasma (3 MK) in the decay phase of a large stellar flare on the evolved giant star HR 9024. They explained this cool upflow as a CME candidate. At optical wavelengths, spectral line asymmetries are more frequently reported (e.g., Leitzinger et al. 2014; Honda et al. 2018; Koller et al. 2021), although their explanations sometimes are not straightforward. For instance, Honda et al. (2018) found that the Hα line in a stellar flare on the dMe star EV Lac exhibited a long-lived blueward asymmetry, with a velocity of ∼100 km s⁻¹, even throughout the whole duration of the flare. They speculated that the blue asymmetries might be caused by possible filament activation/eruption and evaporation, but little is yet known about the definite origin of such long-duration blue asymmetries. Koller et al. (2021) have searched for asymmetries of the Balmer lines in time-resolved spectroscopy provided by the Sloan Digital Sky Survey data release 14 and reported six CME candidates with excess flux in line wings, with the highest projected velocities of 300–700 km s⁻¹. In fact, most (5/6) of these asymmetries appear in the red wings. Koller et al. (2021) suggested that these red-wing asymmetries may be caused by several possible processes, i.e., backward-directed CMEs occurring close to the limb, falling cool material during the eruption, or severe chromospheric condensation.

In our case, apart from obvious line broadening, red-wing asymmetries also appear during the flare process. From Figure 5, it is clear that the redward asymmetry becomes more evident after the impulsive increase in the Hα intensity. This strongly suggests that such line asymmetries are a flare-induced phenomenon. The line asymmetries we detected last toward the end of our observations, when the Hα intensity starts to decrease. This is timed to coincide with the line broadening after the flare peak, namely the continuing nonthermal heating. Moreover, the line asymmetries demonstrate an increasing amplitude until the end of our observations, with a redshifted velocity at 100–200 km s⁻¹. According to our current knowledge of solar flares, the long-lived red-wing enhancement may originate from three possible scenarios. (1) A flare-driven coronal rain process. The typical downward velocity of solar coronal rain is more than 100 km s⁻¹; as revealed in many observations of solar flares, coronal rain results from the radiative cooling process of hot flaring plasma that is trapped in post-flare loops, thus these loops often become more prominent in the decay phase of flares. These observed features are basically consistent with the properties of red-wing asymmetries we detected here. (2) A flare-induced chromospheric condensation process. In solar flares, chromospheric condensation results from nonthermal heating processes in flares and it usually happens in the impulsive phase, with a typical downflow velocity of tens of kilometers per second (e.g., Ding et al. 1995; Tian et al. 2015; Yu et al. 2020) but sometimes even up to or above 100 km s⁻¹ as well (e.g., Ichimoto & Kurokawa 1984; Tei et al. 2018). Here, the red-wing asymmetry is accompanied by obvious line broadening after the flare peak, thus suggesting that continuing nonthermal heating might still be responsible for the appearance of condensation downflows in the decay phase. (3) A flare-associated filament/prominence eruption either with a nonradial backward propagation or with strong magnetic suppression. In the case of a backward ejected filament, the flare must occur near or at the limb of the stellar disk, where its footpoint-produced Hα emission (especially line broadening) at least can be visible from the Earth. This special requirement regarding spatial location suggests that a backward filament eruption is rare. Another situation is a failed filament eruption, which was first reported on the Sun (Ji et al. 2003) and recently supposed to be common on M-type stars due to their strong magnetic suppression (Alvarado-Gómez et al. 2018; Li et al. 2020b). In this case, the ejected filament material may be quickly confined by an overlying magnetic field after the flare peak and may then demonstrate more and more prominent downward mass falling to the stellar disk, which thus might cause the increasing amplitude of red-wing asymmetries in our observation. Meanwhile, the ejecting filament at the flare onset is expected to cause a short-lived blueshift before we detected the red-wing asymmetries. But this blueshift was absent from the LAMOST observation, possibly due to its poor temporal resolution (around 23 mins). We caution that, by only using chromospheric spectral observations, it is very difficult to distinguish from our current observation between these three possible scenarios or their combinations, since the observed Hα radiation from the star could be a superposition of flare, post-flare loop, and prominence eruption. In future, to better understand these physical processes in stellar flares, new time-resolved spectroscopic observations using multi-temperature lines are needed.

The speed of the moving plasma is 100–200 km s⁻¹, far below the star's escape velocity. However, we should remember that this speed only represents the lower limit of the true velocity of the moving material, because the Doppler shift can only be used to measure the velocity component along the line of sight, and ejected filament material originating from the limb very likely has a small velocity component along the line of sight. Thus, we cannot exclude the possibility that a backward filament eruption might drive a CME if the red-wing asymmetries were to be an eruption. Moreover, as shown in Figure 4(e), the moving plasma responsible for the redward asymmetry likely decelerates. The average deceleration rate is computed as 7.6 m s⁻², which is much smaller than the gravity of the star. In the filament/prominence eruption scenario, this very small deceleration rate might support the idea that the ejecting filament material is flying in a nonradial direction that is nearly perpendicular to the line of sight, suffering a severe projection effect.

On the other hand, the acceleration of plasma during solar coronal rain is about g_eff/3 (Antolin 2020). If we assume that the falling materials move downward with a constant acceleration of g_eff/3 and the propagation distance scales with the radius of the star, the maximum velocity of the falling plasma on a star can be estimated as follows:

$$\begin{eqnarray}&&v={v}_{\odot }\sqrt{\displaystyle \frac{{gR}}{{g}_{\odot }{R}_{\odot }}},\end{eqnarray} \tag{ 8 }$$

where v_☉ is the maximum velocity of plasma during coronal rain on the Sun, which is about 100 km s⁻¹ according to Antolin (2020). g and g_☉ are the gravitational accelerations on the star and the Sun, respectively, and the stellar and solar radii are represented by R and R_☉, respectively. Our calculation yielded a value of 126 km s⁻¹ for v, which is consistent with the velocities inferred from our RB asymmetry analysis, supporting the possible coronal rain scenario mentioned above.