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

Young solar-type stars frequently produce superflares, serving as a unique window into the young Sun-Earth environments. Large solar flares are closely linked to coronal mass ejections (CMEs) associated with filament/prominence eruptions, but observational evidence for stellar superflares remains scarce. Here, we present a 12-day, multiwavelength campaign observation of young solar-type star EK Draconis (G1.5V, 50–120 Myr age) utilizing the Transiting Exoplanet Survey Satellite, the Neutron star Interior Composition ExploreR, and the Seimei telescope. The star has previously exhibited blueshifted Hα absorptions as evidence for a filament eruption associated with a superflare. Our simultaneous optical and X-ray observations identified three superflares of 1.5 × 1033–1.2 × 1034 erg. We report the first discovery of two prominence eruptions on a solar-type star, observed as blueshifted Hα emissions at speeds of 690 and 430 km s−1 and masses of 1.1 × 1019 and 3.2 × 1017 g, respectively. The faster, massive event shows a candidate of post-flare X-ray dimming with the amplitude of up to ∼10%. Several observational aspects consistently point to the occurrence of a fast CME associated with this event. The comparative analysis of the estimated length scales of flare loops, prominences, possible dimming region, and starspots provides the overall picture of the eruptive phenomena. Furthermore, the energy partition of the observed superflares in the optical and X-ray bands is consistent with flares from the Sun, M-dwarfs, and close binaries, yielding the unified empirical relations. These discoveries provide profound implications of the impact of these eruptive events on early Venus, Earth, and Mars and young exoplanets.


INTRODUCTION
Solar and stellar flares are the most powerful explosive phenomena in atmospheres of cool stars observed in the wide range of wavelengths, from radio to X-ray bands (see Shibata & Magara 2011;Benz 2017, for reviews).They are thought to be caused by the conversion of magnetic energy stored in stellar coronae into kinetic and thermal energy via magnetic reconnection.These events on the Sun are frequently accompanied by coronal mass ejections (CMEs) of hundreds of millions of tons of hot coronal material propagating into the interplanetary space at a few hundreds to a few thousands of km s −1 .X-rays and Extreme Ultraviolet (EUV) emission from large solar flares and associated fast CMEs often have severe impacts on the planetary magnetospheres and ionospheres (see Airapetian et al. 2020;Temmer 2021;Cliver et al. 2022, for reviews).
The largest flare energy observed on our Sun is about 10 32 erg (e.g., Emslie et al. 2012), while larger flares with the energy over 10 33 erg referred to as "superflares" have never been reported in the modern solar observations (e.g., Aulanier et al. 2013).However, recent studies have identified multiple possible extreme solar energetic events with energies up to 10 34 ergs have occurred in the past millennia using cosmogenic radionuclides ( 14 C, 10 Be, 36 Cl) recorded in tree rings and ice cores (e.g., Miyake et al. 2019;Usoskin 2023).Furthermore, recent observations of solar-type (G-type main-sequence) stars have provided important insights into the rate of occurrence of superflares from the young solar-type stars and the analogs of the present-day Sun.Vigorous observational studies of stellar flares over the last 30 yrs revealed that young rapidly-rotating solar-type stars show frequent superflares at a rate of 1 event per day (age of ∼100 Myr; Audard et al. 1999Audard et al. , 2000;;Namekata et al. 2022c;Yamashita et al. 2022;Colombo et al. 2022), while mature slowly-rotating solar-type stars show superflares at much lower occurrence rates of about 1 event per a few thousand years (age of several Gyr; Maehara et al. 2012;Shibayama et al. 2013;Notsu et al. 2019;Okamoto et al. 2021).These data suggest that superflares from the young Sun could have had significant impact on the magnetospheric and atmospheric environments of the young Venus, Earth, and Mars, on a time scale of 300 Myr that is controlled the stellar rotation rate (Tu et al. 2015).Also, the detection of superflares on the slowly-rotating stars imply a possibility that superflares can occur on the present-day, mature Sun.Such studies are critical for understanding the frequency and impact of these events on our civilization and its infrastructures.Thus, superflares on solar-type stars have been highlighted from the context of solar physics (e.g., Shibata et al. 2013;Aulanier et al. 2013;Cliver et al. 2022), planetary science (e.g., Fujii et al. 2018;Linsky 2019;Airapetian et al. 2020), biochemistry (e.g., Kobayashi et al. 2023), and historical and geophysical studies (e.g., Miyake et al. 2012;Usoskin 2023;Hayakawa et al. 2017).
High magnetic activity of young solar-type stars is characterized not only by stellar superflares but also by strong average surface magnetic fields (e.g., Vidotto et al. 2014;Kochukhov et al. 2020) and giant starspots covering a substantial portion of the hemisphere (up to tens of percent) (e.g., Waite et al. 2017;Maehara et al. 2017;Namekata et al. 2019Namekata et al. , 2020a;;Notsu et al. 2019;Herbst et al. 2021;Okamoto et al. 2021;Yamashita et al. 2022).The magnetic energy stored in active regions associated with starspots power X-ray and EUV bright (up to 10 30 erg s −1 ), hot (up to 10 MK) and dense (10 10−11 cm −3 ) coronae and massive winds (e.g., Güdel et al. 1995Güdel et al. , 1997a;;Pevtsov et al. 2003;Güdel 2004;Wright et al. 2011;Takasao et al. 2020;Shoda & Takasao 2021;Airapetian et al. 2021;Toriumi & Airapetian 2022;Namekata et al. 2023).Observations of the nearby young solar-type stars in various phases of evolution have provided insights into the time history of our Sun and its heliosphere, which have had a profound influence on the young Earth over the span of millions to billions of years (e.g., Güdel et al. 1997b;Güdel 2007;Johnstone et al. 2019;Okamoto et al. 2021).According to recent geophysical and fossil data, life on Earth is thought to have originated between 4.4 and 3.8 billion years ago (Westall et al. 2023).This time period is known as the "Hadean Eon", during which the Earth would have been subjected to intense XUV radiation and frequent solar flares and CMEs from the young Sun.The XUV radiation from the young Sun could have triggered photoionization processes and contributed to the atmospheric escape of the early Earth (Airapetian et al. 2017;Yoshida & Kuramoto 2021).Moreover, the powerful energetic particles associated with flares/CMEs could have induced chemical reactions, resulting in the formation of greenhouse gases and prebiotic chemistry in the Earth's primitive atmosphere (Airapetian et al. 2016;Kobayashi et al. 2023).These processes may intimately contribute to the formation and disruption of habitable worlds, conditions crucial for the emergence and sustainment of life.As there is little direct evidence of the strong XUV flare radiation and energetic CMEs of the young Sun from the present-day Solar system, we need to infer them from observations of young solar-type stars at ages of a few tens to hundreds of millions of years characteristic of the Hadean period.
Several nearby solar-type stars have been considered as likely proxies for the young Sun (e.g., Güdel et al. 1997b;Güdel 2007;Johnstone et al. 2019).These include κ 1 Ceti (G5V, ∼ 600 Myr old), EK Draconis (EK Dra, G1.5V, ∼ 125 Myr, Waite et al. 2017;S ¸enavcı et al. 2021), and DS Tucanae A (DS Tuc A, G6V, ∼ 45 Myr old).Audard et al. (1999Audard et al. ( , 2000) ) reported on X-ray flare frequencies for κ 1 Ceti and EK Dra, using data from the Extreme Ultraviolet Explorer (EUVE ).They examined the contribution of flares to coronal heating and investigated the relationship between flares and quiescent X-ray emissions.Later, Hamaguchi et al. (2023) reported X-ray superflares on κ 1 Ceti using the Neutron star Interior Composition ExploreR (NICER).Ayres (2015) detailed UV spectroscopic observations of the decay phase of a superflare on EK Dra, captured by the Hubble Space Telescope.They concluded that Si IV and O IV UV emission lines forming at 8×10 5 K redshifted at 30-40 km s −1 are the signatures of the post-flare loop downflows.Guarcello et al. (2019) conducted multi-wavelength observations of a superflare on the active solar-type star HII 345 (G8V), using the Kepler Space Telescope (K2 mission, 30 min cadence) and XMM-Newton.They estimated the energy partition between X-ray and white-light flares and derived the length scale of the flare loop from X-ray data analysis, although its time delay between the peaks in the optical and X-ray bands was not clearly resolved.More recently, Pillitteri et al. (2022) reported two superflares on DS Tuc A, an exoplanet-hosting star, in X-ray and Near-Ultraviolet (NUV) wavelengths using XMM-Newton.The X-ray emission peaks are notably delayed compared to NUV peaks, showing the relationship between heating processes associated with the Neupert effect as observed in solar flare events (Neupert 1968;Dennis & Zarro 1993).Specifically, these delays between peaks inform us about the impulsive sources that transport energetic particle energy flux from the solar/stellar coronae and deposit heating in lower chromospheric layers, which respond back in the form of evaporated flows filling coronal loops with hot and dense X-ray emitting plasma.While these studies have uncovered many critical aspects of flaring atmospheres of magnetically active stars, multi-wavelength observations remain relatively scarce.The Kepler Space Telescope established a robust relationship between rotation (age) and the frequency of white-light flares (WLFs) for solar-type stars (Maehara et al. 2012;Shibayama et al. 2013;Notsu et al. 2019;Okamoto et al. 2021).Application of the extensive Kepler 's statistics to assess impact of flares on exoplanetary environments requires characterizion of the radiative flux from superflares on young solar-type stars in the wide range of wavelengths: from XUV to optical.
Recent optical spectroscopic observations of superflares on young solar-type stars provided indications of stellar CMEs.Namekata et al. (2022d) reported the first detection of the optical Hα spectra of a superflare of ∼ 10 33 erg on the young solar-type star EK Dra (the flare is labeled as "E4" in the following sections).The Hα spectra show a blueshifted absorption as evidence of a gigantic stellar filament eruption (10 4 K plasma) with a velocity of 510 km s −1 and mass of 10 18 g (approximately 10 times larger than the largest solar CMEs).In the standard solar flare model, photospheric motions produce a stressed magnetic flux rope connected to a filament, which becomes unstable and is ejected into the solar corona, initiating a flare through magnetic reconnection (Shibata & Magara 2011;Lynch et al. 2019).These prominence signatures relate to the formation of magnetic flux ropes in the solar and stellar corona, a phenomenon widely researched in solar studies.On the Sun, filament or prominence eruptions (cores of CMEs) are often associated with interplanetary CMEs, the sources of space weather (Cliver et al. 2022).By directly comparing with solar and Sun-as-a-star observations, Namekata et al. (2022d) suggested that a massive high-energy stellar CME is probably associated with the eruption of an extended filament, affecting the young planetary systems and stellar mass loss.On the other hand, Namekata et al. (2022c) discovered another 10 34 erg-class superflare on EK Dra (labeled as "E5") which did not show significant blueshifted Hα spectra, casting a doubt in the one-to-one relationship between filament eruptions/CMEs and superflares.Alternatively, Namekata et al. (2022c) suggested that the superflare could have occurred near the stellar limb and possibly is radiated from flare loops (Heinzel & Shibata 2018), resulting in no blue/redshifts due to the location.The above observational studies demonstrate that time-resolved optical spectroscopic observations are the promising -and currently the only approach to studying CMEs on young solar-type stars, which are also beneficial for understanding other chromospheric heating phenomena (cf.Švestka et al. 1962;Ichimoto & Kurokawa 1984;Canfield & Gayley 1987;Kowalski et al. 2013;Namekata et al. 2020b;Graham et al. 2020;Namekata et al. 2022a;Otsu et al. 2022;Namizaki et al. 2023, for solar and M-dwarf flares).However, to date, only two such observational events have been reported in the literature.To further our understanding of the impact on young planetary environments, we need to significantly increase the number of samples.
The search for signatures of stellar CMEs on other types of stars, such as M/K-dwarfs and active close binaries, has been relatively successful, by using the strategy based on solar CME observations (see, recent reviews by Osten & Wolk 2017;Leitzinger & Odert 2022;Namekata et al. 2022b).This is thanks to the long history of ground-and space-based observations, which is a result of their frequent flares and faint stellar quiescent background (e.g., Pettersen et al. 1984;Pettersen 1989;Hawley & Fisher 1992).Representatively, spectroscopic observations revealed many flare events associated with blueshifted emission line profiles, called "blue asymmetries", in optical Hydrogen Balmer lines (Houdebine et al. 1990;Eason et al. 1992;Gunn et al. 1994;Crespo-Chacón et al. 2006;Fuhrmeister et al. 2008Fuhrmeister et al. , 2011Fuhrmeister et al. , 2018;;Vida et al. 2016Vida et al. , 2019;;Flores Soriano & Strassmeier 2017;Honda et al. 2018;Koller et al. 2020;Muheki et al. 2020a,b;Maehara et al. 2021;Inoue et al. 2023), UV lines (Hawley et al. 2007;Leitzinger et al. 2014), and low-temperature X-ray lines (Argiroffi et al. 2019) during flares on M dwarfs, close binaries, and subgiant/giant stars.These signatures represent the direct evidence of plasma motions in line-of-sight (LOS) direction and the detected number is very large (>100).Although their relatively low velocities of 100-300 km s −1 in most cases make it difficult to distinguish whether they originate from prominence eruptions/CMEs or flare-related surface phenomena, some of them are characterized as prominence eruptions leading to CMEs by the faster velocities than the stellar surface escape velocities (Houdebine et al. 1990;Vida et al. 2016Vida et al. , 2019;;Inoue et al. 2023).Another promising method to indentify stellar CMEs is based on the detection of the XUV dimming observed during solar CMEs forming as a direct result of evacuation and outward propagation of materials from the solar corona, which can last for hours (Mason et al. 2016;Jin et al. 2022).Veronig et al. (2021) reported on possible post-flare X-ray and EUV dimming in cool dwarfs by using archive dataset, and later Loyd et al. (2022) found post-flare dimming in coronal Far Ultraviolet (FUV) lines from ϵ Eridani (K2V).Furthermore, Zic et al. (2020) reported on a possible type-IV radio burst after a flare on Proxima Centauri (M5.5V), although no detections of type-II burst have been reported to date (Crosley & Osten 2018a,b;Villadsen & Hallinan 2019).These eruptive events on cool dwarfs have been in general highlighted as a possible driver of mass and angular momentum loss evolution of the low-mass stars (e.g., Osten & Wolk 2015;Cranmer 2017;Wood et al. 2021) and as a threat to the exoplanet habitability (e.g., Segura et al. 2010;Tilley et al. 2019;Yamashiki et al. 2019;Airapetian et al. 2020;Chen et al. 2021).Although there is still a lack of conclusive observational evidence, these observations provide a hint for potential approaches that could further be extended to the search of eruptions from solar-type stars, as discussed in our study.Given the above context, the aim of this study is threefold: (1) to understand the radiative energy distribution and evolution of flares among white-light, Hα, and X-ray wavelength; (2) to study the diversity of Hα spectra from superflares on young solar-type stars by increasing the sample size; and (3) to provide a comprehensive understanding of stellar CMEs through multi-wavelength observations.To achieve these goals, we conducted the optical and X-ray monitoring observations of the young solar-type star, EK Dra, over a 12-day period from 2022 April 10 to 2022 April 21 (refer to Table 1).This period coincides with the optical photometric observations performed by the Transiting Exoplanet Survey Satellite (TESS, Ricker et al. 2015, Section 2.2).The high-time-cadence optical spectroscopic observations were conducted with 3.8m Seimei telescope in Okayama Observatory, Japan (Section 2.3).The X-ray monitoring observations were also performed by NICER onboard International Space Station (ISS) (Section 2.4).This paper is organized as follows.Section 2 outlines the observational details, while Section 3 describes the analyses and results.In Section 4, we investigate the radiative properties of superflares and their length scales.Section 5 examines the characteristics of prominence eruptions and a potential X-ray dimming.Section 6 integrates the information from Sections 4 and 5 to address our main questions and present our conclusions.In Section 6.4, we outline our future studies, which include modeling prominence eruptions and starspots, as discussed in our forthcoming paper II (Namekata et al. in preparation).

EK Draconis
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 Hadean period on Earth (Waite et al. 2017;S ¸enavcı et al. 2021).It is a rapidly rotating star with the 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 the 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. 1999Audard et al. , 2000;;Ayres 2015;Namekata et al. 2022d,c).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 (cf.Section "Low mass companion" in Supplementary Information of Namekata et al. 2022d).

TESS
The TESS observed EK Dra (TIC 159613900) in its Sector 50 from 26 March to 22 April 2022 1 .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 2 .We used eleanor for photometric analysis to extract light curves from TESS full-frame image (FFI) data (Feinstein et al. 2019).The eleanor yields "raw" flux, "corrected" flux that is regressed against instrumental effects, and "PCA (Principal Component Analysis)" flux that 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.

3.8-m Seimei Telescope
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-8000 Å.The campaign ran from 2022 April 10 to 2022 April 20, with the exposure time set at either 60 or 120 seconds, depending on weather conditions.We followed the data reduction procedure outlined in Namekata et al. (2020bNamekata et al. ( , 2022d,c) ,c) using IRAF 3 and PyRAF 4 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 −1 ; 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 pre-flare subtracted equivalent width as ∆EW = -(EW -EW preflare ), where an emission is positive value.
We processed the datasets 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 six 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-12 keV (approximately 1900 cm −2 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.

Light Curves
Figure 2 overviews the observed multi-wavelength 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 Note-"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.detected in X-ray.Furthermore, for Flare E3, an increase has been detected in Hα, WL, and X-rays altogether.In the following, 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.

TESS White-light Flare Analysis
We quantified the white-light flare (WLF) emission in the TESS -band, following the procedures as taken in Namekata et al. (2022d,c).The detection of flares first requires the removal of background rotational brightness variations.We removed it by using the Fast Fourier Transformation numpy.fft.fft and a low-pass filter with the cut-off frequency at 3.0 d −1 .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 .The same as Figures 3 and 4, but for the superflare observed on 2022 April 17 (E3).The basal level of Hα EW is a median value of -25∼-5 min, while the basal level of soft X-ray is a mean value of block 1.
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.0119for flares E1, E2, and E3, respectively).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 33 erg, flare E2 to be 1.22 ±0.02 × 10 34 erg, and flare E3 to be 3.4 ±0.1 × 10 33 erg; hence, all of these are classified as "superflares".Namekata et al. (2022c) reported the largest energy scale of superflares on EK Dra is ∼ 5 × 10 34 erg, so these are relatively smaller than the upper limit.Also, we derived the full-width at half maximum (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.curvefit.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 of 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 ∼9,000 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 broad-band stellar spectra with the effective temperature of 9,000-10,000 K (e.g., Hawley & Fisher 1992), and thus many stellar flare studies use this temperature as a template.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 3,570-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 4 ) and the increase in emission area (A ∝ T −4 ), both of which are cancelled by each other when using the Equation (1) of Shibayama et al. (2013) (i.e., E = σT 4 A).Therefore, our assumption may provide an uncertainty of the flare energy by a factor of few.We need to perform multi-band observations of superflares in U, B, V and near-UV bands to derive better constraints on the flare energy (Brasseur et al. 2023  (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.

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.The method (1): the first method simply fits the emission component with a single Gaussian function (referred to as "one component" fit).The method (2): the second method first fitted 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 v=-5±18 km/s t=58.9 min 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.leastsquares.The 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 flare 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 min, where method 1 fits successfully.Regarding method 2, we show the fitting results up to t=25 min, 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 −1 at t = -4.8min for both methods.However, the period from t = -9.To improve the signalto-noise ratio for fitting, the data has been binned into 5-minute intervals.The magenta and green lines represent the results of fitting with a two-component Gaussian, and the dashed lines correspond to the sum of the two Gaussians.Please see Figure 7 for the other spectra of flare E2 without significant blueshifted components.
Table 3. Properties of the prominence and filament eruptions on EK Dra.

Hα Asymmetry
Note-Vmax is the maximum blueshifted velocity (max(−V )).V typical is the typical blueshifted velocity around when the EW takes a maximum value.(two component) to include a conservative estimate (where the conservative value is the maximum velocity after t > -2.5 min).
(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 −1 These values are summarized in Tables 3 and 4.
Figure 11 shows the temporal evolution of the EW, velocity, and velocity dispersion of the fitted parameters for the superflares E1 and E2.For the 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.
Figure 12 shows the relation between velocity and velocity dispersion for blueshifted velocity of the flare 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.
Note-tWLF 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 min to t = 25.0 min.For E2, the duration from t = 15±2.5min to t = 25±2.5min.For E3, the lower limit is the duration where the velocity < 0 km s −1 and the upper limit is the duration where the velocity is < -100 km s −1 ."Timing BlueAsym." is the timing when the blue shift components appear relative to WLFs."Timing X.dim." is the timing when the possible X-ray dimming appears relative to WLFs.t X.dim. is the duration of the possible X-ray dimming.amp 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 backgroundsubtracted X-ray count rates.∆EM X.dim. is the change of emission measure due to the possible X-ray dimming.L X.dim. is the length scale of the possible X-ray dimming region.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−t0)/τ .Also, the duration and timing of the apperance of the blueshifted components are summarized in Table 4.

Estimation of Energy and Timescale
The Hα radiated energy is estimated from the light curve of ∆EW whose emission is larger than the one fifths of the peak values.It is calculated by multiplying the enhanced Hα ∆EW by the continuum flux and integrating in time.The continuum flux of EK Dra around Hα is derived as 1.57W m −2 nm −1 at 1 AU (Namekata et al. 2022d) from flux-calibrated EK Dra's spectrum and the stellar distance given by Gaia Data Release 2 (Lindegren et al. 2018).As a result, the Hα flare radiation energy were estimated as 4.9 ±0.2 × 10 32 erg, 1.8 ±0.1 × 10 32 erg, and 2.9 ±0.3 × 10 31 erg for flare E1, E2, and E3, respectively (Table 2).
Note that the above radiation energies of the flares E1 and E2 include blueshifted components.For the flare E1, the radiation energy of the possible central component is estimated as 9.1×10 31 erg, which is only 19 % of the total radiation energy, meaning that most of the radiations come from blueshifted component.Unlike the flare E1, for the flare E2, the radiation energy of the possible central component is estimated as 1.7×10 31 erg, which is 97 % of the total radiation energy, meaning that most of the radiations come from the central component.
The duration of the Hα flares are estimated as the one at FWHM of the smoothed flaring light curves of Hα ∆EW, summarized in Table 2.Note that the most of the radiation of the flare E1 comes from blueshifted component, so the duration of the flaring emission from footpoints or flare loops would be different from the obtained here.

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)7 and The Tuebingen-Boulder ISM absorption model (TBAbs)8 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 18 cm −2 based on the previous study (Güdel et al. 1995) since it is too small to consistently determine for each frame.Linsky et al. (2022) reported a similar N H value of 1.7 × 10 18 cm −2 .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.5 and 6 for the fitting parameters.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, 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 since it was not determined consistently when it was set as a free parameter.The results of the spectral fits are summarized in Figure 13 and 14, while the fitting parameters are compiled in Table 6.The timeresolved temperature (T ) and emission measure (EM ) are plotted in Figure 15.Since 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.
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 "bolometric X-ray" in this study; Kawai et al. 2022), 0.5-7.9keV (Guarcello et al. 2019), 0.2-12 keV (Kuznetsov & Kolotkov 2021;Stelzer et al. 2022), and 1.55-12.4keV (=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.

Light Curve Fitting
NICER's orbit did not enough cover the total flare evolution, 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 For flare E2, two patterns of light curve extrapolation are presented: assuming the peak time is the same as that of the white light (blue), and assuming the peak time is delayed by 10 minutes compared to the white light.The X-ray peak of flare E2 could possibly be higher by a factor of 4.2 (green) and 5.7 (blue) than the observed peak.The X-ray energy of flare E2 could possibly be higher by a factor of 2.7 (green) and 3.2 (blue) than the observed energy.6.9 +0.3  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 ∝ (EM 0.5 ) ξ EM,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.
just an upper limit.By assuming that the flare peak time is 0 min or 10 min 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.2than 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 min later than the peak time of the 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's orbit.Here we discuss the time evolution of flaring emissions in multi-wavelength, particularly focusing on flare E3. Figure 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 +3 −2 min, and the time delay of X-ray against WL is 7 +3 −2 min.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 non-thermal radiation (dL X /dt, Veronig et al. 2002).The time delay of time-derivative X-ray against Hα is -3 +1 −2 min, and the time delay of differential X-ray against WL is -3 +2 −4 min.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.
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 (non-thermal radiation) (cf.Veronig et al. 2002), which is the evidence of the heating by non-thermal electrons followed by the chromospheric evaporation increasing the coronal density of the flaring loops (Shibata & Magara 2011).Since 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" The time lags of maximum correlation coefficients, rmax, is defined as the most likely value.The 90% percentile gives the lag time uncertainty.The time delay of X-ray against Hα is 6 +3 −2 min, the time delay of X-ray against WL is 7 +3 −2 min, the time delay of differential X-ray against Hα is -3 +1 −2 min, and the time delay of differential X-ray against WL is -3 +2 −4 min.See Tristan et al. (2023) for the method.
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 non-thermal 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. 2017bNamekata et al. , 2020bNamekata et al. , 2022c)).These multi-wavelength light curve variations should be compared with the future numerical calculations.

Flare Radiative Energy Partition
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 that 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 multi-wavelength emission mechanisms of the flares of EK Dra we detected are equivalent to those occurring in flares on other types of stars.
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.Equations 1 and 2 are plotted in 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 E WL,bol 10 33 erg = 10 0.48±0.19• E X,bol 10 33 erg 0.94±0.07

E
. (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;Notsu et al. 2019;Feinstein et al. 2019;Okamoto et al. 2021;Tu et al. 2020Tu 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: From Equations 4 and 5, another empirical relation can be deduced as Note that the scaling factor in Equation 6 is so different from that in Equation 3 because Equation 6 is a linear relation while Equation 3 is a non-linear 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 an M-dwarf AD Leo (∼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 32 erg).In particular, the GOES band is located at a higher energy compared to other energy bands compared in Table 8, and therefore is sensitive to the flare's temperature increase around 10 7 K.Since it is widely known that the flare temperature increases as the energy scale increases (Tsuboi et al. 1998(Tsuboi et al. , 2016;;Shibata & Yokoyama 1999, 2002;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 35 •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 well-known Neupert effect, by assuming that WLF is a good proxy of hard X-ray radiation.On the other hand, a non-linear 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 TSI = 0.32×10 32 •(F GOES /(10 −4 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 (cf.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 The Sun (G2V) (5), ( 6) Note-E WL,bol is the bolometric white-light flare energy under the assumption of 10,000 K blackbody radiation.E X is the X-ray flare energy in the energy band of the first culmn."Sp.Type" is the spectral type of the star.* The Kepler-band flare energy is reported.§ The TESS-band flare energy is reported.The authors also mentioned that the E WL,bol /E X values are lower limits because all of the X-ray flare evolution were not covered by NICER.References: (1) Guarcello et al. ( 2019  Note-"BVAmp" is the brightness variation amplitude relative to average flux in TESS band due to the stellar rotation."Rot.Phase" is the timing when the superflares occur relative to stellar rotational phase."Loc.min" and "Loc.max" is the local minimum and local maximum of the light curve, respectively, and "Rise" and "Decline" is the rising and decline phase of the light curve, respectively.As a reference, stellar radius is 0.94 times solar radius ∼6.55×10 10 cm, and stellar disk area is ∼13.5×10 21cm 2 .§ The loop length is estimated with the method by Shibata & Yokoyama (1999, 2002) under the assumption of coronal electron density ne = 10 11 cm −3 .% The loop length is estimated with the method by Namekata et al. (2017b) under the assumption of coronal density n = 10 2 n ⊙ ∼ 10 11 cm −3 .ß The estimation from rise timescale of Reale (2007).† The data is taken from Namekata et al. (2022d).* The data is taken from Namekata et al. (2022c).
scaling laws of Namekata et al. (2017b) and Cliver et al. (2022).In the future, it is expected that multi-wavelength observations of superflares with 10 37−39 ergs in RS CVn type binary stars, for example, could lead to an understanding of the extension capability of these scaling laws.

Length Scale of Flare Loop
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 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 assume 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 where EM peak is the peak EM, 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 10 cm −3 , and Ness et al. ( 2004) also reported the consistent value (e.g., Güdel 2004).Therefore, the assumption of coronal density of ∼ 10 11 cm −3 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 11 cm −3 , the loop length of flare E2 and flare E3 are derived to be 1.0×10 11 cm and 2.1×10 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 bracket).Reale (2007) applied their one-dimensional hydrodynamic flare model to the rise phase of the flare, suggesting an equation to estimate the flare loop (half-)length: 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 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's parameters are available for now.But 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 WLFs timescale is proportional to reconnection timescale (∝Alfven timescale) and flare energy E flare is determined by the magnetic energy.The derived relation is in the formula of 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 (cf., Shibata & Yokoyama 1999, 2002's model assumed the pre-flare coronal density).The normalized factor of n is assumed as solar value n ⊙ , which is usually considered as ∼ 10 9 cm −3 .Under the assumption of n ∼ 100n ⊙ ∼ 10 11 cm −3 (Güdel et al. 1995), the loop length of flare E1, E2, and flare E3 are derived to be 0.73×10 10 cm, 1.16×10 10 cm, and 0.67×10 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 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 WLFs energy-duration relationship proposed by Namekata et al. (2017b).The method proposed by Shibata & Yokoyama (1999, 2002) also yields a close value that are approximately twice as large as those from other methods for flare E3 (excluding flare E2 where the flare peak was not captured).Thus, these 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 to compare between flares The length of the flare loops derived here will be compared with the length scale of prominences and starspots in the following 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).

Discovery of Blueshifted "Emission" Spectra as Evidence of Stellar Prominence Eruptions
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 solar-type star (G dwarf), while similar emission profile has 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 −1 after a superflare (labeled as flare "E4" in Table 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 a 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 that the observed blueshifted emission component 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.15km s −2 for the one component fit (or 0.23±0.14km s −2 for the two component fit).This corresponds remarkably well with EK Dra's surface gravity of 0.30±0.05km s −2 , strengthening the picture that a prominence-like cool plasma is the 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.04km s −2 (Namekata et al. 2022d), indicating the 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 −2 .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 due to stellar rotation (cf.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 −1 , 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 −1 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 −1 are observed in spatially resolved flare ribbon Hα spectra (Tei et al. 2018;Švestka et al. 1962).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 the 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 Figure 20.Velocity of solar and stellar prominence/filament eruptions (blueshift events) as a function of flare GOES flux (bottom axis) or bolometric WLF energy (top axis).Data of filament/prominence eruptions on the G-dwarf EK Dra are derived from this study and Namekata et al. (2022d).Data of solar filament eruptions are taken from Seki et al. (2021).Data for M-dwarfs, RS CVn, and YSO are taken from Moschou et al. (2019), Inoue et al. (2023), andMaehara et al. (2021).The cyan solid line is a power-law fitting line for all the solar and stellar data.The black dashed line is a line with a power-law index of 1/6 which represents the upper edge of the solar and stellar data Takahashi et al. (2016).
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.

Velocity, Velocity Dispersion, and Their Time Evolution
The observational properties of the blueshifted component are estimated in Section 3.3.1 and 3.3.2.Here, we review the obtained results in terms of the blushifted velocity, 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 min, appearing at the start of the white-light flare or possibly earlier (the duration of the white light was 20 min).The maximum blueshifted velocity is 330 ±35 -690 ±92 km s −1 for one component fit (490 ±33 -690 ±93 km s −1 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 −2 (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 −2 , exhibiting a pattern of decay synchronous with the velocity.
(E2) This flare was the weak, slowest, short-lived, prominence eruption.The ∆EW of the blueshifted emission component was 0.034 Å.The duration was 10 ±5 min, appearing at the peak of the WLF.The maximum velocity is 430 ±19 km s −1 No significant acceleration and deceleration phase was observed, although possible deceleration was estimated as 0.12 ±0.20 km s −2 .The typical velocity dispersion was 55 km s −2 .
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 looks that there is a weak positive trends if we include solar and stellar data as 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.
Given the fact that 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 is 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 CME,upper = V 0 •(F GOES ) 1/6 , based on the assumption that the flare radiation energy E WL,bol ∼ 10 35 F GOES should be roughly less than CME kinetic energy 0.5 Kotani et al. (2023) suggested 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 ∝ F 2/3 GOES , therefore a similar scaling law for the upper limit of the prominence velociy V p can be also derived by following the Takahashi et al. (2016).
where F GOES can be replaced with 10 −35 E WL,bol , V 0 ∼5200 km s −1 is roughly determined by observational upper limit 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 −1 , 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 −1 of flare E1 and E4, even considering the velocity dispersion of the instrumental broadenings (∼75 km s −1 ; 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 −1 (or even up to 880-1080 km s −1 ), considering the velocity dispersion.
These larger velocities exceeds the escape velocity of EK Dra ∼670 km s −1 , 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 other indicator in X-ray in Section 6.1.

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 28 erg s −1 for flare E1 and 0.47 ±1.6 × 10 28 erg s −1 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 The green and red points corresponds to the stellar prominence and filament eruptions on EK Dra, respectively (this study, Namekata et al. 2022d).The stellar data are the maximum line-of-sight velocities while the solar data are maximum radial velocity.The solid black indicates an empirical solar threshold that can roughly distinguish whether a filament eruption is associated with CMEs or not in the formula of (Vr max/100 km s −1 )(L/100 Mm) 0.96 = 0.8 (Seki et al. 2021).
optical depth τ = 1 ∼ 100 gives the theoretical value of surface flux of the prominence F blue as F blue = 10 log 10 (τ )+5 (log 10 (τ ) < −0.22) 10 0.55•log 10 (τ )+4.9 (log 10 (τ ) > −0.22) Using these values, the projected surface area of the prominence A p can be expressed as This equation gives the surface area of the prominence as 1.8×10 22 -6.3×10 22 cm 2 for flare E1 and 5.3×10 20 -1.9×10 21 cm 2 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 p = L 2 p , 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): The estimated length scale of the prominence is (8.3 − 30) × 10 10 cm (1.3-4.5 R star ) ∼ (83 − 300) × 10 10 cm (13-45 R star ) for flare E1 and (1.4 − 5.1) × 10 10 cm (0.2-0.8 R star ) ∼ (14 − 51) × 10 10 cm (2.2-7.8R star ) for flare E2.Since 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 the assumption of aspect ratio is?The correction of the prominence aspect ratio with Equation 16tend to give extremely large length scale.This may indicate that the prominences could be not necessarily a filamentary structure after the eruptions and could look like a halo-type prominence eruptions (cf.Parenti 2014).Anyway, we do not know at all 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 will be revisited when we discuss the general picture of the relationship between 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 (V r max /100 kms −1 )(L p /100 Mm) 0.96 = 0.8, ( to determine whether a solar prominence eruption is associated with a CME, where V 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.

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 p ∼ A p Dn H m H , where 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: 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 Here we roughly assume the prominence ionization fraction n e /n H = 0.17 ∼ 0.47 (cf.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 prominences are 1.3 +2.9 −0.9 × 10 20 g (i.e., 4.0 × 10 19 g -4.2 × 10 20 g) and 1.7 +3.7 −1.1 × 10 18 g (i.e., 5.2 × 10 17 g -5.3 × 10 18 g) for flare 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 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 event are 5.8 +12.8 −4.0 × 10 34 erg (i.e., 1.8 × 10 34 erg -1.9 × 10 35 erg) and 1.2 +2.7 −0.8 × 10 33 erg (i.e., 3.8 × 10 32 erg -3.9 × 10 33 erg), respectively.Thus, the kinetic energy of the prominence is about three times larger for the E1 event (E WL,bol = 1.5 × 10 33 erg) and about fifty times smaller for the E2 event (E WL,bol = 12.2 × 10 33 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 on diversity in Section 6.3).
Figure 22 shows the mass and kinetic energy as a function of flare radiated energy for flare E1, E2, and E4 in comparison with solar filament eruptions/CMEs (Namekata et al. 2022d;Kotani et al. 2023) and their empirical Figure 22.Mass and kinetic energy of solar/stellar filament/prominence eruptions and CMEs as a function of bolometric WLF energy or GOES X-ray energy.The green and red points correspond to the stellar prominence and filament eruptions on EK Dra, respectively (this study, Namekata et al. 2022d).The black crosses are solar CME data (Yashiro & Gopalswamy 2009;Drake et al. 2013) while the green pluses are data for solar prominence/filament eruptions taken from Namekata et al. (2022d) and Kotani et al. (2023).The cyan dashed and magenta dotted lines are observational fits for solar CMEs expressed as MCME ∝ E 0.59 bol and E kin ∝ E 1.05 bol (cyan color; Drake et al. 2013) and MCME ∝ E 0.7 bol (magenta color; Aarnio et al. 2012).The relations of E kin = E bol , E kin = E bol /10, and E kin = E bol /100 are also plotted with gray color in the right panel.Here, GOES X-ray energy/flux is related to bolometric energy with the relations of E bol = 100 EGOES (Emslie et al. 2012) and E bol [erg] = 10 35 FGOES [W m −2 ] ( Shibata et al. 2013) to connect solar data to stellar data.Also, we assume bolometric energy ≈ bolometric WLF energy (Kretzschmar et al. 2010;Kretzschmar 2011;Osten & Wolk 2015).
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 much different factors such as surface gravity, atmospheric density, etc, while other studies has 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 small energy regime (10 28−30 erg).
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 local magnetic field strength of the star.
Also, the kinetic energy of E1 event aligns with or exceeds the solar flare energy-CME kinetic energy scaling relation.In contrast, the kinetic energies of E2, 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 expanding dense and cool lower portion of the self-similarly expanding flux rope.Consequently, the velocity of prominence or filament Note-The lower value is no assumption of the prominence's aspect ratio while the upper value is derived with the assumption of the aspect ratio of the prominence.The left y-axis is the case when the prominence is cubic structure, while right y-axis is the cause when the prominence is cylinder (i.e., filamentary structure).The blue colored region corresponds to the assumed optical depth.
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 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.

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 10000 as references.Solar prominence eruptions often shows 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 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 −1 can reach up to only 2.4 stellar radius within 1 hour (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.
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 × 109 -6.3 × 10 10 cm −3 .This is the actually consistent with the solar prominence observations of 10 9−11 cm −3 (Parenti 2014).This consistency in electron density supports the assumed range of the optical depth of this study.
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 10 cm −3 with a temperature of ∼ 2 × 10 6 K (Güdel et al. 1995).This means that the gas pressure in outside corona is ∼40 times larger than solar ones because solar coronal electron density is ∼10 9 cm −3 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 −3 ), i.e., ∼ 4 × 10 10−12 cm −3 .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 9 -6.3 × 10 10 cm −3 ).However, since 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 −1 and the timescale of ∼10 minutes, the prominence could be expanded into ∼ 3 × 10 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 −3 .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) have been 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 cover the realistic prominence electron density.Table 10 includes estimates with much more extreme optical depth τ =10000 as a reference.Since Equation 13shows 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 τ =10000.

Candidate of Possible X-ray Dimming after Flares and Prominence Eruptions
We investigated whether the prominence eruptions in flares E1 and E2 are associated with X-ray dimmings, as another evidence of CMEs in multi-wavelength (e.g., Veronig et al. 2021).Figures 24(b) and 25 show the ISS-orbitalbinned X-ray light curves (see also Appendix A for more details of flare E1).We found that the X-ray count rates in post flare phase of flare E1 (i.e., block 5 and 6) decrease compared to the possible pre-flare level (block 3) by -5.8 ±5.1 % and -8.5 ±5.6 % for block 5 and 6, respectively (Figures 24(b)).The error bars are calculated by considering those in 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 's WLF and that in block 4 can include the decay phase in flare E1, so we interpret block 3 (and possibly block 1) is representative of the pre-flare level.To remove the possible contamination of background X-ray flux (from the NICER's non-imaging field-of-view of 5 arcmin diameter 9 ; Gendreau et al. 2016), we also calculated a background-subtracted X-ray count rates in Figure 24(b) using 3C50 background model (nibackgen3c50; Remillard et al. 2022), which produce the dimming amplitude of -8.8 ±6.2 % and -14.7 ±7.1 % for each block.On the other hands, Figure 25 shows that there were no significant dimming events for flare E2 and E3.
Then, is this X-ray flux decrease a CME-related X-ray dimming?Indeed, the NICER's X-ray data is not completely continuous one 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 to identify whether it could be a dimming.Also, the net exposure time of this date is relatively short compared to other dates, the data quality is also lower than usual.In addition, the definitions of "per-flare level" is very difficult for finding XUV dimming in active stars (Veronig et al. 2021;Loyd et al. 2022) since their "pre-flare level" 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 is supportive of the possible association of the X-ray dimming, that could be indicative of the occurrence of unknown CMEs.Therefore, in the following, we call this an X-ray decrease of flare E1 as "possible     Note-The hydrogen column density N H in TBAbs is fixed as 3 × 10 18 cm −2 .The abundance in apec is fixed as typical value of 0.23.
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 database10 .The hydrogen column density is fixed as N H = 3×10 18 cm −2 , which was the same assumption as Güdel et al. (1997a).Abundance relative to solar one is fixed as 0.23 in this analysis, which is a median values for two component fitting results for all orbit 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 component.
The decrease in total emission measure (EM tot = EM 1 + EM 2 ) of block 5 and 6 compared to the pre-flare orbit (block 3) is 7.2 ±10.7 ×10 51 cm −3 (7.8 ±11.7 %) and 12.7 ±9.2 ×10 51 cm −3 (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 rate.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 2 e V = n 2 e L 3 dim .
V ≈ n −2 e ∆EM tot = (7.The estimated mass that escape can be estimated in the formula of = (1.67 × 10 −24 g) • (10 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 −3 would be a reasonable assumption for EK Dra's active region (see, Section 4.3; Güdel et al. 1995;Ness et al. 2004;Güdel 2004).
The assumption of a coronal density n e of 10 11 cm −3 gives the estimates of the volume V = (7.2±11.7 − 12.7 ±9.2 ) × 10 29 cm 3 , the length scale L dim = (0.90 +0.3 −0.9 −1.1 +0.2 −0.4 )×10 10 cm, and the mass M dim = (1.2±2.0 −2.1 ±1.5 )×10 17 g.Likewise, the assumption of a coronal density n e of 10 10 cm −3 gives different estimations as V = (7.2±11.7 − 12.7 ±9.2 ) × 10 31 cm 3 , L dim = (4.2+1.6 −4.2 − 5.0 +1.0 −1.7 ) × 10 10 cm, and M dim = (1.2±2.0 − 2.1 ±1.5 ) × 10 18 g.The length scale of the coronal source of the possible X-ray dimming region ((0.90 +0.3  −0.9 − 1.1 +0.2 −0.4 ) × 10 10 cm or (4.2 +1.6  −4.2 − 5.0 +1.0 −1.7 ) × 10 10 cm) is roughly consistent with -or slightly larger than -the flaring loop length (0.73 × 10 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 compare 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 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 enough smaller than the estimated mass of prominence (4.0 × 10 19 g -4.2 × 10 20 g).Therefore, it is not surprising even if an X-ray dimming associated with the gigantic prominence eruption in E1 is detected.These comparison in length scale and mass, however, largely depends on the assumed physical quantities of stellar corona and prominences (cf.Section 5.2).More detailed comparisons to ensure that there is no inconsistency between the possible dimming and prominences requires further models (cf., Loyd et al. 2022;Leitzinger et al. 2022).

DISCUSSION AND CONCLUSION
We report on three superflares (labeled as E1-E3) from the young solar-type star, EK Dra, using our simultaneous optical and X-ray observations taken with Seimei, TESS, and NICER.The time-resolved Hα spectra were characterized in the relation with simultaneous WL and X-ray flares for the first time on solar-type stars.A short summary of the preceding sections is as follows: Section 4 describes the multi-wavelength time evolution, and energy partitions, and length scales of superflares.Section 4.1 suggests that occurrence mechanism of superflares on solar-type stars could be common as those of solar and M-dwarf flares, as time delay of X-ray against WL and Hα suggest the chromospheric heating followed by chromospheric evaporations.Section 4.2 provide meaningful empirical laws for flare radiative energy partition into X-ray and WL, linking the well established Kepler /TESS photometry to assessing their impact on exoplanetary environment.Section 5 reports on the first discovery and properties of gigantic prominence eruptions (E1 and E2) on a solar-type star.Section 5.1 gives a conclusion that the detected blueshifted emissions are stellar prominence eruptions from several observational considerations.Section 5.3 reports and discusses a possible candidate of the associated post-flare X-ray dimming (E1), which enables the first simultaneous multi-wavelength diagnostics of stellar CMEs on solar-type stars in different signature (i.e., the Balmer line blueshifts and the possible candidate of X-ray dimming).
Here, bringing all information together, we will discuss the following remaining topics entitled in Section 6.1-6.3 and give a conclusion for each topic.In the following, we also mention previous studies reporting two Hα superflares on the EK Dra, where one of them shows filament eruptions (labeled as "E4", Namekata et al. 2022d) and the other do not (labeled as "E5", Namekata et al. 2022c), as summarized in Tables 3, 4, and 9. 6.1.Can the Prominence Eruptions be Associated with Stellar CMEs?
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 −1 (conservatively, 330-490 km s −1 ), whereas those in flare E2 was 430 km s −1 (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 −1 , it was gradually decelerated and eventually showed almost zero velocity.This may suggest that the 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 are indicative of the occurrence of the associated CMEs which consist of outer coronal mass and fast prominence components, especially in the flare E1: • The large velocity dispersion of 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 E1 event could reach 990-1080 km s −1 (conservatively, 630-880 km s −1 ), 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 lead to CMEs (e.g., Wood et al. 2016, Appendix B).These suggest that a part of the fast prominence components might have escaped.
• 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 flare E1 and E2 could be associated with CMEs, considering their gigantic length scales and velocities (Section 5.2.2).
• 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 32 erg.Since 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 (cf.Namekata et al. 2022d, Supplementary Information).The initial kinetic energy of the prominence eruptions from events E1 and E2 are 5.8 +12.8 −4.0 × 10 34 erg and 1.2 +2.7 −0.8 × 10 33 erg, respectively (Section 5.2.4), which are 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 multi-polar 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).
• After the flare E1, a possible candidate of X-ray dimming were detected, as an indication of escaping the surrounding hot coronal plasma associated with CMEs (Section 5.3).Although the determination of the preflare X-ray level is an issue in detecting stellar XUV dimming, the simultaneous detection with Hα blueshift event is supportive of the detection of X-ray dimming.The estimated upper limit of length scale of the dimming region of up to ∼ (9 − 11) × 10 9 cm is roughly consistent with the length scale of the flares ∼ 7.3 × 10 9 cm.The estimated mass of escaping coronal plasma of 1.2 × 10 17 -2.1 × 10 18 g is enough small that the erupted prominence mass of 4.0 × 10 19 -4.2 × 10 20 g.Therefore, it is not surprising that the detected X-ray dimming amplitudes is likely to be interpreted as CME-related dimming (Section 5.3).However, again we should note that the NICER's Here we uses the L 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 27data 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 the flare loop length and producing eruptive prominences/filaments at greater scales (at least 1-10 times of 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).

Energy Scales
Our observations provides 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: where f E mag is the available magnetic energy stored in active regions associated with starspots, f is the fraction of the non-potential magnetic energy (here, f = 0.1), E mag is calculated as B 2 L 3 spot /8π, 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 is much larger than X-rays and WL radiative energies for flare E1, but the relation is opposite for flare E2 and E4.This diversity may require some theoretical explanations, which is beyoud the scope of this paper.
Our flare observations of EK Dra shows the energy partition into X-rays and WL appeared 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 Xrays (Section 4.2).The derived Equation 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 emits the Hα flare radiation energy of ∼30 % of WLF energy, which is an extremely large portion compared to usual EK Dra's and M-dwarf superflares of 1 % (see, Section 5.1).Hence, E1 is one of the outliers in the derived flare energy scaling relationship 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 emission from prominence eruptions.

Time Evolution
Table 4 catalogs the timescales and timing of blueshift events and potential X-ray dimming events relative to WLFs.E1 flare event offers a probable sequence of events, as discussed in Section 6.1, whereby a stellar prominence begins to erupt almost concurrently with the brightenings of WLFs, potentially leading to the appearance of X-ray dimming approximately 2 hours after WLFs.On the other hand, E2/E4 events have different timescales and timings for the emergence of blueshift events, as seen in Table 4, though they lack detections or observations of X-ray dimming.This suggests that the appearances and formation/acceleration processes of stellar prominences/filaments can have a diverse nature (cf.Section 6.3).

Rotational Phase
Additionally, all five detected superflares (E1∼E5) on EK Dra have simultaneous TESS 's long-term photometric data, allowing us to explore the relationship between flares and the stellar rotational phase.The rotational phase at which the flares occur is summarized in Table 9.By assuming the existence of a single dominant active region causing the rotational variation, we can infer that the flare E1 with prominence eruptions occurred when the dominant active region was near the disk center due to rotation.This seems inconsistent with the occurrence of prominence eruptions outside the stellar limb, suggesting that the source of the prominence eruption is not the dominant active region or that there may be multiple sizeble active regions forming the complex rotational modulation.The E2 flare occurred as the dominant active regions neared the limb, which aligns with the occurrence of a prominence eruption.Finally, the E4 flare event occurred when the dominant spots were on the opposite side of the star, suggesting that the associated filament eruption was produced over the relatively small active regions (Namekata et al. 2022d).These conclusions are not yet definitive, given that observations of other solar-type stars suggest there may be multiple starspots on the surface even in single-spot-like optical light curves (Morris et al. 2017;Basri 2018;Namekata et al. 2020a).The Doppler Imaging of EK Dra, performed in the different epoch from our campaign, shows the gigantic starspots, changing its latitude in time ranging from the equator up to high latitude close to the pole (e.g., Strassmeier & Rice 1998;Rosén et al. 2016;Waite et al. 2017;Järvinen et al. 2018).The near-pole spots usually do not appear in optical light curves but could be the sites of the filament/prominence eruptions.Furthermore, the dark spot causing photometric variability may be not necessarily co-located with regions of strong large-scale magnetic fields (Waite et al. 2017), so we also requires a magnetic-field mapping for further understandings.In our future work, we plan further investigations into the origins of these eruptive phenomena by performing multi-spot and/or magnetic-field mapping using the same observing periods (see Section 6.4).

Diversity and Frequency of Prominence/Filament Eruptions from Young Solar-type Stars
Tables 3, 4, and 9 provide preliminary insights into the diversity and frequency of eruptive flares on young solartype 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 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 does 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 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 field of active regions (Sun et al. 2022;Toriumi et al. 2017).Our limited sample of eruption events shows that the association rate of stellar prominence/filament eruptions against superflares is 3/5 (60 %) on the young solar-type stars.Blueshifts can be detected with our spectrograph (R ∼ 2000, δv ∼ 150 km s −1 ) only when the eruption has relatively high (>150 km s −1 ) LOS velocity.Considering this limitation on the LOS velocity, we note that the absence of blueshifts in flare 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 region 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-sclae photospheric motion driven magnetic flux rope energization suggest an ejection of a CME from a young solar-type star, κ 1 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 young solar-type stars are relatively non-axisymmetric 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), and, 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 recent M-dwarf flare survey shows that M dwarfs blueshift association rate is not so much high (e.g., 7 events among 41 flares on M-dwarfs, YZ CMi, EV Lac, and AD Leo; 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 Table 3, 4, and 9, and also mentioned in Section 6.1 and 6.2.This indicates that there could be even more diversity if we increase a sample size.We will present and discuss further observational results in our forthcoming papers (see Section 6.4).

Future Studies
Our series of studies, including Namekata et al. (2022c,d), 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 the 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(Ikuta et al. , 2023)).Additionally, at the almost similar period of the campaign periods in April 2022, we conducted high-dispersion optical spectroscopy of EK Dra using the High Dispersion Echelle Spectrograph (HIDES; Izumiura 1999) onboard 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 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 multi-wavelength rotational modulations observed in Hα, X-ray, and ground-based B-band photometric observations (Ikuta et al. in prep.).
We will also perform further modeling of the prominence eruptions from flare E1.Specifically, the Hα data offers several key constraints on the locations and trajectories of the prominence eruptions (cf.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 (cf.Leitzinger et al. 2022).Ultimately, we will carry out 3D numerical simulations of prominence eruptions/CMEs, incorporating the data from surface starspots and Zeeman Doppler Imaging (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.8m Seimei telescope, 2m Nayuta telescope, Okayama 1.88m telescope, and Bulgarian 2m 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 (c.f., Getman et al. 2021;Getman & Feigelson 2021;Bouvier et al. 2023) to main sequences, e.g., DS Tuc A and κ 1 Ceti, to eventually understand the Sun in time.
Bulgarian NSF grant No.KP-06-N58/3 (2021).S.V.J. acknowledges the support of the DFG priority program SPP 1992 "Exploring the Diversity of Extrasolar Planets" (JE 701/5-1).A.A.V. acknowledges funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No 817540, ASTROFLOW).The spectroscopic data used in this paper were obtained through the program 22A-N-CN06 (PI: K.N.) with the 3.8m Seimei telescope, which is located at Okayama Observatory of Kyoto University.This paper includes data collected with the TESS mission, obtained from the MAST data archive at the Space Telescope Science Institute (STScI).Funding for the TESS mission is provided by the NASA Explorer Program.STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555.Some of the data presented in this paper were obtained from the Mikulski Archive for Space Telescopes (MAST) at the Space Telescope Science Institute.The specific observations analyzed can be accessed via 10.17909/xgds-j146.NICER analysis software and data calibration were provided by the NASA NICER mission and the Astrophysics Explorers Program. Figure 27 presents the expanded X-ray count rates and exposure fractions of NICER for the ISS-orbital block 3 to 6 before and after the flare E1.The nicerl2, xspec, and lcurve in HEASoft was used for extracting the data in good time intervals of NICER observations, which results in the exposure fractions of each exposure bin not necessarily being unity.The nicerl2 task functions screens out exposures under unfavorable conditions, notably in situations with a high level of particle background or optical light contamination.When NICER encounters these suboptimal conditions, the Good Time Intervals (GTI) are fragmented.Consequently, this fragmentation leads to a decrease in the fractional exposures of the corresponding light curve bins.From an analysis of our dataset mkf file, it is observed that the TOT OVER COUNT is relatively high, nearing a value of ∼200.This raises the possibility that the nicerl2 task may have excised certain segments from the exposures in each light curve bin, especially those not meeting the standard screening criteria.This implies that the relatively high noise level of NICER X-ray count rates in short time cadence for the data on 2022 April 10 is attributable to the shortened exposure time resulting from the bad time corrections.

Figure 1 .
Figure 1.Examples of Hα spectra of EK Dra on 10 Apr 2022.A black solid line is the pre-flare template for flare E1, and a red dashed line is the flaring spectrum of flare E1.

Figure 2 .
Figure2.Overall light curves of EK Dra for (a) Hα EW, (b) TESS 's optical white light, and (c) NICER's X-ray count rate (0.5-3 keV), taken from 2022 April 9 to 2022 April 21.Hα EW is calculated within the wavelength range of 6562.8±10Å.Note that a negative (or negative-direction) EW value corresponds to emission by definition.TESS 's optical white light fluxes are normalized to the average flux.The vertical dashed lines represent detected superflares, labeled as "E1", "E2", and "E3".
tWL and tHα is the FWHM duration of the WL and Hα flare.τWL is the e-folding decay time of the flare.ae The energy and flux values are already corrected for light curve extrapolations.* The value in bracket is the radiation energy of the central component without the blueshifted component.‡ 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).

Figure 3 .
Figure 3.Light curves for the superflare observed on 2022 April 10 (E1).(a) The red line is pre-flare subtracted Hα equivalent width (∆EW = − (EW − EW preflare )) with the grayed error bars.Note that positive ∆EW values correspond to emissions.The Hα EW is calculated within a range of 6562.8±20Å, which differs from the definition applied to flares E2 and E3 (i.e., 6562.8±10Å).(b) Background-subtracted TESS WL (= (F − F bkg )/Fav ).(c) Soft X-ray count rate (0.5 -3 keV).In panel (c), each ISS orbit on this day is labeled as block 1-6.(d) Hα ∆EW and WL, each normalized by their respective peak values.In each panel, the vertical dashed line marks the TESS flare onset time, and the horizontal dashed line depicts the basal level for each wavelength (Hα: median value of -80∼-20 min, soft X-ray: mean value of block 3).

Figure 4 .
Figure4.The same as Figure3, but for the superflare observed on 2022 April 16 (E2).In panel (a), the Hα EW is calculated within a range of 6562.8±10Å. Panel (d) displays Hα ∆EW, WL, and soft X-ray, each normalized by their respective peak values.The basal level of Hα ∆EW is a median value of -45∼-15 min, while the basal level of the soft X-ray is the mean value of block 2.

Figure 6 .
Figure 6.Time evolution of pre-flare subtracted Hα spectra of the superflare observed on 2022 April 10 (E1).(a) Reference light curves in Hα ∆EW (scaled on the left vertical axis) and white light (scaled on the right vertical axis).(b) A dynamic spectrum of the pre-flare subtracted Hα spectrum, plotted on the time-wavelength plane.The red (blue) color represents emission (absorption) signals, with the upper part indicating shorter wavelengths (i.e., blueshift), following the unified format used in our series of papers (Namekata et al. 2022d,c).The color bar shows the non-dimensional intensity which is normalize by the continuum level.The dashed lines indicate the stellar surface gravity (g * ) and half of the surface gravity (0.5 • g * ).(c) Time-binned, pre-flare subtracted Hα spectra during the superflare.Detailed spectra can be seen in Figure 9.The 1-σ errors in continuum level is plotted as error bars for reference.

Figure 6
Figure 6, 7, and 8 shows pre-flare subtracted Hα spectra of flare 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 min (see, Section 3.3.2).

Figure 9 .
Figure9.Temporal evolution of pre-flare-subtracted spectra during superflare on 2022 April 10 (E1).Observations were carried out with an exposure time of 2 (or 1) minute(s) and a readout time of 17 seconds.The 1σ continuum-level scattering in each spectrum is represented by error bars.The blue line corresponds to the results of fitting with a single-component Gaussian, while the red and green lines denote the results of fitting with a two-component Gaussian.The fitting process was conducted using scipy.optimize.leastsquares, taking into account the error bars.

Figure 10 .
Figure10.Temporal evolution of pre-flare-subtracted spectra during superflare on 2022 April 16 (E2).To improve the signalto-noise ratio for fitting, the data has been binned into 5-minute intervals.The magenta and green lines represent the results of fitting with a two-component Gaussian, and the dashed lines correspond to the sum of the two Gaussians.Please see Figure7for the other spectra of flare E2 without significant blueshifted components.
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.05km s −2 .L Hα,blue is the luminosity of the blueshifted component.Mp is the mass of the erupted prominences/filament. E kin is the kinetic energy of the prominence/filament eruptions.§ 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).

]Figure 11 .
Figure 11.Temporal evolution of the blueshifted emission line components for (Upper) flare E1 on 2022 April 10 and (Lower) superlare E2 on 2022 April 16.(Left) Temporal variation in EW (note that positive values mean emissions), (Middle) Temporal variation in velocity, (Right) Temporal variation in velocity dispersion (Gaussian's standard deviation σ).The blue squares represent the result of fitting with a one-component Gaussian, while the red and green lines represent the blueshifted component and the central component, respectively, from the two-component Gaussian fit.

Figure 12 .
Figure 12.Relationship between the velocity and velocity dispersion of the blueshifted component.The data for E1 and E2 are from Figure 11.The filament eruption (blueshifted absorptions) on EK Dra on 2020 April 5 (labeled as "E4") was taken from Namekata et al. (2022d) (labeled as "N22a").

Figure 13 .
Figure13.NICER X-ray spectrum of the superflare on 2022 April 16 (E2) and the corresponding fitting results (applied to the range of 0.3-8.0keV).(a) Fitting results for the pre-flare spectrum of block 2 of E2, fitted with four apec components (each represented by a black line).The blue line represents the total model spectrum, which is used as the pre-flare spectrum in panels (b) through (f).(b-f) Observed and model spectra when flare block 3 is divided into blocks 3a through 3e.Refer to Figure15for the division times and Tables5 and 6for the fitting parameters.

NICERFigure 15 .Figure 16 .
Figure 15.Temporal evolution of the X-ray spectral fitting results for superflares E2 (A) and E3 (B).(a) Light curves of X-ray count rates with a time resolution of 2 minutes in the 0.5-3 keV range.(b) Hardness ratio, which is obtained by dividing the flux in the 0.5-1 keV range by the flux in the 1-3 keV range.(c) Temporal change in temperature, with different colors indicating different time resolutions.(d) Temporal change in the Emission Measure (EM).

4.
TIME EVOLUTION, ENERGETICS, AND LENGTH SCALES OF SUPERFLARES 4.1.Time Evolution of Optical and X-ray Flares

Figure 17 .
Figure 17.The cross-correlation coefficients for different lag times among light curves for the superflare on 2022 April 17 (E3).The time lags of maximum correlation coefficients, rmax, is defined as the most likely value.The 90% percentile gives the lag time uncertainty.The time delay of X-ray against Hα is 6 +3 −2 min, the time delay of X-ray against WL is 7 +3 −2 min, the time delay of differential X-ray against Hα is -3 +1 −2 min, and the time delay of differential X-ray against WL is -3 +2 −4 min.See Tristan  et al. (2023)  for the method.

,Figure 18 .
Figure18.Relations between Hα, X-ray, and white-light energies for solar and stellar flares.(left)The relation between Hα and Bolometric X-ray energies (0.1-100 keV = 0.124-124 Å).Data of superflares on the G-dwarf EK Dra detected in this study are represented by red squares.Other solar and stellar data are taken fromKawai et al. (2022).The function fitted to the observational data with a power-law is shown by the cyan line.(right) The relation between bolometric white-light and Hα energies.Data of superflares on the G-dwarf EK Dra, which are associated with filament/prominence eruptions, are represented by filled red squares, whereas those without eruptions are represented by open red squares (this study,Namekata et al. 2022d, Namekata et al. 2022c).For E1, the total energy including the blueshift component is plotted as the upper limit, with the radiative energy of the central component from the two-component fit being plotted as the lower limit.Data of M-dwarf flares are taken fromMaehara et al. (2021),Namizaki et al. (2023), and(Notsu et al. 2023).Solar flare data are taken fromNamekata et al. (2022a).The function fitted to the observational data, excluding flare E1, with a power-law is shown by the cyan line.

Figure 19 .
Figure 19.Relation between the bolometric white-light flare energy and the peak GOES-band soft X-ray flux (1-8 Å = 1.55-12.4keV) for solar and stellar flares.Data for the superflares (E2, E3) from the G-dwarf EK Dra derived in this paper are plotted as red squares.Other data for superflares from G/F-dwarfs (orange circles) are taken from Guarcello et al. (2019).Data for superflares for M/K-dwarfs are taken from Guarcello et al. (2019), Kuznetsov & Kolotkov (2021), and Stelzer et al. (2022).Note that data from Guarcello et al. (2019) have been converted to GOES band flux and bolometric WLF flare energy.The solar observational scaling laws from Cliver et al. (2022) (dashed line) and Namekata et al. (2017b) (solid line), as well as the empirical law from Shibata et al. (2013) (dash dotted line), are also plotted.

Figure 21 .
Figure 21.Comparison between the velocity (Vr max) and length scale (L) of filament/prominence eruptions on the Sun and G-dwarf EK Dra.The blue circles and orange crosses indicate the solar filament eruptions with and without CMEs, respectively.The green and red points corresponds to the stellar prominence and filament eruptions on EK Dra, respectively (this study,Namekata et al. 2022d).The stellar data are the maximum line-of-sight velocities while the solar data are maximum radial velocity.The solid black indicates an empirical solar threshold that can roughly distinguish whether a filament eruption is associated with CMEs or not in the formula of (Vr max/100 km s −1 )(L/100 Mm) 0.96 = 0.8(Seki et al. 2021).

Figure 23 .
Figure23.Estimated electron density as a function of optical depth.The left y-axis is the case when the prominence is cubic structure, while right y-axis is the cause when the prominence is cylinder (i.e., filamentary structure).The blue colored region corresponds to the assumed optical depth.

Figure 24 .
Figure 24.Possible X-ray (0.5-3 keV) dimming after the gigantic prominence eruptions associated with the superflare on 2022 April 10 (E1).(a) Reference TESS light curve.(b) The X-ray light curve.The blue points represent the average of the gray data with 2-minute resolution for each orbital period.The red points denote the count rate after NICER's background subtraction.The area filled in blue represents the time interval believed to show dimming relative to the pre-flare (block 3).The area filled with gray and red represent the pre-flare level with error bars for blue and red points data, respectively.(c), (d) The temporal evolution of the temperature and emission measure of two components.The area filled with gray in panel (d) represents the pre-flare level with error bars.

Figure 25
Figure25.X-ray (0.5-3 keV) light curves binned for each ISS's orbit around superflares on 2022 April 16 (E2) and 2022 April 17 (E3).The gray points represent the data with a time resolution of 2 minutes, while the blue points correspond to the average of the gray data for each ISS orbit (∼ 90 min period).The horizontal lines denote the pre-flare values.For flare E3, the average value of the previous orbit and the pre-flare value just before the flare are plotted.

Figure 27 .Figure 28 .Figure 29 .
Figure27.(upper) NICER 0.5-3 keV X-ray count rates and (lower) fractional exposure rate for each orbit of block 3-6 on 2022 April 10.Fractional exposure rate is the "FRACEXP" column of the lcurve outputs.The value is the live-time fraction for each time-bin and takes values between zero and one after the deadtime correction with "LDEADTIME".The corrected rate and errors in upper panels are calculated by dividing original values by "FRACEXP".

Table 1 .
Observation logNote- † These are calculated by multiplying the number of frames times the each exposure time (plus readout time ∼17 sec).
PyRAF is part of the stscipython package of astronomical data analysis tools and is a product of the Science Software Branch at the Space Telescope Science Institute.
2 https://mast.stsci.edu/portal/Mashup/Clients/Mast/Portal.html 3 IRAF is distributed by the National Optical Astronomy Observatories, which are operated by the Association of Universities for Research in Astronomy, Inc., under cooperate agreement with the National Science Foundation. 4

Table 2 .
The list of superflares on EK Dra in this study and past observations.

Table 4 .
Properties of the prominence and filament eruptions (continued from Table3) and the possible X-ray dimming.tWLFtblueTiming BlueAsym.X-ray dimming Timing X.dim.t X.dim.amp X.dim.∆EM X.dim.L X.dim.

Table 5 .
Best-fit values of the pre-flare NICER X-ray spectra for E2 and E3.Note-The hydrogen column density N H in TBAbs is fixed as 3 × 10 18 cm −2 .

Table 6 .
Results of NICER X-ray spectral fit for flares E2 and E3.
Note-The hydrogen column density N H in TBAbs is fixed as 3 × 10 18 cm −2 .The abundance in apec is fixed as a pre-flare value of each flare in Table5(0.28for flare E2 and 0.38 for flare E3).

Table 7 .
Summary of NICER X-ray flare parameters.

Table 8 .
Summary of X-ray and bolometric white-light energy of stellar flares in this study and other studies.

Table 9 .
Length scale and area of the superflares, prominence/filament eruptions, and starspots on EK Dra in this study and past observations.

Table 10 .
Length scale of prominences in different prominence parameters for flare E1 and E2.

Table 11 .
Results of NICER X-ray spectral fit for flares E1 X-ray dimming.