A Possible Mechanism for the “Late Phase” in Stellar White-light Flares

M dwarf flares observed by the Transiting Exoplanet Survey Satellite (TESS) sometimes exhibit a peak-bump light-curve morphology, characterized by a secondary, gradual peak well after the main, impulsive peak. A similar late phase is frequently detected in solar flares observed in the extreme ultraviolet from longer hot coronal loops distinct from the impulsive flare structures. White-light emission has also been observed in off-limb solar flare loops. Here, we perform a suite of one-dimensional hydrodynamic loop simulations for M dwarf flares inspired by these solar examples. Our results suggest that coronal plasma condensation following impulsive flare heating can yield high electron number density in the loop, allowing it to contribute significantly to the optical light curves via free-bound and free–free emission mechanisms. Our simulation results qualitatively agree with TESS observations: the longer evolutionary timescale of coronal loops produces a distinct, secondary emission peak; its intensity increases with the injected flare energy. We argue that coronal plasma condensation is a possible mechanism for the TESS late-phase flares.


Introduction
The number of observed stellar flares has increased dramatically since the launch of the Kepler mission (Maehara et al. 2012;Balona 2015;Davenport 2016;Van Doorsselaere et al. 2017;Notsu et al. 2019;Yang & Liu 2019) and the Transiting Exoplanet Survey Satellite (TESS; Günther et al. 2020;Tu et al. 2020Tu et al. , 2021;;Crowley et al. 2022;Pietras et al. 2022).Both stellar and solar flares are thought to be driven by the release of magnetic energy from magnetic reconnection (Priest & Forbes 2002;Benz & Güdel 2010), and while they share many common characteristics, stellar flares can be significantly more energetic than their solar cousins.Stellar superflares can release up to 10 33 -10 36 erg bolometric energy (e.g., Maehara et al. 2012;Günther et al. 2020) compared to solar flares, the greatest of which release less than or equal to about 10 32 -10 33 erg (see discussions of the Carrington flare level in Hayakawa et al. 2023, and references therein).This makes space weather conditions around these stars more hazardous than in the heliosphere.It is worth noting that some of the stellar flares have also been recorded with energy below 10 30 erg (Howard & MacGregor 2022;Pietras et al. 2022).
The majority of optical stellar flare light curves exhibit an exponential rise and a gradual decay (e.g., Davenport et al. 2014).Such an impulsive peak is thought to originate from the low stellar atmosphere (e.g., Hudson et al. 2006).A fraction of flare light curves also display finer features such as quasiperiodic pulsations (QPPs; Rodono 1974;Zimovets et al. 2021).Based on certain solar observations and models, studies have suggested that some QPPs may result from the dynamic response of the stellar atmosphere to impulsive heating (Nakariakov & Melnikov 2009;Van Doorsselaere et al. 2016;Zimovets et al. 2021).
In a recent survey, Howard & MacGregor (2022) found that 42.3% of 3792 M dwarf superflares (energy above 10 32 erg) observed by TESS showed complex morphology in their light curves (Howard & MacGregor 2022).Furthermore, 17% (31 total) of these complex flares exhibited a peak-bump morphology, with a large, highly impulsive peak followed by a second, more gradual Gaussian peak.In another survey of ∼140,000 flares from over 25,000 late-type stars, Pietras et al. (2022) found that approximately 40% of intense long-duration flares can be described by double flare profiles, with the second component dominating the decay phase.As shown in Figure 1, we choose one peak-bump flare from Howard & MacGregor (2022) as an example, and fit the light curve using a double flare profile (see Equation (3) in Pietras et al. 2022).The secondary, late phase has a distinct timescale of tens of minutes, and a significant amplitude, typically a few percent of the stellar flux.Its occurrence rate was too high to be explained by chance occurrence.
The peak-bump morphology is in fact frequently detected in solar flares, typically in the extreme ultraviolet (EUV) and soft X-ray (SXR) wavelengths (Woods et al. 2011;Qiu et al. 2012;Dai et al. 2013Dai et al. , 2018;;Liu et al. 2013Liu et al. , 2015;;Sun et al. 2013;Dai & Ding 2018;Chen et al. 2020b).The EUV/SXR late phase originates from long coronal loops that are magnetically connected to the main flare site.During the impulsive phase, these loops are rapidly filled with high-temperature plasma, which requires a longer time to cool via conduction compared to shorter loops (Reale 2014).A secondary peak will be delayed in the spectral irradiance of cooler lines, with a timescale ranging from tens of minutes to hours.Additional heating and further magnetic reconnection may also contribute to the delay, and produce late-phase hot emission at somewhat lower temperatures than in the primary event.
What is the cause of the late phase in TESS optical light curves?Does it originate from the coronal loops, similar to the solar counterpart observed in hotter EUV/SXR emitting plasma?Before answering these questions, we note that the solar flare energy content is on the lower end of the known stellar flare population.Its impact on total solar irradiance is at most a few 10 −5 (Kretzschmar et al. 2010), orders of magnitude less than the detected stellar counterparts.Spatially resolved solar white-light flare observations suggest that the enhancements of the optical continuum originate from the lower atmosphere. 6The locally integrated light curves are mostly single peaked, with rare exceptions (e.g., Matthews et al. 2003;Kerr & Fletcher 2014;Hao et al. 2017).
Off-limb white-light structures have been observed in the gradual phase of some solar flares, implying a significant enhancement of the density in the coronal loops (e.g., Hiei et al. 1992;Fremstad et al. 2023, and references therein).More examples of such observations (Martínez Oliveros et al. 2014;Fremstad et al. 2023) have been obtained from the Helioseismic and Magnetic Imager (Schou et al. 2012), which provides pseudo-continuum images derived from the Fe I 6173.3Å line observations.The most famous recent example is from the bright solar flare SOL2017-09-10T15:35. Situated on the west limb, it hosted post-reconnection coronal loops emitting brilliantly in the continuum against the background (Jejčič et al. 2018;Zhao et al. 2021;Fleishman et al. 2022;Martínez Oliveros et al. 2022).The loop apexes reached about 18 Mm (25″); the continuum emission in those loops peaked about 30 minutes after the SXR peak.
Based on such observations, Heinzel et al. (2017) and Heinzel & Shibata (2018) proposed that optically thin emission from the corona can meaningfully contribute to the white-light stellar flare emission.Their calculations showed that the intensity of the observations can constrain the electron number density (n e ) and temperature (T) within a specific parameter regime.A high electron density n e  10 12 -10 13 cm −3 is unequivocally required, leading to free-free and free-bound emission mechanisms dominating over Thomson scattering.While such high coronal densities have indeed been reported for some solar flares (Hiei et al. 1992;Heinzel et al. 2017;Jejčič et al. 2018), it remains unclear what physical mechanism can generate the specific (n e , T) values sufficient to produce intense enhancements of the continuum in the coronal loops (see, for example, the discussion in Section 2.3 of Kerr 2023).
We note that another solar phenomenon, known as coronal rain (e.g., Scullion et al. 2016;Antolin & Froment 2022;Mason & Kniezewski 2022, and references therein), results from the rapid increase in local coronal density by orders of magnitude after the flare impulsive phase.It is generally explained by thermal instability (Parker 1953;Field 1965;Claes & Keppens 2019; Figure 1.An example of a TESS flare light curve with a peak-bump morphology from the TESS Input Catalog (TIC) ID 272232401.The x-axis is a fraction of a day starting from Barycentric TESS Julian Date (BTJD) 2238.Black dots represent data outside the fitting range, red circles indicate data used for fitting, and gray circles denote excluded data around the time of 2238.23.These data are excluded due to a small local maximum that may come from another region of the star and is not well described by the peak-bump morphology.The light curve is fitted using an empirical double flare template f (t) from Pietras et al. (2022): , where t represents time.The fitted impulsive phase (first portion of f (t)), late phase (second portion of f (t)), and total flare model are shown as gold dashed-dotted, cyan dashed, and green solid lines, respectively.The best-fit parameters are hr, and D 2 = 1.12 hr −1 .The χ 2 of the fit is 2.72.Claes et al. 2020), where the plasma cooling rate increases drastically due to the density enhancement related to radiative energy loss.In a detailed magnetohydrodynamic (MHD) simulation, condensed plasma in the form of coronal rain appeared tens of minutes after the onset of magnetic reconnection (Ruan et al. 2021).This mechanism has been evoked to explain the in situ formation of solar prominences (Xia & Keppens 2016;Kaneko & Yokoyama 2017;Zhou et al. 2020), and can potentially produce dense recombining plasma capable of white-light flare emission.
In this work, we propose coronal plasma condensation7 as a possible mechanism for the observed peak-bump morphology in TESS stellar flare light curves.To probe the underlying mechanism, we perform a suite of one-dimensional (1D) hydrodynamic (HD) loop simulations for typical M dwarf flare parameters.We show that the simulated, optically thin whitelight emission, in terms of the evolution timescale and the total flux, is qualitatively in agreement with TESS observations.Below, Section 2 describes our 1D HD simulations and forward synthesis of the TESS light curve.Section 3 presents our modeling results.Section 4 presents a summary and discussion.

HD Simulation
The solar flare occurs in a 3D volume, but performing parameter studies efficiently and simulating a realistic chromosphere exceeds the current computational capabilities of 3D models.However, solar flares consistently occur in loop structures that follow the geometry of magnetic field lines in the corona (Benz 2017).These magnetic field lines guide the flow of hot plasma and produce coronal loops with enhanced emission.The standard solar flare model suggests that magnetic reconnection heats the plasma in the loop, converting magnetic energy into heat, bulk motions, and radiation emissions such as X-rays.Due to the confined and elongated nature of the loop, it can be approximated as a 1D structure for modeling purposes, which simplifies and facilitates simulations using 1D HD models (e.g., McClymont & Canfield 1983).This approach has been successfully used in many solar flare simulations, including the physics of energy and radiation transport, such as nonthermal particles, Alfvén waves, thermal conduction, and radiation (Kerr 2022, 2023, andreferences therein).Onedimensional models have also been employed to simulate M dwarf flares (Allred et al. 2015;Kowalski 2023).
We use the open-source Message Passing Interface Adaptive Mesh Refinement Versatile Advection Code 2.0 (MPI-AMRVAC; Xia et al. 2018;Keppens et al. 2021) to simulate a flaring loop via a 1D HD model.The MPI-AMRVAC code is a parallized partial differential equation solver framework that contains many different numerical schemes in multi-dimension with multiple physics modules.This code has been successfully applied to simulating solar flare/eruption (Fang et al. 2016;Ruan et al. 2020Ruan et al. , 2023;;Guo et al. 2021Guo et al. , 2023;;Zhong et al. 2021Zhong et al. , 2023;;Druett et al. 2023) as well as the formation and evolution of cold material in the solar corona (Xia & Keppens 2016;Xia et al. 2017;Zhou et al. 2018Zhou et al. , 2020Zhou et al. , 2023;;Hermans & Keppens 2021;Jenkins & Keppens 2022).For the simulation in this paper, we employ the HLLC flux scheme (Toro et al. 1994) and a five-order weighted essentially nonoscillatory slope limiter (Liu et al. 1994).We consider optical thin radiation loss from the radiation loss curve by Colgan et al. (2008) and include thermal conduction using Spitzer conductivity.The saturation of thermal conduction is implemented when the electron is as fast as the speed of sound (Cowie & McKee 1977).We also implement a transition region adaptive conduction method to capture mass evaporation and energy exchange more accurately (Johnston & Bradshaw 2019;Zhou et al. 2021).The 1D HD loop simulation based on MPI-AMRVAC has been validated in prior research (Xia et al. 2011;Zhou et al. 2014Zhou et al. , 2020Zhou et al. , 2021)), which demonstrated its capability and flexibility required for this study.We note also that since flare plasma is generally confined by the magnetic field of the loops, a 1D field-aligned approximation is a reasonable assumption (for an extensive discussion of the utility of 1D flare loop models to understand the physics of flares see Kerr 2022).Due to the absence of treatment of optically thick radiation, nonthermal electrons, waves, and turbulence, this code is less suitable for studies of the lower atmosphere than other (radiation-) HD codes such as RADYN (Allred et al. 2015).Those aspects, however, are not essential to our goals and are beyond the scope of this paper, which focuses on the flaring coronal dynamics.
As in most previous work, we model an individual flare loop as a semi-circular tube with a uniform cross section.The loop coordinate l ranges from −πR L /2 to πR L /2 from one end to the other, where R L is the major loop radius, and represents the loop height as well.For the initial condition, we adopt a simplified atmospheric model with a temperature profile T defined by where h is the height.Here, T b = 10 kK is the temperature at the footpoint, h = 0.This value is based on the observed chromospheric temperatures in M dwarfs, which can range from several thousand to tens of thousands of Kelvin (Cram & Mullan 1979;Mauas 2000).We note that RADYN radiative HD simulations by Kowalski et al. (2016) revealed that a temperature of approximately 10,000 K is required to accurately replicate flare emission on M stars in white light.At the loop top (apex), the loop coordinate is l t = 0; the temperature is T t = 6 MK, typical for M dwarf corona (Allred et al. 2015).The transition region has a height of = h 5 tr Mm, and a thickness of = w 250 tr km.These values are on the higher end of the typical solar values (Zhang et al. 1998;Tian et al. 2008).Such a piecewise-like temperature profile has been widely used for simplified solar atmospheres in previous studies (Bradshaw & Mason 2003;Fang et al. 2016;Johnston & Bradshaw 2019).The hot coronal portion of the loop is maintained using a uniform background heating rate of Q bg = 2 × 10 −2 erg cm −3 s −1 , which is in the typical range of solar coronal heating power (Sakurai 2017).
We set the number density at the loop apex as n e = 10 10 cm −3 .The initial density and pressure along the loop are computed using ideal gas law assuming a hydrostatic atmosphere.The stellar parameters used here are the median of those M dwarfs with peak-bump light curves from Howard & MacGregor (2022): effective temperature T eff = 3323 K, radius R å = 0.48R ☉ , and stellar surface gravity acceleration g = 5.65 × 10 4 cm s −2 .The model first runs with only background heating for 2 × 10 4 s, allowing the atmosphere to relax into equilibrium.
Following an empirical scaling law based on stellar flare observations (Shibata & Yokoyama 1999, 2002), we set the nominal flare heating rate as Q = B 2 V a /(4πL), where B is the characteristic magnetic field strength, L = πR L is the loop length, and is the Alfvén speed.The mass density ρ corresponds to a nominal electron number density 5 × 10 9 cm −3 , which is the value at the loop top after the relaxation (see 10 10 cm −3 initial value).
In this work, we investigate two loop radii (apex heights), 10 and 30 Mm, which are typical for average and large solar active regions (Benz 2017) and are well consistent with stellar flare statistics (Shibata & Yokoyama 1999, 2002;Namekata et al. 2017).We further consider three different magnetic field strengths, 400, 600, and 800 G.The choice of 400 G is based on the microwave observations of a coronal current sheet during a large solar flare SOL2017-09-10T15:35 (Chen et al. 2020a).The other two values are meant for the M dwarfs, where the mean photospheric field can reach the kilogauss range with much larger starspots (Berdyugina 2005;Shulyak et al. 2019;Reiners et al. 2022).
As illustrated in Figure 2(a), we simulate the flare energy injection as time-dependent, localized heating sources at the loop apex (H t ) and the footpoints (H b ).The former mimics the impulsive heating directly caused by magnetic reconnection, whereas the latter mimics a secondary, gradual heating of the lower atmosphere.The approach is widely used in 1D flare simulations (see discussion in Section 4).The two heating terms are formulated as Here, l t = 0 is the loop coordinate of the apex, and = h 5 tr Mm is the height of the transition region as defined earlier.We use λ = 1 Mm as the typical spatial scale of the HD flare heating model (Bradshaw & Mason 2003), which is compatible with the calculation based on the nonthermal electrons' energy deposition layer (Radziszewski et al. 2020).This value for the loop-top source is justified by the bow shock length along loops in 2D simulations (Chen et al. 2015;Ruan et al. 2020).We note that the footpoint heating is applied at the transition region height, = = h h 5 tr Mm, rather than the actual footpoint h = 0 Mm.The total energy input E is estimated as the temporally and spatially integrated heating H t + 2H b (note the factor of 2 is to account for both footpoints; see the notes for Table 1 for more detail).The piecewise function f i (t), where i stands for loop top (t) or footpoint (b), is defined as  Normalized temporal profiles of the heating rate for the apex f t (t) (blue solid line) and footpoints f b (t) (dotted lines) for the three different cases listed in Table 1 (orange, green, red).The decay timescales of the footpoint heating rate are 14.5, 29.5, and 59.5 minutes, respectively.The left and right portions depict the initial impulsive loop-top heating and the late-time gradual footpoint heating processes, respectively.They are discontinuous and have different scales on the horizontal axis.
definition.Finally, a scaling factor k i , defined as makes the time-integrated energy from one footpoint and the loop-top heating sources to be approximately equal.We consider this equal-partition assumption to be appropriate given the lack of knowledge of the heating mechanism.
The heating time profiles used in this study are shown in Figure 2(b).We assume the loop top and the footpoint heating start and reach their peak at the same time with identical rising timescales.Specifically, we use t s = 0 s, t p = 30 s, and τ r = 30 s for both sources.For the more impulsive loop-top source, the decay timescale is fixed at τ d,t = 30 s.For the more gradual footpoint sources, we consider three different decay timescales, τ d,b = 14.5, 29.5, and 59.5 minutes.The corresponding end times are t e,b = 15, 30, and 60 minutes, respectively.We note that the 60 s impulsive heating profile is typical for a single loop in the solar cases (Bradshaw & Cargill 2013;Qiu & Longcope 2016).The gradual heating durations, 15 and 30 minutes, are typical for the solar cases (Qiu & Longcope 2016;Zhu et al. 2018).The 60 minute cases are longer than what has been reported for solar flares, but may well be within the range of more energetic stellar flares.
For this study, we perform a total of 10 simulations with various combinations of free parameters, as summarized in Table 1.We divide these simulations into four groups.For the first group (Cases 1-6), we fix the decay time τ d,b = 29.5 minutes, and explore the effect of changing loop radius R L and magnetic field B. For the second group (Cases 7 and 8), we fix R L = 10 Mm and B = 800 G, and explore the effect of the changing decay timescale of footpoint heating.For the third group (Cases 9-10), we experiment with having only the footpoint or loop-top heating source.In the fourth group (Cases 11-13), we examine the impact of asymmetric footpoint heating by scaling the footpoint heating H b (t) at l = πR L /2 by factors of 0.9, 0.5, and 0.1, respectively.All other parameters remain the same as those in Case 1.

Emission Synthesis
To investigate gradual phase emission, we synthesize the optically thin continuum emission I ν assuming the loop is filled with a completely ionized, hydrogen plasma (Heinzel et al. 2017;Heinzel & Shibata 2018;Jejčič et al. 2018).We consider only contribution from the coronal portion of the loop that is above the critical height where T(h) > 0.2 MK in the initial, relaxed atmosphere.We consider only the wavelength range between 534 and 1060 nm used by the TESS filter.
We assume an off-limb flare loop in local thermal equilibrium and ignore the background intensity.The optically thin hydrogen emission from free-free ( n I ff ) and free-bound ( n I fb ) mechanism can be expressed as where D is the geometric depth of the emitting plasma along the line of sight, assumed to be equal to the loop length D = L for simplicity.Additionally, B ν (T) is the Planck function, κ ff and κ bf are the hydrogen free-free and bound-free opacity, given by 1.166 10 , 1 .7 Here, n e is the electron density; i and ν i represent the principle quantum number and the continuum limit frequency of the specific spectral series, respectively; g ff (ν, T) and g bf (i, ν) are Notes.
a The decay timescale of the footpoint sources, which provides the gradual heating.In Case 10, there is no footpoint source.
b The total energy, calculated as ò , represents the sum of the heating at the loop top and two footpoints.The temporal profile of the heating has a triangle shape with durations of τ r,t + τ d,t and τ r,b + τ d,b for H t and H b , respectively (Equation ( 3)).The first term is estimated to be 0.5πλL 2 Q(τ r,t + τ d,t ), where we use λ = 1 Mm as the heating scale, L = πR L as the loop length, and πL 2 as the flaring area.The second term is the energy from two footpoints, estimated to be ).The associated total energies are calculated accordingly.c Cases 9 and 10 only contain the footpoint and loop-top sources, respectively.d Cases 11-13 contain asymmetric heat source settings.The heat sources at footpoint l = πR L /2 are reduced by multiplying a factor of 0.9, 0.5, and 0.1, respectively.
the corresponding Gaunt factors, assumed to be unity, and k B is the Boltzmann constant.The free-bound emissions are calculated from Paschen to Humphreys continua.Furthermore, the Thomson scattering by electrons is formulated as where σ T = 6.65 × 10 −25 cm −2 is the cross section for Thomson scattering, and n J inc is the incident intensity, which equals B ν (T eff ) times a dilution factor of 0.4 (Heinzel & Shibata 2018).Finally, the total flare loop emission is The observed TESS stellar flare light curve is generally normalized by the quiescent flux from the stellar disk.To compare our model with observations, we estimate the relative flare flux by calculating the ratio between the emission from the coronal flare loop and the background emission.The latter is approximated by a uniform disk of blackbody radiation B ν (T eff ).As a first-order estimate, we ignore the limb-darkening effect.To determine the emission from the flare loop, we integrate I ν f (ν) across frequency ν and along the loop l, multiply it by the spatial length scale πL to emulate the effect of multiple loops.Finally, both the flare and the stellar disk flux are scaled by the TESS filter response function f (ν) (Ricker et al. 2015):

Thermodynamic Evolution
The modeled loop thermodynamic evolution is shown in the top three rows of Figures 3-7.We use colors to illustrate the evolution of the physical parameters in the loop as a function of time.In all cases, the atmosphere reached a steady state after a relaxation period of 2 × 10 4 s before the heating source turns on at t = 0.The two dense and cool layers near the footpoints correspond to the chromosphere and transition region, while the remaining region in the central part is the hot and tenuous corona.
Shortly after the heating onset, the temperature near the loop apex l = 0 rapidly increases to nearly 108 K, whereas the density decreases briefly.This is owing to the impulsive energy injection at the loop top, and is most pronounced in Cases 4-6.In the meantime, the top of the chromosphere is heated and generates a strong evaporation flow.Shortly after, the temperature of the coronal portion of the loop rises to several tens of million Kelvin.As strong upflows toward the loop top develop due to the expansion of the chromospheric plasma, the number density of the loop gradually increases to about 10 11 -10 12 cm −3 in all cases.
As the plasma cools, the increased density in the coronal loop triggers thermal instability.For shorter loops (Cases 1-3), a dramatic increase in density (to n e > 5 × 10 12 cm −3 ) occurs near l = 0 at t = 7 minutes, which is accompanied by a rapid drop in temperature to around 10 4 K.This cool plasma remains at the loop top for an extended period of time.In contrast, cool plasma in longer loops (Cases 4-6) first appears around l ≈ ± 35 Mm, corresponding to a height of 12 Mm above the footpoints, approximately 16 minutes after the onset of heating.
The condensed plasma (white contours in the number density diagrams) is co-spatial with the front edge of the rapid upflow (same locations in the velocity diagrams).The pattern suggests that the material is propelled by an evaporation acoustic shock with speeds ranging from 50-170 km s −1 (see the local speed of sound of 10-170 km s −1 ).
As the footpoint heating diminishes, the condensed plasma starts descending toward the chromosphere 8 around t ≈ 30 minutes in Cases 1-6.During this phase, footpoint velocities reach approximately the freefall value, about 100 km s −1 for the shorter loops (Cases 1-3), and 200 km s −1 for longer loops (Cases 4-6).Meanwhile, global catastrophic cooling can occur elsewhere in the loop without condensation (Cargill & Bradshaw 2013).This is evidenced by the rapid decrease in temperature in Case 10 (white arrow).The temperature and density of the cool coronal plasma can also be modulated by acoustic shock waves that are being reflected in the loop.The oscillating patterns are clearly visible in the density, temperature, and velocity profiles from 20-30 minutes in Cases 4-6, which will result in oscillations in the synthesized TESS light curves.
We find that the thermodynamic evolution in Cases 7 and 8 is similar to that in Case 3. The shorter (longer) footpoint heating duration yields an earlier (later) onset of the falling plasma.The values are about 16 and 45 minutes for Cases 7 and 8, respectively.Sustained footpoint heating in Case 8 also delays global catastrophic cooling from happening in that model.
We find that the gradual footpoint heating is crucial to coronal condensation in our model, based on Cases 9 and 10.With only footpoint heating, the evolution of Case 9 is similar to Case 3. The difference appears to be a delayed initial coronal condensation at about 15 minutes.With only a loop-top source, Case 10 does not produce any coronal plasma condensation.A catastrophic cooling phase is accompanied by a gradual decrease in density.This result is consistent with previous findings (Antiochos 1980a;Reep et al. 2020;Antolin & Froment 2022).
Cases 11-13 demonstrate that with increased asymmetry of footpoint heating, both coronal condensation and flare emission diminish.Case 11, differing minimally from Case 1, shows a slightly reduced duration of condensed gas in the corona.In Case 12, secondary and tertiary condensation episodes occur around the 20 and 30 minute marks, respectively.These cycles of evaporation and condensation may be attributed to thermal nonequilibrium (Froment et al. 2015(Froment et al. , 2017(Froment et al. , 2020)).Subsequent cycles exhibit decreasing amounts of condensed plasma due to reduced footpoint heating.In the highly asymmetric Case 13, the negligible plasma condensation in the corona is accompanied by a substantial flow between footpoints.This finding, characterized by a large heating ratio of 10 between the two footpoints, agrees well with the results in Klimchuk & Luna (2019).

Signatures of a Late-phase Flare
Our synthesized light curves (Cases 1-9 and 11) qualitatively reproduce the two main features of the TESS late-phase light curves reported in Howard & MacGregor (2022).First, the coronal emission from the plasma condensations exhibits a gradual peak, about 20 minutes after the onset of impulsive heating.Second, the maximum relative fluxes of the simulated late-phase range from 10 −3 to 10 −1 .These features are quantitatively summarized in Table 2.
We list below several observations of the modeled latephase flare.
1.For a fixed loop length, greater energy injection leads to greater peak flux.For the longer loops (Cases 4-6), the peak value is reached during rapid oscillations driven by reflecting shocks.2. Longer loop lengths lead to delayed coronal condensation and shorter late-phase duration.The peak times of the late phase do not exhibit an obvious trend.3. Longer footpoint heating leads to longer late-phase duration.4. The absence of gradual and/or strong asymmetric footpoint heating leads to no coronal condensation, and therefore no late-phase flare emission.
We note that our models only apply to the optically thin corona.The optically thick emission from the lower atmosphere, which is believed to dominate the impulsive phase, is naturally missing from the model.The only exception is from the first few minutes in Cases 4-6 (Figure 4).The impulsive peak comes from the sudden, large amount of plasma evaporated into the lower corona (38  |l|  42 Mm), which is associated with the onset of heating.This is most visible in panel (a) of Figure 10 in the Appendix.
In order to account for the optical emission from dense layers (which may very well be optically thick), one needs to properly model the flare chromosphere with sufficient spatial resolution and including non-LTE and nonequilibrium radiation transfer.That is, coupling the radiative transfer equation within the HD equations self-consistently (e.g., Carlsson & Stein 1992, 1995, 1997;Heinzel et al. 2016), or employing the flux-limited diffusion approximation method (Levermore & Pomraning 1981;  Moens et al. 2022).However, our focus is on the corona, whose thermodynamics with the optically thin radiation loss have been validated by prior research (Antiochos 1980a(Antiochos , 1980b;;Zhou et al. 2020).Assessing the contribution from the denser lower atmosphere representing the impulsive phase is beyond the scope of this work, which focuses on the late-phase emission that we speculate originates in the corona.
In the Appendix, we estimate the contribution of various emission mechanisms during the late phase.We find that the hydrogen free-free continuum dominates the white-light emission during the early stages, whereas the hydrogen recombination continua dominate the late stages.Thomson scattering is negligible in all cases.
While our synthesis assumes simultaneous energization of all loops, it is important to acknowledge that the timing and relative flux of the late phase may vary if there are time delays between successive loop activations.Staggering the activations within a short time span, ranging from seconds to several minutes, could slightly extend the duration and magnitude of the late phase.This adjustment represents a slight change in the superposition of individual loops but is not expected to significantly affect the main findings of our study.
One point of interest is the life cycle of plasma and its emission signature, from its initial evaporation in the chromosphere, through condensation in the corona, to its ultimate return to the chromosphere.We investigate this by tracking a single plasma parcel in Case 1.The location of the parcel is determined by integrating the velocity in time from its initial position; its volume is simply - n e 1 . The results are depicted in Figure 8.As the parcel initially ascends (8-11 minutes), its density decreases sharply due to significant expansion.The temperature, however, drastically increases to above 10 MK coronal value due to footpoint heating.The density (2 × 10 11 cm −3 ) and temperature (10 7 K) remain relatively stable as the ascension continues.Upon reaching the loop apex (18 minutes), thermal instability sets in, causing the temperature to plummet from 10 7 to 10 4 K, and the number density to increase by two orders of magnitude.The subsequent descent is characterized by gradual temperature and density decrease.
It is evident from Figure 8 that the synthetic emission correlates with the number density.This is because both freefree and bound-free opacities are proportional to the square of the number density, k µ n n e 2 , as shown in Equations ( 6) and (7).The total emission varies as I ν ∝ n e after the parcel's volume - n e 1 is factored in.Furthermore, the dominating emission mechanism changes over time due to the ).After the thermal instability occurs, the temperature decrease results in a reduction of overall emission.The free-bound emission declines less significantly than the freefree emission, so it dominates the late phase, especially between 30 and 35 minutes.In fact, we show that the late-phase emission predominantly stems from the cold, condensed coronal plasma cooler than 0.158 MK, with contributions from the hotter plasma being negligible (Appendix and Figure 13).

Summary and Discussion
In this paper, we perform a suite of 1D HD simulations for M dwarf flares with energy from about 10 33 -10 34 erg.We assess the optically thin emission from the coronal plasma (Heinzel et al. 2017;Heinzel & Shibata 2018), and find that coronal plasma condensation can lead to significant emission in the optical continuum.The synthetic light curves exhibit a pronounced secondary peak, whose delay time (from the initial flare heating) and the relative magnitude are qualitatively consistent with the observations from Howard & MacGregor (2022).We thus propose coronal plasma condensation as a Our simulations suggest that the late-phase magnitude increases with flare energy.A substantial amount of heating is required to bring the chromospheric plasma into the corona before condensation can take place, and to raise the electron density sufficiently high to produce meaningful white-light emission.The least energetic event (Case 11) here has a late phase of merely 0.5% of the stellar flux, but requires 1.35 × 10 33 erg of energy input, about 10 times more energetic than the most intense solar flares such as SOL2017-09-06T12:00.This may explain why the white-light coronal loops are rarely observed in solar eruptions.Furthermore, the Sun is brighter than the M dwarf we studied, and the lack of systematic observation of these factors might also contribute to the rare sightings of these phenomena.
Our findings, as outlined in Table 2, demonstrate that the energy emitted through the free-bound mechanism within TESSʼs optical range is about twice as large as that emitted through the free-free mechanism.Notably, the contribution of Thomson scattering to the total flare energy is minimal.On the other hand, the coronal plasma condensation process gives rise to a substantial amount of compressed cold plasma, predominantly situated at the loop's apex, creating a density that is two to three orders denser than the hot counterpart, as demonstrated in Figures 3-7.Under the theoretical framework provided by Heinzel & Shibata (2018) and Jejčič et al. (2018), the collision rate is shown to be proportional to the number density, resulting in the free-free and free-bound emission being directly proportional to the square of the density (see Equations ( 6) and ( 7)).The dependence of emission on temperature, on the other hand, is moderate.This is evidenced by a maximum tripling in emission when the temperature declines from 1 to 0.01 million Kelvin (refer to Figures 4 and 5 in Jejčič et al. 2018).We can intuitively understand the importance of the cold plasma through the following simple calculation.Assuming the condensed plasma accounts for a mere 1% volume at the total loop, owing to its high density, the aggregate emission from these condensed materials would be between 100 and 10,000 times greater than that from the rest of the loop's hot portion, a difference starkly portrayed in Figure 12 in the Appendix.
We note that the emission synthesis is carried out separately from the simulation, making it not fully self-consistent with the thermal instability process in the simulation.Detailed exploration of the consequences is beyond the scope of this work.Still, it is important to acknowledge the critical role played by the condensed cold plasma component in generating the observed TESS flux throughout the entire late phase.
The simulated late-phase flare is accompanied by fast plasma evaporation and draining, whose velocities range from 50-150 km s −1 .These motions may be probed with optical and UV spectral lines that are sensitive to relevant temperatures.Fast Doppler velocities have been reported for stellar flare observations, from several tens to hundreds km s −1 , that might have resulted from the chromospheric evaporation and/ or coronal rain (Argiroffi et al. 2019;Wu et al. 2022;Namizaki et al. 2023).
The reflecting acoustic shocks in the modeled flare loops modulate the plasma density and temperature, causing the synthetic light curves to exhibit oscillations reminiscent of those in TESS flare observations.The oscillation periods in Case 4-6 are approximately 1.5-2.5 minutes, which falls within the observed range (2-36 minutes; Howard & Mac-Gregor 2022).This may serve as a possible mechanism for the QPPs observed in some stellar flares (Nakariakov & Melnikov 2009;Van Doorsselaere et al. 2016;Zimovets et al. 2021).Our simplified model of course cannot account for the complex   MHD wave modes that are known to be important to flare dynamics.
In this study, we use loop-top and footpoint heating sources to mimic the flare energy injection.Physically, the loop-top heating is directly related to the outflow of the magnetic field reconnection, either from the bow shock on the loop top or the collision of energetic particles from the reconnection (Masuda et al. 1994;Shibata et al. 1995;Forbes & Acton 1996;Guidoni & Longcope 2010;Unverferth & Longcope 2020, 2021;Fleishman et al. 2022).The footpoint heating, on the other hand, can be attributed to turbulence dissipation of Alfvénic waves following loop retraction (e.g., Ashfield & Longcope 2023), the thermalization of energetic particles, which is nonthermal electrons or protons, (e.g., Brown 1971;Emslie 1978;Holman et al. 2011), Alfvén wave dissipation (e.g., Emslie & Sturrock 1982;Fletcher & Hudson 2008;Kerr et al. 2016;Reep et al. 2016), thermal conduction, and other wave-particle processes (Kowalski 2023).Many recent flare models, such as RADYN (Carlsson & Stein 1992, 1995, 1997;Allred et al. 2005Allred et al. , 2015) ) and HYDRAD (Reep et al. 2013), typically adopt nonthermal particles as the energy injection mechanism (though they have explored alternatives).Thus, they tend to focus on the initial impulsive footpoint heating, and it has been recognized that they do not capture the longerduration gradual phase (see discussions in Reep et al. 2020;Allred et al. 2022;Kerr 2022Kerr , 2023)).Investigating the physical nature of these heating sources is beyond the scope of this paper.Instead, we take a mechanism-agnostic approach with our experiments, being concerned primarily with the magnitude of heating and the resulting effects.
The extended, gradual footpoint heating appears to be crucial for triggering coronal condensation.9Reep et al. (2020) found that footpoint heating by electron beams alone cannot produce coronal condensations in flare simulations.Some other mechanism must act alongside impulsive footpoint heating, further motivating our deposition of energy directly into the footpoint portion of our flare loop that continues energy transport through the gradual phase.
Evidence supporting the gradual heating phase is abundant in solar observations (Qiu & Longcope 2016;Zhu et al. 2018).It could be due to wave turbulence (Ashfield & Longcope 2023), turbulent suppression of thermal conduction (Allred et al. 2022), or long-lasting magnetic reconnection (Ruan et al. 2021).For long-lasting magnetic reconnection, an extended current sheet trailing a coronal mass ejection may be required (Chen et al. 2020a).
For the footpoint sources, pronounced heating asymmetry diminishes the likelihood of thermal nonequilibrium, as evidenced by our Cases 11-13 and Klimchuk & Luna (2019).As the energy input is not expected to be symmetric in most flare loops, this could account for the infrequent occurrence of late-phase activity in many superflares.
It is worth noting that we do not consider the variation in the locations of the footpoint sources, which has been reported in solar flares (Radziszewski et al. 2020).Our justifications are threefold.First, these variations are due to nonthermal electrons propagating from the corona to the lower atmosphere, losing energy in the upper and middle chromosphere via thick target collisions.Deposition of energy in deeper layers is also possible via the beam's interaction with Langmuir and ionacoustic waves (Kowalski 2023).However, the electron beam alone does not appear to be sufficient to induce coronal condensation in HD flare modeling (Reep et al. 2020).Second, while there is strong evidence of nonthermal electrons during the impulsive phase in solar flares, our focus here is rather on the gradual phase.There is no compelling evidence for substantial amounts of nonthermal electrons in these later stages.Third, since the time variation for footpoint source locations is not well constrained, we opt to use a constant height as in this exploratory study.The location at = h 5 tr Mm and the thickness λ = 1 Mm aligns with the height from turbulence heating in Ashfield & Longcope (2023) and the thickness in Radziszewski et al. (2020).
Our emission synthesis assumes a local thermal equilibrium condition and optically thin continuum radiation, which can be valid in a range of coronal plasma parameters (Heinzel et al. 2017;Heinzel & Shibata 2018).It lacks proper treatment of the dense lower layers, which has been touched upon by prior studies of impulsive chromospheric flare sources using radiation HD simulations (e.g., Kowalski et al. 2015Kowalski et al. , 2017;;Kleint et al. 2016;Kowalski 2022).We therefore do not model the impulsive phase of the flare, and focus on the late phase in the corona instead.We find that the maximum optical depth of the condensed plasma is around unity. 10n reality, the condensed plasma in the flare loops will not be perfectly aligned along the line of sight; they will also not evolve synchronously.The total emission will then come from a superposition of many segments of optically thin plasma: our conclusions are thus expected to hold qualitatively.
To further improve our understanding of stellar flares, additional observations are needed.Combining spectroscopic observations with the TESS white-light data could provide valuable insights.Spectroscopy can offer detailed information about the temperature, density, and composition of the flaring plasma, complementing broadband photometric observations.Moreover, EUV observations could help determine if the EUV late phase observed in solar flares is also present in stellar flares.By studying these aspects in greater detail, we can gain a more comprehensive understanding of the physical processes underlying stellar flares.

Appendix Contribution of Various Emission Mechanisms
The relative contribution from various emission mechanisms along the flaring loop are calculated as ò ò To differentiate the contribution from the hot and cold materials in the corona, we apply a temperature threshold of 0.158 MK, equivalent to the hydrogen ionization energy of 13.6 eV.The emission calculations are limited to the corona region displayed in Figure 9, as illustrated in Figure 12.    Figure 13.The emission contribution from the hot and cold materials (above or below 0.158 MK) in the corona for all cases.

Figure 2 .
Figure 2. (a) Illustration of our 1D loop model.The locations of the heating sources are colored in brown for the loop apex and footpoints as in Equation (2).(b)Normalized temporal profiles of the heating rate for the apex f t (t) (blue solid line) and footpoints f b (t) (dotted lines) for the three different cases listed in Table1(orange, green, red).The decay timescales of the footpoint heating rate are 14.5, 29.5, and 59.5 minutes, respectively.The left and right portions depict the initial impulsive loop-top heating and the late-time gradual footpoint heating processes, respectively.They are discontinuous and have different scales on the horizontal axis.

Figure 3 .
Figure 3. Modeled flare atmosphere and synthesized TESS light curves.Cases 1-3 are shown in columns (a)-(c), respectively (fixed loop radius R L = 10 Mm, nominal footpoint heating timescale τ d,b = 29.5 minutes, and varying magnetic fields B = 400, 600, and 800 G).The top three rows show the evolution of electron number density n e , temperature T, and velocity v along the flare loop.The vertical axis is the loop coordinate l, and the horizontal axis is time t.The loop apex is located at the center l = 0; the footpoints are at the two ends.Positive (negative) velocities at positive (negative) loop coordinates indicate downflows, i.e., flows from the loop top to the footpoints, and are shown as red (blue).Time t = 0 marks the beginning of the heating input.The contours are for high-density region in the corona n e = 5 × 10 12 cm −3 .The dashed and dashed-dotted vertical lines in the second row indicate the end time of the loop-top and footpoint sources, respectively.The bottom row shows the synthesized TESS flare emission from the coronal plasma, normalized to the stellar disk emission.The total loop emission is shown in red, and the contributions from Thomson scattering, free-free emission, and free-bound emission are shown in blue, orange, and green, respectively.

Figure 6 .
Figure 6.Same as Figure 3, but for Cases 9 and 10 (R L = 10 Mm) with only footpoint and loop-top heating sources, respectively.The white arrow points out the region with catastrophic cooling and no plasma condensation.

Figure 8 .
Figure 8. Tracing a plasma parcel (initial position l = −13.65Mm at time t = 8 minutes) in Case 1 by integrating its velocity over time.Top the trajectory of the parcel overlaid on the velocity stack plot.Top right: evolution of temperature.Bottom left: evolution number density.Bottom right: the corresponding relative TESS flux, and contribution from free-free, free-bound, and Thomson scattering mechanisms.
The total relative emission δF/F is defined as δF Th /F + δF bf /F + δF ff /F.Figures 9-12 show the relative emission evolution defined in the above Equations (A1)-(A3) along the flare loop for all cases.The synthesized TESS light curves in Figures3-5are computed by integrating the aforementioned emissions within the corona region.The corona region is defined as the area between the two red dashed lines, which indicate the location where the temperature reaches 1 MK prior to the onset of heating.

Figure 9 .
Figure 9.The relative emission contribution along the loop and its evolution in time, which is evaluated from Equations (A1)-(A3).The area between the two red dashed lines indicates the corona portion.

Table 2
Late-phase Features from Simulations a Peak time of the late phase relative to the start time of flare heating.bDuration of the late phase where magnitude is greater than 10 −3 (approximate TESS error level).cTherelative energy contribution in TESSʼs optical range (corresponding to total flare energy in each case) of free-bound, free-free, and Thomson scattering emissions from the corona during the flare.