Multiwavelength Campaign Observations of a Young Solar-type Star, EK Draconis. I. Discovery of Prominence Eruptions Associated with Superflares

EK Dra (HD 129333, G1.5V) is a young solar-type star at an age of 50–125 Myr. Its effective temperature of 5560–5700 K, radius of 0.94 R_⊙, and mass of 0.95 M_⊙ make this star one of the best proxies for an infant Sun at the time of the Hadean period on Earth (Waite et al. 2017; Senavcıet al. 2021). It is a rapidly rotating star, with a rotation period of 2.77 days (Audard et al. 1999; Ayres 2015; Namekata et al. 2022d), and it exhibits a high level of magnetic activity in the form of a hot and dense corona (e.g., Scelsi et al. 2005; Linsky et al. 2012; Namekata et al. 2023), strong surface magnetic field and large starspots (e.g., Berdyugina & Järvinen 2005; Waite et al. 2017; Järvinen et al. 2018), and frequent superflares (e.g., Audard et al. 1999, 2000; Ayres 2015; Namekata et al. 2022c; Namekata et al. 2022d). Although a faint low-mass companion is reported ∼20 au away in the projected distance from the G-type primary star (here we call the primary G-dwarf EK Dra), these two stars are not believed to be magnetically connected (Waite et al. 2017), and the G-dwarf is the main source of stellar superflares and associated eruptive events (see Section "Low mass companion" in Supplementary Information of Namekata et al. 2022d).

TESS observed EK Dra (TIC 159613900) in its Sector 50 from 2022 March 26 to 2022 April 22. ²¹ TESS conducts photometric observations using four cameras with the TESS filter in the optical band of 6000–10000 Å. Each sector is observed for approximately 27 days (Ricker et al. 2015). During this period, EK Dra was observed with a 10-minute time cadence mode. The data were obtained from the Multimission Archive at the Space Telescope (MAST) archive. ²² We used eleanor for photometric analysis to extract light curves from TESS full-frame image (FFI) data (Feinstein et al. 2019). Eleanor yields "raw" flux; "corrected" flux, which is regressed against instrumental effects; and "principal component analysis" ("PCA") flux, which has common systematics between targets on the same camera removed. We employed the "raw" flux for flare analysis due to the relatively high noise level of the "corrected" and "PCA" fluxes. For the analysis of rotational modulations, only the long-term rotational modulations in the "raw" flux were calibrated using a smoothed "PCA" flux.

We conducted optical spectroscopic monitoring observations using the Kyoto Okayama Optical Low-dispersion Spectrograph with optical-fiber Integral Field Unit (KOOLS-IFU; Matsubayashi et al. 2019), installed on the 3.8 m Seimei Telescope (Kurita et al. 2020) at Okayama Observatory of Kyoto University. The KOOLS-IFU is an optical spectrograph that provides a spectral resolution of R ∼2000 and covers a wavelength range from 5800 to 8000 Å. The campaign ran from 2022 April 10 to 2022 April 20, with the exposure time set at either 60 or 120 s, depending on weather conditions. We followed the data reduction procedure outlined in Namekata et al. (2020b), Namekata et al. (2022c), Namekata et al. (2022d) using IRAF ²³ and PyRAF ²⁴ packages. This produced a wavelength-calibrated and continuum-normalized 1D spectrum for each frame. The wavelength of the spectra was calibrated to account for the radial velocity shift caused by Earth's motion and the radial velocity of EK Dra (–20.7 km s⁻¹; Lindegren et al. 2018). Figure 1 shows examples of the obtained Hα spectrum of EK Dra, which is typical G-dwarf spectra with absorption (even during flares). We computed the Hα equivalent width (hereafter "EW," which refers to Hα emission integrated for 6562.8–10 Å ∼6562.8+10 Å), and used this value to plot the Hα light curves. Note that we defined an EW in terms of an emission profile to be negative. In this literature, on the other hand, we often refer to pre-flare subtracted equivalent width as ΔEW = —(EW–EW${}_{\mathrm{preflare}}$), where an emission is a positive value.

Zoom In Zoom Out Reset image size

NICER is a soft X-ray instrument (0.2–12 keV) on board the ISS (Gendreau et al. 2016). NICER monitored EK Dra from 2022 April 10 to 2022 April 21 (Observation ID: 5202490102–5202490113). ²⁵ These observations were scheduled concurrently with those from the 3.8 m Seimei telescope and TESS. We downloaded the data from NASA's website. ²⁶ The exposure time for each ISS orbit (about 90 minutes) varied from a few minutes to a few tens of minutes, depending on the dates.

We processed the data sets using the NICER calibration version CALDB XTI(20221001) and the nicerl2 tool within HEASoft version 6.31 (and NICERDAS version V010c). NICER's X-ray Timing Instrument (XTI) consists of 56 X-ray modules. In our data processing with nicerl2, we eliminated 6 modules (id = 11, 14, 20, 22, 34, 60) that were either inactive or noisy. We performed barycentric correction for each event file using barycor.

The XTI has a large collection area between 0.2 and 12 keV (approximately 1900 cm⁻² at 1.5 keV). We used xselect and lcurve to produce a light curve with a 2-minute time cadence and an energy band of 0.5–3 keV. As NICER is not an imaging instrument, there is some background contamination in the higher and lower energy bands. Thus, we chose 0.5–3 keV as a conservative energy range where the variation in background is relatively negligible.

We extracted the time-resolved X-ray spectra using xselect. The background files were created with nibackgen3c50 (Remillard et al. 2022). We performed spectral model fitting with xspec. We will explain the details of the model in Section 3.4.

Figure 2 overviews the observed multiwavelength light curves of EK Dra. As shown in the figure, we have labeled the superflares, targeted in this paper, as "E1" to "E3" (Table 2). We are targeting the events where there is an increase in Hα or X-ray that occurs simultaneously with the flares from TESS. Figures 3 to 5 show the enlarged light curves of flares E1 to E3. Flare E1 has shown the largest Hα emission increase (ΔEW of ∼0.7 Å) ever reported for solar-type stars, and an increase has been confirmed in both Hα and white light (WL). Simultaneous observation data in X-rays do not exist for flare E1 due to the observation gap of the ISS, but data before and after the flare do exist also in X-ray. Flare E2 shows a simultaneous increase in Hα and WL, and an increase during its decay phase has been detected in X-ray. Furthermore, for Flare E3, an increase has been detected in Hα, WL, and X-rays altogether. We will describe the TESS data analyses and results in Section 3.2, Hα spectral analyses and results in Section 3.3, and X-ray, light-curve, and spectral analyses and results in Section 3.4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 2. The List of Superflares on EK Dra in This Study and Past Observations

	GOES	EX,bol	EWL,bol	EHα	tWL	τWL	tHα	Ref.
		(0.1–100 keV)			(FWHM)	(e-decay)	(FWHM)
		[1033 erg ]	[1033 erg ]	[1031 erg ]	(min)	(min)	(min)
2022 Apr 10 (E1)	⋯	⋯	1.5±0.1	48.7±1.7(9.1±2.6) a	20	18	48	(1)
2022 Apr 16 (E2)	X4000–5400 b	${12.3}_{-0.5}^{+0.3}$-${16.7}_{-0.7}^{+0.4}$ b	12.2±0.2	17.7±0.7	40	20	54	(1)
2022 Apr 17 (E3)	X3200	${1.14}_{-0.05}^{+0.06}$	3.4±0.1	2.9±0.3	30	9.6	23	(1)
2020 Apr 05 (E4)	⋯	⋯	2.0±0.1	1.7±0.1 c	6.0	5.4	7.8 c	(2)
2020 Mar 14 (E5)	⋯	⋯	26±3	40±4	130	26	130	(3)

Notes. "GOES" is the goes X-ray class. E_X,bol is the bolometric X-ray flare energy in 0.1–100 keV. Note that the value of "GOES" and E_X,bol are calibrated for the light-curve extrapolation. t_WL and t_Hα are the FWHM duration of the WL and Hα flare. τ_WL is the e-folding decay time of the flare.

^a The value in bracket is the radiation energy of the central component without the blueshifted component. ^b The energy and flux values are already corrected for light-curve extrapolations. ^c The flare emission on April 5 in Namekata et al. (2022d) is very short-lived, and its light curve is a combination of blueshift absorption and flaring emission, so the Hα flare duration is expected to be underestimated. References: (1) This study, (2) Namekata et al. (2022d), (3) Namekata et al. (2022c).

Download table as:  ASCIITypeset image

We quantified the white-light flare (WLF) emission in the TESS band, following the procedures as taken in Namekata et al. (2022c) and Namekata et al. (2022d). The detection of flares first requires the removal of background rotational brightness variations. We removed them by using the fast Fourier transform numpy.fft.fft and a low-pass filter with the cutoff frequency at 3.0 day⁻¹. The detrended light curves are shown in panels (a) of Figures 3, 4, and 5. This study did not perform an automatic flare detection, but simply identified the presence of significant WLFs in the TESS light curve, specifically those occurring simultaneously with Hα and X-ray flares. In this context, a "significant" flare refers to an event where more than two consecutive data points exceed the photometric error (approximately 100 ppm) as computed by eleanor. As noted in Section 3.1, we confirmed the presence of significant flare increases corresponding to each of the Hα and X-ray flares E1, E2, and E3, as in the vertical lines in panels (a) of Figures 3, 4, and 5, respectively. Throughout this paper, we define the first flare point in TESS light curve as the flare standard time (BJD-2459680.0329, 2459685.9980, 2459687.0119 for flares E1, E2, and E3, respectively).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The calculation of bolometric WL flare energy follows the conventional technique proposed by Shibayama et al. 2013. We assumed a 10,000 K blackbody radiation spectrum for the flares and computed the bolometric radiation energy. As a result, we derived the WLF energy (E_WL,bol) of flare E1 to be 1.5_±0.1 × 10³³ erg, flare E2 to be 1.22_±0.02 × 10³⁴ erg, and flare E3 to be 3.4_±0.1 × 10³³ erg; hence, all of these are classified as "superflares." Namekata et al. 2022c reported that the largest energy scale of superflares on EK Dra is ∼5 × 10³⁴ erg, so these are relatively smaller than the upper limit. Also, we derived the FWHM duration (t_WL) and e-folding decay time of WLFs (τ_WL). The e-folding decay times (τ_WL) were obtained by fitting the decay phase of WLFs with an exponential function and scipy.optimize.curve_fit. These values are presented in Table 2.

To better constrain the white-light flare energy, we need to discuss the assumption of 10,000 K blackbody temperature. Most solar flare studies that employ spatially resolved data reported the emission temperature of 5,000–7,000 K (e.g., Watanabe et al. 2013; Kerr & Fletcher 2014; Kleint et al. 2016; Namekata et al. 2017b), while only Sun-as-a-star observations suggest the hotter temperature of the blackbody spectrum ∼9000 K (Kretzschmar et al. 2010; Kretzschmar 2011). On the other hand, a historically well-known stellar superflare on an M-dwarf, AD Leonis, shows the broadband stellar spectra with the effective temperature of 9000–10,000 K (e.g., Hawley & Fisher 1992), and thus many stellar flare studies use this temperature as a template. A recent statistical survey shows that the flare temperatures vary depending on flare energies for M-dwarf flares (Howard et al. 2020). More recently, COROT's two-band photometric observations constrained the temperatures of G-dwarfs superflares as 3570–24,300 K (Rabello Soares et al. 2022). Hence, no unified model for radiation temperature of solar and stellar WLFs is currently available. Namekata et al. 2017b concluded that the uncertainty in the emission temperature does not largely affect the estimation of the WLF energy because the decrease in the temperature leads to the decrease in surface luminosity (∝T⁴) and the increase in emission area (A ∝ T⁻⁴), both of which are canceled by each other when using Equation (1) of Shibayama et al. (2013; i.e., E = σ T⁴ A). Therefore, our assumption may provide an uncertainty of the flare energy by a factor of few. We need to perform multiband observations of superflares in U, B, V, and near-UV bands to derive better constraints on the flare energy (Brasseur et al. 2023).

3.3.1. Overview of Hα Flare Spectra

Figures 6, 7, and 8 shows pre-flare subtracted Hα spectra of flares E1, E2, and E3. Each flare show different flare spectra as follows:

• (E1)  
  The flare E1 exhibits significant blueshifted emission profiles during the flare. During the impulsive phase of WLF, it does not demonstrate strong emission in the Hα line center and most of the radiation comes from the blueshifted component (see Section 3.3.3). Following the noticeable blueshift profiles, the blueshift wavelength gradually reverts to the Hα central wavelength, disappearing within approximately 34 minutes (see Section 3.3.2).
• (E2)  
  The superflare E2 primarily manifests emission in the Hα line center. However, the period spanning roughly 10 to 15 minutes from the TESS flare onset displays significant blueshifted emission profiles that exceed the 3σ noise level of the far-wing continuum.
• (E3)  
  The superflare E3 exhibits the weakest Hα EW enhancement and does not display any significant Hα line asymmetry.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.3.2. Spectral Fitting

In regard to flares E1 and E2, which demonstrated explicit blueshift phenomena, we performed a spectral fit using Gaussian functions to derive the velocity, velocity dispersion, and EW of the blueshifted component. In this process, we conducted two types of fits based on whether or not to isolate and fit the central flare component. Method (1) simply fits the emission component with a single Gaussian function (referred to as "one component" fit). Method (2) first fits the central component by using the spectra only the longer wavelength of Hα line center (>6562.8 Å) with a Gaussian function to obtain the central Hα component. Following that, we subtracted this fit function from the original spectrum and then fitted the residual spectrum again with a Gaussian function, which we defined as the blueshifted component (referred to as "two components" fit). In both methods, we evaluated the errors on the y-axis based on the data scattering of the far-wing continuum level and performed the fitting using scipy.optimize.least_squares. Method (2) has also been employed in studies such as Maehara et al. (2021) and Namizaki et al. (2023). We use a Gaussian function, which involves fewer parameters than a Voigt function, because of the low signal-to-noise ratio (S/N) of our data. Figures 9 and 10 show the fitting results for the flares E1 and E2, respectively. The following is the summary of the obtained results of fitting:

• (E1)  
  For flare E1, we utilized both methods (1) and (2). This is because the profile seems to be totally blueshifted in its initial phase. Figure 9 displays the fitting results up to t = 63 minutes, where method 1 fits successfully. Regarding method (2), we show the fitting results up to t = 25 minutes, where the significant blueshifted component still remained after removing the central component. The maximum blueshifted velocity V_max, defined as $| \max \{V(t)\}|$, was 690_±82 km s⁻¹ at t = −4.8 minutes for both methods. However, the period from t = −9.3 minutes to −4.8 minutes, where relatively high velocities are present, does not yet have a high signal compared to noise. We therefore sometimes describe the maximum velocity as V_max=330_±35−690_±92 km s⁻¹ (one-component) and 490_±33−690_±93 km s⁻¹ (two-component) to include a conservative estimate (where the conservative value is the maximum velocity after t > −2.5 minutes).
• (E2)  
  For superflare E2, considering that the blueshifted component is clearly separated from the central component, we only applied method (2). The maximum velocity V_max is 430_±19 km s⁻¹.

These values are summarized in Tables 3 and 4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 3. Properties of the Prominence and Filament Eruptions on EK Dra

	Hα asymmetry	${V}_{\max }$	Vtypical	Vdisp,typ	−dV/dt	LHα,blue	Mp	Ekin	Ref.
		(km s−1)	(km s−1)	(km s−1)	[km s−2]	[1028 erg s−1]	[1018 g]	[1032 erg ]
2022 Apr 10 (E1)	Blue emission	330-690	300	300	0.34±0.15	16.6±1.7	${130}_{-90}^{+290}$	${580}_{-400}^{+1280}$	(1)
	(two comp.) a	(490-690)	(390)	(230)	(0.23±0.14)	(9.8±1.3)	(–)	(–)
2022 Apr 16 (E2)	Blue emission	430	380	55	0.12±0.20	0.47±0.16	${1.7}_{-1.1}^{+3.7}$	${12}_{-8}^{+27}$	(1)
2022 Apr 17 (E3)	No	⋯	⋯	⋯	⋯	⋯	⋯	⋯	(1)
2020 Apr 05 (E4)	Blue absorption	510	260	220	0.34±0.04	⋯	${1.1}_{-0.9}^{+4.2}$	${3.5}_{-3.0}^{+14.0}$	(2)
2020 Mar 14 (E5)	No (Red?)	⋯	⋯	⋯	⋯	⋯	⋯	⋯	(3)

Note. ${V}_{\max }$ is the maximum blueshifted velocity ($\max (-V)$). V_typical is the typical blueshifted velocity around when the EW takes a maximum value. V_typical is the typical velocity dispersion (σ) in the fitted Gaussian function. −dV/dt is the deceleration of the blueshifted component. EK Dra's surface gravity is 0.30_±0.05 km s⁻². L_Hα,blue is the luminosity of the blueshifted component. M_p is the mass of the erupted prominences/filament. E_kin is the kinetic energy of the prominence/filament eruptions.

^a The values in parentheses "()" are the values of the two-component fitting of the flare E1. References: (1) This study, (2) Namekata et al. (2022d), (3) Namekata et al. (2022c).

Download table as:  ASCIITypeset image

Table 4. Properties of the Prominence and Filament Eruptions (continued from Table 3) and the Possible X-Ray Dimming

	tWLF	tblue	TimingBlue asym.	X-ray dimming	Timing${}_{{\rm{X}}.\dim .}$	${t}_{{\rm{X}}.\dim .}$	amp${}_{{\rm{X}}.\dim .}$	${\rm{\Delta }}{\mathrm{EM}}_{{\rm{X}}.\dim .}$	${L}_{{\rm{X}}.\dim .}$
	(min)	(min)				(min)	(%)	(1051 cm−3)	(1010 cm)
(E1)	20	34±2	∼start of WLF	Yes?(block5,6)	133 minutes After WLF	>91	5.8-8.5	(–)	(–)
	(–)	(–)	(–)	(bkg subtracted)	(–)	(–)	8.8–14.7	7.2–12.7	0.90–5.0
(E2)	40	10±5	∼peak of WLF	No	⋯	⋯	⋯	⋯	⋯
(E3)	30	⋯	⋯	No	⋯	⋯	⋯	⋯	⋯
(E4)	6.0	58-86	∼end of WLF	(No obs.)	⋯	⋯	⋯	⋯	⋯
(E5)	130	⋯	⋯	(No obs.)	⋯	⋯	⋯	⋯	⋯

Note. t_WLF is the FWHM duration of WLFs as a reference. t_blue is the duration of blueshifted components. For E1, the duration from t = −7.1 minutes to t = 25.0 minutes. For E2, the duration from t = 15 ± 2.5 minutes to t = 25 ± 2.5 minutes. For E3, the lower limit is the duration where the velocity < 0 km s⁻¹, and the upper limit is the duration where the velocity is < −100 km s⁻¹. "Timing_{Blue Asym.}" is the timing when the blueshift components appear relative to WLFs. "Timing${}_{{\rm{X}}.\dim .}$" is the timing when the possible X-ray dimming appears relative to WLFs. ${t}_{{\rm{X}}.\dim .}$ is the duration of the possible X-ray dimming. ${\mathrm{amp}}_{{\rm{X}}.\dim .}$ is the X-ray dimming amplitude for the X-ray count rates (0.5–3 keV). The value in parentheses "()" is the dimming amplitude for the background-subtracted X-ray count rates. ${\rm{\Delta }}{\mathrm{EM}}_{{\rm{X}}.\dim .}$ is the change of emission measure due to the possible X-ray dimming. ${L}_{{\rm{X}}.\dim .}$ is the length scale of the possible X-ray dimming region.

Download table as:  ASCIITypeset image

Figure 11 shows the temporal evolution of the EW, velocity, and velocity dispersion of the fitted parameters for the superflares E1 and E2. For flare E1, both results for one-component and two-component fitting are plotted. The flare E1 shows clear temporal decrease in the blueshifted EW, velocity, and velocity dispersion for both fitting results. The superflare E2 shows a weak decreasing trend of the blueshifted velocity, but no clear evolution for blueshifted EW and velocity dispersion.

Zoom In Zoom Out Reset image size

Figure 12 shows the relation between velocity and velocity dispersion for blueshifted velocity of the flares E1 and E2. As a reference, the filament eruption event on EK Dra reported by Namekata et al. (2022d; "E4" events) is added in Figure 12. The superflare E2 and E4 show clear positive correlation between these two parameters.

Zoom In Zoom Out Reset image size

Based on the temporal evolution of the blueshifted velocity, some basic values of blueshifted components are estimated and summarized in Tables 3 and 4. We estimated the typical blueshifted velocity V_typical and velocity dispersion V_disp,typ as follows: For flare E1, the values are estimated when the EWs of the blueshifted component are most prominent, while for flare E2, the values are estimated as a median value. The deceleration of the blueshifted components −dV/dt (=−A/τ) are obtained by fitting the time evolution of velocity with an exponential function=${{Ae}}^{-(t-{t}_{0})/\tau }$. Also, the duration and timing of the appearance of the blueshifted components are summarized in Table 4.

3.3.3. Estimation of Energy and Timescale

3.4.1. Spectral Fitting for X-Ray Spectra

For the flares E2 and E3, an increase in intensity was detected in NICER's X-ray data, simultaneously with TESS WL and Hα flares. We performed a spectral analysis on the increased intensity during the flares using an emission spectrum from collisionally ionized diffuse gas calculated from the AtomDB atomic database (apec) ²⁷ and the Tuebingen-Boulder ISM absorption model (TBAbs) ²⁸ in xspec. First, we defined the data from the ISS orbit prior to the flares as "pre-flare orbits" (block 2 for flare E2 in Figure 4(c), block 3 for flare E3 in Figure 5(c)) and carried out a spectral fit of the pre-flare spectra. The objective is to remove the pre-flare component from the flare component, thus we utilized as many as four apec models and TBAbs. Three apec models did not adequately reproduce the overall spectral shape. Only the ISM hydrogen column density N_H is fixed as 3 × 10¹⁸ cm⁻² based on the previous study (Güdel et al. 1995) because it is too small to consistently determine for each frame. Linsky et al. (2022) reported a similar N_H value of 1.7 × 10¹⁸ cm⁻². In each case, the N_H value is so small that the factor 2 difference does not significantly affect the fitting results. The fitted spectra are presented in Figures 13(a) and 14(a), while the parameters are summarized in Table 5. The error bars are obtained as a result of the xspec fitting.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 5. Best-fit Values of the Pre-flare NICER X-Ray Spectra for E2 and E3

Component	Flare E2 (block 2)				Flare E3 (block 1)
	T (keV)	(T [107 K])	EM (1052 cm−3)	Abund.	T (keV)	(T [107 K])	EM (1052 cm−3)	Abund.
1	2.06±0.53	(2.39±0.62)	1.58±0.52	0.28±0.04	2.77±0.62	(3.21±0.71)	2.12±0.38	0.38±0.06
2	0.94±0.02	(1.09±0.03)	5.28±0.74	⋯	0.96±0.04	(1.11±0.05)	3.30±0.42	⋯
3	0.40±0.08	(0.47±0.09)	2.63±0.77	⋯	0.50±0.09	(0.58±0.10)	1.94±0.33	⋯
4	0.19±0.38	(0.22±0.44)	0.18±0.62	⋯	0.17±0.03	(0.20±0.04)	0.99±0.15	⋯
χ2(dof)	184(168)				190(177)
reduced χ2	1.10				1.07

Note. The hydrogen column density N_H in TBAbs is fixed as 3 × 10¹⁸ cm⁻².

Download table as:  ASCIITypeset image

We fixed the parameters of the obtained pre-flare spectrum and added a new apec model to perform a spectral fit of the flare component. The flare duration was divided into either two or five segments (named "block"). Initially, a 1-sec-cadence light curve was created in 0.5–3 keV. Based on this, time bins were generated to divide the total count value of the orbit during the flare by half or into fifths. For flare E2 (or E3), the two segmented bins are labeled as blocks 2A–2B (or blocks 3A–3B for flare E3) and the five segmented bins are labeled as blocks 2a–2e (or blocks 3a–3e for flare E3). See Figure 15 for more details of the time bins. A spectral fit was conducted for all these bins. We estimated temperature (T) and emission measure (EM) as free parameters. The coronal abundance was fixed as the pre-flare values because it was not determined consistently when it was set as a free parameter. The results of the spectral fits are summarized in Figures 13 and 14, while the fitting parameters are compiled in Table 6. The time-resolved temperature (T) and emission measure (EM) are plotted in Figure 15. Because the procedure "nicerarf" does not work for the barycentric corrected event data, the above spectral fit was performed for non-barycentric-corrected event data. After the spectral analysis, we corrected the times for a barycentric correction.

Zoom In Zoom Out Reset image size

Table 6. Results of NICER X-Ray Spectral Fit for Flares E2 and E3

Flare E2's block	3	3A	3B	3a	3b	3c	3d	3e	–
T (keV)	1.14±0.01	1.27±0.02	1.01±0.02	1.35±0.03	1.19±0.03	1.20±0.03	1.05±0.03	0.96±0.03	⋯
(T [107 K])	1.32±0.01	1.48±0.02	1.17±0.02	1.57±0.04	1.38±0.03	1.39±0.03	1.22±0.04	1.11±0.03	⋯
EM [1052 cm−3]	8.02±0.14	11.09±0.24	5.48±0.16	13.13±0.40	9.52±0.34	8.03±0.32	5.36±0.27	4.31±0.23	⋯
χ2(dof)	414(220)	250(168)	232(177)	141(126)	166(125)	185(129)	152(127)	112(118)	⋯
reduced χ2	1.88	1.49	1.31	1.12	1.33	1.44	1.19	0.95	⋯

Flare E3's block	2	2A	2B	2a	2b	2c	2d	2e	3
T (keV)	2.67±0.16	3.22±0.28	1.92±0.10	4.50±0.90	3.13±0.38	2.41±0.24	1.44±0.05	5.30±1.42	1.00±0.06
(T [107 K])	3.09±0.18	3.74±0.33	2.22±0.11	5.23±1.04	3.63±0.44	2.80±0.27	1.66±0.06	6.15±1.64	1.16±0.07
EM [1052 cm−3]	6.09±0.12	8.11±0.22	6.92±0.20	5.66±0.29	10.35±0.38	9.90±0.36	7.23±0.35	4.66±0.28	0.70±0.07
χ2(dof)	388(243)	246(184)	278(177)	125(128)	151(127)	145(128)	177(123)	138(120)	279(203)
reduced χ2	1.60	1.34	1.57	0.97	1.19	1.13	1.44	1.15	1.37

Note. The hydrogen column density N_H in TBAbs is fixed as 3 × 10¹⁸ cm⁻². The abundance in apec is fixed as a pre-flare value of each flare in Table 5 (0.28 for flare E2 and 0.38 for flare E3).

Download table as:  ASCIITypeset image

The flare radiation fluxes and energies are estimated in the energy band of 0.5–3 keV (an observed band), 0.1–100 keV (referred to as a "bolometric X-ray" in this study; Kawai et al. 2022), 0.5–7.9 keV (Guarcello et al. 2019), 0.2–12 keV (Kuznetsov & Kolotkov 2021; Stelzer et al. 2022), and 1.55–12.4 keV (=1–8 Å, GOES band, Stelzer et al. 2022) by using flux command in xspec. We calculated the flare energies in these different bands to compare with previous studies. Table 7 shows the peak X-ray flux and total X-ray energy in each band.

Table 7. Summary of NICER X-Ray Flare Parameters

Flare	EMpeak	TEM,peak	kTslope	τrise	τdecay	ξEM,T	LX	EX
	[1052 cm−3]	[107 K]	[keV ksec−1]	(ks)	(ks)		[1029 erg s−1]	[1032 erg ]
				(0.5–3 keV)			(0.5–3 keV)
E2	13.1±0.4	1.57±0.04	-0.26±0.91	(<3.96)	1.95±0.19	0.603	${6.3}_{-0.4}^{+0.2}$	${22.9}_{-1.0}^{+0.5}$
(ext.)							(26.6${}_{-1.7}^{+1.0}$)	(61.9${}_{-2.6}^{+1.4}$)
E3	10.4±0.4	3.63±0.44	-3.32±1.49	0.7±0.2	0.63±0.10	4.28	${4.9}_{-0.2}^{+0.3}$	${6.9}_{-0.1}^{+0.3}$
Flare	LX	EX	LX	EX	LX	EX	LX	EX
	[1029 erg s−1]	[1032 erg ]	[1029 erg s−1]	[1032 erg ]	[1029 erg s−1]	[1032 erg ]	[1029 erg s−1]	[1032 erg ]
	(0.1–100 keV)		(0.5–7.9 keV)		(0.2–12 keV)		(1.55–12.4 keV/GOES)
E2	${8.2}_{-0.5}^{+0.4}$	${29.3}_{-1.2}^{+0.7}$	${5.9}_{-0.3}^{+0.2}$	${21.0}_{-0.7}^{+0.5}$	7.4${}_{-0.5}^{+0.3}$	${26.1}_{-1.1}^{+0.7}$	1.9${}_{-0.1}^{+0.2}$	${5.3}_{-0.3}^{+0.4}$
(ext.)	(34.5${}_{-2.1}^{+1.5}$)	(79.1${}_{-3.4}^{+2.0}$)	(24.8${}_{-1.2}^{+1.0}$)	(56.6${}_{-2.0}^{+1.4}$)	(31.1${}_{-2.1}^{+1.4}$)	(70.6${}_{-2.9}^{+2.0}$)	(8.0${}_{-0.6}^{+0.8}$)	(14.3${}_{-0.8}^{+1.1}$)
E3	${8.1}_{-0.6}^{+0.7}$	${11.4}_{-0.5}^{+0.6}$	${6.2}_{-0.4}^{+0.5}$	${8.6}_{-0.3}^{+0.4}$	7.6${}_{-0.6}^{+0.7}$	${10.5}_{-0.5}^{+0.5}$	4.0${}_{-0.5}^{+0.6}$	${5.4}_{-0.5}^{+0.4}$

Note. "ext." means the extrapolated value by the light-curve fits. EM_peak is the peak emission measure in the observational window. T_EM,peak is the temperature when the emission measure (EM) is maximum value in the observational window. kT_slope is the slope of the temperature (kT) evolution. τ_rise is the rise time of 0.5–3 keV light curve from its onset to peak. τ_decay is the e-fodling decay time of 0.5–3 keV light curve. ξ_EM,T is the power-law index of the relation between temperature (T) and emission measure (EM) in decay phase in the formula of $T\propto {({\mathrm{EM}}^{0.5})}^{{\xi }_{\mathrm{EM},{\rm{T}}}}$. L_X is the X-ray luminosity and E_X is the X-ray energy. The energy and flux values are not corrected for light-curve extrapolations in this table.

Download table as:  ASCIITypeset image

3.4.2. Light-curve Fitting

NICER's orbit did not cover the total flare evolution enough, so we have extrapolated the light curve by fitting the decay phase with exponential function and fitting the rise phase with linear function (see Figure 16). The rise time and decay time are summarized in Table 7. For flare E2, the flare rise phase was not observed, so the rise timescale is just an upper limit. By assuming that the flare peak time is 0 minutes or 10 minutes later than the peak time of the TESS's WLF, we extrapolated the possible peak X-ray flux. As a result of the extrapolation, the X-ray peak of flare E2 is possibly higher by a factor of 4.2–5.7 than the observed peak. The X-ray energy of flare E2 is possibly higher by a factor of 2.7–3.2 than the energy calculated from the observed light curve. In this study, we use the possible X-ray flare peak value of flare E2 as the flux when the flare peak time is 10 minutes later than the peak time of TESS's WLF. For flare E3, we did not perform any extrapolation for the flux and energy because most of the flare evolution was covered by the one NICER orbit.

Zoom In Zoom Out Reset image size

Here we discuss the time evolution of flaring emissions in multiwavelength, particularly focusing on flare E3. Figures 5 reveals that the rise of the soft X-rays appears to be clearly delayed in comparison to the white-light and Hα emissions. We used the method employed by Mitra-Kraev et al. (2005; later, Tristan et al. 2023) to estimate the most likely delay time of soft X-ray flare L_X(t) and time-derivative soft X-ray dL_X(t)/dt (>0). Figure 17 shows the cross-correlation coefficient as a function of time lags. The time lags of maximum of the Pearson's correlation coefficients, ${r}_{\max }$, is defined as the most likely value. The time delay of X-ray against Hα is ${6}_{-2}^{+3}$ minutes, and the time delay of X-ray against WL is ${7}_{-2}^{+3}$ minutes. The error bar is obtained the same way as Tristan et al. (2023). Also, we calculated the time-derivative soft X-ray light curve, as a good proxy of the nonthermal radiation (dL_X/dt, Veronig et al. 2002). The time delay of time-derivative X-ray against Hα is −3${}_{-2}^{+1}$ minutes, and the time delay of differential X-ray against WL is −3${}_{-4}^{+2}$ minutes. The time-derivative soft X-ray light curve (dL_X/dt) shows an increase almost simultaneous with the WL and Hα, with an overall delay of a few minutes for the entire light curve, but peaking at nearly the same time.

Zoom In Zoom Out Reset image size

The close correspondence between the impulsive phase of WL and Hα with the time-derivative X-ray light curve may be an indication of the Neupert effect (Neupert 1968; Tristan et al. 2023). The theoretical Neupert effect is usually expressed as a delay of soft X-rays (thermal emission) against hard X-rays (nonthermal radiation; see Veronig et al. 2002), which is the evidence of the heating by nonthermal electrons followed by the chromospheric evaporation increasing the coronal density of the flaring loops (Shibata & Magara 2011). Because the WL emissions (and some of the Hα emissions) are known to be a good proxy of hard X-ray emission in the Sun (Krucker et al. 2011; Namekata et al. 2017b), the time delay of soft X-ray against WL (and possibly Hα) for flare E3 is interpreted as the "empirical" Neupert effect (e.g., Tristan et al. 2023). This interpretation suggests that flare E3 occurs through the same physical mechanism as solar flares (the same conclusion has been obtained for M-dwarf flares; Hawley et al. 1995; Tristan et al. 2023). Pillitteri et al. (2022) found a clear signature of the "empirical" Neupert effect in 200–300 nm near-UV band and soft X-ray relation for superflares on a young solar-type star, DS Tuc A (age of 40 Myr, spectral type of G6). Pillitteri et al. (2022) and this study have the same conclusion about the flare heating mechanism on young solar-type stars.

Furthermore, Hα and white light potentially appear to lag slightly behind the time-derivative soft X-rays, or seem to have a longer decay. This suggests that in addition to the heating caused by nonthermal electrons corresponding to hard X-rays, heating due to thermal conduction and/or radiative backwarming also contribute to Hα and WL emissions in decay (Machado et al. 1989; Isobe et al. 2007; Namekata et al. 2017b, 2020b, 2022c). These multiwavelength light-curve variations should be compared with the future numerical calculations.

Figure 18(a) shows the relation between Hα and the bolometric soft X-ray flare radiation energy for solar and stellar flares, while Figure 18(b) shows the relation between the bolometric WL and Hα flare radiation energy. We found that the data of the superflares on EK Dra are almost consistent with a trend from solar flares and stellar flares on M-dwarf and RS CVn stars. This suggests that the multiwavelength emission mechanisms of the flares of EK Dra we detected are equivalent to those occurring in flares on other types of stars.

Zoom In Zoom Out Reset image size

In Figure 18(b), the Hα energy of flare E1 is plotted (labeled as "E1") with the upper limit as the total emission energy including the blueshifted spectrum, and the lower limit as the emission energy assuming the presence of a central component. The upper limit (and even possibly the lower limit) is clearly 1–1.5 orders of magnitude larger in Hα emission energy relative to WL emission energy than the other EK Dra flares. This can be attributed to the fact that the emission source of the blueshifted Hα line in flare E1 is different from the other EK Dra flares (not the flare loop or footpoints), suggesting a radiation contribution from stellar prominence eruption, as discussed in Section 5.1.

Incorporating the data from EK Dra (except for flare E1), we derived empirical rules regarding the relationship between Hα and X-ray, as well as WL and Hα. By linking these, we deduced an empirical scaling relation for energy distributions at the bolometric X-ray (0.1–100 keV) and bolometric WL wavelengths, which are releasing a large amount of energy of flares.

$$\begin{eqnarray}&&\displaystyle \frac{{E}_{{\rm{H}}\alpha }}{{10}^{33}\,\mathrm{erg}}={10}^{-1.15\pm 0.08}\cdot {\left(\displaystyle \frac{{E}_{{\rm{X}},\mathrm{bol}}}{{10}^{33}\,\mathrm{erg}}\right)}^{0.90\pm 0.03},\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{E}_{\mathrm{WL},\mathrm{bol}}}{{10}^{33}\,\mathrm{erg}}={10}^{1.68\pm 0.15}\cdot {\left(\displaystyle \frac{{E}_{{\rm{H}}\alpha }}{{10}^{33}\,\mathrm{erg}}\right)}^{1.04\pm 0.07},\end{eqnarray} \tag{ 2 }$$

Equations (1) and (2) are plotted on cyan lines in Figure 18. From Equations (1) and (2), an empirical relation between the X-ray (E_X,bol) and WL (E_WL,bol) flare energy can be deduced as

$$\begin{eqnarray}&&\displaystyle \frac{{E}_{\mathrm{WL},\mathrm{bol}}}{{10}^{33}\,\mathrm{erg}}={10}^{0.48\pm 0.19}\cdot {\left(\displaystyle \frac{{E}_{{\rm{X}},\mathrm{bol}}}{{10}^{33}\,\mathrm{erg}}\right)}^{0.94\pm 0.07}.\end{eqnarray} \tag{ 3 }$$

While this relationship appears to be almost linear, there may be a tendency for X-ray emission energy to slightly dominate as the flare scale increases. This relationship would be useful for statistically investigating the high-energy radiation environment around host stars by connecting to the large statistics of WLFs established by TESS and Kepler (Maehara et al. 2012; Davenport 2016; Feinstein et al. 2019; Notsu et al. 2019; Tu et al. 2020; Okamoto et al. 2021; Tu et al. 2021; Vida et al. 2021).

Furthermore, assuming that the emission energy at different wavelengths can be described by a linear relationship, a corresponding relationship is derived below for reference:

$$\begin{eqnarray}&&{E}_{{\rm{H}}\alpha }={10}^{-1.06\pm 0.09}\cdot {E}_{{\rm{X}},\mathrm{bol}},\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{E}_{\mathrm{WL},\mathrm{bol}}={10}^{1.59\pm 0.06}\cdot {E}_{{\rm{H}}\alpha },\end{eqnarray} \tag{ 5 }$$

From Equations (4) and (5), another empirical relation can be deduced as

$$\begin{eqnarray}&&{E}_{\mathrm{WL},\mathrm{bol}}={10}^{0.54\pm 0.11}\cdot {E}_{{\rm{X}},\mathrm{bol}}\end{eqnarray} \tag{ 6 }$$

Note that the scaling factor in Equation (6) is so different from that in Equation (3) because Equation (6) is a linear relation, whereas Equation (3) is a nonlinear one. In the above, we collected and compared X-ray data in the 0.1–100 keV energy range.

Many other studies derive emission energy in different energy bands, causing inconsistency in comparison (see Section 3.4). In this paper, to compare with other previous studies, we recalculated the flare energy of EK Dra in all energy bands shown in previous studies. The results are summarized in Table 8. As a result, it was found that the white-light-to-X-ray energy ratios (E_WL/E_X) were comparable to those of previous studies in all energy bands. However, in the data compared with solar flares in the GOES band, it was found that solar flares show relatively larger E_WL/E_X ratios (∼70–100) than superflares on EK Dra (∼6.2 and 7.0–22.6) and M-dwarf AD Leonis (∼13). This may suggest, as mentioned above, that in large-scale stellar superflares (10^33-34 erg), X-ray emission energy may become stronger than small-scale solar flares (<10³² erg). In particular, the GOES band is at a higher energy compared with other energy bands in Table 8, and therefore is sensitive to the flare's temperature increase around 10⁷ K. Because it is widely known that the flare temperature increases as the energy scale increases (Tsuboi et al. 1998; Shibata & Yokoyama 1999, 2002; Tsuboi et al. 2016; Getman et al. 2021), the energy in the GOES band may tend to increase more than proportionally.

Lastly, Figure 19 shows the relation between the GOES X-ray peak flux of the flare (F_GOES) and the energy of the white-light flare (E_WLF,Bol). Conventionally, some studies have used an empirical relationship E_WLF,Bol = 10³⁵ · F_GOES (Shibata et al. 2013) to compare the energy scale of solar and stellar flares because most solar flares are evaluated by F_GOES while most stellar flares are by E_WLF,Bol. Later, Namekata et al. (2017b) confirmed observationally that this power-law index becomes linear (∼1) for solar flare observations, but its absolute value is slightly different from Shibata et al. (2013)'s law, and the relationship was E_WLF,Bol = 10^34.3 · F_GOES. Namekata et al. 2017b suggest that this linear relationship can be explained by the well-known Neupert effect, by assuming that WLF is a good proxy of hard X-ray radiation. On the other hand, a nonlinear scaling relation for total solar irradiance (TSI) and GOES X-ray peak flux of solar flares has been reviewed by Cliver et al. (2022), and the relationship ${E}_{\mathrm{TSI}}=0.32\times {10}^{32}\cdot {({F}_{\mathrm{GOES}}/({10}^{-4}{\mathrm{Wm}}^{-2}))}^{0.72}$ has been derived. While all of these relationships are scaling laws describing solar flares, they differ slightly, and a consistent view has not yet been established (see Cliver et al. 2022). Then, which scaling law can describe superflares on stars? Figure 19 shows the comparison between the scaling laws and stellar superflares on G/M-dwarfs, including EK Dra's superflares with red color. As a result, we found that all the scaling laws intersect without significant contradiction with the stellar data. In detail, it can be said for now that the two superflares of EK Dra prefer the scaling laws of Namekata et al. (2017b) and Cliver et al. (2022). In the future, it is expected that multiwavelength observations of superflares with 10^37–39 erg in RS CVn–type binary stars, for example, could lead to an understanding of the extension capability of these scaling laws.

Zoom In Zoom Out Reset image size

The length scales of the coronal flare loops are estimated by using (1) Shibata & Yokoyama (1999, 2002)'s model (see also Namekata et al. 2017a), (2) Reale (2007)'s model, and (3) Namekata et al. (2017b)'s method. See Table 9 for the result summary.

Shibata & Yokoyama (1999, 2002)'s model is based on the two-dimensional magnetohydrodynamic (MHD) simulations of flares performed in Yokoyama & Shibata (1998). They assumed the balance between heating by magnetic reconnection and cooling by thermal conduction in the flare loop and derived the scaling relation to estimate the flare loop length in the following equation

$$\begin{eqnarray}&&{L}_{\mathrm{loop}}={10}^{9}{\left(\displaystyle \frac{{\mathrm{EM}}_{\mathrm{peak}}}{{10}^{48}{\mathrm{cm}}^{-3}}\right)}^{3/5}{\left(\displaystyle \frac{{n}_{{\rm{e}}}}{{10}^{9}{\mathrm{cm}}^{-3}}\right)}^{-2/5}{\left(\displaystyle \frac{{T}_{\mathrm{EM},\mathrm{peak}}}{{10}^{7}{\rm{K}}}\right)}^{-8/5}\mathrm{cm},\end{eqnarray} \tag{ 7 }$$

where EM_peak is the peak EM, and T_EM,peak is the temperature at the EM peak. This reliable method is later validated by solar flare observations (Namekata et al. 2017a). Güdel et al. 1995 reported that the coronal density of EK Dra is >4 × 10¹⁰ cm⁻³, and Ness et al. (2004) also reported the consistent value (e.g., Güdel 2004). Therefore, the assumption of coronal density of ∼10¹¹ cm⁻³ would be a reasonable value for EK Dra's coronal active region. Under the assumption of the pre-flare coronal density in active region n_e=10¹¹ cm⁻³, the loop length of flare E2 and flare E3 are derived to be 1.0 × 10¹¹ cm and 2.1 × 10¹⁰ cm, respectively (see Table 9). The advantage of this approach lies in the fact that information on temporal variation is not required as long as the flare peak is captured. However, as in the case of flare E2, where the flare peak is not captured, it should be noted that accurate values may not be derived (therefore, the value in Table 9 is in brackets).

Reale (2007) applied a one-dimensional hydrodynamic flare model to the rise phase of the flare, suggesting an equation to estimate the flare loop (half-)length:

$$\begin{eqnarray}&&{L}_{\mathrm{loop}}=950\sqrt{{T}_{\mathrm{peak}}}{\tau }_{\mathrm{rise}}\psi \ \mathrm{cm},\end{eqnarray} \tag{ 8 }$$

where T_peak is the maximum temperature and ψ = T_peak/T_EM,peak. This gives the estimate of the loop length of the flare E3 as 0.68 × 10¹⁰ cm. Note that Reale et al. (1997) initially derived the scaling relation to derive the flare loop length from the decay phase of the X-ray flare, but it requires the instrument-calibrated parameters, and no NICER parameters are available for now. But as future development may enable us to estimate the loop length from the decay phase, we derived the required parameters in Table 5.

Namekata et al. (2017b) derived the scaling relation to derive the flare loop length from the WLF's energy and duration. They used the relationship that the WLF's timescale is proportional to the reconnection timescale (∝Alfvén timescale), and flare energy ${E}_{\mathrm{flare}}$ is determined by the magnetic energy. The derived relation is in the formula of

$$\begin{eqnarray}&&{L}_{\mathrm{loop}}\sim 1.64\times {10}^{9}{\left(\displaystyle \frac{{\tau }_{\mathrm{WL}}}{100{\rm{s}}}\right)}^{2/5}{\left(\displaystyle \frac{{E}_{\mathrm{WLF}}}{{10}^{30}\mathrm{erg}}\right)}^{1/5}{\left(\displaystyle \frac{n}{{n}_{\odot }}\right)}^{-1/5}\mathrm{cm},\end{eqnarray} \tag{ 9 }$$

where the coefficient was calibrated by solar coronal loop observation, τ_WL is the e-folding decay time of WLFs, E_WLF is the bolometric WLF energy, and n is the pre-flare coronal number density. Note that the dependence on n is not explicitly included in some previous studies (e.g., Namekata et al. 2022c), but we included it because we need to consider it to discuss the results consistently with other methods (see Shibata & Yokoyama 1999, 2002's model, which assumed the pre-flare coronal density). The normalized factor of n is assumed as solar value n_⊙, which is usually considered as ∼10⁹ cm⁻³ . Under the assumption of n ∼ 100n_⊙ ∼ 10¹¹ cm⁻³ (Güdel et al. 1995), the loop length of flares E1, E2, and E3 are derived to be 0.73 × 10¹⁰ cm, 1.16 × 10¹⁰ cm, and 0.67 × 10¹⁰ cm, respectively (see Table 9).

In summary, the flare loop length estimated from X-ray and WL flares are listed in Table 9. The estimated loop length is roughly ∼0.5 − 2 × 10¹⁰ cm, which corresponds to ∼8 − 30% of stellar radius. We found that the loop lengths derived from the X-ray rise time by Reale (2007)'s methods are consistent with those derived from WLF's energy-duration relationship proposed by Namekata et al. (2017b). The method proposed by Shibata & Yokoyama (1999, 2002) also yields a close value that is approximately twice as large as those from other methods for flare E3 (excluding flare E2, where the flare peak was not captured). Thus, this consistency supports the prediction capability of these different methods, which are widely used but not often simultaneously applied. Particularly, in this campaign, as all flares have been observed by TESS, the estimations from WLFs by Namekata et al. (2017b) appear to be very useful for comparing among flares.

The length of the flare loops derived here will be compared with the length scale of prominences and starspots in Section 6.2. This will enable new discussions on the spatial scale ratios of different aspects of eruptive events. Moreover, the length of the flare loops is important for constraining the initial height of prominence eruptions, and is expected to be useful in future comparisons with numerical calculations in our paper II (see Section 6.4).

As revealed in Section 3, among the three Hα superflares we detected (labeled E1, E2, E3), two superflares (E1, E2) exhibited a significant blueshifted Hα "emission" profile. This is the first discovery of the blueshifted Hα emission profile on a solar-type star (G-dwarf), while similar emission profiles have been frequently reported on M-dwarfs. The profile exhibits a more prominent blueshifted component compared to those observed in M-type stars: In flare E1, the entire spectrum is blueshifted in the initial phase, and in flare E2, a bump-like structure is observed at a wavelength different from the Hα center. This suggests that something different from the usual phenomena seen at the Hα center (i.e., flare ribbons) is occurring and moving in the LOS direction. As mentioned in the introduction, Namekata et al. (2022d) detected a blueshifted "absorption" component moving at a maximum velocity of approximately −510 km s⁻¹ after a superflare (labeled as flare "E4" in Tables 3 and 4), and comparison with solar observations confirmed that this was due to the filament eruption occurring inside the stellar disk. Our new discovery in this study proposed a picture similar to the Sun, where on the solar-type star EK Dra, the Hα line appears as an absorption component when cool plasma are present within the stellar disk (called filament) and as an emission component when they extend beyond the disk (called prominence). Given this context, the blueshifted emission component detected this time could likely be the first evidence of prominence eruption on a solar-type star (i.e., a G-dwarf).

The temporal variation of the velocity strongly supports the idea that the observed blueshifted emission components are stellar prominence eruptions. From the temporal evolution of the shifted velocity, it is possible to derive the initial deceleration rate (discussed in more detail in Section 5.2). For flare E1, the deceleration rate is calculated to be 0.34 ± 0.15 km s⁻² for the one-component fit (or 0.23 ± 0.14 km s⁻² for the two-component fit). This corresponds remarkably well with EK Dra's surface gravity of 0.30 ± 0.05 km s⁻², strengthening the picture that a prominence-like cool plasma is ejected and decelerated by the surface gravity of the star. Interestingly, this value is very similar to those of the filament eruption events E4 of 0.34±0.04 km s⁻² (Namekata et al. 2022d), indicating that a similar phenomenon is occurring except for their radiation signature. As for flare E2, the deceleration rate comes out as 0.12 ± 0.20 km s⁻². Although the error is large due to the time resolution of the spectral fit and the weakness of the signal, the order of magnitude matches the surface gravity.

What are other possible causes for the observed blueshift emission profile? First, the active regions/quiescent prominences moving in LOS because of stellar rotation (see Jardine et al. 2020) can potentially contribute to the blueshift, but given that the stellar rotational velocity (vsini) is 16.4 ± 0.1 km s⁻¹, a possibility of the sudden appearance of active regions/quiescent prominences due to stellar rotation can be rejected. For instance, the co-rotation radius at a velocity of 500 km s⁻¹ is 30 stellar radius. Given the sudden occurrence of blueshift associated with the WLF originating from the stellar surface, events at 30 stellar radii are likely negligible. Second, the presence of a cool plasma above the chromospheric evaporation flows can cause blueshifts as discussed in flare studies from M-type stars (Honda et al. 2018). Very rarely in solar flares, blueshifted emission line profiles of <100 km s⁻¹ are observed in spatially resolved flare ribbon Hα spectra (Švestka et al. 1962; Tei et al. 2018). The solar observations are interpreted by a process where high-energy particles penetrate deep into the chromosphere, causing evaporation, and cool plasma in the upper chromosphere is accelerated before being heated (Tei et al. 2018). However, it is difficult to explain the phenomena we observed in EK Dra (especially E1) with this process. The occurrence of evaporation implies simultaneous chromospheric heating and condensation flows, with chromospheric radiation from denser plasma expected in the Hα core or redshifted components. However, in flare E1, initially almost all Hα emission components are blueshifted, and even when it decelerates and overlaps with the line core, the blueshifted component's EW is a dominant feature, which contradicts the expected time evolution. Therefore, it is difficult to explain the observations (especially E1) caused by evaporated flows.

Furthermore, as mentioned in Section 4.2, for flare E1, the total radiation energy of the Hα line (including the blueshifted component and line center component) is exceptionally large, around ∼30% relative to the WL energy, which is unusual considering the values of other EK Dra flares E2 ∼ E5 (∼1%). This suggests that most of the Hα radiation for flare E1 comes from a source spatially different from the WLF source, i.e., an erupted prominence.

Based on the above discussion, we concluded that we discovered prominence eruptions on a solar-type star for the first time (undoubtedly for flare E1). Like the Sun, our study has made it clear that on solar-type stars, it is possible to determine whether the prominence/filament has erupted beyond/within the stellar visible disk based on whether it appears as an emission/absorption (see summary in Section 6.3). In the following, we will conduct a more detailed analysis of the prominence eruptions observed in this study.

5.2.1. Velocity, Velocity Dispersion, and Their Time Evolution

The observational properties of the blueshifted component are estimated in Sections 3.3.1 and 3.3.2. Here, we review the obtained results in terms of the blushifted velocity, the velocity dispersion, and their time evolution.

• (E1)  
  This flare was the most remarkable, fastest, and long-lived prominence eruption. The ΔEW of the blueshifted component was 0.45–0.70 Å. The duration was 34_±2 minutes, appearing at the start of the white-light flare or possibly earlier (the duration of the white light was 20 minutes). The maximum blueshifted velocity is 330_±35–690_±92 km s⁻¹ for one-component fit (490_±33−690_±93 km s⁻¹ for two-component). No significant acceleration time was observed; instead, it showed an exponential decay of velocity, with an initial deceleration rate of 0.34_±0.15 km s⁻² (already discussed in Section 5.1 in comparison wih surface gravity). The velocity eventually approach to zero velocity, and emission component disappeared without showing significant redshifts. The typical velocity dispersion was 300 km s⁻², exhibiting a pattern of decay synchronous with the velocity.
• (E2)  
  This flare was the weak, slow, short-lived prominence eruption. The ΔEW of the blueshifted emission component was 0.034 Å. The duration was 10_±5 minutes, appearing at the peak of the WLF. The maximum velocity is 430_±19 km s⁻¹. No significant acceleration and deceleration phases were observed, although possible deceleration was estimated as 0.12_±0.20 km s⁻². The typical velocity dispersion was 55 km s⁻².

How can we understand the velocity in the relation of solar flares and blueshift events on other types of stars? Figure 20 shows the velocity of prominence/filament eruptions (blueshift events) as a function of flare's GOES-band flux (bottom axis) or bolometric WLF energy (top axis) for solar and stellar flares. The velocity of prominence/filament eruptions on EK Dra is comparable to the largest values of solar prominence eruptions, while some stellar prominence eruptions from "megaflares" (10^35-36 erg) on RS CVn–type stars/M-dwarfs seems to have a much larger velocity than solar flares (Houdebine et al. 1990; Inoue et al. 2023). It appears that there is a weak positive trend if we include solar and stellar data as

$$\begin{eqnarray}&&V=590\cdot {({F}_{\mathrm{GOES}})}^{0.10\pm 0.02}.\end{eqnarray} \tag{ 10 }$$

Note that this correlation may be influenced by observational bias because it is likely that only those with sufficiently high velocities are being detected for stellar flares (due to spectral resolution, flare contamination, and LOS uncertainty). Therefore, this trend may be an apparent one, and there can be a large diversity in velocity.

Zoom In Zoom Out Reset image size

Because some stellar superflares exhibit prominence eruptions at much higher velocities than the Sun, it can be anticipated that the upper limits of stellar observational data are of significance. Takahashi et al. (2016) suggested that the upper limit of the CME speed can be theoretically described by the scaling law of ${V}_{\mathrm{CME},\mathrm{upper}}={V}_{0}\cdot {({F}_{\mathrm{GOES}})}^{1/6}$, based on the assumption that the flare radiation energy E_WL,bol ∼ 10³⁵ F_GOES should be roughly less than CME kinetic energy $0.5\cdot {M}_{\mathrm{CME}}\cdot {V}_{\mathrm{CME},\mathrm{upper}}^{2}$. Recently, Kotani et al. (2023) suggested that the power-law index for the relation between erupted mass and flare energy scale is the same for CMEs and prominence eruptions, i.e., $M\propto {F}_{\mathrm{GOES}}^{2/3}$; therefore, a similar scaling law for the upper limit of the prominence velociy V_p can be also derived by following Takahashi et al. (2016):

$$\begin{eqnarray}&&{F}_{\mathrm{GOES}}\propto \displaystyle \frac{1}{2}{M}_{{\rm{p}}}{V}_{{\rm{p}},\mathrm{upper}}^{2}\propto {F}_{\mathrm{GOES}}^{2/3}{V}_{{\rm{p}},\mathrm{upper}}^{2},\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{V}_{{\rm{p}},\mathrm{upper}}={V}_{0}\cdot {({F}_{\mathrm{GOES}})}^{1/6}.\end{eqnarray} \tag{ 12 }$$

Where F_GOES can be replaced with 10⁻³⁵ E_WL,bol, V₀ ∼ 5200 km s⁻¹ is roughly determined by observational upper limits for solar flares/filament eruptions. In stellar flares, there have been no observed events exceeding this empirical scaling in velocity (Equation (12)), implying that this established solar scaling law poses no contradictions in describing the upper limits of velocities for filament/prominence eruptions observed in the case of superflares.

Lastly, we would like to focus on a new aspect, the velocity dispersion of the blueshifted components. For flare E1, we discovered a clear decay in velocity dispersion (Figure 11). Furthermore, it was found to have a positive correlation with the absolute value of the velocity (Figure 12). Given that the spatial distribution of thermal and turbulent velocities of prominences are on the order of several tens of km s⁻¹, similar to the Sun (e.g., Sakaue et al. 2018; Namekata et al. 2022d), we cannot explain a velocity dispersion of approximately 300 km s⁻¹ of flares E1 and E4, even considering the velocity dispersion of the instrumental broadenings (∼75 km s⁻¹; Matsubayashi et al. 2019). Actually, this broad velocity dispersion has also been detected in M-dwarf flares (Vida et al. 2016; Leitzinger et al. 2022). Here, by assuming that prominence eruptions are a phenomenon in which loop structures expand self-similarly, we propose that this velocity dispersion may be observing components from spatially different regions. In other words, the slower parts could be near the foot of the loop-like prominence, and the faster parts near the top. Considering this self-similarly expanding loop structure, the fact that velocity dispersion has a positive correlation with the average speed is also reasonable. If this is true, the obtained central velocity of the prominence could be just an average one, with the velocity dispersion indicating the existence of both faster and slower plasma. Therefore, we suggest that the maximum speed of flare E1 could reach up to 630–990 km s⁻¹ (or even up to 880–1080 km s⁻¹), considering the velocity dispersion.

These larger velocities exceed the escape velocity of EK Dra, ∼670 km s⁻¹, and it can be suggested from the Hα line observation alone that flare E1 (and even flare E2) could be connected to a CME. The possibility that these are connected to CMEs will be discussed in further detail combined with another indicator in X-ray in Section 6.1.

5.2.2. Area and Length Scale

Here, we estimate the area and length scale of the prominence. We basically follow the method described by Inoue et al. (2023), but the assumption of optical depth and electron density is original based on the discussion in Section 5.2.4. The luminosity of the blueshift component (L_Hα,blue) at the time of maximum EW is derived as 16.6_±1.7 × 10²⁸ erg s⁻¹ for flare E1 and 0.47_±1.6 × 10²⁸ erg s⁻¹ for flare E2 as shown in Table 3, based on the results of the Gaussian spectral fit. Assuming an optical depth τ = 10 ∼ 1000 for our prominences (see Section 5.2.4 for the assumption), we derive the area, length scale, mass, and kinetic energy based on the theoretical correlations between prominence plasma parameters and the emitted radiation derived by Heinzel et al. (1994). Note that Heinzel et al. (1994)'s work is based on the model of solar quiescent prominence, and the parameters in eruptive promiences on active stars could be different, although further radiative transfer modeling is beyond the scope of this paper. The assumption of the optical depth τ = 1 ∼ 100 gives the theoretical value of surface flux of the prominence F_blue as

$$\begin{eqnarray}{F}_{\mathrm{blue}}=\left\{\begin{array}{ll}{10}^{{\mathrm{log}}_{10}(\tau )+5} & ({\mathrm{log}}_{10}(\tau )\lt -0.22)\\ {10}^{0.55\cdot {\mathrm{log}}_{10}(\tau )+4.9} & ({\mathrm{log}}_{10}(\tau )\gt -0.22)\end{array}\right.\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}=\ \left\{\begin{array}{ll}2.8\times {10}^{5} & (\tau =10)\\ 3.5\times {10}^{6} & (\tau =1000)\end{array}\right.[\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{cm}}^{-2}\,{\mathrm{sr}}^{-1}].\end{eqnarray} \tag{ 14 }$$

Using these values, the projected surface area of the prominence A_p can be expressed as

$$\begin{eqnarray}&&{A}_{{\rm{p}}}={L}_{{\rm{H}}\alpha ,\mathrm{blue}}/(2\pi {F}_{\mathrm{blue}}).\end{eqnarray} \tag{ 15 }$$

This equation gives the surface area of the prominence as 1.8 × 10²²–6.3 × 10²² cm² for flare E1 and 5.3 × 10²⁰–1.9 × 10²¹ cm² for flare E2. Furthermore, the length scale of the prominence L_p is derived as follows, assuming a simple hypothesis that the prominence is roughly square (${A}_{{\rm{p}}}={L}_{{\rm{p}}}^{2}$, the assumption of "cubic" erupted prominences) for the lower limit, and an aspect ratio of the length and width of the prominence of 1:10 for the upper limit (A_p = L_p · L_p/10, the assumption of "cylinder"-like erupted prominences):

$$\begin{eqnarray}&&{L}_{{\rm{p}}}={({A}_{{\rm{p}}})}^{0.5}\sim 10\cdot {({A}_{{\rm{p}}})}^{0.5}.\end{eqnarray} \tag{ 16 }$$

The estimated length scale of the prominence is (8.3 − 30) × 10¹⁰ cm (1.3–4.5 R_star) ∼(83 − 300) × 10¹⁰ cm (13–45 R_star) for flare E1 and (1.4 − 5.1) × 10¹⁰ cm (0.2–0.8 R_star) ∼(14 − 51) × 10¹⁰ cm (2.2–7.8 R_star) for flare E2. Because we used the maximum EW value for calculations, the observed length scale would be the one after the prominence is erupted and expands, so it would be possible that the estimated length scale of prominence is larger than the stellar radius.

How significant is the assumption of aspect ratio? The correction of the prominence aspect ratio with Equation (16) tends to give an extremely large length scale. This may indicate that the prominences could be not necessarily a filamentary structure after the eruptions and could look like halo-type prominence eruptions (see Parenti 2014). Anyway, we do not know the shape of the stellar prominence for now, so to include all the possibilities, we use the values including the aspect ratio correction as an upper limit of the prominence length scale, following the methods in Namekata et al. (2022d). Actually, this assumption does not significantly affect the following discussion of length scale because the lower limit value is important, but it will be revisited when we discuss the general picture of the relationship among flares, prominences, and starspots in Section 6.2.

Figure 21 shows the plot of the prominence length scale and the maximum radial velocity. In the figure, the prominence/filament eruptions of EK Dra are compared with solar filament eruptions (Seki et al. 2021). Seki et al. (2021) proposed a criterion, namely

$$\begin{eqnarray}&&{({V}_{{\rm{r}}\_\max }/100\,\mathrm{km}\ {{\rm{s}}}^{-1})({L}_{{\rm{p}}}/100\,\mathrm{Mm})}^{0.96}=0.8,\end{eqnarray} \tag{ 17 }$$

to determine whether a solar prominence eruption is associated with a CME, where ${V}_{{\rm{r}}\_\max }$ is the maximum radial velocity of filament eruption. In other words, the longer the prominence and the higher the velocity (i.e., the larger the kinetic energy), the higher the likelihood of it leading to a CME. This threshold was first applied to stellar filament eruptions on EK Dra by Namekata et al. (2022d). While this threshold is fundamentally applicable to the Sun, EK Dra has similar surface gravity, surface temperature, and atmospheric stratification as the Sun, suggesting the potential applicability of this threshold. As can be seen in the figure, even the lower limit values of prominence eruptions E1 and E2 are positioned well above this criterion. This supports the possibility that the prominence eruptions associated with flares E1 and E2 may ultimately lead to stellar CMEs. Note, however, that the different activity levels between the Sun and EK Dra could potentially affect relative change in this criterion. Future theoretical investigations are required for this threshold.

Zoom In Zoom Out Reset image size

5.2.3. Mass and Kinetic Energy

Here, under the assumption mentioned above, with an optical depth of τ = 10 ∼ 1000 for the prominences, we determine the prominence mass M_p. The mass can generally be calculated as ${M}_{{\rm{p}}}\sim {A}_{{\rm{p}}}{{Dn}}_{{\rm{H}}}{m}_{{\rm{H}}}$, with the spatial LOS depth of the prominence as D, the number density of hydrogen as n_H, and the mass of hydrogen as m_H. Here, based on the numerical calculations by Heinzel et al. (1994), a relationship has been derived between the optical depth of the prominence and the emission measure EM_blue, as follows:

$$\begin{eqnarray}{\mathrm{EM}}_{\mathrm{blue}}={{Dn}}_{{\rm{e}}}^{2}=\left\{\begin{array}{ll}8.1\times {10}^{29} & (\tau =10)\\ 3.3\times {10}^{32} & (\tau =1000)\end{array}\right.[{\mathrm{cm}}^{-5}]\end{eqnarray} \tag{ 18 }$$

Here, based on the assumption of "cubic" or "cylinder"-like structure introduced in Section 5.2.2, we assume that the depth of the prominence is approximately equal to L_p (cubic) or 0.1 × L_p (cylinder). By using the emission measure EM_blue value under the assumed optical depth, we can describe the M_p as

$$\begin{eqnarray}&&{M}_{{\rm{p}}}\sim {A}_{{\rm{p}}}{{Dn}}_{{\rm{H}}}{m}_{{\rm{H}}}\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&\ =\,{A}_{{\rm{p}}}D{\left(\displaystyle \frac{{\mathrm{EM}}_{\mathrm{blue}}}{D}\right)}^{0.5}{\left(\displaystyle \frac{{n}_{{\rm{e}}}}{{n}_{{\rm{H}}}}\right)}^{-1}{m}_{{\rm{H}}},\end{eqnarray} \tag{ 20 }$$

Here, we roughly assume the prominence ionization fraction n_e/n_H = 0.17 ∼ 0.47 (see Table 1 of Labrosse et al. 2010). Note that this value is modified in Notsu et al. (2023) because the original values used in Inoue et al. (2023) was a mistake. As a result, the mass of prominence is ${1.3}_{-0.9}^{+2.9}\times {10}^{20}$ g (i.e., 4.0 × 10¹⁹ g—4.2 × 10²⁰ g) and ${1.7}_{-1.1}^{+3.7}\times {10}^{18}$ g (i.e., 5.2 × 10¹⁷ g—5.3 × 10¹⁸ g) for flares E1 and E2, respectively. Here, the average of the upper and lower limits in the log scale is denoted as the representative value (the same applies to kinetic energy).

Also, the prominence kinetic energy can be derived from

$$\begin{eqnarray}&&{E}_{\mathrm{kin}}=\displaystyle \frac{1}{2}{M}_{p}{V}_{\mathrm{typical}}^{2},\end{eqnarray} \tag{ 21 }$$

where V_typical is the blueshifted velocity at the maximum Hα EW. As a result, the kinetic energy of prominence for E1 and E2 eruptive events is ${5.8}_{-4.0}^{+12.8}\times {10}^{34}$ erg (i.e., 1.8 × 10³⁴ erg —1.9 × 10³⁵ erg) and ${1.2}_{-0.8}^{+2.7}\times {10}^{33}$ erg (i.e., 3.8 × 10³² erg —3.9 × 10³³ erg), respectively. Thus, the kinetic energy of the prominence is about three times larger for the E1 event (E_WL,bol = 1.5 × 10³³ erg) and about 50 times smaller for the E2 event (E_WL,bol = 12.2 × 10³³ erg) than the radiative energy of the WLF (E_WL,bol) (see Figure 21(b)). This indicates a diversity depending on the flares; flare E1, which has a smaller energy scale for the WLF than flare E2, is larger in all aspects such as mass, speed, and kinetic energy (see the discussion of diversity in Section 6.3).

Figure 22 shows the mass and kinetic energy as a function of flare-radiated energy for flares E1, E2, and E4 in comparison with solar filament eruptions/CMEs (Namekata et al. 2022d; Kotani et al. 2023) and their empirical scaling relations (Aarnio et al. 2012; Drake et al. 2013). Here, we focus comparisons among solar-type stars, the Sun and EK Dra, to show a solar-stellar comparison without many different factors, such as surface gravity, atmospheric density, etc., while other studies have already focused on the comparison among different spectral types. The original figures are taken from Namekata et al. (2022d), and here we added two newly discovered stellar prominence eruptions on EK Dra in green-filled symbols. In addition, small-scale solar filament eruptions/surges investigated by Kotani et al. (2023) were newly added for increasing samples in a small energy regime (10^28-30 erg ).

Zoom In Zoom Out Reset image size

We found that the mass of the gigantic prominence eruption in flare E1 exceeds the predictions of solar empirical scaling relations between flares and CMEs, while the measurements for E2 and E4 are consistent with these relations. Our results suggest that, like solar CMEs/filament eruptions, the ejected mass against a given flare energy scale exhibits considerable scattering in the case of young Sun-like stars, spanning one to two orders of magnitude. Kotani et al. 2023 theoretically suggested that the scattering may be explained by difference in magnetic field strength supporting the prominence. Hence, our findings in diversity in mass (and kinetic energy below) may be related to the diversity in the local magnetic field strength of the star.

Also, the kinetic energy of the E1 event aligns with or exceeds the solar scaling relation between flare energy and CME kinetic energy. In contrast, the kinetic energies of the E2 and E4 events and solar filament eruptions/surges are smaller by one to two orders of magnitude than the scaling relation. The smaller kinetic energies associated with solar-stellar prominence/filament eruptions have been explained in the context of two scenarios, as discussed in Namekata et al. (2022d). The first explanation posits that the velocity of stellar CMEs and prominence/filament eruptions can be suppressed by an overlying robust magnetic field, leading to a smaller allocation of kinetic energy (Alvarado-Gómez et al. 2018). The second explanation attributes it to the different velocity between CMEs and prominence or filament eruptions (Maehara et al. 2021; Namekata et al. 2022d). Solar prominence or filament eruptions are represented by the expanding dense and cool lower portion of the self-similarly expanding flux rope. Consequently, the velocity of prominence or filament eruption is 4–8 times less than that of overlying CMEs (Gopalswamy et al. 2003). This perspective naturally accounts for the smaller kinetic energy observed in solar and stellar prominence or filament eruptions than the solar CME's scaling relation. However, given the significant scattering in solar data, the broad range observed in the kinetic energy of stellar prominence/filament eruptions on EK Dra appears plausible by solar analogy. The cause of this wide variation of solar and stellar prominences/CMEs warrants further investigation in the future.

5.2.4. Optical Depth and Electron Density of Prominences

Then, is the assumption of optical depth τ = 10–1000 appropriate? Table 10 summarizes the estimated prominence length scales for different values on the optical depth τ of 0.1, 1, 10, 100, 1000, and 10,000 as references. Solar prominence eruptions often show optical depth of τ ∼ 1 in Hα line (e.g., Sakaue et al. 2018; Namekata et al. 2022d). The optical depth τ = 0.1 and 1 gives an extremely large value of the length scale of at least >24 R_star and >8.5 R_star, respectively, for flare E1. These values look to be unrealistic considering that the prominence expanding velocity of ∼500 km s⁻¹ can reach up to only 2.4 stellar radius within 1 hr (although plane-of-sky prominence velocities might be larger than the LOS velocities). Therefore, it would be reasonable to assume the optical depth of τ = 10 as a lower limit.

Table 10. Length Scale of Prominences in Different Prominence Parameters for Flares E1 and E2

	Prominence length
	E1		E2
(optical depth)	(Rstar)	(1010 cm)	(Rstar)	(1010 cm)
τ = 0.1	24–240	160–1600	4.1–41	27–270
τ = 1	8.5–85	56–560	1.5–15	9.7–97
τ = 10	4.5–45	30–300	0.78–7.8	5.1–51
τ = 100	2.4–24	16–160	0.41–4.1	2.7–27
τ = 1000	1.3–13	8.3–83	0.22–2.2	1.4–14
τ = 10, 000	0.67–6.7	4.4–44	0.12–1.2	0.77–7.7

Note. The lower value is not an assumption of the prominence's aspect ratio, while the upper value is derived with the assumption of the aspect ratio of the prominence.

Download table as:  ASCIITypeset image

Next, let us discuss the validity of the assumed optical depth from the consistensy in the electron density. Figure 23 shows the estimated electron density as a function of assumed optical depth. The electron densities are estimated from Equation (18). As you can see, the optical depth of 10–1000 gives the electron densities of 1.7 × 10⁹–6.3 × 10¹⁰ cm⁻³. This is actually consistent with the solar prominence observations of 10^9–11 cm⁻³ (Parenti 2014). This consistency in electron density supports the assumed range of the optical depth of this study.

Zoom In Zoom Out Reset image size

Then, what about the electron densities of stellar erupting prominences on active stars? Past X-ray observations indicate that the electron density in the quiescent corona of EK Dra is at least 4 × 10¹⁰ cm⁻³ with a temperature of ∼2 × 10⁶ K (Güdel et al. 1995). This means that the gas pressure in outside corona is ∼40 times larger than in solar ones because solar coronal electron density is ∼10⁹ cm⁻³ with the similar temperature. Considering the gas pressure balance between outside corona and prominences, we can estimate that the prominence density in EK Dra can be ∼40 times larger than the solar values (10^9-11 cm⁻³), i.e., ∼4 × 10^10–12 cm⁻³. The plasma β of solar prominences are roughly ∼1, so this could be the upper limit value of the quiescent prominence electron density. This value (∼4×10^10-12) could be a little bit larger than the estimated prominence electron density (1.7 × 10⁹–6.3 × 10¹⁰ cm⁻³). However, because our prominence parameters are estimated from the Hα luminosity of ∼10 minutes after the prominence initiation for flare E1, the prominence could be more dilute than the original value. Under the assumption of the expansion velocity of 500 km s⁻¹ and the timescale of ∼10 minutes, the prominence could be expanded into ∼3 × 10¹⁰ cm, which is ∼4 times larger than the flare loop size. If we assume that the prominence expands in a cubic way, the volume could expand into ∼64 larger than the original volume, resulting in ∼1/64 of the original electron density, i.e., ∼6 × 10^8–10 cm⁻³. This expanded prominence case in active stars is approximately consistent with the electron density estimated from the optical depth and cubic shape (Figure 23).

In summary, the optical depth of 0.1–1 (∼solar values) is doubtful considering the unrealistically large prominence size. The upper limit of the optical depth could be difficult to be constrained only from the size, but the assumed optical depth of <1000 well covers the realistic prominence electron density. Table 10 includes estimates with much more extreme optical depth τ = 10,000 as a reference. Because Equation (13) shows that the dependence of optical depth on the length scale L_p ∝ (τ)^∼−0.25 is not so strong, the discussion will not significantly change even though we assume τ = 10,000.

We investigated whether the prominence eruptions in flares E1 and E2 are associated with X-ray dimmings, as other evidence of CMEs in multiwavelength (e.g., Veronig et al. 2021). Figures 24(b) and 25 show the ISS-orbital-binned X-ray light curves (see also Appendix A for more details on flare E1). We found that the X-ray count rates in the post-flare phase of flare E1 (i.e., blocks 5 and 6) decrease compared to the possible pre-flare level (block 3) by −5.8_±5.1 % and −8.5_±5.6 % for blocks 5 and 6, respectively (Figure 24(b)). The error bars are calculated by considering those in the pre-flare phase (block 3) and those in corresponding phases. The tiny X-ray enhancement in block 2 would be related to a simultaneous TESS WLF, and in block 4 can include the decay phase in flare E1, so we interpret block 3 (and possibly block 1) as representative of the pre-flare level. To remove the possible contamination of background X-ray flux (from the NICER's nonimaging field of view of 5' diameter ²⁹ ; Gendreau et al. 2016), we also calculated a background-subtracted X-ray count rate in Figure 24(b) using a 3C50 background model (nibackgen3c50; Remillard et al. 2022), which produced the dimming amplitude of −8.8_±6.2 % and −14.7_±7.1 % for each block. On the other hand, Figure 25 shows that there were no significant dimming events for flare E2 and E3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Then, is this X-ray flux decrease a CME-related X-ray dimming? Indeed, the NICER X-ray data are not completely continuous, because NICER is on the ISS, so the full variation of the X-ray light curve is not obtained. This is one of the difficulties of identifying whether it could be a dimming. Also, the net exposure time of this date is relatively short compared to other dates, and the data quality is lower than usual. In addition, the definitions of "pre-flare level" is very difficult for finding XUV dimming in active stars (Veronig et al. 2021; Loyd et al. 2022) because their "pre-flare levels" are essentially a superposition of frequent small (but solar M/X-class) flares and therefore variable. Considering these, it may be not possible to conclusively suggest that an X-ray dimming is detected. However, simultaneous detection with Hα blueshifts in flare E1 supports the possible association of the X-ray dimming, which could indicate the occurrence of unknown CMEs. Therefore, in the following, we call this an X-ray decrease of flare E1 as "possible X-ray dimming" or just a "candidate" and investigate its properties and whether they are consistent with flare and prominence eruption sizes, following the previous studies on M/K-dwarfs (Veronig et al. 2021; Loyd et al. 2022).

For estimating the emission measure of the possible X-ray dimming, we performed the spectral fitting of a background-subtracted X-ray spectrum for each ISS orbit on 2022 April 10 with two-component apec models, which are emission spectra from collisionally ionized diffuse gas calculated from the AtomDB atomic database. ³⁰ The hydrogen column density is fixed as N_H = 3 × 10¹⁸ cm⁻², which was the same assumption made by Güdel et al. (1997a). Abundance relative to solar one is fixed as 0.23 in this analysis, which is a median value for two-component fitting results for all orbits in this campaign. Figure 26 shows the X-ray spectra for each orbit and the fitting results, whose parameters are summarized in Table 11. Figure 24(c) and (d) show the time evolution of the temperature and emission measure of two plasma components.

Zoom In Zoom Out Reset image size

The decrease in total emission measure (EM_tot = EM₁ + EM₂) of blocks 5 and 6 compared to the pre-flare orbit (block 3) is 7.2_±10.7 × 10⁵¹ cm⁻³ (7.8_±11.7 %) and 12.7_±9.2 × 10⁵¹ cm⁻³ (13.8_±10.0 %), respectively. Here, the error bars are a bit large but this is because the EM was obtained by the spectral fitting while the error bars in the light curves were obtained by energy-integrated count rates. This is roughly consistent with the dimming amplitude in their light curve [5.8_±5.1 (block 5)–8.5_±5.6 (block 6) % for non-background-subtracted count rates or 8.8_±6.2–14.7_±7.1 % for background-subtracted count rates]. We can estimate the volume and the length scale of the corona from the equation of ΔEM_tot = n_e n_H V ≈${n}_{{\rm{e}}}^{2}V$ = ${n}_{{\rm{e}}}^{2}{L}_{\dim }^{3}$.

$$\begin{eqnarray}\begin{array}{rcl}V & \approx & {n}_{{\rm{e}}}^{-2}{\rm{\Delta }}{\mathrm{EM}}_{\mathrm{tot}}=({7.2}_{\pm 11.7}-{12.7}_{\pm 9.2})\\ & & \times \,{10}^{31}\cdot {\left(\displaystyle \frac{{n}_{{\rm{e}}}}{{10}^{10}{\mathrm{cm}}^{-3}}\right)}^{-2}{\mathrm{cm}}^{3}\end{array}\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}\begin{array}{rcl}{L}_{\dim } & \approx & {n}_{{\rm{e}}}^{-2/3}{\rm{\Delta }}{\mathrm{EM}}_{\mathrm{tot}}^{1/3}=({4.2}_{-4.2}^{+1.6}-{5.0}_{-1.7}^{+1.0})\\ & & \times \,{10}^{10}\cdot {\left(\displaystyle \frac{{n}_{{\rm{e}}}}{{10}^{10}{\mathrm{cm}}^{-3}}\right)}^{-2/3}\mathrm{cm}\end{array}\end{eqnarray} \tag{ 23 }$$

The mass that escapes can be estimated in the formula of

$$\begin{eqnarray}&&\begin{array}{l}{M}_{\dim }={m}_{{\rm{H}}}{n}_{{\rm{H}}}V\approx {m}_{{\rm{H}}}{n}_{{\rm{e}}}V={m}_{{\rm{H}}}{n}_{{\rm{e}}}^{-1}{\rm{\Delta }}{\mathrm{EM}}_{\mathrm{tot}}\end{array}\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}&&\begin{array}{l}={(1.67\times {10}^{-24}\,{\rm{g}})\cdot ({10}^{10}{\mathrm{cm}}^{-3})}^{-1}\cdot {\left(\displaystyle \frac{{n}_{{\rm{e}}}}{{10}^{10}{\mathrm{cm}}^{-3}}\right)}^{-1}\\ \ \cdot \,(({7.2}_{\pm 11.7}-{12.7}_{\pm 9.2})\times {10}^{51}\,{\mathrm{cm}}^{-3})\end{array}\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}&&=\ ({1.2}_{\pm 2.0}-{2.1}_{\pm 1.5})\times {10}^{18}\cdot {\left(\displaystyle \frac{{n}_{{\rm{e}}}}{{10}^{10}{\mathrm{cm}}^{-3}}\right)}^{-1}\,{\rm{g}}\end{eqnarray} \tag{ 26 }$$

It is reasonable to assume that this dimming region possesses the typical density of the active regions, as numerical simulations by Shiota et al. (2005) show that the dimming region can be considered as the immediate adjacent space to the magnetic field lines undergoing reconnection. Therefore, the assumption of coronal density of 10^10–11 cm⁻³ would be reasonable for EK Dra's active region (see Section 4.3; Güdel et al. 1995; Güdel 2004; Ness et al. 2004). The assumption of a coronal density n_e of 10¹¹ cm⁻³ gives the estimates of the volume V = (7.2_±11.7 − 12.7_±9.2) × 10²⁹ cm³, the length scale ${L}_{\dim }\,=({0.90}_{-0.9}^{+0.3}-{1.1}_{-0.4}^{+0.2})\times {10}^{10}$ cm, and the mass ${M}_{\dim }\,=({1.2}_{\pm 2.0}-{2.1}_{\pm 1.5})\times {10}^{17}$ g. Likewise, the assumption of a coronal density n_e of 10¹⁰ cm⁻³ gives different estimations as V = (7.2_±11.7 − 12.7_±9.2) × 10³¹ cm³, ${L}_{\dim }=({4.2}_{-4.2}^{+1.6}-{5.0}_{-1.7}^{+1.0})\times {10}^{10}$ cm, and ${M}_{\dim }=({1.2}_{\pm 2.0}-{2.1}_{\pm 1.5})\times {10}^{18}$ g.

The length scale of the coronal source of the possible X-ray dimming region ($({0.90}_{-0.9}^{+0.3}-{1.1}_{-0.4}^{+0.2})\times {10}^{10}$ cm or $({4.2}_{-4.2}^{+1.6}-{5.0}_{-1.7}^{+1.0})\times {10}^{10}$ cm) is roughly consistent with—or slightly larger than—the flaring loop length (0.73 × 10¹⁰ cm) estimated in Section 4.3. This supports that the observed post-flare decrease in soft X-ray count rate can be related to the escape of hot coronal plasma associated with a CME. Now, if this X-ray dimming were indeed the so-called dimming, how would it be compared to the properties of the prominence? The length scale of the dimming region is actually much smaller by one to two orders of magnitude than the estimated length scale of prominence eruption ((8.3–300) × 10¹⁰ cm). It is speculated that the prominence gradually expands and its area increases as it erupts, as in the case of solar prominence/filament eruptions. The length scale of the prominence estimated here corresponds to the value when its surface area (i.e., EW) was the largest. Therefore, it is reasonable that estimated prominence length scale is larger than the length scale of the flare region and dimming region. Moreover, the escaped coronal mass estimated from the possible X-ray dimming is smaller than the estimated mass of prominence (4.0 × 10¹⁹ g—4.2 × 10²⁰ g). Therefore, it is not surprising even if an X-ray dimming associated with the gigantic prominence eruption in E1 is detected. These comparisons in length scale and mass, however, largely depend on the assumed physical quantities of stellar corona and prominences (see Section 5.2). More detailed comparisons to ensure that there is no inconsistency between the possible dimming and prominences require further models (see Leitzinger et al. 2022; Loyd et al. 2022).

In Section 5, we report the discovery of gigantic stellar prominence eruptions associated with superflares (E1 and E2). The occurrence of the prominence eruption is conclusively suggested in Section 5.1. The maximum velocity of the prominence eruption in flare E1 was 690 km s⁻¹ (conservatively, 330–490 km s⁻¹), and 430 km s⁻¹ in flare E2 (Section 5.2.1). While the maximum velocity in E1 event is slightly larger than the escape velocity of EK Dra of ∼670 km s⁻¹, it was gradually decelerated and eventually showed almost zero velocity. This may suggest that most of the emitting mass of the prominence stalled, and it is unclear if the prominence itself ultimately escaped. However, we consider that the following points indicate the occurrence of the associated CMEs which consist of outer coronal mass and fast prominence components, especially in flare E1:

• 1.  
  The large velocity dispersion of a blueshifted Hα component can be interpreted as spatial velocity distribution in the expanding 3D loop structure (Section 5.2.1), which is also apparent in the Sun-as-a-star analysis (Appendix B). Based on this interpretation, the maximum velocity component in the E1 event could reach 990–1080 km s⁻¹ (conservatively, 630–880 km s⁻¹), which is significantly greater than the escape velocity. In the case of the Sun, prominence eruptions represent the lower part of CMEs, and some solar events show that the majority of prominence falls back to the Sun, while the other part of fast prominence components and outer coronal masses actually leads to CMEs (e.g., Wood et al. 2016, Appendix B). These suggest that a part of the fast prominence components might have escaped.
• 2.  
  Following the scheme suggested by Namekata et al. (2022d), comparing solar and stellar eruptions in terms of the length scale and velocity strongly suggests that the prominence eruptions in flares E1 and E2 could be associated with CMEs, considering their gigantic length scales and velocities (Section 5.2.2).
• 3.  
  The numerical simulation by Alvarado-Gómez et al. (2018) suggests that the large-scale magnetic field could suppress CMEs under the average magnetic field strength of 75 G if the kinetic energy of CMEs is below a threshold around ∼3 × 10³² erg . Because the EK Dra's averaged magnetic field strength of 66–89 G is similar to this numerical setup (e.g., Waite et al. 2017), we may be able to apply this threshold to our superflares (see Namekata et al. 2022d, Supplementary Information). The initial kinetic energy of the prominence eruptions from events E1 and E2 is ${5.8}_{-4.0}^{+12.8}\times {10}^{34}$ erg and ${1.2}_{-0.8}^{+2.7}\times {10}^{33}$ erg, respectively (Section 5.2.4), which is significantly larger than this threshold. This indicates that the magnetic suppression might not be an efficient process in these cases. Indeed, the mentioned model employs a dipolar magnetic field, while spectropolarimetric observations of this star show that large-scale magnetic fields have multipolar components (e.g., Rosén et al. 2016), and we will perform further numerical modelings with the observed data in the future (Section 6.4).
• 4.  
  After flare E1, a possible candidate of X-ray dimming was detected, as an indication of escaping the surrounding hot coronal plasma associated with CMEs (Section 5.3). Although the determination of the pre-flare X-ray level is an issue in detecting stellar XUV dimming, the simultaneous detection with the Hα blueshift event supports the detection of X-ray dimming. The estimated upper limit of length scale of the dimming region of up to ∼(9–11) × 10⁹ cm is roughly consistent with the length scale of the flares ∼7.3 × 10⁹ cm. The estimated mass of escaping coronal plasma of 1.2 × 10¹⁷–2.1 × 10¹⁸ g is enough small that the erupted prominence mass is 4.0 × 10¹⁹–4.2 × 10²⁰ g. Therefore, it is not surprising that the detected X-ray dimming amplitudes are likely to be interpreted as CME-related dimming (Section 5.3). However, again we should note that the NICER X-ray observation is usually sparse due to the ISS orbit, and we need to be careful with the treatment of these data (see also Appendix A). More convincing detections should be required for further discussions.

Given these without any suggestive evidence that is against the occurrence of CMEs in flare E1, we concluded that the flare E1 is associated with the occurrence a stellar CME. This conclusion is stronger than that drawn in Namekata et al. (2022d). Furthermore, solar studies indicate that velocities of outer CMEs are typically 4–8 times larger than the velocities of prominences (Gopalswamy et al. 2003), suggesting that the associated CME velocity from flare E1 could reach up to 2800–5500 km s⁻¹.

It is important to mention that the center-of-mass velocity of the giant prominence in E1 evolved from 690 km s⁻¹ to near zero in about 20 minutes, and then disappeared without showing a redshift. This means that the center-of-mass of the prominence stopped midway and disappeared for some reason. The disappearance and terminal zero velocity are puzzling because solar prominence eruptions show redshifted spectra after impulsive blueshifts, as presented in Appendix B. By integrating the blueshift velocity over evolution time (up to t < 50 minutes), the line-of-sight displacement of the center-of-mass is calculated as 6.0_±0.3 × 10¹⁰ cm (0.91_±0.04 R_star) for one-component fitting and 7.5_±0.2 × 10¹⁰ cm (1.14_±0.03 R_star) for two-component fitting. The gravity is expected to be 3.6 times weaker at the height of 0.91 R_star than the surface gravity. The timescale of plasma acceleration from rest to > 50 km s⁻¹ from this height is approximately 10 minutes. It is likely that before it began falling back to the star, the plasma was heated or became more diffuse and thus became invisible in the Hα line. Another possible explanation is that in the later phase, the emissions in line center come from flare ribbons. Also, one may speculate that this zero-velocity erupted prominence could be linked to "slingshot" prominences (e.g., Jardine & Collier Cameron 2019). Although there is no observation of slingshot prominences for EK Dra, the E1 event may provide a possiblity that the flare-related prominence could be one source of the slingshot prominences seen in rapidly rotating stars.

The prominence eruption in flare E2 (this study) and filament eruptions in flare E4 (Namekata et al. 2022d) offer less robust evidence of CMEs than flare E1. The prominence eruption associated with the E2 flare event was relatively small and short-lived, and it did not exhibit any substantial X-ray dimming. The lack of significant X-ray dimming following the E2 event could be attributed to its small dimming amplitude, which could be caused by the relatively small prominence area, an order of magnitude smaller than that observed in the E1 flare event. The filament eruptions associated with the E4 flare event (Namekata et al. 2022d) exhibited a clear blueshifted absorption, but its maximum velocity of 510 km s⁻¹ was lower than the escape velocity, and it eventually showed redshifted absorption in Hα at velocities of a few tens of km s⁻¹. However, these low velocities associated with E2 and E4 flare events do not necessarily rule out the presence of CMEs. These values could represent a lower limit of velocity along the LOS. Also, large velocity dispersion in blueshifted components indicates the existence of faster components than the fitting center, which exceeds the escape velocity (see Section 5.2.1). Furthermore, solar observations show that the velocities of CMEs are typically 4–8 times larger than those of prominence eruptions (Gopalswamy et al. 2003), as also discussed in Section 5.2.1. Consequently, the low velocity (less than the escape velocity) blueshift plsama motions in E2 and E4 may suggest that the associated CME velocities could be as high as ∼1700–3400 km s⁻¹ (E2) and 2000–4100 km s⁻¹ (E4), respectively. Therefore, following the solar analogy, we speculate that both E2 and E4 eruptive events are likely associated with fast CMEs, as suggested by Namekata et al. (2022d).

In this section, we will present an overall picture of the eruptive phenomena observed in the young solar-type star EK Dra. This spans from the occurrence of starspots to the manifestation of superflares and prominence eruptions/CMEs, being compared in terms of length scale, energetics, time evolution, and rotational phase.

6.2.1. Length Scale

Table 9 compiles the length scale of starspots L_spot, superflares (coronal loop length) ${L}_{\mathrm{flare}}$ (Section 4.3), and erupting prominences L_p (Section 5.2.2). The size and area of the starspots are calculated from the brightness variation amplitude (BVAmp) of TESS light curves using the methods of Maehara et al. (2017). Table 4 also compiles the length scale of possible dimming regions (${L}_{\dim }$) for the E1 flare event (see Section 5.3). Consequently, the following summarizes the ratios of these length scales to flare loop length:

$$\begin{eqnarray}&&\begin{array}{l}{L}_{\mathrm{spot}}:{L}_{\mathrm{flare}}:{L}_{{\rm{p}}}:{L}_{\dim }\\ =\,\left\{\begin{array}{ll}2.9:1:(11\sim 410):(1.5\sim 6.8), & \mathrm{for}\,({\rm{E}}1)\\ 2.1:1:(1.2\sim 44):(-), & \mathrm{for}\,({\rm{E}}2)\\ 5.0:1:(2.9\sim 260):(-), & \mathrm{for}\,({\rm{E}}4)\end{array}\right.\end{array}\end{eqnarray} \tag{ 27 }$$

Here we use the ${L}_{\mathrm{flare}}$ estimated by the method of Namekata et al. (2017b) since all of these flares have white-light flare observations (see Section 4.3). Equation (27) data provide an initial picture of eruptive phenomena on young solar-type stars: superflares originating from active regions associated with starspots at length scale of at least 2–5 times larger than the flare loop length and producing eruptive prominences/filaments at greater scales (at least 1–10 times the flare loop length). This picture is not unique but looks different depending on events. These relations are useful when comparing observational data with numerical simulations and solar observations. The uncertainty in the length scale of prominences/filaments is currently high as we consider a wide conservative range of optical depths (τ = 10–1000). The prominence length scale would change depending on the phase of the prominence/filament eruption, as it is expected to expand over time, so we assumed a conservative range of optical depth. Also, ionization equilibrium in moving plasma would be different from quiescent prominence (Heinzel & Rompolt 1987), which can significantly change the relation in Equations (13) and (18). Further constraints could be made by combining several spectral lines (e.g., Hα and Hβ) in future observations or radiative transfer modeling (Okada et al. 2020; Leitzinger et al. 2022).

6.2.2. Energy Scales

Our observations provide an initial picture of flare energetics in eruptive superflares from a young solar-type star. The energetics of the eruptive E1, E2, and E4 events can be compiled as follows:

$$\begin{eqnarray}&&\begin{array}{l}{{fE}}_{\mathrm{mag}}:{E}_{\mathrm{WL},\mathrm{bol}}:{E}_{{\rm{X}},\mathrm{bol}}:{E}_{\mathrm{kin}}\\ =\,\left\{\begin{array}{ll}230:1:(-):{39}_{-27}^{+85}, & \mathrm{for}\ ({\rm{E}}1)\\ 44:1:({1.2}_{\pm 0.2}):{0.10}_{-0.07}^{+0.22}, & \mathrm{for}\ ({\rm{E}}2)\\ 180:1:(-):{0.17}_{-0.02}^{+0.70}, & \mathrm{for}\ ({\rm{E}}4)\end{array}\right.\end{array}\end{eqnarray} \tag{ 28 }$$

where fE_mag is the available magnetic energy stored in active regions associated with starspots, f is the fraction of the nonpotential magnetic energy (here, f = 0.1), E_mag is calculated as ${B}^{2}{L}_{\mathrm{spot}}^{3}/8\pi$, and B is assumed as 1000 G. Both the kinetic and radiation energies are much smaller than the stored magnetic energy of the active regions/starspots. This aligns well with our general expectations that giant starspots can produce superflares and prominence/filament eruptions. The kinetic energies of prominence/filament eruptions are much larger than X-rays and WL radiative energies for flare E1, but the relation is opposite for flares E2 and E4. This diversity may require some theoretical explanations, which are beyoud the scope of this paper.

Our flare observations of EK Dra show that the energy partition into X-rays and WL appear to be roughly equal. However, when incorporating data from solar flares and flares from other stars, and deriving the energy distribution rule to white light and X-rays via Hα, we noticed a trend where higher-energy flares tend to be more dominated by X-rays (Section 4.2). The derived Equations (3) or (6) could be useful when estimating the total amount of X-rays irradiated onto the atmospheres of young exoplanets by superflares. Also, we found that the large prominence eruptions of flare E1 emit the Hα flare radiation energy of ∼30% of WLF energy, which is an extremely large portion compared to the usual EK Dra and M-dwarf superflares of 1% (see Section 5.1). Hence, E1 is one of the outliers in the derived flare energy-scaling relationship. This may indicate that the flare statistics based on WLFs with Kepler and TESS can be affected to some extent by the contribution of Hα line emissions from prominence eruptions.

6.2.3. Time Evolution

6.2.4. Rotational Phase

Tables 3, 4, and 9 provide preliminary insights into the diversity and frequency of eruptive flares on young solar-type stars. First, we demonstrate that blueshifted Hα components can appear both in emission and absorption from the young solar-type star EK Dra. This discovery reveals that the well-known signatures of filament (absorption) and prominence (emission) observed on the Sun can also hold true for young solar-type stars. While it might be expected that young solar-type stars emit intense UV radiation, including Lyα, that can alter the excitation states of hydrogen within the prominence (Leitzinger et al. 2022), our observations indicate that the radiative signatures of filaments and prominences do not significantly change even in young, active stars. In other words, the Hα source function, S_λ , of the prominence relative to the background intensity I_λ seems to remain below unity (S_λ /I_λ < 1). This could be a very strong constraint in stellar prominence whose emissivity is unknown.

While the occurrence frequency of prominence/filament eruptions (and CMEs) is a crucial factor for habitability of (exo-)planetary environments, its value is poorly understood and is usually estimated based on flare frequency by assuming a one-to-one relationship between stellar flares and CMEs. For example, it is known that large M/X-class solar flares do not necessarily produce CMEs as they could be supressed by overlying coronal magnetic fields of active regions (Toriumi et al. 2017; Sun et al. 2022). Our limited sample of eruption events shows that the association rate of stellar prominence/filament eruptions against superflares is 3/5 (60%) on young solar-type stars. Blueshifts can be detected with our spectrograph (R ∼ 2000, δ v ∼ 150 km s⁻¹) only when the eruption has relatively high (>150 km s⁻¹) LOS velocity. Considering this limitation on the LOS velocity, we note that the absence of blueshifts in flares E3 and E5 does not necessarily imply the absence of eruptions. Therefore, the expected association rate of prominence/filament eruptions (and possibly CMEs) could be higher than 3/5. This evidence indicates that the strong suppression by local magnetic fields of active regions may not play a significant role in the case of a gigantic eruption associated with the superflare releasing large portions of the magnetic energy stored in an active region (see Section 6.2).

Recent numerical simulations of large-scale photospheric motion driven magnetic flux rope energization suggest an ejection of a CME from a young solar-type star, κ¹ Ceti (Lynch et al. 2019). Such scenarios can be applicable to global multipole magnetic-field environments of young solar-type stars. Indeed, the magnetic fields observed on EK Dra and other young solar-type stars are relatively nonaxisymmetric compared to active M-dwarfs (Donati & Landstreet 2009), and the magnetic fields on solar-type stars are organized on smaller spatial scales than fully convective M-dwarfs (See et al. 2019); therefore, the fall-off with radius is much faster. Thus, we speculate that young solar-type stars can be relatively more CME-productive, at least than M-dwarfs. On the other hand, it is expected that ejection of CMEs from M-dwarfs with a dipolar magnetic field could be suppressed by the strong global dipole magnetic fields of the star (Alvarado-Gómez et al. 2018). This may explain why a recent M-dwarf flare survey shows that M-dwarfs' blueshift association rate is not so high (e.g., seven events among 41 flares on M-dwarfs YZ CMi, EV Lac, and AD Leonis; Notsu et al. 2023), compared to our estimate of expected CMEs from EK Dra.

Presently, all flares E1–E5 exhibit different properties in prominence/filament properties in terms of the radiation signatures, mass, duration, timing, and X-ray counterpart, as summarized in Tables 3, 4, and 9, and also mentioned in Sections 6.1 and 6.2. This indicates that there could be even more diversity if we increase the sample size. We will present and discuss further observational results in our forthcoming papers (see Section 6.4).

Our series of studies, including Namekata et al. (2022c) and Namekata et al. (2022d), has offered numerous observational constraints on the eruptions of prominences and filaments from young solar-type stars. This section outlines our future research plans for refining our understanding of these eruptive phenomena.

First, we plan to determine the maps of starspots and large-scale magnetic fields to address questions such as "Where do prominence and filament eruptions originate?" (refer to Section 6.2). We will model the rotational modulation in the optical light curves from TESS using our previously established starspot-mapping method (Ikuta et al. 2020, 2023). Additionally, at the almost similar period of the campaign periods in 2022 April, we conducted high-dispersion optical spectroscopy of EK Dra using the High Dispersion Echelle Spectrograph (HIDES; Izumiura 1999) on board the 1.88 m reflector at Okayama Observatory in Japan. We have also utilized the high-resolution echelle spectropolarimeter NARVAL, fiber-fed from the Cassegrain focus of the 2 m Bernard Lyot Telescope at the Pic du Midi observatory in France. These data will be analyzed using Doppler imaging and Zeeman Doppler imaging (ZDI) techniques to estimate the surface maps on EK Dra. The resulting maps will be compared with each other, with the occurrence of prominence eruptions (discussed in paper II), and with multiwavelength rotational modulations observed in Hα, X-ray, and ground-based B-band photometric observations (K. Ikuta et al. 2023, in preparation).

We will also perform further modeling of the prominence eruptions from flare E1. Specifically, the Hα data offer several key constraints on the locations and trajectories of the prominence eruptions (see Section 6.2). For instance, an emission Hα profile without any absorption signal implies that a prominence's location is consistently outside the stellar limb from the beginning to the end of the event. The occurrence of a WLF indicates that the eruptive flare's footpoints are not—or at least not entirely—rooted behind the stellar disk. Moreover, the nearly consistent deceleration with surface gravity suggests the eruptions could occur in directions close to LOS. These constraints limit the conditions under which prominence eruptions can occur, thus enabling us to estimate the eruptions' kinematics and evolution. In paper II, we will model these processes using a simple model and reveal the probable inclination angle relative to the LoS and its physical mechanism. By using the model, we plan to apply the radiative transfer to estimate the changes in the surface emission of the prominence as its size/height increases (see Leitzinger et al. 2022). Ultimately, we will carry out 3D numerical simulations of prominence eruptions/CMEs, incorporating the data from surface starspots and ZDI maps (e.g., Lynch et al. 2019; Airapetian et al. 2021; Jin et al. 2022).

Moreover, our observational campaign of two nearby young solar-type stars, including EK Dra, has been in progress since 2020 with the 3.8 m Seimei telescope, 2 m Nayuta telescope, Okayama 1.88 m telescope, and Bulgarian 2 m RCC telescope at Rozhen National Astronomical Observatory. We aim to provide further statistical results of the Hα eruptive phenomena by increasing the superflare sample size to more than 10. This expanded statistical study will further constrain the frequency and diversity discussed in Section 6.3. The likelihood of viewing the eruptive events as prominences or filaments will be modeled and compared with observed diversity with more samples. We will later extend this approach to solar-mass stars of various ages—from weak-line/classical T-tauri stars (see Getman & Feigelson 2021; Getman et al. 2021; Bouvier et al. 2023) to main sequences, e.g., DS Tuc A and κ¹ Ceti—to eventually understand the Sun.