# The Accretion History of EX Lup:
A Century of Bursts, Outbursts, and Quiescence

Mu-Tian Wang (王牧天) School of Astronomy and Space Science, Nanjing University,
Nanjing, 210023, People’s Republic of China.  Key Laboratory of Modern
Astronomy and Astrophysics, Ministry of Education, Nanjing, 210023, People’s
Republic of China  Gregory J. Herczeg (沈雷歌) Kavli Institute for Astronomy and
Astrophysics, Peking University, No.5 Yiheyuan Road, Haidian District, Beijing
100871, People’s Republic of China  Department of Astronomy, Peking
University, No.5 Yiheyuan Road, Haidian District, Beijing 100871, People’s
Republic of China  Hui-Gen Liu (刘慧根) School of Astronomy and Space Science,
Nanjing University, Nanjing, 210023, People’s Republic of China.  Key
Laboratory of Modern Astronomy and Astrophysics, Ministry of Education,
Nanjing, 210023, People’s Republic of China  Min Fang (房敏) Purple Mountain
Observatory, Chinese Academy of Sciences, No.10 Yuanhua Road, Qixia District,
Nanjing 210023, People’s Republic of China  School of Astronomy and Space
Science, University of Science and Technology of China, Hefei, Anhui 230026,
People’s Republic of China  Doug Johnstone NRC Herzberg Astronomy and
Astrophysics, 5071 West Saanich Rd, Victoria, BC, V9E 2E7, Canada  Department
of Physics and Astronomy, University of Victoria, Victoria, BC, V8P 5C2,
Canada  Ho-Gyu Lee  Korea Astronomy and Space Science Institute, 776, Daedeok-
daero, Yuseong-gu, Daejeon, 34055, Republic of Korea  Frederick M. Walter
Department of Physics & Astronomy, Stony Brook University, Stony Brook NY
11794-3800, USA  Franz-Josef Hambsch Groupe Européen d’Observations Stellaires
(GEOS), 23 Parc de Levesville, 28300 Bailleau l’Evêque, France  Bundesdeutsche
Arbeitsgemeinschaft für Veränderliche Sterne (BAV), Munsterdamm 90, 12169
Berlin, Germany  Vereniging Voor Sterrenkunde (VVS), Oostmeers 122 C, 8000
Brugge, Belgium  Carlos Contreras Peña  Department of Physics and Astronomy,
Seoul National University, 1 Gwanak-ro, Gwanak-gu, Seoul 08826, Korea  Jeong-
Eun Lee Department of Physics and Astronomy, Seoul National University, 1
Gwanak-ro, Gwanak-gu, Seoul 08826, Korea  Mervyn Millward  York Creek
Observatory, George Town, Tasmania, Australia  Andrew Pearce  Skynet Robotic
Telescope Network, Perth Observatory, Western Australia  Berto Monard
Kleinkaroo Observatory, Calitzdorp, Western Cape, South Africa  Lihang Zhou
(周立杭)  Department of Astronomy, Peking University, No.5 Yiheyuan Road, Haidian
District, Beijing 100871, People’s Republic of China Mu-Tian Wang, Gregory
Herczeg, Hui-Gen Liu<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
(Received 2023 June 15; Revised 2023 August 3; Accepted 2023 August 16)

###### Abstract

EX Lup is the archetype for the class of young stars that undergoes repeated
accretion outbursts of $\sim 5$ mag at optical wavelengths and that last for
months. Despite extensive monitoring that dates back 130 years, the accretion
history of EX Lup remains mostly qualitative and has large uncertainties. We
assess historical accretion rates of EX Lup by applying correlations between
optical brightness and accretion, developed on multi-band magnitude photometry
of the $\sim 2$ mag optical burst in 2022. Two distinct classes of bursts
occur: major outbursts ($\Delta V\sim 5$ mag) have year-long durations, are
rare, reach accretion rates of $\dot{M}_{\rm acc}\sim 10^{-7}$ $M_{\odot}$
yr-1 at peak, and have a total accreted mass of around 0.1 Earth masses. The
characteristic bursts ($\Delta V\sim 2$ mag) have durations of $\sim 2-3$
months, are more common, reach accretion rates of $\dot{M}_{\rm acc}\sim
10^{-8}$ $M_{\odot}$ yr-1 at peak, and have a total accreted mass of around
$10^{-3}$ Earth masses. The distribution of total accreted mass in the full
set of bursts is poorly described by a power law, which suggests different
driving causes behind the major outburst and characteristic bursts. The total
mass accreted during two classes of bursts is around two times the masses
accreted during quiescence. Our analysis of the light curves reveals a color-
dependent time lag in the 2022 post-burst light curve, attributed to the
presence of both hot and cool spots on the stellar surface.

Light curves(918), Classical T Tauri stars(252), Stellar accretion(1578),
Eruptive variable stars(476)

##  1 Introduction

Outbursts of accretion onto young stars play a significant role in the
assembly of stellar mass (Fischer et al., 2022) and the chemistry of any
envelope and the protoplanetary disk (e.g. Lee, 2007; Jørgensen et al., 2015;
Hsieh et al., 2019). Eruptive young stellar objects (YSOs) are conventionally
classified by their photometric and spectroscopic characteristics into two
categories, FU Ori-type and EX Lup-type objects, after their namesakes. The FU
Ori-type outbursts reach accretion rates of $10^{-6}-10^{-4}$ $M_{\odot}$
yr-1, a factor of $10^{3}-10^{5}$ higher than the quiescent accretion rates of
$\sim 10^{-9}-10^{-8}$ $M_{\odot}$ yr-1, and can last for decades or even
centuries, while the smaller EX Lup-type outbursts have accretion rates of
$\sim 10^{-7}$ $M_{\odot}$ yr-1 that last for months (see reviews by Audard et
al., 2014; Hartmann et al., 2016). Follow-up spectroscopy is often used to
classify EX Lup-type outbursts from emission lines (e.g. Connelley & Reipurth,
2018), though modern searches reveal that many accretion outbursts have
characteristics intermediate between the EX Lup and FU Ori classes (e.g
Contreras Peña et al., 2017; Guo et al., 2021).

We focus on EX Lup, the namesake of a class of outbursts and a member of the
Lupus star-forming region (e.g. Alcalá et al., 2017) that is known for
repeated bursts. Monitoring since 1893 revealed large outbursts in 1945, 1955,
and 2008 (e.g. McLaughlin, 1946; Herbig et al., 1992; Aspin et al., 2010),
marking bookends for the studies of EX Lup by George Herbig (Herbig, 1950,
2008). Interspersed between these outbursts were many smaller ($\Delta V\sim
2-3$ mag) bursts that lasted for at most a few months (Herbig et al., 2001;
Herbig, 2007). The accretion rates have been estimated to range from
$10^{-8}-10^{-10}$ $M_{\odot}$ yr-1 in quiescence (e.g. Sipos et al., 2009;
Sicilia-Aguilar et al., 2015; Alcalá et al., 2017) and from $10^{-7}-10^{-6}$
$M_{\odot}$ yr-1 for the large outbursts (e.g., Aspin et al. 2010; Juhász et
al. 2012, see also Giannini et al. (2022) for other sources), but with
differences between studies of 1–2 orders of magnitude.

Measurements of the accretion rates during both quiescence and outburst
provide a quantified mass measurement that can be used to constrain the
instability physics. For EX Lup, the outbursts are likely caused by
instabilities in the innermost disk or in the magnetic star-disk connection
(e.g. D’Angelo & Spruit, 2010; Armitage, 2016). The largest outbursts can also
have a significant impact on the disk structure and chemistry. The 2008
outburst led to the depletion of gas in the inner disk (Banzatti et al.,
2015), moved snowlines to larger radii (White et al., 2020), and created
comet-forming materials (Ábrahám et al., 2009). Organic molecules reformed
after destruction during the outburst, while crystalline grains that formed
during the burst were transported outward (Kóspál et al., 2023).

The most recent burst of EX Lup, in 2022 and initially identified by Zhou &
Herczeg (2022) and Kóspál et al. (2022), was intensely monitored with multi-
band photometry, including by astronomers who contributed to the database of
the American Association of Variable Star Observers (AAVSO). The exquisite
light curve provides an opportunity to study the changes in brightness and
mass accretion rate. Correlations between brightness and accretion rate during
the 2022 burst aid further estimates of historical accretion rates from the
long-term photometric light curve. Our results are mostly similar to the
contemporaneous and independent analysis of the 2022 burst and historical
light curve by Cruz Sáenz de Miera et al. 2023.

In this paper, we re-evaluate the accretion history of EX Lup based on
correlations developed from the burst in 2022. As the archetype of its class
and with historical monitoring that dates back to 1893, the accretion history
can be reconstructed for over a century. With the improved knowledge about the
variability link between mass accretion rate and photometry, we return to the
century-long photometry data of EX Lup to measure the mass accreted in bursts
and compare that to the mass accreted during quiescence. In Section 2 we
describe the historical and modern photometry and spectroscopy of EX Lup, and
review the stellar properties of EX Lup in Section 3. In Section 4 we provide
a photometric analysis of EX Lup’s historical and modern light curves. In
Section 5 we describe our method to translate multi-wavelength photometry to
mass accretion rates. In Section 6 we estimate the accreted masses during
burst and in quiescent.

##  2 Historical and Recent Light Curves and Spectra

This section summarizes the photometric and spectroscopic data we obtained for
analysis in this paper. The photometric data is shown in Figure 1.

Figure 1: The full light curve collections of EX Lup. The photometry data
covered in this paper come from DASCH (blue petagons), AAVSO (black dots),
ASAS-SN (red and blue diamonds), NEOWISE (black and orange triangles), TESS
(orange squares), and SMARTS/1.3-m telescope (red and green crosses). The
timestamps of several outbursts and bursts from EX Lup are also indicated in
the figure.

###  2.1 DASCH light curves

We obtain the EX Lup $B$-band light curves available between 1893 and 1954
from the DASCH database. The Digital Access to a Sky Century@Harvard (DASCH,
Laycock et al., 2010; Tang et al., 2013) is a project to digitize historical
data from a collection of photographic plates hosted by the Harvard College
Observatory between 1885 to 1992. The plate-derived photometry was further
calibrated by AAVSO Photometric All-Sky Survey (APASS Zacharias et al., 2013)
and is converted into Johnson B magnitudes.
111http://dasch.rc.fas.harvard.edu/lightcurve.php

###  2.2 AAVSO Light Curves

American Association of Variable Star Observers (AAVSO) database has a wealth
of photometry of EX Lup from 1955 to 2022. We obtained the full collection of
EX Lup photometry from the AAVSO website222https://www.aavso.org/ on Nov. 10
2022. The majority of AAVSO photometry is in the format of Visual magnitude
before 2008, sourced from the Royal Astronomical Society of New Zealand
(RASNZ) and in the non-modern format of the visible band, as observed by
different observers, most especially Albert Jones in the Visual band before
1970 and between 1987-2007.

Herbig et al. (1992) transformed the 1955-1956 Visual photometry from Jones to
the Johnson V band by

$V=0.895\times Vis.+0.63.$ (1)

We extrapolate this relation to photometry from Jones between 1955-1993. After
1994, we infer that Jones reported photometry in the standard Johnson V
magnitude, because the minimum visual magnitude changed from $\sim$14 to
$\sim$13, close to the minimum V band magnitude (13.25) mag observed by ASAS-
SN.

The recent 2022 burst was also covered with frequent full-length multi-color
$UBV(RI)_{c}$ photometric monitoring. Based on the record from AAVSO, the
campaign starts on 2022-03-13 and ends on 2022-10-17. The U-band photometry,
particularly powerful for accretion rate measurements (Gullbring et al.,
1998), was reported along with $BV(RI)_{c}$ by Franz-Josef Hambsch; 14 other
observers also contributed to this light curve. Thanks to contributions from
all the observers, the seven-month multi-color photometry monitoring is kept
on a daily cadence, with very few one-day gaps.

###  2.3 ASAS-SN Light Curves

The All-Sky Automated Survey for Supernovae (ASAS-SN) catalog (Shappee et al.,
2014; Kochanek et al., 2017) provides $V$-band photometry from 2016-3-10 to
2018-9-17 and $g$ band photometry from 2018-6-13 to present, downloaded on
2023-07-17 from the ASAS-SN website333https://asas-sn.osu.edu/. On a typical
night, ASAS-SN will observe EX Lup three times with 90-second exposures, but
sometimes there are cases where only one or two exposures are available
because of the weather. These differences have a negligible effect on the
photometric analysis presented in Section 4.2.

The ASAS-SN camera has a field of view of $4.5\times 4.5^{\circ}$ and a pixel
size of 8′′. The contamination by other stars is negligible within the large
aperture. One source separated by $8\farcs 7$ has Gaia $G=16.4$ mag, roughly 4
mag fainter than EX Lup, while all others are fainter than 18 mag (Gaia
Collaboration et al., 2021).

###  2.4 TESS Photometry

EX Lup was observed by TESS in Sector 12 and 39 with a 20-second cadence. We
downloaded the SPOC-processed TESS light curves of EX Lup in these two
sectors, and convert the TESS flux in units of electron per second $f_{e}$ to
TESS magnitude $T=-2.5\log(f_{e})+20.44$ from the formula in TESS Instrument
Handbook. Although the pixel scale of TESS is 21′′, EX Lup dominates the total
fluxes within the one-pixel range. The TESS filter extends from $\sim
6000-10000$ Å.

###  2.5 NEOWISE Photometry

The mid-IR photometry arises from all-sky observations of the WISE telescope.
WISE surveyed the entire sky from 2010 to 2011 in bands W1 (centered at 3.4
$\mu$m) and W2 (centered at 4.6 $\mu$m), with spatial resolutions of 6.1″ and
6.4″, respectively (Wright et al., 2010; Mainzer et al., 2011). In 2013
September, WISE was reactivated as the NEOWISE-reactivation mission
(NEOWISE-R, Mainzer et al., 2014). NEOWISE-R is still operating and the latest
released data set contains observations until mid-December 2022. For each
visit to a particular area of the sky, WISE performs several photometric
observations over a period of $\sim$few days. Each area of the sky is observed
in a similar way every $\sim$ 6 months. The data presented were obtained
between 2010 and 2022 and downloaded from the single-epoch catalogues stored
at the NASA/IPAC Infrared Science Archive (IRSA) catalogues (NEOWISE Team,
2020), using a 3″ radius from the coordinates of the YSO.

###  2.6 Optical SMARTS/Andicam Photometric Monitoring

EX Lup was observed in Johnson $B$ and $V$ band by SMARTS/Andicam on the 1.3-m
telescope (Subasavage et al., 2010) from 2004-01-27 to 2005-06-29 in
quiescence on a weekly cadence, and from 2008-01-22 to 2009-07-30 during
outburst on nightly cadence. In total, SMARTS/Andicam obtained 297 $B$-band
observations on 292 nights and 286 $V$-band observations on 286 nights.

###  2.7 Optical-near-IR VLT/X-Shooter spectrum

EX Lup was observed by VLT/X-Shooter on UTC 2010-05-04 03:44:51.779 as part of
program 085.C-0764 (PI H.M. Günther) and was initially published by Alcalá et
al. (2017). The X-Shooter spectrum is flux-calibrated and covers from
3000–25000 Å. The synthetic magnitudes derived for the X-Shooter spectrum are
Johnson U=13.927, V=12.956, B=14.087, SDSS g = 13.633 mag, 2MASS J = 9.80 mag,
after applying filter curves from the Spanish Virtual Observatory (Rodrigo et
al., 2012; Rodrigo & Solano, 2020). The estimated extinction-corrected and
accretion-subtracted magnitude of the photosphere is U = 15.474, V=13.055,
B=14.535, SDSS g = 13.892 mag, 2MASS J = 9.89 mag with typical uncertainties
of $\sim 0.05$ mag. At $J$-band, this indicates a veiling of $\sim 8$% during
quiescence, similar to $J$-band veiling measurements for other accreting low-
mass stars (Fischer et al., 2011; Sousa et al., 2023).

###  2.8 Optical Keck/HIRES spectrum

We obtained a flux-calibrated high-resolution spectrum of EX Lup with Keck
I/HIRES on MJD 54609.416 (2008 May 23), during the large burst in 2008,
originally published in Aspin et al. (2010). Fluxes were calibrated using two
spectra of the spectrophotometric standard EG 274 (Hamuy et al., 1992), both
obtained within $1$ hr of the EX Lup spectra and bracketing the airmass of the
EX Lup observation.

###  2.9 Near-IR Gemini South/IGRINS spectrum

We obtained a high-resolution ($R\sim 45,000$) spectrum of EX Lup from 1.4–2.5
$\mu$m using Immersion Grating Infrared Spectrograph (IGRINS, Park et al.,
2014) at Gemini South on 2022-03-14 (MJD 59652.36), obtained in GS-2022A-Q-109
(PI H-G Lee) through the Korean partnership in Gemini Observatory. The
observations were obtained in 2$\times$ABBA pattern with 131-second
integrations at each step. The A0 star HIP78445 was observed as a standard
star to correct for the telluric absorption. The data were reduced using the
standard IGRINS pipeline (Lee et al., 2017) with a modified version of optimal
extraction (Horne, 1986).

The flux in the EX Lup spectrum was calibrated from the relative count rate in
the EX Lup and telluric standard star. For both EX Lup and HIP 78445, the slit
viewer indicates that the stars were well centered. We also corrected for
differences in slot losses by measuring the spatial profile in the cross-
dispersion direction and assuming a Gaussian point-spread function. From this
process, we measure $K=8.07$ with an uncertainty of $\sim 0.1$ mag.

Figure 2: The best-fit (red) of EX Lup X-Shooter spectrum (black) taken on May
4, 2010. The spectrum is dereddened with $A_{v}=$ 0.1 mag, fitted with a
photosphere template (blue) and a slab model (grey).

Figure 3: The Balmer Jump (top) and U-band spectra (bottom) for EX Lup in
quiescence (red) and outburst (black). The accretion spectrum of EX Lup during
outburst has a bluer slope and maintains a Balmer Jump. The U-band spectrum
during the outburst is highly structured, indicating that the accretion flow
has an optically thick photosphere. In quiescence, the U-band spectrum of EX
Lup (and most other accreting young stars) is mostly featureless.

##  3 Properties of EX Lup

In this section, we re-evaluate the properties of EX Lup, including revising
the extinction from $A_{V}=1.1$ mag in Alcalá et al. (2017) to $0.1$ mag based
on an analysis of the same X-Shooter spectrum. A more complete analysis of the
X-Shooter spectrum and our revised parameters is presented in Appendix A.

The optical spectrum of EX Lup in quiescence is the combination of a cool
stellar photosphere, including TiO bands in absorption, and excess accretion
continuum and line emission. To determine the accretion components in the
spectrum, we test several non-accreting, low-mass YSO spectra to serve as
photospheric templates of EX Lup (Manara et al., 2013). We find the most
appropriate photospheric template from Manara et al. (2013) and extinction
that yields an accretion continuum spectrum that is most free of photospheric
features (see Appendix for details).

Our final measurements are a spectral type of M0.5 (corresponding to $3810$ K
using the conversion of Herczeg & Hillenbrand 2014) and extinction $A_{V}=0.1$
mag, determined by extinction curve of $R_{V}=3.1$ (Cardelli et al., 1989).
This fit yields a residual (accretion) spectrum with no evident molecular band
absorption features but with fluctuations that are larger than seen for many
other accreting young stars. The dereddened EX Lup spectrum, the photospheric
spectrum template (TWA 14), and the subtracted accretion spectrum are
displayed in Figure 2.

In contrast to our spectral fits, Alcalá et al. (2017) measured M0 and an
extinction $A_{V}=1.1$ mag, estimated using a K7 template. The lower
extinction from our work is consistent with the expectation of low extinction
from the non-detection of DIBs (see correlations in Carvalho & Hillenbrand,
2022) in both in-burst and quiescent spectra. The spectral type is roughly
consistent with the effective temperature measurement of $3750$ K by Sipos et
al. (2009).

When combined with the distance of $154.72\pm 0.36$ pc (Gaia Collaboration et
al., 2022), the extinction- and veiling-corrected photospheric $J$-band
magnitude of 9.89, and bolometric correction from Pecaut & Mamajek (2013) for
an M0.5 star, we obtain a luminosity of $0.43$ L⊙ and radius $1.50$ R⊙. In the
Somers et al. (2020) pre-main sequence models, these parameters correspond to
a mass of 0.83 M⊙ and age of 4.6 Myr for 50% spot coverage, which is adopted
here. The 0% spot models correspond to 0.55 M⊙ and 2.3 Myr.

###  3.1 The spectroscopic accretion rate during quiescence

From the X-Shooter spectrum, the accretion luminosity is $0.0256$ L⊙, as
measured from the excess continuum from 3500–5000 Å. The bolometric correction
is calculated from a plane-parallel slab model of Valenti et al. (1993), as
implemented by Herczeg & Hillenbrand (2008). This accretion luminosity leads
to an accretion rate of $1.77\times 10^{-9}$ M⊙ yr-1.

The accretion luminosity measured here is 8 times smaller than the 0.2 L⊙ and
the accretion rate is 20 times smaller than the $3.9\times 10^{-8}$ from
Alcalá et al. (2017), using the exact same spectrum. However, most of this
difference is attributed to the higher extinction of $A_{V}=1.1$ mag in Alcalá
et al. (2017), which changes the accretion luminosity by a factor of $\sim
4.3$ and the radius estimate from red-optical emission by $\sim 1.5$. The
larger pre-Gaia distance of 200 pc, adopted by Alcalá et al. (2017), accounts
for a factor of 1.7 in accretion luminosity and 1.3 in stellar radius.

###  3.2 The accretion rate during the 2008 outburst

The optical spectrum of EX Lup in late May 2008, near the end of the large
outburst (Figure 3), consists of a continuum with a forest of emission lines.
The emission line spectrum, also reported by Aspin et al. (2010) and Sicilia-
Aguilar et al. (2012), is similar to the optical spectra of other heavily
veiled, high accretion rate young stars (e.g. Hamann & Persson, 1992; Gahm et
al., 2008). During quiescence, these lines are weak but detectable, after
subtracting photospheric templates (Sicilia-Aguilar et al., 2015). The
accretion continuum during the burst spectrum of 2008 is much bluer than the
accretion continuum during quiescence, suggestive of the outburst producing a
hotter temperature in the accretion shock region where the optical emission is
produced.

The continuum emission at $<3700$ Å is brighter than that at longer
wavelengths, consistent with expectations for the Balmer jump that
characterizes magnetospheric accretion (e.g. Calvet & Gullbring, 1998). When
modeled as H continuum emission, the accretion luminosity is $2.00$ L⊙,
leading to a mass accretion rate of $1.5\times 10^{-7}$ M⊙ yr-1. The stellar
photosphere is no longer detectable, so all of the optical emission is now
attributed to accretion.

While typical accreting young stars, including EX Lup in quiescence, have
Balmer continuum spectra that are mostly featureless, the 2008 burst spectrum
features a complex array of absorption features short of the Balmer jump. The
accretion flow seems to be sufficiently optically thick for photosphere-like
absorption lines to appear. An analysis of the detailed spectral features is
beyond the scope of this paper.

###  3.3 Spectroscopic classification of bursts and outbursts

Spectroscopy is historically used to discriminate between types of bursts.
Bona fide EX Lup-type outbursts show strong CO overtone emission, indicating a
warm disk surface, along with strong Br$\gamma$ emission produced by
magnetospheric accretion (see, e.g., Lorenzetti et al., 2012; Connelley &
Reipurth, 2018; Contreras Peña et al., 2023; Ghosh et al., 2023). Weaker
accretion bursts should produce Br$\gamma$ emission but are not expected to
show CO emission; in stronger bursts the magnetospheric accretion paradigm may
break down (e.g. Hartmann & Kenyon, 1996).

The major outburst of EX Lup in 2008, the archetype of EX Lup-type outbursts,
featured strong CO overtone emission and Br$\gamma$ emission (e.g. Aspin et
al., 2010; Kóspál et al., 2011). The strength of the CO emission was much
weaker at the end of the burst than at the beginning (Aspin et al. 2010, see
also changes in fundamental emission in Goto et al. 2011), either because of a
decrease in the gas in the inner disk (Banzatti et al., 2015) or a decrease in
heating of the disk surface.

The 2022 burst shows CO in absorption, consistent with the classification as a
characteristic burst rather than a major EX Lup-type outburst. While the depth
of CO bands is shallower during the 2022 burst than during quiescence in 2010,
this change is explained by a veiling increasing by 0.65 at 2.3 $\mu$m. This
increase in veiling is more than expected from the synthetic K-band brightness
(8.46 in 2010 compared with 8.31 during the 2022 burst), even accounting for
stronger dust emission and 10% fainter photospheric emission at 2.3 $\mu$m
than at 2.15 $\mu$m. Although no CO emission is detected, we cannot rule out
the presence of weak CO emission on top of the photospheric absorption.

##  4 Analyzing Photometry from Old and New Bursts

In this section, we begin by describing the historical brightness of EX Lup,
following previous descriptions (e.g. McLaughlin, 1946; Herbig et al., 1992;
Herbig, 2007). We then update this description for modern monitoring campaigns
by ASAS-SN and by other astronomers, followed by the detailed analysis of a
burst in 2022.

###  4.1 Overview of the Historical Lightcurve of EX Lup

The photometric monitoring of EX Lup began with plate photometry by Harvard
Observatory in 1890, which is displayed in the top panel of Figure 1. From
1890-1920, EX Lup had a minimum brightness of 14 mag in the B-band,
accompanied by irregular fluctuations between 14-11 mag (McLaughlin, 1946). As
the observational cadence increased between 1920 and 1955, the photometric
behavior of EX Lup became more apparent. A long-term brightening event
($\Delta B\sim 2$ mag) began in 1937 and lasted for 3-4 years, followed by a
year-long outburst with $\Delta B\sim 4$ in 1945. This photometric variability
pattern repeated again between 1947-1956, with similar magnitude change and
durations of long-term brightening and outburst compared to the previous
cycle, including another large outburst in 1955 (Herbig et al., 1992).

In the decades after 1955, EX Lup was more frequently observed by the AAVSO.
From 1955–1975, EX Lup had a base level of $V=13.5$ mag, with sporadic bursts
with $\Delta V\sim 2$ mag and durations of less than a year (Bateson et al.,
1991). The bursts were then rare from 1975-1993 as seen in third panel of
Figure 1. Activity resumed with seven bursts between 1993–2002, each with a
duration and photometric amplitude similar to the bursts in the 1960s and
intervals of around 1-2 years (Lehmann et al., 1995; Herbig et al., 2001;
Herbig, 2007). Another large outburst occurred in 2008, with a 300-day
duration and 5 magnitudes change in the V-band (Aspin et al., 2010; Juhász et
al., 2012).

The historical photometry from DASCH and AAVSO reveals a complex ensemble of
photometric variability from EX Lup. When in quiescence, EX Lup has a night-
to-night range in the variability of $<0.6$ mag in $g$-band, typical of other
accreting young stars (e.g. Venuti et al., 2021; Wang et al., 2023). The
bursts with amplitudes of $\sim 2$ magnitudes and durations of a few months
are more common than the outbursts with amplitudes of 4-5 magnitude and
durations of $\sim 6$ months.

###  4.2 Photometric Variability Analysis From Recent Light Curves

Figure 4: Top: NEOWISE $W_{1}$ (orange) and $W_{2}$ (black) photometry.
Middle: The ASAS-SN $g$\- (red) and $V$-band (blue) light curves of EX Lup
between 2015 and 2022, shown in blue and red dots, respectively. The bursts in
2020 and 2022 are indicated with black sticks. The two gray vertical lines
indicate the epochs when two sectors of TESS observations are taken. Bottom:
the Lomb-Scargle periodograms for $V$-band (left) and $g$-band (right) light
curves. The $g$-band light curves in 2020 and 2022 are excluded from
periodogram analysis due to the presence of bursts. Figure 5: Top: the TESS
light curves in 2019 and 2021. Contemporaneous ASAS-SN $g$-band light curves
are also plotted with offset. For the clarity of display, $g$-band photometry
taken within one day is binned as one, and TESS light curves are binned to a
10-minute cadence. Bottom: the Lomb-Scargle Periodogram for TESS during two
sectors (each $\sim$ 24 days) and ASAS-SN $g$-band light curves in 2019-2020
and 2021-2022.

The ASAS-SN light curve taken between 2015 and 2023 traces almost a decade of
the behavior of EX Lup with a near-nightly cadence. The $V$-band light curve
between 2015 and 2018 shows a slow brightening of 0.1 mag in mean brightness
and no bursts. Two consecutive bursts are also observed in 2020, though with a
sparse cadence. The larger 2022 burst is well-monitored in $g$-band. The
$g$-band light curve during quiescent periods of 2021 has a median value that
is 0.2 magnitude brighter compared to the light curve in 2019.

TESS photometry allows us to compare the light curves before and after the
2020 bursts in greater detail (Figure 5). The TESS light curves in 2019 show
periodic modulations with a peak-to-peak amplitude of 0.4 mag. The $g$-band
light curves also follow a periodic modulation with an amplitude of 0.5 mags.
After the 2020 burst, the periodic modulation is still present in the TESS
light curves, but additional stochastic bursts emerge in the brightness
minima. The $g$-band photometry is brighter in 2021 than in the 2019 TESS
visit, indicating a higher accretion rate. As a result, the 2021 TESS
lightcurve has a higher minimum brightness, with accretion filling in the
minimum brightness from the photosphere, and shows several small, readily
detected bursts.

Figure 6: Left: the ASAS-SN $g$ and AAVSO multi-color light curves of 2022
burst. The light curves are offset with different values for clarity. Right:
Lomb-Scargle periodogram of post-burst (after MJD=59725.0) light curves, the
color scheme is the same as the left plot. The 7.417-day periodicity
discovered in radial velocity measurements is also indicated (Kóspál et al.,
2014).

The top panel of Figure 4 shows the contemporaneous NEOWISE W1 and W2 mid-
infrared photometry. The mid-IR emission brightened with the 2020 and 2022
bursts, although the mid-IR emission was also brighter by 0.3 mag in epochs in
2017 and 2018 without any evident bursts in ASAS-SN light curves. We further
discuss the connection between mid-IR emission and mass accretion rates in
Section 5.2.

A Lomb-Scargle Periodogram analysis of the ASAS-SN lightcurve shows a peak of
7.415$\pm$0.024 days for the $V$-band and 7.411$\pm$0.051 days for the
$g$-band (bottom panels of Fig. 5). We found no evidence of a change in
periodicity in the light curves taken each year, though the power of the peak
periodicity depends on the level of burst activities in each year. For 2019,
the Lomb-Scargle periodogram reveals a periodicity of 7.420$\pm$0.103 for the
$g$-band and 7.394$\pm$0.826 day for the TESS light curve. In 2021, the peak
periodicity of TESS light curves is $7.321\pm 0.836$ days, similar to that in
2019. The $g$-band peak periodicity remains around $7.332\pm 0.172$ days, but
the Lomb-Scargle power for this peak decreases, and a new peak emerges at
around 10 days, which is not seen from the periodograms of earlier $g$-band
light curves. This change might be the result of elevated accretion activity
in 2021.

Figure 7 shows the time evolution of spot modulation light curves. The ASAS-SN
and TESS light curves between 2016 and 2023 are folded to $P=7.42$ day and
JD=2457000. We also include a segment of $V$-band light curves taken between
2004 and 2005. The spot modulation profiles remain steady and with a small
scatter between 2017 and 2019. After the 2020 burst, the $g$-band photometry
brightens by 0.2 and the brightness minima has a larger scatter. This larger
scatter is consistent with the stochastic bursts seen at the minimum of
modulation during the TESS observations in 2021 (Figure 5). After the 2022
burst, the modulation phase changed, and the original phases of brightness
minima turned into maxima. We will further discuss this in Section 6.4.

Figure 7: The phase-folded light curves of ASAS-SN $V-$ (red), $g-$band
(blue), and TESS (black) light curves, all folded to $P=7.42$ day and
JD=2457200.

###  4.3 The 2022 Burst

Figure 8: Phase-folded post-burst light curves in UBVRI and $g$ band from
AAVSO. The multi-color light curves are phase-folded with P$=7.385$ day and to
the U-band best-fit reference date JD$=2459716.96$. The red vertical lines
indicate the location of brightness maxima for different wavelengths to show
the phase lags. The blue errorbar at the lower left of each panel shows the
median photometric errors in each band. Figure 9: 2022 burst on the color-
magnitude diagram (CMD). The color is calculated using photometry taken at
epochs separated by no more than one day. The data points are color-coded by
the time since JD=2459600.0, the time when the burst reached its brightness
maximum. The fainting of burst proceeded until JD=2459725.0, after which the
variability is dominated by the spot modulations. The dotted black line is the
color trend due to accretion enhancement, on the scale of one to fifteen times
the accretion level inferred from the X-Shooter spectrum. The red line is
interstellar reddening law of $A_{V}=0.5$ mag and $R_{V}=3.1$.

EX Lup was reported to undergo a brightness rise from ASAS-SN $g$-band light
curves early in 2022 (Zhou & Herczeg, 2022; Kóspál et al., 2022), triggering
intense follow-up by astronomers who reported results to AAVSO. The collection
of multi-color light curves of the 2022 burst is displayed in the left panel
of Figure 6. On 2022-03-25 (MJD=59664.0), EX Lup reached its maximum, from
where it faded until June 2022, when the $g$ photometry returned to pre-burst
brightness. A tentative dip around 2022-08-08 (MJD=59800.0) and a slow fade
after 2022-09-17 (MJD 59840.0) are also seen in light curves.

The following two subsubsections analyze the variability through periodicity
seen through the burst and the changes in color versus brightness. The
periodicity during the burst is an accretion hot spot (see also Sicilia-
Aguilar et al., 2015), further modeled and discussed in Section 6.4. The color
changes versus brightness establish correlations with accretion rate, which we
then apply to historical light curves.

Figure 10: Top: The ASAS-SN $g$ and NEOWISE W1/W2 photometry before and after
the 2022 burst. Bottom: the color-magnitude diagram of W1 and W1-W2. The color
trend predicted by the extinction of $A_{V}=4$ mag is indicated with a black
arrow.

####  4.3.1 Spot Modulation of the lightcurves

Spot modulation is detected during and after the 2022 burst, as first reported
during the rise by Kóspál et al. (2022). The post-burst light curves after
2022-05-25 (MJD 59725.0) show a strong periodicity peak in Lomb-Scargle
Periodograms (right panel of Figure 6). The in-burst light curves before
2022-05-25 (MJD 59725.0) also display spot modulations.

We fit the in-burst and post-burst light curves of each band with a sinusoidal
function to determine the spot period ($P_{\rm sin}$), amplitude, and
reference date ($T_{0}$). Before this fit, we removed the fading trend with a
linear fit (in magnitudes). For the post-burst light curve, the fades in the
light curves between MJD of 59780.0 and 59810.0 and after MJD=59840.0 are
excluded.

The least-square fits of a sinusoid to the in-burst light curve yield a mean
period of 7.55 days, which then decreases back to a mean of 7.385 days in
post-burst light curves, both from multi-band fits (summarized in Table 1). We
also include the peak periods $P_{\rm LSP}$ from the Lomb-Scargle periodogram
(LSP) in Table 1. The mean periods derived from in-burst light curves are 0.1
day longer than the periods measured during quiescence but within
uncertainties; the longer measured periods are probably due to the larger
photometric scatter and shorter observation span of in-burst photometry.

The spot amplitudes derived from in-burst light curves range from 0.11 to 0.27
mag, while the amplitudes from post-burst light curves range from 0.1 to 0.33
mag. The in-burst spot amplitudes are generally smaller than those from post-
burst. The amplitudes derived from shorter wavelengths are larger, as expected
from the presence of cool or hot spots on the stellar surface (Herbst et al.,
1994). However, we measure differences in the reference date of the sinusoidal
model for multi-color light curves, which suggests color-dependent phase lag
in spot modulation between different wavelengths. To remove the bias of
different periods, we re-fit the multi-color light curves but fix the periods
at 7.55 days for in-burst light curves and 7.385 days for post-burst light
curves, with reference dates listed in Table 1. Both in-burst and post-burst
light curves show a longer lag towards longer wavelengths. For example, the
reference date difference between the U band and the I band in post-burst
light curves is $\sim$ 1.5 days with a significance of 12.6 $\sigma$, while in
the in-burst light curve, the $U-I$ phase lag is 1.25 days and has a lower
significance of 5.3 $\sigma$.

To better visualize this phase lag, we present the phase-folded post-burst
light curves in Figure 8. The phase lag between light curve peaks from
different observing wavelengths is also reported by Kóspál et al. (2014), who
used photometry from the optical V-band to near-infrared _Spitzer_ band,
though that was only a tentative signature owing to the sparse photometry
coverage. We will further discuss this phase lag feature in Section 6.4.

Table 1: Sinusoidal Fits to Spot Modulations in 2022 Burst

Band | PLSP | Psin | Amplitude | $T_{0}$
---|---|---|---|---
| (day) | (day) | (mag) | (JD-2459000.0)
|  | In-BurstaaBetween MJD-59670.0 and MJD-59725.0 |
U | $7.53\pm 0.58$ | $7.54\pm 0.14$ | $0.27\pm 0.05$ | $717.08\pm 0.19$
B | $7.52\pm 0.53$ | $7.52\pm 0.11$ | $0.25\pm 0.04$ | $717.30\pm 0.16$
g | $7.57\pm 0.51$ | $7.51\pm 0.11$ | $0.27\pm 0.04$ | $717.29\pm 0.18$
V | $7.56\pm 0.48$ | $7.56\pm 0.08$ | $0.26\pm 0.03$ | $717.48\pm 0.11$
R | $7.62\pm 0.51$ | $7.52\pm 0.08$ | $0.18\pm 0.02$ | $717.71\pm 0.13$
I | $7.58\pm 0.54$ | $7.61\pm 0.07$ | $0.12\pm 0.01$ | $718.34\pm 0.13$
|  | Post-BurstbbBetween MJD-59725.0 and MJD-59875.0, but the segments between MJD-59780.0 and MJD-59810.0 is excluded. |
U | $7.39\pm 0.15$ | $7.39\pm 0.02$ | $0.33\pm 0.03$ | $717.08\pm 0.10$
B | $7.40\pm 0.16$ | $7.40\pm 0.02$ | $0.25\pm 0.02$ | $717.34\pm 0.08$
g | $7.42\pm 0.17$ | $7.39\pm 0.03$ | $0.18\pm 0.03$ | $717.45\pm 0.16$
V | $7.39\pm 0.44$ | $7.39\pm 0.02$ | $0.14\pm 0.01$ | $717.89\pm 0.08$
R | $7.37\pm 0.17$ | $7.36\pm 0.02$ | $0.12\pm 0.01$ | $718.40\pm 0.09$
I | $7.38\pm 0.52$ | $7.38\pm 0.02$ | $0.10\pm 0.01$ | $718.63\pm 0.07$

Note. — $P_{\rm sin}$ is the period determined from sinusoidal fits, and
$P_{\rm LSP}$ is from Lomb-Scargle periodogram analysis. The reference date
$T_{0}$ is obtained from sinusoidal modeling by fixing P = 7.385 days for
post-burst and P = 7.550 days for in-burst light curves.

####  4.3.2 Color-Magnitude Diagrams

The multi-color photometry during the fade of the 2022 burst is shown as
color-magnitude diagrams in Figure 9. In general, EX Lup gets fainter because
the accretion continuum is weaker while the spot-averaged photosphere does not
change. Since the accretion continuum is hotter and has a bluer color than the
photosphere, EX Lup appears brighter at bluer colors. However, during
accretion bursts, the photosphere is $10-40$ times fainter than the accretion
continuum in $U$ and $B$ bands and thus contributes little to the $U-B$ color.
The behavior during such bursts depends primarily on the shape of the
accretion continuum, and thus change in $U-B$ is small as accretion decreases.

With the optical spectroscopy described in Section 3.1, we compare the
observed color-magnitude trends with the ones synthesized from the accretion
model fitted from X-Shooter spectra (see Section 3.1), marked by the black
dotted line in Figure 9. The predicted color trend is obtained by increasing
the accretion fluxes from one to fifteen times that of the base, so that the
predicted magnitude ranges are comparable to the observations. The predicted
trends due to increasing accretion agree well with observed $B-V$, $V-R$, and
$R-I$ color. However, the observed $U-B$ colors are redder than predicted
during the burst, constant in color at $U<13$, indicating that the in-burst
color of accretion at short wavelengths is redder than that measured from the
quiescent emission. This $U-B$ behavior indicates that the Balmer jump is
small, as expected when the accretion flow is optically thick (see for example
the models of Calvet & Gullbring 1998, where the optically thin pre-shock gas
has a larger Balmer jump than the optically thick post-shock gas).

We also analyze the mid-infrared color behavior from WISE/NEOWISE, including
epochs that cover before, during, and after the 2022 burst. The ASAS-SN $g$,
the NEOWISE W1/W2 photometry, and the CMD are shown in Figure 10. In the first
two epochs in 2021, EX Lup is quiescent and its location on the color-
magnitude diagram is localized. In the epoch 20 days before the 2022 burst
peak, W1 is 0.4 mag brighter with a similar W1-W2 color. Fifty days after the
end of the burst, EX Lup is fainter in both W1 and W2 bands and slightly bluer
than the pre-burst color, indicating some minor change in either the inner
disk or in contributions from the accretion flow. The trace of color
variability does not follow the prediction of extinction events.

##  5 Accretion Rates from Light Curves

In this section, we develop methods to convert photometry into accretion
rates. We then calculate accretion rates from the past 10 years of high-
quality monitoring, including the 2022 burst, and then assess accretion rates
during the historical major outbursts and characteristic bursts. Finally, we
compare our accretion rates, obtained from excess accretion continuum
emission, to accretion rates using other diagnostics and past accretion rates
from the literature.

###  5.1 Conversion Between Excess Emission and Accretion Luminosity

For EX Lup, as with most other accreting young stars, the accretion emission
dominates over the photosphere at short wavelengths, especially at
$\lambda<4200$ Å (see Figure 2 and also spectroscopic fits in, e.g., Pittman
et al. 2022). Among the most common photometric bands, the $U$-band best
diagnoses the relative change in accretion level, although some correction for
photospheric contribution is still needed. We first calculate the accretion
rate from $U$-band photometry and then use the $U$-band measurements to
establish conversions for $B$, $g$, and $V$ bands.

The excess $U$-band flux due to accretion is converted to accretion luminosity
using the relationship developed by Gullbring et al. (1998),

$\log(L_{\rm acc}/L_{\odot})=1.09^{+0.04}_{-0.18}~{}\log(L_{\rm
U}/L_{\odot})+0.98^{+0.02}_{-0.07}.$ (2)

The U-band luminosity of EX Lup, $L_{\rm U}$, is calculated as

$L_{\rm U}=4\pi d^{2}F_{\rm U,zero}\Delta\lambda_{U}10^{-0.4U_{\rm ex}},$ (3)

where $d=154.72$ pc is the distance of EX Lup, $F_{\rm U,zero}$ is the $U$
band zero magnitude flux, $\Delta\lambda_{U}=680$ Å is the approximate width
of the $U$ band filter, and $U_{\rm ex}$ is the $U$ band magnitude, corrected
for extinction and after subtracting the underlying photospheric flux.

Leveraging the 2022 burst and post-burst multi-color photometry monitoring, we
investigate how the accretion-induced excess emissions in the $U$ band are
related to those in other wavelengths. In principle, spots, accretion
variability, and dust extinction can all contribute to observed flux
variations, so genuine excess flux from accretion should be isolated from
other causes of variability. We subtract the in-burst and post-burst light
curves with the best-fit spot models in Table 1. The two segments during the
2022 burst with signs of a slow fade and are also excluded from this
calculation.

We assume the observed fluxes consist only of contributions from accretion and
the photosphere, so any observed flux above the photospheric level is caused
by accretion alone. For the adopted data, the excess emission is calculated in
magnitudes for different bands. For example, the excess emission in $U$ band
$U_{\rm ex}$ is

$U_{\rm ex}=-2.5\log_{10}(10^{-0.4(U_{\rm obs}-A_{U})}-10^{-0.4U_{\rm pho}}),$
(4)

where $U_{\rm obs}$ is the observed $U$ band magnitudes, $A_{U}=1.56~{}A_{V}$
is the extinction in Johnson $U$ band, and $U_{\rm pho}$ is the photosphere
fluxes in the U band, synthesized from X-Shooter spectrum analysis (see Sec
3.1).

Following the procedures above, we calculate the excess emission in U, B,
ASASN-SN $g$, and V band during the 2022 burst. The multi-color photometry is
contemporaneous, with excess emission measurements separated by less than 0.5
days. Figure 11 shows that the $B_{\rm ex}$, $g_{\rm ex}$, and $V_{\rm ex}$
are well correlated to the near-simultaneous $U_{\rm ex}$, with
relationships444 The slopes $\Delta U_{\rm ex}/\Delta(BgV)_{\rm ex}$ fit to
the post-burst data are smaller than those from in-burst data. However, the
measurements from the in-burst light curves span a much larger range ($\Delta
m\sim 2-3$ mag, compared to ($\Delta m\sim 1$ mag) during quiescence, so the
combined fit is adopted here as more reliable. to the post-burst, in-burst,
and combined light curves described by best-fit lines in Table 2.

A comparison between our spectroscopic accretion rate measurements and the
accretion rates independently calculated from photometry establishes that the
photometry-based accretion rates may be extrapolated to the large outbursts.
Figure 11 shows that excess emissions from our photometric method are similar
to those obtained from synthetic photometry from the HIRES spectrum in May
2008, from UBV photometry from WFI (Juhász et al., 2012), and from quiescence
with the X-shooter spectrum, with an average difference of 0.2 dex in the
accretion rates (see also Section 5.4). In general, these three sets of excess
emissions are beyond the calibrated range from the 2022 burst but are in
agreement with the extrapolation of the calibrated relation. Therefore, we
conclude that the excess emission relations we derived from the 2022 burst are
applicable to a wide range, from the quiescent accretion state to the far more
luminous outburst events.

Figure 11: The relations of excess emission due to accretion in the U, B, g,
and V band, derived from photometry and spectroscopy. Open black and blue dots
are the excess emissions measured from 2022 in-burst and post-burst light
curves, respectively. The gray line is the best-fit linear model to the 2022
burst excess emission measurements. For comparison, we also add excess
emission measurements during the 2008 outburst from photometry (yellow, Juhász
et al., 2012) and spectrum (green, Aspin et al., 2010), and during the
quiescent state from X-Shooter spectrum (red, Figure 2). Table 2: Linear Fits
to the Excess Emission

|  | Slope | Intercept | rms [mag]
---|---|---|---|---
| In-burst | $0.93\pm 0.02$ | $0.42\pm 0.25$ | 0.11
$U_{\rm ex}-B_{\rm ex}$ | Post-burst | $0.80\pm 0.06$ | $2.26\pm 0.79$ | 0.09
| Joint-Fit | $0.94\pm 0.01$ | $0.36\pm 0.14$ | 0.10
| In-burst | $1.00\pm 0.05$ | $-0.24\pm 0.65$ | 0.18
$U_{\rm ex}-g_{\rm ex}$ | Post-burst | $0.57\pm 0.13$ | $5.56\pm 1.76$ | 0.16
| Joint-Fit | $0.99\pm 0.03$ | $-0.09\pm 0.38$ | 0.18
| In-burst | $0.81\pm 0.03$ | $2.03\pm 0.39$ | 0.21
$U_{\rm ex}-V_{\rm ex}$ | Post-burst | $0.50\pm 0.04$ | $6.37\pm 0.61$ | 0.11
| Joint-Fit | $0.82\pm 0.02$ | $1.89\pm 0.22$ | 0.18

Note. — The lower limit of excess emission relation: $B_{\rm ex,low}$=14.9
mag, $g_{\rm ex,low}$=14.79 mag, $V_{\rm ex,low}$=14.83 mag.

###  5.2 Calculating accretion Rates from 2022 Burst

Mass accretion rates $\dot{M}_{\rm acc}$ can be calculated as

$\dot{M}_{\rm acc}=(1-\frac{R_{*}}{R_{t}})^{-1}\frac{L_{\rm
acc}R_{*}}{GM_{*}}$ (5)

where $M_{*}=0.83~{}M_{\odot}$ and $R_{*}=1.50~{}R_{\odot}$ are the adopted
mass and radius of EX Lup, the truncation radius $R_{t}$ is assumed to be $5$
$R_{*}$ (adopted for consistency with other accretion rate measurements, see
Hartmann et al., 2016)555The truncation radius is typically assumed $5$
$R_{*}$, based on the balance between the gas pressure and magnetic pressure
(e.g. Johnstone et al., 2014). However, the truncation radius should depend on
$\dot{M}^{2/7}$ (e.g. Johns-Krull & Herczeg, 2007), so epochs with larger
accretion rates may have a smaller truncation radius., and the accretion
luminosity $L_{\rm acc}$ is calculated from the U-band excess luminosity from
Eq. 2 or from the excess fluxes in B, g, and V band through the relations in
Table 2.

Figure 12 presents the histograms of accretion rates measured from the excess
emission of $U$, $B$, $g$, and $V$ band during the 2022 post-burst light
curves. The two dips in brightness are excluded. The photometric variations
induced by spot modulations are removed by the sinusoidal models (as listed in
Table 1). Observations separated by less than two hours are binned into a
single epoch. We also fit the Gaussian probability density functions to the
histograms and show the best-fit mean value and 1$\sigma$ variance of the
Gaussian function on the upper right of each panel.

After the 2022 burst finished, the excess emission measurements from $g$, $B$,
and $V$ all reveal a similar mean accretion rate of $3.33\times 10^{-9}$
$M_{\odot}$ yr-1, two times stronger than the accretion rate measurement of
$1.77\times 10^{-9}$$M_{\odot}$ yr-1 in the X-Shooter spectrum in 2010. This
increase in accretion rate quantifies what is seen visually, the steady
brightening in the $g$-band of 0.6 mag from 2019–2023 (see Fig. 4). The peak
accretion rate of the 2022 burst is $1.74\times 10^{-8}$ $M_{\odot}$ yr-1, as
measured from U-band excess, which is $\sim 5$ times higher than the post-
burst accretion rates and $\sim 10$ times higher than the accretion rate prior
to the burst.

The Gaussian sigma spread of post-burst accretion rates is around 0.08 dex,
lower than the stochastic accretion variability measured for other accreting
young stars (for example, see the analysis of a large sample by Venuti et al.
2021 and in-depth analyses of individual stars by, e.g., Zsidi et al. 2022 and
Herczeg et al. 2023). This distribution should be attributed to the day-to-
month intrinsic accretion variability, due to the fact that the median cadence
of post-burst observations is 1-2 days and the timespan of light curves is
around 4 months, and we have removed the photometric modulations caused by
spots. In the post-burst light curves, the rotational modulations induce
larger photometric variations than the accretion variability.

If the rotational modulations were not corrected, the spread of the accretion
rate derived from apparent photometry would increase to 0.14-0.22 dex.
Assuming the intrinsic accretion variability and modulation are independent,
the intrinsic accretion variability only makes up 36-66% of the apparent
accretion rate spread. This agrees with the finding of Venuti et al. (2014)
that rotational modulation of spots usually drives more significant
photometric variations than intrinsic accretion variability on the timescales
of the rotation period, at least for some objects.

Figure 13 shows the relation between NEOWISE W1 photometry and the mass
accretion rates of EX Lup. The brightening of W1 tends to correlate with the
increase in mass accretion rates. The brightest mid-IR points correspond to
2022-02, the rising stage of the 2022 burst, when the mid-IR emission is 0.35
mag brighter than quiescence while the accretion rate is 0.4 dex brighter than
the pre-burst epoch, 2021-02. However, in some epochs the $W1$ and $W2$ bands
brightened without an accompanying change in accretion rate. For example,
albeit the $g$-band light curves show a 0.1 dex increase in mass accretion
rates between epochs in 2018-08 and 2019-02, the W1 band faded by 0.3 mag from
2018 to 2019. Mid-IR emission is produced by the warm dust in the disk,
although the exact interpretation of this emission and variability is
challenging (e.g. McGinnis et al., 2015; Rebull et al., 2015). The dominant
emission should take the form of the irradiated disk surface (except in more
extreme cases where viscous heating dominates, e.g., Rodriguez & Hillenbrand
2022). The correlations between mass accretion rates and mid-infrared emission
exhibited by EX Lup should help efforts (e.g. Contreras Peña et al., 2020,
2023) to calibrate extensive analyses of long-term mid-IR variability (e.g.
Park et al., 2021; Zakri et al., 2022).

Figure 12: The histograms of mass accretion rates measured from the excess U,
B, ASAS-SN g, and V band emission in the post-burst light curves of the 2022
burst. Gaussian fits the $\log\dot{M}_{\rm acc}$ histograms are plotted with
red dashed lines. Figure 13: The relations between NEOWISE W1 photometry and
mass accretion rates $\dot{M}_{\rm acc}$. The median mass accretion rates are
derived from ASAS-SN $g$(blue) and $V$(red) band photometry taken within 4
days apart from the NEOWISE observation epochs, and the errorbars are the
variability range of mass accretion rates in the 8-day time bins. The filled
areas are the extrapolation ranges of our excess emission relations for $g$
and $V$ fluxes.

###  5.3 Accretion Rates During Major Outbursts and Characteristic Bursts

A visual inspection is performed to identify the bursts from EX Lup historical
light curves between 1945 to 2022. We first select a subset of ten bursts from
AAVSO light curves that are well characterized from the photometric data
(Figure 14). These bursts have near-complete photometric coverage, with well-
determined burst durations set by the time when EX Lup is continuously
brighter than its quiescent state. We also include an outburst event at 1945,
though with a poorly estimated burst duration, identified from DASCH light
curves in Figure 14 to compare with the other two outburst events. For
completeness, we also include some less characterized bursts interrupted by
long data gaps or inconsistent photometry data from different AAVSO observers,
for which we cannot derive reliable $M_{\rm acc}$ from photometry. The light
curves for this set of poorly characterized bursts are displayed in Figure 15.

The complete information of each burst, including the $V$ band magnitude and
accretion rates at peak, burst duration, and total accreted mass, is tabulated
in Table 5.3. We measure the mass accretion rates and accreted masses during
the burst for these ten bursts by converting the observed visible magnitude to
the Johnson V band using Eq. 1. Then accretion rates are calculated via the
excess emission relations in Table 2. The total accreted mass in the burst
$M_{\rm acc,~{}b}$ is obtained by integrating all the photometry-derived mass
accretion rates over the burst duration. We also calculate the ratio of
$M_{\rm acc,~{}b}$ to quiescently accreted mass assuming the same duration
$M_{\rm acc,~{}q}$ to provide a sense of the extra masses accreted during the
burst.

The major outbursts in 1955 and 2008 have much longer durations (hundreds of
days) and reach $m_{V}\sim 8$ mag at peak. The peak mass accretion rate during
these major eruptions is $\sim 3-5\times 10^{-7}$ $M_{\odot}$ yr-1, hundreds
of times larger than the average accretion rate during quiescence. Around
$0.1$ Earth masses can be accreted during these outbursts, which would
otherwise take several decades to consume in the quiescent state. The
characteristic bursts have accretion rates of $2-5\times 10^{-8}$ $M_{\odot}$
yr-1, ten times lower than the major outbursts.

The 1945 outburst does not have as complete photometric coverage as the other
two outbursts, with photometry during only the fading stage spanning $\sim$200
days. Nevertheless, we put constraints on the mass accretion parameters by
using a spline model fitted to longer segments of light curves to simulate the
whole profile of the outburst, and by only using the in-fade photometry. These
two methods estimate the upper and lower limit of accreted mass of the 1945
burst, respectively. The lower limit of peak mass accretion rate of the 1945
outburst is $~{}2.3\times 10^{-7}$ $M_{\odot}$ yr-1, the magnitude of which is
consistent with the 1955 and 2008 outbursts. The accreted mass is 6.2 to
8.9$\times 10^{-8}~{}M_{\odot}$, similar to the 1955 outburst, but still lower
than the 2008 outburst.

We also identify two long-term brightening events in DASCH light curves in
1937 and 1947 (see top panel of Figure 1), each with a duration of 3-4 years,
much longer than the typical duration of bursts (a few months). Taking the B
magnitude calibrated by the APASS catalog and the excess emission relation in
Table 2, these two events brightened $\sim 1$ mag in the B band, translating
to a relative change in $\dot{M}_{\rm acc}$ of 0.47 dex and a peak accretion
rate of $1.8\times 10^{-8}$ $M_{\odot}$ yr-1. This accretion rate is 5 times
larger than the mean accretion rates derived from the recent light curves and
is close to the peak accretion rate in the 2022 burst.

Figure 14: A selection of historical bursts in EX Lup. The three figures in
the first row show EX Lup’s historical outbursts in 1945, 1955, and 2008. The
1945 outburst is additionally fitted with a spline in red. The figures in the
second and third rows show the characteristic bursts. The blue points are
$B$-band magnitudes from DASCH light curves. The open circles are visual
photometry from Albert Jones and are transformed to the Johnson $V$ band via
Eq. 1. The black points are transformed AAVSO visual magnitude, and the black
crosses are photometry made in the Johnson $V$ band from AAVSO. Red points are
SMARTS $V$-band light curves. The equivalent mass accretion rates derived from
excess emission relation are displayed on the right axes. The estimated
accreted mass $M_{\rm acc,b}$ during each burst and its ratio to the quiescent
accretion $M_{\rm acc,q}$ with the same duration is displayed on the upper
right in each panel. The red dashed line stands for the V magnitude baseline
for the mean quiescent rate ($1.7\times 10^{-9}$$M_{\odot}$ yr-1). The blue
dashed line is the $V$-band fluxes of EX Lup’s photosphere.  Figure 15: An
overview of less characterized bursts seen from EX Lup. The photometry symbol
scheme is the same as Figure 14. These bursts are hard to have robust accreted
masses estimates, either due to incomplete photometry coverage caused by
seasonal observational gaps or inconsistent data reported by different
observers.  Table 3: The Historical Burst Properties of EX Lup during
1945-2022

Burst Epoch | Duration | Peak Magnitude1 | Peak $\dot{M}_{\rm acc}$ | $M_{\rm acc,b}$ | $M_{\rm acc,b}/M_{\rm acc,q}$ | Reference
---|---|---|---|---|---|---
| (day) | (V band) | ($\times 10^{-9}\hbox{$M_{\odot}$ yr${}^{-1}$}$) | ($\times 10^{-9}M_{\odot}$) |  |
1945 | 250-500 | 9.5($B$) | >232.94 | 62-89 | 43-54 | -
1955 Apr. | 350 | 8.3 | 358.84 | 72.32 | 53.34 | a,f
1956 Mar. | 400 | 11.2 | 28.83 | 9.69 | 4.45 | f
1973 Aug. | 80 | 11.3 | 26.34 | 1.92 | 6.77 | e
1993 Apr. | 60 | 10.7 | 44.39 | 2.73 | 8.33 | -
1994 Jan. | 80 | 11.2 | 28.27 | 3.09 | 7.32 | b
1995 Apr. | 90 | 11.4 | 23.02 | 2.31 | 13.80 | a,b,d
1998 Jun. | 110 | 11.0 | 38.02 | 6.22 | 14.74 | c
1999 Jun. | 80 | 10.9 | 37.92 | 5.81 | 9.39 | c
2008 Jan. | 200 | 7.8 | 551.37 | 125.14 | 120.87 | c
2022 Mar. | 70 | 10.9 | 45.67 | 2.71 | 7.88 | a
1959 Sep. | $>60^{3}$ | 11.1 | - | - | - | -
1962 Aug. | $>80^{3}$ | 11.2 | - | - | - | -
1964 May | 403 | 11.64 | - | - | - | -
1966 May | 1004 | 11.7 | - | - | - | -
1968 Jan. | 503 | 11.9 | - | - | - | -
1972 Jan. | $<100^{3}$ | 11.6 | - | - | - | -
1977 Feb. | $<120^{3}$ | 10.65 | - | - | - | -
2002 Jul. | $<210^{3}$ | 10.3 | - | - | - | c
2013 | $<500$ | $>11.0$ | - | - | - | -
2020 | - | $>11.0(g)$ | - | - | - | -

11footnotetext: Converted from Visual magnitude via Eq. 1.

22footnotetext: Measured from Johnson V band.

33footnotetext: Sparse photometry coverage.

44footnotetext: Inconsistent visual magnitude between different observers.

Note. — $M_{\rm acc,b}$ is the accreted mass during each burst, and $M_{\rm
acc,q}$ is the accreted mass in the same duration of each burst and assuming
quiescent accretion rate of $\dot{M}_{\rm acc,q}=1.7\times 10^{-9}$
$M_{\odot}$ yr-1 (see Section 3.1).

References. — a: Aspin et al. (2010); b: Herbig et al. (2001); c: Herbig
(2007); d: Lehmann et al. (1995); e: Kóspál et al. (2022); f: Bateson et al.
(1991)

###  5.4 Comparisons to other accretion rate measurements for EX Lup

The accretion luminosities presented in this work are obtained by directly
measuring the excess accretion continuum emission. Table 4 compares our
accretion rates with selected accretion rates from the literature and from
other spectroscopic diagnostics. The accretion rates presented in this work
are much more likely to be underestimated than overestimated due to the
exclusion of low-density flows (see Espaillat et al. 2021, Pittman et al. 2022
and references therein) and the exclusion of line emission (see discussions in
Alcalá et al. 2014 and Herczeg et al. 2023).

The Br$\gamma$ emission is detected in emission, weakly during quiescence and
then stronger as the accretion rate increases. The accretion rates inferred
from the Br$\gamma$ luminosity are well correlated with the accretion rates
from the continuum, lending some confidence to the application of global
correlations between Br$\gamma$ line and continuum luminosity (e.g. Alcalá et
al., 2017) to EX Lup bursts and outbursts. However, the Pa$\beta$ emission
measured by Sipos et al. (2009) leads to an accretion rate that is a factor of
10 lower than the continuum-based accretion rates.

The quiescent accretion rates are systematically larger than the accretion
rates of $\sim 10^{-10}$ $M_{\odot}$ yr-1 that are estimated from H$\alpha$
line widths (e.g. Sicilia-Aguilar et al., 2012, 2015). The difference is
likely attributed to different methodologies. Excess continuum emission should
be considered more reliable, since multiple factors affect the H$\alpha$ line
profile (see also, e.g. Manara et al., 2017).

Table 4: Comparison of Accretion Rates Between Spectra and Photometry  |  |  |  | $\dot{M}_{\rm acc}$ |
---|---|---|---|---|---
Instrument | Date | Accretion | $\dot{M}_{\rm acc}$ | From Excess Emission7 | Ref.
|  | Diagnostics | ($\times 10^{-9}~{}M_{\odot}$ yr-1) | ($\times 10^{-9}~{}M_{\odot}$ yr-1) |
|  |  |  | $U$ | $B$ | $g$ | $V$ |
Keck/HIRES | 2008-05-23 | Continuum | 150 | 124.02 | 79.50 | 72.45 | 69.40 | This work1
Keck/NIRSPEC | 2008-05-23 | Br $\gamma$ | 140 | 124.02 | 79.50 | 72.45 | 69.40 | This work2,3
NTT/SOFI | 2008-04-20 | Br $\gamma$ | 200 | 84.57 | 60.85 | 136.35 | 132.58 | a4,6
MPG/FEROS | 2008-05-08 | H$\alpha$ 10% EW | 30 | - | - | - | 120 | b5
GeminiSouth/IGRINS | 2022-03-14 | Br $\gamma$ | 8.9 | - | - | 9.4 | 10.5 | This work2
VLT/X-Shooter | 2010-05-04 | Continuum | 1.77 | 1.74 | 1.22 | 0.97 | 1.76 | This work1
VLT/X-Shooter | 2010-05-04 | Br $\gamma$ | 1.34 | 1.74 | 1.22 | 0.97 | 1.76 | This work1,2
MPG/FEROS | 2007-07-29 | H$\alpha$ 10% EW | 0.3 | - | - | - | 4.1 | b5
NTT/SOFI | 2001-05-04 | Pa$\beta$ | 0.4 | - | - | - | 5.13 | c6
a: Juhász et al. (2012); b: Sicilia-Aguilar et al. (2012); c: Sipos et al.
(2009).
1\. Synthetic magnitudes from the spectrum.
2\. $L_{\rm line}-L_{\rm acc}$ from Alcalá et al. (2017).
3\. Equivalent width from Aspin et al. (2010), scaled with K-band photometry
from (Juhász et al., 2012).
4\. Photometry from WFI observation (Juhász et al., 2012).
5\. $\dot{M}_{\rm acc}-EW_{H\alpha}$ from Natta et al. (2004)
6\. $L_{\rm line}-L_{\rm acc}$ from Muzerolle et al. (1998)
7\. Assuming $A_{V}=0.1$ mag, $A_{U}=1.5A_{V}$, $A_{B}=1.3A_{V}$,
$A_{g}=1.2A_{V}$.

##  6 Discussion

In this Discussion, we first provide an overview of the major outbursts and
minor bursts detected over the past 130 years (Section 6.1). This accretion
history tells us that over the past 70 years, about two times more mass has
accreted onto EX Lup during outbursts than in quiescence, as described in
Section 6.2. The description of the bursts and quiescence is used in Section
6.3 to offer constraints and challenges to the two hypothetical explanations
for EX Lup-type bursts. We then describe in Section 6.4 the spotted surface of
EX Lup and consequences for photometric variability, which needs to be
considered for the accretion rates during quiescence.

###  6.1 A Historical Overview of Accretion Bursts onto EX Lup

Figure 16 shows the historical light curve in accretion rate from 1893–2022.
The accretion rates are derived from mean magnitudes in half-year bins for
DASCH light curves, 100-day bins for AAVSO light curves, and 60-day bins for
ASAS-SN light curves. The century-long accretion history shows a gradual
decrease in accretion rate up to 1 dex from the 1940s to the 1980s. The large
spikes in the accretion rate evolution are attributed to accretion bursts,
which we further characterize in this subsection and evaluate their roles in
EX Lup’s historical accretion in Section 6.2.

The historical bursts from EX Lup are divided here into two categories:
characteristic bursts and major outbursts. The characteristic bursts have a
typical peak magnitude of $m_{V}\sim 11.5-10.5$ and last for $\lesssim$ 100
days. The mass accretion rate during these characteristic bursts is elevated
above the quiescent level by 10–30 at their peak. These characteristic bursts
accrete masses $\sim 10^{-3}$ Earth masses, or around ten times the
quiescently accreted mass during the same time period. The recent 2022 burst
also falls into the category of a characteristic burst.

Figure 17 shows the histogram of $M_{\rm acc}$ of ten well-characterized
bursts. Six equal-width bins are used to plot the histogram. Admittedly, the
number of bursts is small and the choice of the bin will affect the
statistics, but it generally reflects a trend that the characteristic bursts
are more common in EX Lup. No bursts in the sample have $M_{\rm acc}$ between
$10^{-8.0}\sim 10^{-7.25}~{}M_{\odot}$.

To obtain the statistical significance for the lack of bursts in this
intermediate $M_{\rm acc}$ bin, we fit a power law to the number of bursts
versus the accreted masses. Then we draw samples from the fitted power law and
perform the Kolmogorov-Smirnov test with our sample. The $p-$value is
$1.08\times 10^{-5}$, suggesting the accreted mass of bursts seen from EX Lup
are inconsistent with a single power-law distribution. The two distinctive
classes of bursts seen from EX Lup suggest that the driving causes of these
bursts are not regulated by self-organized criticality (Bak et al., 1987). In
other words, the distribution of mass accreted in bursts form two distinct
populations, the characteristic bursts, and outbursts, which might not share
common driving mechanisms.

The adopted uncertainty of extinction and the excess emission relation
introduces 15% and 40% uncertainties to the accretion rates and total accreted
mass estimate, which is not sufficient to reconcile the gaps in accreted mass
histogram. If the truncation radius shrinks during bursts, the accretion rate
during bursts would be expected to be even higher than measured here. The
efforts of compiling a complete EX Lup bursts record would also fail due to
the sporadic data gaps, so we cannot completely rule out bursts of
intermediate size. However, the major reason for incompleteness is the $\sim
2-3$-month seasonal data gaps every year, comparable to the typical duration
of characteristic bursts ($\sim 100$ days).

If the extinction is $A_{V}=1.1$ mag, as measured by Alcalá et al. (2017) and
adopted by Cruz Sáenz de Miera et al. (2023), all accretion rates would be
roughly a factor of $\sim 6.5$ times higher than estimated here (see change in
accretion luminosity and radius described in §3.1).

Figure 16: The mass accretion rates of EX Lup between 1890 and 2022. DASCH
$B$-band, AAVSO $V$-band, and ASAS-SN $g$\- and $V$\- band light curves are
used to derive the accretion rates. Filled and hollow points are within and
beyond the calibrated ranges of excess emission relations, respectively.
Figure 17: The histogram of accreted masses $M_{\rm acc}$ of ten bursts from
Figure 14. The black sticks at the bottom mark the value of $M_{\rm acc}$ for
each burst. The error bars of the number of bursts in each $M_{\rm acc}$ bin
are the Poisson error due to the small sample size. The dashed line is a power
law relation fitted to the burst number in the three accreted masses bins.
However, the KS test suggests the accreted mass of 10 bursts is not likely
sampled from such power law distribution (with a p-value of $1.08\times
10^{-5}$), due to the non-detections between $10^{-8.0}\sim
10^{-7.25}~{}M_{\odot}$.

###  6.2 The Ratio of Accreted Masses During Bursts and Quiescent State

We have now tabulated the history of accretion onto EX Lup over the past
century, including a tabulation of bursts and accretion rates and an
assessment of accretion during quiescent periods. These accretion rates are
calculated by converting historical photometry to accretion rate through
correlations developed using the multi-color photometry in the 2022 burst. To
move forward, we ask the question: _is more mass accreted by EX Lup during
bursts or during quiescent?_

An exact comparison is difficult, owing to the non-uniform photometric data
coverage in 130 years, especially because the sparse monitoring before 1955
affects our knowledge about the true frequency of bursts in EX Lup. As a
consequence, we restrict the discussion to the 70-year photometry from AAVSO
and ASAS-SN light curves, during which we have a complete census of major
bursts and a near-complete census of characteristic bursts. For the
characteristic bursts, the photometry may not be as complete due to the
seasonal data gaps (two months per year), therefore we estimate completeness
of 84% to correct for the characteristic burst frequency.

As shown in Table 5.3. The major outbursts in 1955 and another in 2008 have a
total (combined) mass accreted of $1.97\times 10^{-7}~{}M_{\odot}$. The mean
accreted masses in characteristic bursts is $3.68\times 10^{-9}~{}M_{\odot}$.
Assuming the apparent frequency of characteristic bursts is two per decade,
further compensating the 84% completeness and plus two major outbursts, the
accreted masses during bursts $M_{\rm acc,~{}b}$ between 1955 and 2022 is

$\begin{split}M_{\rm acc,~{}b}=\frac{14}{0.84}\times 3.68\times
10^{-9}M_{\odot}+1.97\times 10^{-7}M_{\odot}\\\ =2.6\times
10^{-7}~{}M_{\odot}.\end{split}$ (6)

Despite the higher frequency of characteristic bursts, the major outbursts are
responsible for $\sim 80\%$ of masses accreted during all bursts.

Taking the mass accretion rate during quiescence of $\dot{M}_{\rm
acc,~{}q}=1.7\times 10^{-9}$ $M_{\odot}$ yr-1, as measured in Section 3.1, a
total of $1.1\times 10^{-7}~{}M_{\odot}$ is accreted during quiescent between
1955 and 2022. Therefore, In this case, $M_{\rm acc,~{}q}/M_{\rm acc,~{}b}\sim
1:2$, the amount of mass accreted during the bursts is two times larger than
the mass accreted in quiescence. However, the accretion rate of EX Lup during
the quiescence is variable. For example, in Section 5.2, we found that EX Lup
exhibited a $\sim 0.2$ dex accretion rate increase between 2010 and 2022, as
seen by the gradual brightening in ASAS-SN $g$ light curves (see Figure 4
also). The DASCH light curves also exhibit two long-term brightening events
around 1938 and 1949, both of which lasted 2-3 years. The peak accretion rate
of $1.8\times 10^{-8}$ $M_{\odot}$ yr-1 is similar to those measured from the
AAVSO characteristic bursts, indicating that the quiescent accretion can also
reach the $\dot{M}_{\rm acc}$ value that characteristics burst can attain at
peak, but in a slower way. These changes bring the quiescent and burst
accretion totals closer to equality.

Since the outburst-consumed mass dominates the total accreted mass budget in
all bursts, we caution that the ratio calculated above will be largely
sensitive to the outburst frequency. The apparent outburst frequency in AAVSO
68-year data is one event per 34 years, which is similar to the outburst
frequency between 1890-2022 (one event per 44 years), assuming the 1945
outburst is the only outburst present between 1890-1953. Therefore, we cannot
draw conclusions on whether the outburst frequency is overestimated or not
based on historical data.

We conclude that roughly half of the mass of EX Lup accretes during major
accretion bursts and half accretes during characteristic bursts and
quiescence. The quiescent accretion rate of 1.7-3.3$\times 10^{-9}$
$M_{\odot}$ yr-1 measured across the past ten years is similar to (slightly
lower than) the accretion rate of $3\times 10^{-9}$ $M_{\odot}$ yr-1 expected
for a 0.83 M⊙ star, based on the relationship666scaled for Gaia DR3 distances,
as in Manara et al. (2021), and with masses rescaled to the evolutionary
tracks with 50% spots of Somers et al. (2020). for complete samples of disks
in Lupus and Cham I from Alcalá et al. (2017) and Manara et al. (2017).
Incorporating the full accretion history, including periods of quiescence with
an elevated steady accretion rate, increases the accretion rate of EX Lup by a
factor of $\sim 4$, a factor of a few higher than expected for its mass (from
a sample of single-epoch accretion rates). The ratio of gas accreted in bursts
versus quiescence is robust to extinction uncertainty, since all accretion
rates would be shifted by a similar amount.

###  6.3 Mass loading in the innermost disk

The detailed accretion history of EX Lup should help to motivate improvements
and constrain models to understand these short outbursts. Accretion outbursts
are generically explained by an excess of mass in the inner disk, leading to
excess accretion until that mass has accreted or otherwise dispersed (see,
e.g. Hartmann & Kenyon, 1996). The excess mass may be located only at the
inner disk or may be spread across larger radii.

For EX Lup, the outbursts are strictly an inner disk phenomenon, for two
reasons. First, as seen in Figure 18, the rise times are rapid, with fades
that are slower but still on timescales associated with the inner disk.
Second, the accretion rates of $>10^{-7}$ M⊙ yr-1 during large outbursts are
near the range where the optical and near-IR spectrum should start to show
evidence for accretion heating of the viscous disk, including as CO absorption
(Liu et al., 2022), if the system were in steady state. The lack of these
features indicates that the innermost disk is not being fed at such a large
rate, otherwise the size of the heated disk would be larger and would dominate
the spectrum. Once the excess mass in the innermost disk is accreted, the
outburst fades.

The EX Lup-type outbursts occur as gas builds up at or very near the inner
disk and is then released onto the star. In the framework of D’Angelo & Spruit
(2010, 2012), the magnetic field of the star connects to the disk at a
truncation radius that is larger than the corotation radius, thereby
preventing gas from accreting. Gas piles up at the truncation radius until the
truncation radius moves interior to the corotation radius, producing a large
flow of gas onto the star. An alternative hypothesis proposed by Armitage
(2016) suggests that cycles of the stellar dynamo can trigger a burst of
accretion through a change in polarity, which generates a more efficient Hall
effect.

Our results pose challenges for both of these interpretations. For the
magnetospheric instability model, the accretion should cease until the
truncation radius gets smaller, but EX Lup during quiescent periods accretes
gas at a similar rate as other accreting young stars. In this case, the
truncation radius must be leaky. Moreover, whatever mechanism that causes gas
to build up at the star would probably be common.

For the stellar dynamo triggering hypothesis, large outbursts should occur at
the (semi-) regular intervals expected for changes in polarity, but outbursts
occurred in 1945, 1955, and 2008, not at regular intervals. In this scenario,
the triggering of characteristic and large bursts would have to be distinct,
with the characteristic outbursts unexplained. The accretion burst should
begin where the Hall effect is most efficient, beyond the innermost disk.

Whatever the physics, the phenomenon must be a cyclical build-up and release
of gas at the inner disk. During quiescence, the star-disk accretion rate is
lower than the accretion rate through the disk. This process is similar to the
toy model that Lee et al. (2020) described for the periodic variable EC 53
(V371 Ser), except that outbursts of EX Lup are much larger and less frequent.

Figure 18: The rise and fall slopes in magnitudes per day, measured from ten
bursts in Figure 14 (excluding the 1945 outburst), color-coded by the accreted
mass during the burst.

###  6.4 The Phase Lag of Rotational Modulation

The periods measured between 2004 and 2023 for light curves and also from
absorption lines from Kóspál et al. (2014) remain steady at $\sim 7.415$ days
(Table 5), during periods of weak and strong accretion. The periodicity is
interpreted by Sicilia-Aguilar et al. (2012) and Sicilia-Aguilar et al. (2015)
as a stable and rotating, non-axisymmetric accretion column and hot spot tied
to the stellar surface. In Section 4.3.1 we described the characteristics of
spot modulation in EX Lup multi-color light curves. Here we summarize them as
follows

  * •

The profile of spot modulation follows sinusoidal waves. The peak-to-peak
amplitude of modulation increases to shorter wavelengths, from 0.20 mag ($\sim
10\%$) in the I band to 0.6 mag ($\sim 75\%$ in flux changes) in the U band.

  * •

The light curves exhibit a wavelength-dependent phase-lag, with $I$ band
lagged by $\Delta\phi\sim 0.23$ (1.5 days) compared to the $U$ band in post-
burst light curves. In-burst light curves have phase lags as well, though
smaller than those in post-burst light curves.

Following the interpretation of Sicilia-Aguilar et al. (2015), we develop a
toy model with a cool and a hot spot to explain the phase lag and change in
phase during 2022. Although a detailed description of spot parameters is
beyond the scope of this study, we present a nominal model that qualitatively
reproduces the observed lag in light curves. A single circular spot on the
surface will produce rotational modulation that is the same at all wavelengths
and is therefore insufficient to explain the phase differences. With multiple
spots, the modulations generated by each spot may have different phases,
depending on the longitudes and temperatures. Two circular cool spots (or two
hot spots) would produce entirely anti-phase rotational modulation, when their
longitude difference is $180^{\circ}$, while the longitude difference of one
cool and hot spot should be $0^{\circ}$ to produce the anti-phase modulations.
The phase lag from multi-color photometry would be the result of color-
dependent modulation amplitudes between spots. The spots with different
temperatures, thus different spot contrast to the photosphere, would also
change the overall modulation phase at different wavelengths.

A schematic illustration of the two-component spot model and its synthetic
light curves are summarized in Figure 19. In our proposed model, EX Lup has a
rotation period of 7.42 days and an inclination of 30∘ (as estimated by Sipos
et al. 2009 and Sicilia-Aguilar et al. 2015 from accretion diagnostics,
consistent with the $32.4^{\circ}$ inclination of the cold dust disk imaged
with ALMA by Hales et al. 2018). The photosphere temperature is adopted as
$T_{\rm pho}=3850~{}K$ (Alcalá et al., 2017), with two adjacent hot and cool
spots with large temperature differences (Figure 19 b). For simplicity, we
adopt blackbody emission to calculate the spot-photosphere contrast at
different wavelengths. However, the $\sim 0.4$ mag amplitude in the TESS light
curve indicates that the cool spot must cover $\sim 50$% of the visible
stellar surface (see, e.g., 80% coverage fraction for a similar amplitude
variability of LkCa 4, measured by Gully-Santiago et al. 2017). The hot spot
is related to a confined area at the footpoint of the accretion column (see,
Grosso et al., 2010; Sicilia-Aguilar et al., 2015), so the size of the hot
spot should be smaller than the proposed cool spot, but with much higher
temperature compared to the photosphere to induce significant flux variations
in short wavelengths. Due to the color difference and the large size of the
cool spot, the amplitudes of the combined spot modulations differ between
wavelengths, with the U-band light curves largely modulated by the hot spot
and the I-band light curve more influenced by the cool spot (Figure 19 b).

The simulated phase-lag between $U$ and $BVRI$ band is similar to the observed
ones (Figure 8). The model parameters, listed in Table 6, serve only as a
guide. The phase lag and amplitudes have a degenerate solution in spot
contrast (the spot temperature) and their distribution in longitude, as seen
in Figure 19 c, where we show how the phase and morphology of combined I-band
light curves change with the longitude difference between two spots.

This two-spot model for the periodicity of EX Lup, inferred from the phase
lag, is supported by the change in phase from before and after the 2022 burst
(see Figure 7 in §4.2). Before the 2022 burst, the accretion was weaker, so
the periodicity was modulated by the rotation of the cool spot. After the 2022
burst, the accretion is stronger, so the rotational modulation is determined
by the rotation of the hot spot, at a different phase than the cool spot. The
stochastic bursts gradually develop in the modulation minima between 2019-2021
(as is also seen in TESS light curve Figure 5), so the region hosting these
bursts should be close to the large cool spot in the longitude phase, as
visualized in Figure 19.

A phase lag of $\Delta\phi_{U-I}\sim 0.05$ was measured previously for the
active star HR 1099, interpreted as two cool spots and one hot spot on the
stellar surface (Zhai et al., 1994). More relevant analogs are accreting young
stars, which feature photospheric cool spots and magnetic hot spots. Espaillat
et al. (2021) measured a time-lag of a tenth of the rotational period between
UV and optical peaks for the classical T Tauri star GM Aur. In a follow-up
analysis, Robinson et al. (2022) found time lags in TESS and UBVRI light
curves for 14 classical T Tauri stars, suggesting these longitudinal
stratifications of accretion columns might be common in the CTTSs. The time
lag is interpreted by the non-azimuthally symmetric accretion hot spot region
that has a radial density gradient distribution, leading to different SED
profiles at different rotational phases. Although this description seems
discrepant with our interpretation and model, the formulation and
interpretation are actually quite similar to the presence of hot and cool
spots.

Table 5: Comparisons between EX Lup’s Light Curve Modulation Periods and RV period. Time | Data | Periodogram Period
---|---|---
|  | (day)
2004-2005 | $V$ | 7.406$\pm$0.098
2007-2012 | Absorption lines aaWhile the periodicities from emission lines are more scattered between 6-8 days, see the Table A.2 in Sicilia-Aguilar et al. (2015) for details  | 7.417$\pm 0$.001
2016-2018 | $V$ | 7.419$\pm$0.061
2018-2022 | $g$ | 7.412$\pm$0.042
2022-2023 | $g$ | 7.411$\pm$0.102
Table 6: The two-component spot model of Figure 19

| $T_{\rm spot}$ | Radius | Latitude | Longitude
---|---|---|---|---
Hot spot | 6500 K | 4∘ | 35∘ | 0∘
Cool spot | 3000 K | 30∘ | 35∘ | 55∘

Note. — The effective temperature of photosphere and the inclination of the
star is assumed to be 3850 K and 30∘.

Figure 19: A schematic illustration of a two-component spot model that can
produce the observed multi-color light curves. (a) The multi-color light
curves from the two-component spot model. The predicted phase lags between U
and other bands are indicated. (b) The intensity plot of the stellar surface
for the two-component spot model, with hot (blue) and cool (dark red) spots,
and uniform photosphere background (red). The adopted model parameters are
listed in Table 6. (c) The phase shift and morphological change of synthetic I
band light curves according to different longitude differences of two spots
$\Delta\Omega$. (d) The resulting phase lags of brightness maxima between U
and I band light curves as a function of spots’ longitude differences. This
set of figures is produced by the python package starry (Luger et al., 2019).

##  7 Conclusions

In this work, we have assembled the EX Lup light curve from 1895 to 2023 to
assess the role of accretion bursts and outbursts in this archetype of its
class. The brightness increases from historical photometric data are divided
into two categories: the major outbursts (occurred in 1945, 1955, and 2008)
and characteristic bursts (14 detected in the past 70 years). Three outbursts
occurred during 130 years of observations, each with large photometric
variations up to 5 mags and year-long durations.

We build the historical record of accretion onto EX Lup in Table 5.3 by
leveraging correlations during multi-band photometry of the characteristic
burst in 2022, along with past flux-calibrated spectra. The accretion rate
reaches a peak of $\sim 5\times 10^{-7}$ $M_{\odot}$ yr-1 during large bursts
and $\sim 2-5\times 10^{-8}$ $M_{\odot}$ yr-1 during characteristic bursts,
and wanders from $10^{-9}$ to $10^{-8}$ $M_{\odot}$ yr-1 during quiescent
periods of steady accretion. The two outbursts in 1955 and 2008 accrete
$\sim$0.1 Earth masses, or around 80% of the total burst-consumed masses in
seventy years, while the characteristic bursts are more common than the
outburst, but they typically only consume a few $10^{-3}$ Earth masses during
the burst. The burst-consumed mass is roughly 2 times the mass accreted during
quiescent periods, a comparison limited by uncertainty in the frequency of
large bursts and a more precise long-term baseline for quiescent accretion
rate.

The characteristic burst in 2022 and the quiescent periods before and after
the burst feature spot modulation with a periodicity of $\sim$7.415 days,
revealed through consistent, daily monitoring by ASAS-SN and AAVSO. These
observations are consistent with the scenario that a hot spot originated in
the shock region linking the stellar surface and the non-axisymmetric
accretion column, which is modulated by the stellar rotation (see, e.g.,
Sicilia-Aguilar et al., 2015). Moreover, color-dependent time lags and phase
shifts from 2022 post-burst multi-color light curves indicate the presence of
both a cool and a hot spot.

Our methodology based on photometry to estimate accretion rates provides an
approach to exploit the historical data for long-monitored young variable
stars. New surveys such as LSST, WFST, and TiDO, combined with ASAS-SN and
AAVSO observers will open a new dimension to compare the characteristics of
eruptive young stars (see, e.g., Bonito et al., 2023): do other stars that
historically show repetitive eruptions, especially EXors, also exhibit two-
mode bursts and accrete a significant amount of masses in bursts during their
final stages?

As for EX Lup itself, we expect several more characteristic bursts in the next
decade, based on the clustering of previous bursts between 1993-2002 and the
gradually elevated activities seen from recent light curves. Since the 1955
burst occurred only 10 years after the 1945 burst, we should trigger
additional spectroscopic observations early to follow the accretion and
ejection processes during the burst, since the physics and the consequences
for the disk are studied but still insufficiently explained. Even if we have
to wait until 2060 for the next major outburst, continued monitoring will
allow us to assess why mass builds up in the inner disk and why, eventually,
that dam will break.

##  8 Acknowledgements

We thank the anonymous referee whose comments improve the quality of this
paper. We are extraordinarily grateful for the extensive contributions by many
people, most especially Albert Jones, for historical monitoring of EX Lup. We
thank Kenneth Menzies and other observers for contributions to the AAVSO light
curve during the 2022 burst. We acknowledge with thanks the variable star
observations from the AAVSO International Database contributed by observers
worldwide and used in this research.

We also thank Carlo Manara for discussions and providing the X-Shooter
spectrum obtained in 2010 and Nienke van der Marel for contributions to the
WHT/ISIS spectrum obtained in 2012. GJH thanks Sergei Nayakshin and Guoxingan
Ma for the discussions about this project.

MW and GJH appreciate the Chinese Center for Advanced Science and Technology
for hosting the Protoplanetary Disk and Planet Formation Summer School in
2022, organized by Xue-ning Bai and Ruobing Dong, which facilitated this
research. GJH also appreciates discussions with Will Fischer, Lynne
Hillenbrand, Agnes Kospal, and Mike Dunham during the preparation of
Protostars and Planets VII review, which helped to motivate this work.

This work is supported by the National Natural Science Foundation of China
(grant Nos. 11973028, 11933001, 1803012, 12150009, 12173003) and the National
Key R&D Program of China (2019YFA0706601, 2022YFA1603102). We also acknowledge
the Civil Aerospace Technology Research Project (D050105). D.J. is supported
by NRC Canada and by an NSERC Discovery Grant.

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. The TESS 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 in Sector 12
and 39 can be accessed via MAST (MAST Team, 2021).

## References

  * Ábrahám et al. (2009) Ábrahám, P., Juhász, A., Dullemond, C. P., et al. 2009, Nature, 459, 224, doi: 10.1038/nature08004
  * Alcalá et al. (2014) Alcalá, J. M., Natta, A., Manara, C. F., et al. 2014, A&A, 561, A2, doi: 10.1051/0004-6361/201322254
  * Alcalá et al. (2017) Alcalá, J. M., Manara, C. F., Natta, A., et al. 2017, A&A, 600, A20, doi: 10.1051/0004-6361/201629929
  * Armitage (2016) Armitage, P. J. 2016, ApJ, 833, L15, doi: 10.3847/2041-8213/833/2/L15
  * Aspin et al. (2010) Aspin, C., Reipurth, B., Herczeg, G. J., & Capak, P. 2010, Astrophysical Journal Letters, 719, 50, doi: 10.1088/2041-8205/719/1/L50
  * Audard et al. (2014) Audard, M., Ábrahám, P., Dunham, M. M., et al. 2014, in Protostars and Planets VI, ed. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 387–410, doi: 10.2458/azu_uapress_9780816531240-ch017
  * Bak et al. (1987) Bak, P., Tang, C., & Wiesenfeld, K. 1987, Phys. Rev. Lett., 59, 381, doi: 10.1103/PhysRevLett.59.381
  * Banzatti et al. (2015) Banzatti, A., Pontoppidan, K. M., Bruderer, S., Muzerolle, J., & Meyer, M. R. 2015, Astrophysical Journal Letters, 798, 2, doi: 10.1088/2041-8205/798/1/L16
  * Bateson et al. (1991) Bateson, F. M., McIntosh, R., & Brunt, D. 1991, Royal Astronomical Society of New Zealand Publications of Variable Star Section, 16, 49
  * Bonito et al. (2023) Bonito, R., Venuti, L., Ustamujic, S., et al. 2023, ApJS, 265, 27, doi: 10.3847/1538-4365/acb684
  * Calvet & Gullbring (1998) Calvet, N., & Gullbring, E. 1998, ApJ, 509, 802, doi: 10.1086/306527
  * Cardelli et al. (1989) Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, ApJ, 345, 245, doi: 10.1086/167900
  * Carvalho & Hillenbrand (2022) Carvalho, A. S., & Hillenbrand, L. A. 2022, ApJ, 940, 156, doi: 10.3847/1538-4357/ac9d8e
  * Connelley & Reipurth (2018) Connelley, M. S., & Reipurth, B. 2018, ApJ, 861, 145, doi: 10.3847/1538-4357/aaba7b
  * Contreras Peña et al. (2020) Contreras Peña, C., Johnstone, D., Baek, G., et al. 2020, MNRAS, 495, 3614, doi: 10.1093/mnras/staa1254
  * Contreras Peña et al. (2017) Contreras Peña, C., Lucas, P. W., Minniti, D., et al. 2017, MNRAS, 465, 3011, doi: 10.1093/mnras/stw2801
  * Contreras Peña et al. (2023) Contreras Peña, C., Herczeg, G. J., Ashraf, M., et al. 2023, MNRAS, 521, 5669, doi: 10.1093/mnras/stad820
  * D’Angelo & Spruit (2010) D’Angelo, C. R., & Spruit, H. C. 2010, MNRAS, 406, 1208, doi: 10.1111/j.1365-2966.2010.16749.x
  * D’Angelo & Spruit (2012) —. 2012, MNRAS, 420, 416, doi: 10.1111/j.1365-2966.2011.20046.x
  * Espaillat et al. (2021) Espaillat, C. C., Robinson, C. E., Romanova, M. M., et al. 2021, Nature, 597, doi: 10.1038/s41586-021-03751-5
  * Fischer et al. (2011) Fischer, W., Edwards, S., Hillenbrand, L., & Kwan, J. 2011, ApJ, 730, 73, doi: 10.1088/0004-637X/730/2/73
  * Fischer et al. (2022) Fischer, W. J., Hillenbrand, L. A., Herczeg, G. J., et al. 2022, arXiv e-prints, arXiv:2203.11257, doi: 10.48550/arXiv.2203.11257
  * Gahm et al. (2008) Gahm, G. F., Walter, F. M., Stempels, H. C., Petrov, P. P., & Herczeg, G. J. 2008, A&A, 482, L35, doi: 10.1051/0004-6361:200809488
  * Gaia Collaboration et al. (2021) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2021, A&A, 649, A1, doi: 10.1051/0004-6361/202039657
  * Gaia Collaboration et al. (2022) Gaia Collaboration, Vallenari, A., Brown, A. G. A., et al. 2022, arXiv e-prints, arXiv:2208.00211, doi: 10.48550/arXiv.2208.00211
  * Ghosh et al. (2023) Ghosh, A., Sharma, S., Ninan, J. P., et al. 2023, Journal of Astrophysics and Astronomy, 44, 50, doi: 10.1007/s12036-023-09939-7
  * Giannini et al. (2022) Giannini, T., Giunta, A., Gangi, M., et al. 2022, ApJ, 929, 129, doi: 10.3847/1538-4357/ac5a49
  * Goto et al. (2011) Goto, M., Regály, Z., Dullemond, C. P., et al. 2011, ApJ, 728, 5, doi: 10.1088/0004-637X/728/1/5
  * Grosso et al. (2010) Grosso, N., Hamaguchi, K., Kastner, J. H., Richmond, M. W., & Weintraub, D. A. 2010, Astronomy and Astrophysics, 522, 1, doi: 10.1051/0004-6361/200913850
  * Gullbring et al. (1998) Gullbring, E., Hartmann, L., Briceño, C., & Calvet, N. 1998, ApJ, 492, 323, doi: 10.1086/305032
  * Gully-Santiago et al. (2017) Gully-Santiago, M. A., Herczeg, G. J., Czekala, I., et al. 2017, ApJ, 836, 200, doi: 10.3847/1538-4357/836/2/200
  * Guo et al. (2021) Guo, Z., Lucas, P. W., Contreras Peña, C., et al. 2021, MNRAS, 504, 830, doi: 10.1093/mnras/stab882
  * Hales et al. (2018) Hales, A. S., Pérez, S., Saito, M., et al. 2018, ApJ, 859, 111, doi: 10.3847/1538-4357/aac018
  * Hamann & Persson (1992) Hamann, F., & Persson, S. E. 1992, ApJS, 82, 247, doi: 10.1086/191715
  * Hamuy et al. (1992) Hamuy, M., Walker, A. R., Suntzeff, N. B., et al. 1992, PASP, 104, 533, doi: 10.1086/133028
  * Hartmann et al. (2016) Hartmann, L., Herczeg, G., & Calvet, N. 2016, ARA&A, 54, 135, doi: 10.1146/annurev-astro-081915-023347
  * Hartmann & Kenyon (1996) Hartmann, L., & Kenyon, S. J. 1996, ARA&A, 34, 207, doi: 10.1146/annurev.astro.34.1.207
  * Herbig et al. (2001) Herbig, G., Aspin, C., Gilmore, A., Imhoff, C., & Jones, A. 2001, Publications of the Astronomical Society of the Pacific, 113, 1547, doi: 10.1086/324420
  * Herbig (1950) Herbig, G. H. 1950, PASP, 62, 211, doi: 10.1086/126276
  * Herbig (2007) Herbig, G. H. 2007, The Astronomical Journal, 133, 2679, doi: 10.1086/517494
  * Herbig (2008) Herbig, G. H. 2008, AJ, 135, 637, doi: 10.1088/0004-6256/135/2/637
  * Herbig et al. (1992) Herbig, G. H., Suntzeff, N., & Blanco, B. 1992, Information Bulletin on Variable Stars, 3755, 1
  * Herbst et al. (1994) Herbst, W., Herbst, D. K., Grossman, E. J., & Weinstein, D. 1994, AJ, 108, 1906, doi: 10.1086/117204
  * Herczeg & Hillenbrand (2008) Herczeg, G. J., & Hillenbrand, L. A. 2008, ApJ, 681, 594, doi: 10.1086/586728
  * Herczeg & Hillenbrand (2014) —. 2014, ApJ, 786, 97, doi: 10.1088/0004-637X/786/2/97
  * Horne (1986) Horne, K. 1986, PASP, 98, 609, doi: 10.1086/131801
  * Hsieh et al. (2019) Hsieh, T.-H., Murillo, N. M., Belloche, A., et al. 2019, ApJ, 884, 149, doi: 10.3847/1538-4357/ab425a
  * Johns-Krull & Herczeg (2007) Johns-Krull, C. M., & Herczeg, G. J. 2007, ApJ, 655, 345, doi: 10.1086/508770
  * Johnstone et al. (2014) Johnstone, C. P., Jardine, M., Gregory, S. G., Donati, J. F., & Hussain, G. 2014, MNRAS, 437, 3202, doi: 10.1093/mnras/stt2107
  * Jørgensen et al. (2015) Jørgensen, J. K., Visser, R., Williams, J. P., & Bergin, E. A. 2015, A&A, 579, A23, doi: 10.1051/0004-6361/201425317
  * Juhász et al. (2012) Juhász, A., Dullemond, C. P., Van Boekel, R., et al. 2012, Astrophysical Journal, 744, doi: 10.1088/0004-637X/744/2/118
  * Kochanek et al. (2017) Kochanek, C. S., Shappee, B. J., Stanek, K. Z., et al. 2017, PASP, 129, 104502, doi: 10.1088/1538-3873/aa80d9
  * Kóspál et al. (2022) Kóspál, Á., Fiorellino, E., Ábrahám, P., Giannini, T., & Nisini, B. 2022, Research Notes of the American Astronomical Society, 6, 52, doi: 10.3847/2515-5172/ac5e3c
  * Kóspál et al. (2011) Kóspál, Á., Ábrahám, P., Goto, M., et al. 2011, ApJ, 736, 72, doi: 10.1088/0004-637X/736/1/72
  * Kóspál et al. (2014) Kóspál, A., Mohler-Fischer, M., Sicilia-Aguilar, A., et al. 2014, Astronomy and Astrophysics, 561, 1, doi: 10.1051/0004-6361/201322428
  * Kóspál et al. (2023) Kóspál, Á., Ábrahám, P., Diehl, L., et al. 2023, ApJ, 945, L7, doi: 10.3847/2041-8213/acb58a
  * Laycock et al. (2010) Laycock, S., Tang, S., Grindlay, J., et al. 2010, AJ, 140, 1062, doi: 10.1088/0004-6256/140/4/1062
  * Lee (2007) Lee, J.-E. 2007, Journal of Korean Astronomical Society, 40, 83, doi: 10.5303/JKAS.2007.40.4.083
  * Lee et al. (2017) Lee, J.-J., Gullikson, K., & Kaplan, K. 2017, Igrins/Plp 2.2.0, Zenodo, Zenodo, doi: 10.5281/zenodo.845059
  * Lee et al. (2020) Lee, Y.-H., Johnstone, D., Lee, J.-E., et al. 2020, ApJ, 903, 5, doi: 10.3847/1538-4357/abb6fe
  * Lehmann et al. (1995) Lehmann, T., Reipurth, B., & Brandner, W. 1995, 08, 1.  https://arxiv.org/abs/9506066
  * Liu et al. (2022) Liu, H., Herczeg, G. J., Johnstone, D., et al. 2022, ApJ, 936, 152, doi: 10.3847/1538-4357/ac84d2
  * Lorenzetti et al. (2012) Lorenzetti, D., Antoniucci, S., Giannini, T., et al. 2012, ApJ, 749, 188, doi: 10.1088/0004-637X/749/2/188
  * Luger et al. (2019) Luger, R., Agol, E., Foreman-Mackey, D., et al. 2019, AJ, 157, 64, doi: 10.3847/1538-3881/aae8e5
  * Mainzer et al. (2011) Mainzer, A., Grav, T., Bauer, J., et al. 2011, ApJ, 743, 156, doi: 10.1088/0004-637X/743/2/156
  * Mainzer et al. (2014) Mainzer, A., Bauer, J., Cutri, R. M., et al. 2014, ApJ, 792, 30, doi: 10.1088/0004-637X/792/1/30
  * Manara et al. (2013) Manara, C. F., Testi, L., Rigliaco, E., et al. 2013, A&A, 551, A107, doi: 10.1051/0004-6361/201220921
  * Manara et al. (2013) Manara, C. F., Testi, L., Rigliaco, E., et al. 2013, Astronomy and Astrophysics, 551, doi: 10.1051/0004-6361/201220921
  * Manara et al. (2017) Manara, C. F., Testi, L., Herczeg, G. J., et al. 2017, A&A, 604, A127, doi: 10.1051/0004-6361/201630147
  * Manara et al. (2021) Manara, C. F., Frasca, A., Venuti, L., et al. 2021, arXiv e-prints, arXiv:2103.12446.  https://arxiv.org/abs/2103.12446
  * MAST Team (2021) MAST Team. 2021, TESS Light Curves - All Sectors, STScI/MAST, doi: 10.17909/T9-NMC8-F686
  * McGinnis et al. (2015) McGinnis, P. T., Alencar, S. H. P., Guimarães, M. M., et al. 2015, A&A, 577, A11, doi: 10.1051/0004-6361/201425475
  * McLaughlin (1946) McLaughlin, D. B. 1946, AJ, 52, 109, doi: 10.1086/105935
  * Muzerolle et al. (1998) Muzerolle, J., Calvet, N., & Hartmann, L. 1998, ApJ, 492, 743, doi: 10.1086/305069
  * Natta et al. (2004) Natta, A., Testi, L., Neri, R., Shepherd, D. S., & Wilner, D. J. 2004, A&A, 416, 179, doi: 10.1051/0004-6361:20035620
  * NEOWISE Team (2020) NEOWISE Team. 2020, NEOWISE-R Single Exposure (L1b) Source Table, IPAC, doi: 10.26131/IRSA144
  * Park et al. (2014) Park, C., Jaffe, D. T., Yuk, I.-S., et al. 2014, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9147, Ground-based and Airborne Instrumentation for Astronomy V, ed. S. K. Ramsay, I. S. McLean, & H. Takami, 91471D, doi: 10.1117/12.2056431
  * Park et al. (2021) Park, W., Lee, J.-E., Contreras Peña, C., et al. 2021, ApJ, 920, 132, doi: 10.3847/1538-4357/ac1745
  * Pecaut & Mamajek (2013) Pecaut, M. J., & Mamajek, E. E. 2013, ApJS, 208, 9, doi: 10.1088/0067-0049/208/1/9
  * Pittman et al. (2022) Pittman, C. V., Espaillat, C. C., Robinson, C. E., et al. 2022, AJ, 164, 201, doi: 10.3847/1538-3881/ac898d
  * Rebull et al. (2015) Rebull, L. M., Stauffer, J. R., Cody, A. M., et al. 2015, AJ, 150, 175, doi: 10.1088/0004-6256/150/6/175
  * Robinson et al. (2022) Robinson, C. E., Espaillat, C. C., & Rodriguez, J. E. 2022, ApJ, 935, 54, doi: 10.3847/1538-4357/ac7e51
  * Rodrigo & Solano (2020) Rodrigo, C., & Solano, E. 2020, in XIV.0 Scientific Meeting (virtual) of the Spanish Astronomical Society, 182
  * Rodrigo et al. (2012) Rodrigo, C., Solano, E., & Bayo, A. 2012, SVO Filter Profile Service Version 1.0, IVOA Working Draft 15 October 2012, doi: 10.5479/ADS/bib/2012ivoa.rept.1015R
  * Rodriguez & Hillenbrand (2022) Rodriguez, A. C., & Hillenbrand, L. A. 2022, ApJ, 927, 144, doi: 10.3847/1538-4357/ac496b
  * Shappee et al. (2014) Shappee, B. J., Prieto, J. L., Grupe, D., et al. 2014, ApJ, 788, 48, doi: 10.1088/0004-637X/788/1/48
  * Sicilia-Aguilar et al. (2015) Sicilia-Aguilar, A., Fang, M., Roccatagliata, V., et al. 2015, Astronomy and Astrophysics, 580, doi: 10.1051/0004-6361/201525970
  * Sicilia-Aguilar et al. (2012) Sicilia-Aguilar, A., Kóspál, A., Setiawan, J., et al. 2012, Astronomy and Astrophysics, 544, doi: 10.1051/0004-6361/201118555
  * Sipos et al. (2009) Sipos, N., Ábrahám, P., Acosta-Pulido, J., et al. 2009, Astronomy and Astrophysics, 507, 881, doi: 10.1051/0004-6361/200911641
  * Somers et al. (2020) Somers, G., Cao, L., & Pinsonneault, M. H. 2020, ApJ, 891, 29, doi: 10.3847/1538-4357/ab722e
  * Sousa et al. (2023) Sousa, A. P., Bouvier, J., Alencar, S. H. P., et al. 2023, A&A, 670, A142, doi: 10.1051/0004-6361/202244720
  * Subasavage et al. (2010) Subasavage, J. P., Bailyn, C. D., Smith, R. C., et al. 2010, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 7737, Observatory Operations: Strategies, Processes, and Systems III, ed. D. R. Silva, A. B. Peck, & B. T. Soifer, 77371C, doi: 10.1117/12.859145
  * Tang et al. (2013) Tang, S., Grindlay, J., Los, E., & Servillat, M. 2013, PASP, 125, 857, doi: 10.1086/671760
  * Valenti et al. (1993) Valenti, J. A., Basri, G., & Johns, C. M. 1993, AJ, 106, 2024, doi: 10.1086/116783
  * Venuti et al. (2021) Venuti, L., Cody, A. M., Rebull, L. M., et al. 2021, AJ, 162, 101, doi: 10.3847/1538-3881/ac0536
  * Venuti et al. (2014) Venuti, L., Bouvier, J., Flaccomio, E., et al. 2014, Astronomy and Astrophysics, 570, 1, doi: 10.1051/0004-6361/201423776
  * Wang et al. (2023) Wang, X.-L., Fang, M., Herczeg, G. J., et al. 2023, arXiv e-prints, arXiv:2305.02514, doi: 10.48550/arXiv.2305.02514
  * White et al. (2020) White, J. A., Kóspál, Á., Hughes, A. G., et al. 2020, The Astrophysical Journal, 904, 37, doi: 10.3847/1538-4357/abbb94
  * Wright et al. (2010) Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868, doi: 10.1088/0004-6256/140/6/1868
  * Zacharias et al. (2013) Zacharias, N., Finch, C. T., Girard, T. M., et al. 2013, AJ, 145, 44, doi: 10.1088/0004-6256/145/2/44
  * Zakri et al. (2022) Zakri, W., Megeath, S. T., Fischer, W. J., et al. 2022, ApJ, 924, L23, doi: 10.3847/2041-8213/ac46ae
  * Zhai et al. (1994) Zhai, D. S., Foing, B. H., Cutispoto, G., et al. 1994, A&A, 282, 168
  * Zhou & Herczeg (2022) Zhou, L., & Herczeg, G. J. 2022, The Astronomer’s Telegram, 15271, 1
  * Zsidi et al. (2022) Zsidi, G., Manara, C. F., Kóspál, et al. 2022, Astronomy and Astrophysics, 660, doi: 10.1051/0004-6361/202142203

##  Appendix A Extinction Analysis of EX Lup

In this appendix, we present our combined photosphere and accretion continuum
fits to the VLT/X-Shooter spectrum, with an important downward revision to the
measured extinction. Briefly, we fit the EX Lup X-Shooter spectrum with
several photospheric templates to simultaneously determine the best extinction
and photosphere spectral type for EX Lup.

The first stage of the fitting is carried out individually to three
photospheric templates, TWA 25 (M0), TWA 14 (M0.5), and TWA 13A (M1), obtained
from Manara et al. (2013). For each template, the initial fit has two free
parameters: the extinction of EX Lup $A_{V}$ and a flux scaling factor $s$ of
the photospheric template. We first de-redden the observed EX Lup spectrum
$f_{\rm EX~{}Lup,~{}obs}$ with $A_{V}$, using $R_{V}=3.1$ reddening law
(Cardelli et al., 1989). The photospheric template $f_{\rm pho}\times s$ is
then subtracted from the dereddened spectrum $f_{\rm EX~{}Lup,~{}unred}$ to
yield the accretion continuum, $f_{\rm acc}$, as follows:

$f_{\rm acc}=f_{\rm EX~{}Lup,~{}unred}-f_{\rm pho}\times s$ (A1)

If we assume the accretion fluxes should be constant continuum over 4000-6000
Å, then the $\chi^{2}$ expression is defined as the goodness-of-fit

$\chi^{2}=\sum_{\lambda}\frac{\big{(}f_{\rm
acc,~{}4000-6000}-\texttt{med}(f_{\rm acc,~{}4000-6000}))^{2}}{(0.2\times
f_{\rm EX~{}Lup,~{}obs,~{}4000-6000})^{2}},$ (A2)

with $\texttt{med}(f_{\rm acc,4000-6000})$ as the median value of accretion
fluxes between 4000-6000 Å. The uncertainty in fluxes is treated as the 20% of
fluxes of the observed EX Lup spectrum $f_{\rm obs}$. The best-fit model is
found by minimizing the $\chi^{2}$ defined by Eq. A2 in the range of
$A_{V}=[0.0,1.2]$ mag with steps of 0.1 mag. The range of scaling factor $s$
differs depending on the templates being used.

Figure 20 compares the best-fit $A_{V}$ of EX Lup for three photospheric
templates. The best-fit $A_{V}$ of EX Lup is 0.3, 0.0, and 0.4 mag when using
TWA 25, TWA 14, and TWA 13 as the photospheric templates, respectively. A
change of 10% $\chi^{2}_{min}$ leads to an uncertainty of $A_{v}$ of 0.2-0.3
mag and a 20-30% change of photosphere fluxes. Among the three candidate
templates, TWA 14 is the optimal template that yields the lowest $\chi^{2}$
value. Both the TWA 13 and TWA 25 spectra have a different slope from
dereddened EX Lup spectrum at $\lambda>6500~{}\AA$, resulting in lower and
higher accretion fluxes than the accretion continuum at 4000-6000 $\AA$.

However, TWA 14 is shallower in the depths of several absorption bands between
5500-7500 $\AA$, compared with EX Lup. Since the depth of molecular absorption
bands is inversely proportional to the effective temperature of the
photosphere, and therefore the spectral type, we infer that EX Lup has a
spectral type between TWA 14 and TWA 25.

In the second stage of fitting, we interpolate a suite of synthetic
photospheric templates representative of spectral types between M0 and M0.5,
using a mixture of TWA 14 and TWA 25 spectra with varying weights and search
for a minimum $\chi^{2}$ on a new [$A_{V}$, $s$, $w$] grid.

The new best fit of $A_{V}$ is 0.1 mag with a photosphere template consisting
of 97% TWA 14 fluxes and 3% TWA 25 fluxes. The best model is shown in figure
21. Among the top-ranking models ($\Delta\chi^{2}/\chi^{2}_{min}=10\%$), the
range of $A_{V}$ and weights of TWA 14 spans 0.0-0.2 mag and 0.8-1.0, the
latter implying M0.4-0.5 spectral type of EX Lup. Models with earlier spectral
types favor larger $A_{V}$.

The $\chi^{2}$ of the best-fit model, 56, is not significantly improved
compared to the best model using TWA 14 as the template in Figure 20. In the
new best-fit model, accretion fluxes at $\lambda>7500\AA$ are more aligned
with the fluxes at 4000-6000 Å, though the inconsistencies in molecular band
depths between 5500-7500 $\AA$ still remain. These inconsistencies might come
from causes other than the photospheric spectral types.

We also tried a set of fits using plane parallel slabs to estimate the
accretion (Valenti et al., 1993). These slabs have bluer spectra than the
constant spectrum. We again find a best fit of $A_{V}$ = 0.1 mag with spectral
type M0.45 (Figure 22). The photospheric magnitudes, $A_{V}$, and accretion
spectrum used throughout the paper are taken from this best-fit result.

To summarize the result of our EX Lup extinction analysis, we find that EX Lup
spectra can be best fitted with $A_{V}=0.1$ mag, using a composite
photospheric template with a spectral type of M0.4-0.5. We adopt $A_{V}=0.1$
mag and M0.5 throughout this paper.

Figure 20: Comparisons of EX Lup VLT/X-Shooter spectrum with best-fit
extinction $A_{V}$ for three photosphere templates. The photosphere templates
are taken from the non-accreting YSOs’ VLT/X-Shooter spectra of Manara et al.
(2013), with names and spectral types indicated on the upper left. We assume
$A_{V}=0$ for these photosphere templates. The best-fit $A_{V}$ for each
template is obtained by minimizing the standard deviation (listed lower left)
of accretion fluxes (red) between 4000-6000 $\AA$ and is also displayed on the
upper right. The EX Lup spectrum (gray) is dereddened with best-fit $A_{V}$ in
each panel. The templates are added with constant fluxes (black) to compare
with EX Lup spectrum. Figure 21: The EX Lup spectrum dereddened with best-fit
$A_{V}$ using a mixture of TWA 14 and TWA 25 spectra. The best-fit $A_{V}$ is
0.1 mag with synthetic spectrum made up of 97% TWA 14 and 3% TWA 25. Figure
22: The same as Figure 21 but used an accretion model instead of assuming
constant accretion fluxes. The best-fit $A_{V}$ and spectral type for EX Lup
is 0.1 mag and M0.45. The parameters are adopted throughout the analysis of
this paper.

##  Appendix B Additional optical spectra of EX Lup

We also report several low-resolution optical spectra of EX Lup obtained
during quiescence, with flux calibration that was inadequate for full spectral
fits. A WHT/ISIS spectrum covering 3200-10000 Å at $R\sim 1000$ was obtained
on MJD 56144.869 at an airmass of 2.84, which prevented confidence in our flux
calibration. The spectrum shows a Balmer Jump with an accretion rate of 2–3
$\times 10^{-9}$ M⊙ yr-1 and with a broadband flux that is consistent with the
$\sim M0.5$ spectral type and low extinction.

We also obtained spectra on MJD 57155 and 57195 with the LCO 2m telescope at
Siding Spring Observatory. Neither the spectra could not be flux calibrated
because the observations were separated from the flux calibrator and other
problems with the instrumental setup. However, both show the Balmer Jump,
consistent with ongoing accretion, as well as the molecular features
consistent with the late-K/early-M spectral type.