Accretion variability in RU Lup

The process of accretion in classical T Tauri stars (CTTSs) has been observed to vary on different timescales. Studying this variability is vital to understanding a star's evolution and provides insight into the complex processes at work within. Understanding the dichotomy between continuum veiling and emission line veiling is integral to accurately measuring the amount of veiling present in stellar spectra. Here, 15 roughly consecutive nights of optical spectroscopic data from the spectropolarimeter ESPaDOnS are utilised to characterise the short-term accretion activity in the CTTS, RU Lup, and investigate its relationship with the veiling in the LiI 6707A absorption line. The accretion-tracing HI Balmer series emission lines were studied and used to obtain the accretion luminosity (Lacc) and mass accretion rate (Macc) for each night, which vary by a factor of ~2 between the brightest and dimmest nights. We also measured the veiling using multiple photospheric absorption lines (NaI 5688A, MnI 6021A, and LiI 6707A) for each night. We find the LiI 6707A line provides measurements of veiling that produce a strong, positive correlation with Lacc in the star. When corrected for Li depletion, the average veiling measured in the LiI 6707A line is r_LiI(avg)~3.25+/-0.20, which is consistent with the other photospheric lines studied (r_avg~3.28+/-0.65). We measured short timescale variability in the Lacc and Macc that are intrinsic and not due to geometric effects. Upon comparing the changes in veiling and Lacc, we find a strong, positive correlation. This study provides an example of how this correlation can be used as a tool to determine whether a measured variability is due to extinction or an intrinsic change in accretion. As the determination of veiling is an independent process from measuring Lacc, their relationship allows further exploration of accretion phenomena in young stars.


Introduction
Many processes at work in low-mass (M * ≤ 2 M ) young stars influence their formation and evolution.Two of the most important ones are accretion and ejection of material onto and emanating from young stars, such as T Tauri stars (TTSs).These TTS are low-mass pre-main sequence stars which are commonly divided into two main classes: classical T Tauri stars (CTTSs), which are younger (∼10 6 yr) and more actively accreting than weak-line T Tauri stars (WTTS), which are often older (∼10 7 yr) and at the end of the accretion phase (Walter et al. 1988;Bertout 1989;Hartmann et al. 2016, and references therein).The CTTSs demonstrate a rich emission spectrum mostly due to accretion of material theorised to be structured as funnels, flowing along the stellar magnetic field lines from the inner disk to the stel- The flux calibrated ESPaDOnS data shown in Fig. 3 and Tables 2, A.1, and A.2 are only available at the CDS via anonymous ftp to cdsarc.cds.unistra.fr(130.79.128.5) or via https:// cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/668/A94 Based on observations obtained at the Canada-France-Hawaii Telescope (CFHT) and the Cerro Tololo Inter-American Observatory (CTIO) Small and Medium Aperture Research Telescope Facility (SMARTS).lar surface (e.g.Hartmann et al. 2016, and references therein).Other spectral lines come from the ejection phenomena.They are varied in their structures (e.g.jets and winds) and velocities, depending on the conditions from which they originate; however, it has been shown that a strong link between ejection and accretion exists (e.g.Ray et al. 2007, and references therein).
The accretion luminosity (L acc ) of a young star is produced by the shocked material falling onto the stellar photosphere.The accretion energy released is seen in the form of continuum excess in emission -mostly at ultraviolet (UV) and optical wavelengths -and emission lines.Historically, the L acc of low-mass young stars is estimated by modelling the Paschen and Balmer continua emission as isothermal slabs (Valenti et al. 1993;Gullbring et al. 1998).More often, L acc has been computed from empirical relations between the luminosity of accretion-tracing lines and L acc (see e.g.Gullbring et al. 1998;Muzerolle et al. 1998;Natta et al. 2006).Alcalá et al. (2017) derived more precise relationships between the line luminosity (L line ) of several optical and near-infrared (NIR) accretiontracing emission lines and L acc .
This physical relationship between L line and L acc allows the estimation of the mass accretion rate ( Ṁacc ) Notes. (* ) Calculated from radial velocity measurements. (* * ) Disk inclination with respect to the plane of the sky.(Hartmann et al. 1998).The Ṁacc of CTTSs was calculated by Herczeg & Hillenbrand (2008) by analysing the excess UV and optical emission.A decade later, estimates of the Ṁacc for a large survey of young stars by Alcalá et al. (2017) was derived.
Accretion in young stars is not a steady process, but rather variable (for a review see Fischer et al. 2022, and references therein).Understanding how young stars' accretion activity varies is vital for determining how they gain their mass and what implications that process has on their evolution (Herbst et al. 1994).
Indeed, luminosity amplitude and timescales vary depending on the origin of such variability (see Fig. 3 of Fischer et al. 2022).Temporal changes over timescales of days and amplitudes less than one order of magnitude in the accretion rate (the so-called routine variability) are often attributed to the rotation of the star and are very useful in measuring its rotation period (e.g.Scholz & Jayawardhana 2006;Costigan et al. 2012Costigan et al. , 2014;;Rebull et al. 2014).However, this is not exclusively the case as the rotation contribution might hide any intrinsic changes that come from the stellar-disk magnetospheric interaction.On the other hand, the so-called accretion bursts with larger amplitude (more than one order of magnitude) and timescale variability (from months to years) are more related to the disk's mass transport and refurbishment, and they might originate from disk instabilities and other related mechanisms relevant to EXor and FUor phases (Hartmann et al. 2016).
Therefore, studying the variation of L acc and Ṁacc is important in understanding how TTSs form.Magnetospheric accretion can be divided into two regimes: stable and unstable.These two regimes were predicted by several 3D magnetohydrodynamic (MHD) simulations of disk accretion onto a rotating young star (Romanova et al. 2008;Kulkarni & Romanova 2008, 2009).According to these models, the main factors that determine in which regime a star accretes are the mass accretion rate, the strength of the magnetic field, the stellar rotation rate, and the misalignment between the star's rotation axis and its magnetic poles.
Classical T Tauri star accretion activity and variability can also be studied by using veiling.This is the filling-in of stellar photospheric absorption features such that they appear shallower in observed spectra of an actively accreting young star, due to an accretion-fueled flux excess (both in the continuum and in lines).This phenomenon has been investigated since Joy (1949) first introduced the concept of veiling.With high resolution spectroscopy, it is possible to observe the photospheric absorption lines of the star and this excess due to accretion that effectively veils the spectrum.
This study focusses on the CTTS RU Lup, which is known to have a plethora of strong accretion features (Herczeg et al. 2005;Gahm et al. 2008;Siwak et al. 2016, and references therein).RU Lup is an accreting CTTS of 2−3 Myr that lies in the Lupus star-forming region at a distance of d ≈ 159 pc (Gaia Collaboration 2018).It has a rotational period of 3.71 ± 0.01 (Stempels et al. 2007) and is notable for its strong emission lines and veiling.These accretion-tracing lines are exceptionally broad, for example, the equivalent width (EW) of the Hα emission line is usually found to be greater than 100 Å (Siwak et al. 2016).In their models, this high accretion activity was found by Lamzin et al. (1996) to be responsible for the majority of the observed luminosity, whereby only about 10% is generated by gravitational contraction of the star itself.The Ṁacc was calculated by Herczeg & Hillenbrand (2008) and found to be ∼1.8 × 10 −8 M yr −1 .The study by Alcalá et al. (2017) resulted in values of ∼6.7 × 10 −8 M yr −1 .This spread is indicative of the amount of accretion variability exhibited in RU Lup on the timescales of years.For this study we measured the Ṁacc of RU Lup during 2011 and determined how it changes over the shorter timescale of two weeks.By investigating its short-time stellar-disk magnetospheric interaction, our dataset allowed the exploration of how RU Lup's accretion fits into the category of routine variability, as delineated in Table 1 of Fischer et al. (2022).
Additionally, multiple studies have inferred that RU Lup has a near face-on geometry from varying methods including the analysis of line profiles and measured rotational velocities (see Sect. 2 of Herczeg et al. 2005, and references therein) or the direct measurement of the inner dusty disk geometry and size with NIR interferometry (GRAVITY Collaboration 2021).Because it is seen mostly face-on, RU Lup is an excellent candidate for a study of accretion variability as there are fewer geometrical effects to confound the observations (such as stellar rotation and extinction from the disk).It is also a well-studied and interesting object undergoing strong accretion activity.Stellar properties of RU Lup such as its radius (R ), spectral type (SpT), stellar luminosity (L ), rotational period (P rot ), v sin i, and visual extinction (A V ) can be found in Table 1.CFHT/ESPaDOns is a fibre-fed echelle spectrograph and spectropolarimeter that operates across the entire optical wavelength range (370 nm to 1050 nm).These observations were taken in 'object + sky' mode with a spectral resolution of ∼68 000 and spectral sampling bin size of 1.8 km s −1 .The ESPaDOnS data were reduced with 'Upena', the standard pipeline provided by the CFHT, based on the Libre-ESpRIT package (Donati et al. 1997).Additionally, because no standard star observation was available, the telluric lines in the regions of interest to this study were fitted and removed manually.This included the region in and around the [O i] 6300 Å line, that is the stronger component of the doublet that also includes [O i] 6364 Å, containing strong telluric contamination.This was done by first fitting and subtracting the continuum from the spectra in this region, and then, utilising libraries such as Astropy in Python, each telluric absorption that was blended with the [O i] 6300 Å line emission was identified and fitted with a Gaussian.These Gaussian fits were then subtracted from each night's spectrum.

Observations and data reduction
In addition to telluric contamination, the data contain many photospheric absorption lines, typical of a K7 star.Some of these were blended with the [O i] 6300 Å line and required removal to study its profile.Unlike the telluric absorption lines, it was possible to remove all the photospheric lines in this region (6284−6342 Å) at once.This was done by obtaining a WTTS spectrum of the same spectral type as a template, which contains the same photospheric lines but no emission lines (Martín 1998), and subtracting this from the CTTS spectrum of RU Lup.As WTTSs undergo little-to-no accretion activity (which produces many emission lines), their spectra can be used as templates for the photosphere of CTTSs of the same spectral type.Indeed, the RU Lup spectra was compared with the WTTS spectrum of TAP 45 (same spectral type and similar v sin i = 7.7 ± 0.7 km s −1 ; Nguyen et al. 2012).TAP 45 has a stellar mass of ∼0.79 M , a radius of ∼1.14 R and T eff ∼ 4040 K (Grankin 2013).This template was then re-scaled such that it met the level of the RU Lup continuum; these parameters were adjusted by minimising a χ 2 fit.The resultant broadened and veiled template represents a spectrum containing the fitted photospheric features which can then be used to subtract the contaminating photospheric lines from the regions of interest.This process was completed on a night-by-night basis.
Finally, the radial velocities of the observed lines were calculated using single or multiple Gaussian fits (the detected lines typically show multiple components).These velocities are measured with respect to the systematic velocity of the star (v = +8 ± 2 km s −1 , Takami et al. 2001) in the local standard of rest (LSR) frame.

Photometry and spectroscopic flux calibration
As the ESPaDOnS observations were not flux calibrated, we searched archival data and found simultaneous photometry of RU Lup from the 1.3 m SMARTS/CTIO 1 telescope for the night of 17 June 2011.These ANDICAM photometric data were taken on 17 June 2011, allowing the ESPaDOnS spectra from the same night to be flux calibrated.The data from the B, V, R and I filters were reduced in Python using the 'Astropy' and 'ccdproc' packages to apply biases, darks and flat-fielding (dome flats for V, R and I filters and sky flats for the B filter) to the raw images of RU Lup and the standard stars observed from the same night.The resulting standard stars were flux calibrated using the zero-point fluxes given for the Kitt Peak National Observatory (KPNO) filters used on ANDICAM, and were then used to flux calibrate the science images.The calibrated spectrum of the 17 June 2011 night is shown in Fig. 1, together with the photometric data.
The ESPaDOnS instrument provides spectroscopic data from the B to the I band (3700−10 500 Å).We have computed the synthetic photometry for these bands, using the transmission response of ANDICAM (triangles in Fig. 1).While the B to R synthetic photometry matches the aperture photometry (circles) well, the I band synthetic photometry is about a factor of two higher than the measured value.We believe this to be due to an instrumental issue, and therefore have limited the spectral range we use in this study to the B through R region, where the data are well-calibrated.
To the best of our knowledge, no photometric data are available for the other nights.Therefore, we looked to the from night-to-night, and very small night-to-night flux variations.The profiles of this line show the typical features of forbidden lines in CTTS, namely a high-and a low-velocity component (HVC and LVC, respectively) wherein the HVC is not expected to show any variability on these short timescales (see Sect. is not expected to vary over these timescales either.We estimate that the uncertainty of the calibrated flux is of the order of 5%. In the specific case of RU Lup, our choice is also supported by the results of Gahm et al. (2013), which found that the forbidden line profiles remained remarkably constant over their study's longer timescale of five years.The authors also note that in their flux calibrated spectra, the flux of the [S ii] 6730 Å line remains stable for relatively moderate veiling (r ≤ 2; see Fig. 2 in Gahm et al. 2008).Indeed, this is the case in our data set (see Sect. 3.3).
Our assumptions separately work on both [O i] 6300 Å and [S ii] 6730 Å line flux calibration.Our calibrated spectra are presented in Fig. 3, where the main lines used in this paper are labelled.The [S ii] 6730 Å emission line is the stronger component of the S ii doublet that also includes [S ii] 6716 Å.

Results
Our multi-epoch calibrated spectra are presented in Fig. 3.The spectra, from the B to the R band, show variations in the continuum flux.Many photospheric lines in absorption are present, indicative of a K7 star (White & Basri 2003).Among them, the Li i 6707 Å transition is the most prominent.Additionally, we detect an overwhelming number of emission lines, typically originating from the circumstellar environment of RU Lup such as accretion shocks, winds, jet, and chromosphere.For the purposes of this paper, we select a number of lines that are widely used to trace accretion and ejection activity in CTTS.

Accretion-tracing lines
As previously mentioned, the accretion activity in young stars such as RU Lup produces line emission that allows the study of these processes.There are many known accretion-tracing emission lines at optical wavelengths, and their intensities, profiles, and relative red-and blue-shifts provide information about their origins.
The brightest accretion-tracing lines present in our data are The trend in flux variability seen here is similar to that for the Balmer series.There is an increase in the integrated intensity between the dimmest and brightest nights by a factor of ∼1.9 for He i 5875 Å and ∼2 for He i 6678 Å and He i 7065 Å.
As has been shown by Beristain et al. (2001), the He i 5875 Å line can be decomposed into a narrow component (NC) and a broad component (BC).The NC is generally slightly red-shifted and is believed to originate in the post-shock region near the stellar surface.The BC has been theorised to have multiple origins, with contributions from the accretion columns as well as from a hot wind.We have decomposed the He i 5875 Å and He i 6678 Å lines into two components in our data.Though the components vary from night to night, both profiles indeed show a clear NC and BC every night (Fig. rise to both components correlate well with each other, and therefore with accretion.

Time-dependent accretion luminosity and mass accretion rate
We derive L acc from the luminosity of the accretion-tracing lines, using the well established relations between each line luminosity (L line ) and L acc derived in Alcalá et al. (2017) for CTTSs.L acc is then computed for each night from the average of the L acc values derived from the individual lines (see Fig. E.1).We note that because of the way the L line vs.L acc relation was derived in Alcalá et al. (2017), L acc refers to the continuum excess emission only.
The inferred L acc values found in the He i differ slightly from those of the H i lines.This is likely due to the complexity of processes that the He i lines are tracing, which include accretion, but not exclusively.In particular there can be a strong wind emission (Beristain et al. 2001;Johns-Krull & Basri 2001;McGinnis et al. 2020).Therefore we have decided to estimate L acc using the H i emission lines only.
Our values for L acc range from 0.35 to 0.72 L over the course of 15 days, with an average of 0.5 L (see Table 2).It is worth noting that the observed variability in the accretion luminosity is real, as the line flux variation of the different H i lines is much larger than their flux uncertainties (see Tables A.1 and A.2).
Variations in L acc can in principle be due to variations in extinction.However, there is no evidence in the literature of significant variations in visual extinction which is always very small (e.g.Herczeg et al. 2005;Alcalá et al. 2017) and consistent with this system being oriented face-on.Moreover, Fig. E.1 shows there is no trend of increasing L acc from lines of increasing wavelength.Therefore, we assume in the following that changes in L acc trace physical changes in the accretion rate.
Because this variability cannot be caused by changes in extinction, it must come from accretion variability.The derived range of L acc values is in line with those found in the literature, derived with different methods (e.g.slab modelling of the Balmer jump, luminosity of H i lines, UV excess), namely ∼0.5 L (Herczeg et al. 2005;Antoniucci et al. 2014;Alcalá et al. 2017).
The corresponding value of the mass accretion rate is computed as: where R in is the inner (truncation) radius.The truncation radius is assumed to be R in = 5 R (Gullbring et al. 1998).Values of M and R are as in Table 1.
The Ṁacc ranges from a minimum of ∼3 × 10 −8 M yr −1 to a maximum of ∼7 × 10 −8 M yr −1 .When averaged across the 15 nights, the result is Ṁacc = 4.8 × 10 −8 M yr −1 , in line with the values found in the literature.Small differences are due to slightly different distances (140−160 pc) and masses (0.5−1.1 M ) adopted for RU Lup by different authors.This kind of variability aligns with what Fischer et al. (2022, Table 1) characterise as routine variability wherein there are day-to-day changes of <1−2 mag in the optical and IR.This variability in CTTSs is described in Fischer et al. (2022) as a result of rotational, magnetospheric or inner-disk interactions with the star.In this study, we discuss the impact stellar rotation has (or does not have) on RU Lup's variability in Sect.4.1.
We want to stress here that while our values of Ṁacc depend on the stellar parameters adopted, the time variations of Ṁacc , and those of L acc , do not depend on them.

Time-dependent veiling from metal lines
Additional information on the accretion variability of RU Lup can be obtained by measuring the excess emission in photospheric lines (veiling), as described in Eq. ( 2).The excess emission in a given line is the sum of the continuum emission from the accretion shock at the line wavelength and of the emission of the same line in accreting matter (Hartigan et al. 1989).By measuring the amount of veiling in a CTTS, it is possible to quantify the amount of accretion activity occurring in the star during a given observation.
In Hartigan et al. (1989), a quantitative method was developed to measure the amount of veiling present in a given CTTS.By comparing the spectrum of the accreting CTTS, with a template spectrum that represents only the photosphere (no accretion Notes.The L acc and Ṁacc for each night are determined by averaging the values found from the H i lines, and their uncertainties are estimated from their spread.Columns 5 and 6 contain the veiling measurements for the photospheric lines, Na i 5688 Å and Mn i 6021 Å.The observed veiling measured in the Li i 6707 Å photospheric line is given in Col. 7 (r Li i (obs)) and the abundance corrected veiling is given in Col. 8 (r Li i = r Li i (obs) + ∆r abund.corr. ).The final column contains the χ 2 values of the fits done to calculate the veiling in the Li i 6707 Å absorption line.
( * ) r Li i = r Li i (obs) + ∆r abund.corr., where ∆ r abund.corr.= 1.8, see Sect.present), the dimensionless quantity, veiling, is produced: where r j is the veiling in a specific wavelength interval j, F excess is the excess flux (due to accretion activity) and F photo is the photospheric flux, given from the template star.WTTS are often used for this purpose (e.g.Hartigan et al. 1991;Valenti et al. 1993;Johns-Krull et al. 1999;Dodin & Lamzin 2012;McGinnis et al. 2018).
The veiling present in a specific photospheric line can be similarly derived using the EWs of that line in the CTTS spectrum and the WTTS spectrum: where EW observed is the equivalent width of the observed line in the CTTS spectrum and EW photo is the equivalent width of the same line in the template (WTTS) spectrum.Notably, the numerator in Eq. (3) represents the total flux observed (the sum of the photospheric and excess fluxes), whereas the numerator in Eq. ( 2) represents only the excess flux.
Veiling, being the ratio of the excess flux over the photospheric flux, does not depend on extinction.Although it is not the case for RU Lup, this is very important, as the extinction is sometimes uncertain and it is known to vary with time in a large number of pre-main sequence stars (e.g. the dippers studied by Contreras Peña et al. 2017).
Photospheric lines of different species and properties may exhibit different veiling if accretion-related line-emission is important, as proposed, for example, by Dodin & Lamzin (2012).
This study aims to add to the conversation concerning the distinction between continuum veiling, that is, the contribution to the measurement of veiling due to a continuum excess, and the contribution of specific line emission.We have, therefore, investigated veiling in individual lines in a clean region of our spectra to attain the most accurate measurements possible.This allows  us to delve deeper into how veiling from metal lines relates to the accretion luminosity.
Therefore, we focus our analysis on three intrinsically deep metal lines in the region 5600−6700 Å, namely the Li i 6707 Å, Mn i 6021 Å and Na i 5688 Å lines.The latter were chosen as they are within our wavelength range of interest and have a sufficiently high signal-to-noise ratio (S/N) for veiling measurements.
The Li i 6707 Å line seen in absorption in young K7 stars has the deepest absorption in the optical wavelength range.It has an average EW Li i = 0.17 ± 0.03 Å.We analysed this line by comparing it to a K7 WTTS and estimated the amount of veiling present as described in Sect. 2. An example of our fits is provided in Fig. 4, and in Appendix C our method for calculating veiling is explained in detail.Preliminary analysis of the observed veiling produced an average over the 15 nights of r Li i, avg (obs) = 1.92 ± 0.03.This observed value is affected by there being different abundances of Li in the template star and RU Lup, as discussed in Sect.4.3.The abundance-corrected A94, page 6 of 17 average is r Li i, avg = 3.25 ± 0.20, which is used throughout this study.The nightly results are shown in Col. 8 of Table 2.
Veiling measurements of the Na i 5688 Å and Mn i 6021 Å photospheric lines can be found in Cols. 5 and 6 of Table 2 and their averages over the 15 nights are r Na i, avg = 3.27 ± 0.63 and r Mn i, avg = 3.29 ± 0.67.As these photospheric lines do not have as high an S/N as that of the Li i 6707 Å, and are intrinsically weaker, there is higher uncertainty in the measured veiling.
The additional photospheric metal lines Ti i 6261.10Å, Ca i 6431.10Å, and Ni i 6643.63Å were initially considered as well; however, these lines were too veiled to get accurate measurements.This illustrates the difficulty in measuring the veiling in high accretors as there is a balance wherein the veiling is high enough to be noticeable in the spectrum, but is not so high that the photospheric lines are completely 'washed out'.

Wind and jet lines
In the studied spectral range, there are multiple forbidden emission lines associated with outflowing material from young stars (Edwards et al. 1987;Hamann 1994;Hartigan et al. 1995).
Present in these data with a strong S/N are the [O i] 6300 Å and [S ii] 6730 Å emission lines.The flux calibrated spectra of these lines from the 12 epochs can be seen in Fig. 2 ).These are quite common features in the jets and winds of CTTSs (Rigliaco et al. 2013;Simon et al. 2016;McGinnis et al. 2018).
While the main component of the HVC is tracing the collimated jet, the broad component is often associated with the base of the jet or a disk-wind and the NLVC with a slow-moving compact disk-wind or a photoevaporative wind from the disk (see e.g.Ray et al. 2007;Simon et al. 2016).In the case of RU Lup, the different nature of such components has been recently confirmed by the spectro-astrometric analysis of Whelan et al. (2021).Here the authors showed that the HVC and the broad component come from the extended jet and base of the jet, respectively, whereas the NLVC originates from a very compact ( 8 au) MHD diskwind.
Interestingly, our measurements of the HVC and LVC radial velocities match very well the total jet and wind veloc-ities measured or inferred in other CTTSs of similar mass (see e.g.Ray et al. 2007, and references therein).This is a further indication that the disk geometry is almost face-on (16 • ± 6 • 8 • ; GRAVITY Collaboration 2021), namely the outflow is approaching the observer with a small inclination angle.
Indeed, because of the disk geometry, very little of the redshifted emission from the outflow is observed, as it is located on the side of the star facing away from us, hidden by the disk.This explains why the red-shifted HVC is not detected in our spectra (see Fig. 2).

Discussion
The results presented in Sect. 3 show that accretion in RU Lup varies from night-to-night over a period of 15 nights by about a factor of two.All the analysed accretion tracers (emission line luminosities, profiles, veiling) give consistent results.Notably, the average value of Ṁacc in this period agrees very well with values found in the literature at different epochs (Alcalá et al. 2017).It is therefore likely that the variability pattern we found is typical of RU Lup's accretion at this stage of its evolution.
In addition, as we mentioned in Sect.3.2, the extinction of RU Lup is negligible and has not been taken into consideration in this study.As the EW of each absorption line analysed changes slightly over our short timescales, this indicates that such variability is due to veiling and not to extinction.
As Fig. 2 (Whelan et al. 2021).This star therefore lacks the BLVC which is linked to an MHD disk wind in the innermost part of the disk.

L acc relation to the stellar rotation
Variation in the observed L acc is sometimes associated with the stellar rotation if the line of sight intercepts a single or multiple accretion spots (Costigan et al. 2014), the properties of which are relatively stable over several rotation periods.For instance, Alencar et al. (2018) show that the veiling measured in LkCa15 varies by a factor of ∼5 over the course of a few rotation cycles, and that this variability folds well in phase with the star's rotation period (see their Fig. 5).Similarly, McGinnis et al. (2020) show that the veiling variability and, even more pronounced, the flux  of the He i emission line at 5876 Å of the star IP Tau have a clear modulation at the star's rotation period.However, given the nearly face-on geometry of RU Lup's disk (Gahm et al. 2013; GRAVITY Collaboration 2021), we do not expect that the observed variability is related to stellar rotation.To confirm this hypothesis, we study how L acc varies as a function of the stellar rotation period.In Fig. 5, the variability of L acc is presented across the roughly two weeks of the observations (upper panel) and folded in phase with the rotation of the star (lower panel) which has a period of 3.71 days (Stempels et al. 2007).To quantify the amount of L acc variability potentially due to rotation, we fit a sinusoidal curve to the folded data.The resulting fit (red-dashed curve) is shown in the lower panel of Fig. 5 and the fit parameters are listed in the figure's caption.Similar plots are shown in Fig. 6 for the veiling derived from the Li i 6707 Å line (r Li i ) variability vs. time and folded in phase with stellar rotation (upper and lower panel, respectively).
If the variability seen in these two quantities was solely due to RU Lup's rotation, the lower panels of Figs. 5 and 6 would show a distinct sinusoidal trend.However, the data vary in a more complex way, indicating that an intrinsic variability in the star's accretion exists.Indeed, the curve amplitude in the fit of Fig. 5 is only L acc = 0.04 L implying that, if present, only 7% of the L acc variability can be attributed to stellar rotation.

Intrinsic accretion variability
As the bulk of RU Lup's accretion variability cannot be due to stellar rotation nor to extinction, it must then be related to variations of the mass accretion rate (a factor of 2 or 0.31 dex, see Table 2) onto the star during our 15-day spectroscopic monitoring.On our observational timescale, such variability is similar to what was found in other photometric and spectroscopic studies of CTTSs (Costigan et al. 2012;Venuti et al. 2014;Zsidi et al. 2022).For what concerns RU Lup's accretion variability, long-term photometric monitoring by Siwak et al. (2016) already showed stochastic variability and lack of a single stable periodic pattern in its light-curves.The authors inferred that an unstable accretion regime, typical of CTTSs with high-mass accretion rates, operates in RU Lup.In particular, they argue that RU Lup's photometric variability is due to increased mass accretion rates that can prompt Rayleigh-Taylor (RT) instabilities, moving accretion from a stable to an unstable regime, as predicted by several 3D MHD simulations of disk accretion onto a rotating magnetised star (see Romanova et al. 2008;Kulkarni & Romanova 2008, 2009).
Although several ingredients can influence such a transition (strength of the magnetic field -B, stellar rotation rate, misalignment between the star's rotation axis and its magnetic poles), Romanova et al. (2008) showed that the main player in the development of RT instabilities is the mass accretion rate.
During the stable regime, the accretion disk is blocked by the star's magnetic field, and, at the magnetospheric radius (∼5 R * ), the accreting matter flows around the magnetosphere, forming two funnels that deposit matter near the star's magnetic poles, making accretion relatively stable.An increase in accretion rates and thus in density squeezes the magnetosphere in the equatorial region.As a consequence, it is energetically impossible for the inner disk matter to follow the resulting magnetospheric field lines.RT instabilities develop, causing matter to move towards the star across the field lines, until it reaches field lines that are energetically possible to follow.The matter is wound along these lines, forming miniature funnel flows and accreting at lower latitudes (Romanova et al. 2008;Kulkarni & Romanova 2008, 2009).
A94, page 8 of 17 We can compare our Ṁacc results with the predictions from Romanova et al. (2008), Kulkarni & Romanova (2008, 2009) and verify whether RT instabilities might be the source of the observed variability in RU Lup.The Ṁacc boundary value between the stable and unstable regime in a CTTS is provided by Eq. (1) of Kulkarni & Romanova (2008) and depends on the B, M and R .For our calculations, we use the stellar parameters reported in Table 1.The strength of RU Lup's magnetic field is not well constrained in the literature; however, Johnstone & Penston (1986) provide an upper limit of ∼500 G, which we adopt as B of RU Lup.Equation (1) of Kulkarni & Romanova (2008) with RU Lup's parameters gives Ṁacc ∼ 1.6 × 10 −8 M yr −1 as the critical mass accretion rate whereby the stable regime shifts to an unstable one.This Ṁacc is lower than the observed range of Ṁacc values in our target.Notably, lower B values would provide a lower boundary value for Ṁacc .We can thus confirm the possibility an of unstable accretion regime in RU Lup.Kurosawa & Romanova (2013) employed the results of the 3D simulations of Romanova et al. (2008) to compute timedependent accretion rates and calculate Hydrogen line profiles for a CTTS slightly more massive (0.8 M ) and more magnetised (1000 G) than RU Lup, and with a slightly larger stellar radius (2 R ).As expected from Eq. ( 1) of Kulkarni & Romanova (2008), the accretion rates derived by Kurosawa & Romanova (2013) for such a CTTS in the unstable regime are slightly higher than ours (see lower panel of Fig. 2 of Kurosawa & Romanova 2013).However, the relative Ṁacc variations are comparable to the observed variations in RU Lup.

Veiling in the Li I 6707 Å line
The Li i 6707 Å photospheric line is the deepest absorption line in our spectra and benefits from a high S/N, making it a more enticing photospheric line with which to study the veiling.However, this use of the Li i 6707 Å line potentially suffers from the uncertain Li abundance (A(Li) = log N(Li)) of pre-main sequence stars, due to the onset of Li nuclear burning.The Li abundance decreases as a star gets older, in a way that depends on its stellar properties, that is, mass and temperature (Palla et al. 2007).
It has been hypothesised that WTTSs are simply more evolved CTTSs with dissipated disks (Walter et al. 1988;Bertout 1989).However, there are CTTSs and WTTSs that coexist (e.g.Fig. 15 of Kenyon & Hartmann 1995), and their disk evolution is informed by their star-forming region and circumstellar environment (Bonnell et al. 2006;Rosotti et al. 2014).The details of this complex evolution are beyond the scope of this study, however, the implication is that the amount of Li depletion (a known indicator of stellar age) is often higher in WTTSs than CTTSs, but cannot be assumed (Galli et al. 2015).
Therefore, the template used to measure the veiling needs to take the Li abundance of the target CTTS and the template WTTS into account.For example, if the Li in the template star is depleted, there will be less Li in the photosphere to absorb light and the EW of the Li i line will be less than an un-depleted star.Because the flux of the Li i 6707 Å line in the WTTS is compared with that of the CTTS to evaluate veiling, the abundances of Li must be similar, if not equal, (A(Li) WTTS A(Li) CTTS ) to accurately identify the difference in line depth as a veiling effect.
RU Lup has had its abundance measured as A(Li) > 3.81 dex, which indicates no depletion of Li (Biazzo et al. 2017).On the other hand, TAP 45 has A(Li) values ranging from 2.92 to 3.27, depending on the method used for determining it (see Col. 8 of Table 1 in Sestito et al. 2008, their source #51).This means that our template has an A(Li) from 0.5 to 0.9 dex lower than RU Lup, leading to an underestimate of r Li i .
In order to account for the difference in Li abundance between our WTTS template and RU Lup, we derive a correction factor by comparing the theoretical unveiled EW Li of RU Lup (Biazzo et al. 2017) to the EW Li that we measured from our template star, TAP 45.Biazzo et al. (2017) analyse a sample of YSOs, which includes RU Lup, that were observed with X-shooter (this sample is described in Alcalá et al. 2017;Frasca et al. 2017).They measure the EW Li for each star from the observed data and then unveil these equivalent widths using the values of r 6700 Å estimated in Frasca et al. (2017) from the code ROTFIT (e.g.Frasca et al. 2006Frasca et al. , 2015)).This value of veiling is extrapolated between the values Frasca et al. (2017) found for r 6200 Å and r 7100 Å .The notation EW corr.
Li, CTTS represents the equivalent width of the Li i 6707 Å line in a CTTS with the veiling effect removed (or 'unveiled') (see Sect. 3.1.1 of Biazzo et al. 2017).The value found for RU Lup in Biazzo et al. (2017) is EW corr.Li, RU Lup = 988 ± 10 mÅ.This is illustrated in Fig. 7, where the teal line represents the observed Li i 6707 Å line of RU Lup, the orange line represents the template star, TAP 45, and the dashed light blue line represents the absorption line of RU Lup when unveiled.The difference between the orange line's EW and the light blue line's EW is the difference of veiling measurement between the two templates.We call this the abundance correction factor (∆ r Li ) as it represents the discrepancy between the WTTS template and the template EW Li of RU Lup unveiled: a difference only in Li abundance.
This gives ∆ r Li = 1.8.A resultant total veiling can be found It should be noted that there are large uncertainties in this correction factor, because of the way that the unveiled EW Li of RU Lup was determined.RU Lup is a strongly veiled CTTS in which continuum veiling is likely not the dominant form of veiling (Gahm et al. 2008), meaning that individual photospheric absorption lines suffer from different amounts of veiling due to the filling in of lines by line emission (this will be explained further in the next section Sect.4.4).Therefore using an extrapolation of the veiling from other photospheric lines is not very accurate to determine the amount of veiling in the Li i line.However, the similarities between the spectral types, v sin i, and other stellar properties of RU Lup and the template make the corrected measurements good approximations.In any case, this effect does not affect the relative variation observed in the veiling of the Li i line and the results shown in Fig. 8.
Another method to determine veiling in the Li i line would be to simply compare the EWs of the Li i line in our observed spectra with the EW Li given by Biazzo et al. (2017) using Eq. ( 3).However, we chose instead to use the same method as was used to determine veiling for the Na i and Mn i lines and calculate a correction factor for the Li i abundance to maintain self-consistency.Regarding the measurements of these other two metal lines, it is preferable to compare the observed lines with a WTTS spectrum to determine veiling, rather than use a theoretical unveiled EW of each line in Eq. (3) as was done for Li i by Biazzo et al. (2017).There should not be a measurable difference in the abundance of these elements between our template and RU Lup, and therefore extrapolating values of veiling found from several photospheric lines by Frasca et al. (2017) to get an unveiled EW would only insert additional uncertainties to our measurements.
Even given these caveats, the most interesting finding of this study stands, the relative variation of the veiling with time does not depend on the uncertainties of the Li i 6707 Å abundances in RU Lup itself nor in the template.
Figure 8 shows our result that there is a close proportionality between r Li i and L acc for RU Lup.The absolute value of the accretion luminosity for a given value of r Li i may depend on a number of stellar and accretion properties, and the possibility to derive L acc from the Li i 6707 Å veiling alone needs to be explored further.However, our results suggest that in any individual star the variation of r Li i may be correlated with variabilities in L acc .

Continuum and line veiling
As noted in the previous section, the veiling values of the three metal lines we have studied roughly agree within the uncertainties.However, a closer inspection shows that, while there is a rather tight correlation, the values of the veiling of different lines are not proportional to each other, as shown in Fig. 8.In a system with only continuum excess emission, the veiling in these lines, which spreads over a relatively small interval of wavelengths, should present a 1:1 relationship with one another, as we expect that the continuum will change by a very similar amount over this interval.This, however, is not going to be the case if accretion-related emission in the metal lines is also contributing to the overall veiling values (see Fig. 9).
While the Na i 5688 Å line roughly traces that of the Li i 6707 Å, the Mn i 6021 Å veiling increases sharply as the Li i 6707 Å veiling increases.This implies that the Mn i 6021 Å line emission increases more strongly with an increase in L acc than the Li i 6707 Å and Na i 5688 Å line emission.However, Li i 6707 Å, Na i 5688 Å, and Mn i 6021 Å all show strong reactions to an increase in accretion in RU Lup.
This accretion-related emission in metal lines, 'line veiling', has been observed in a number of accreting TTS (e.g.Bertout 1984;Hartigan et al. 1989;Gahm et al. 2008).In some, the line veiling has been found to be dominant, among them RU Lup (Gahm et al. 2008(Gahm et al. , 2013)).The authors compare the V-band photometric variability of the star with the variability of the veiling and find that the veiling is not correlated, or only weakly related, to the stellar brightness during states of high veiling.They argue that there are two sources of veiling: a continuous excess emission, and narrow emission lines filling-in the photospheric absorption lines, both related to the accretion shocks and foot-prints of the accretion funnels.Line emission becomes the dominant component for high-veiling states (r > 2).In this case, a change in veiling does not lead to a corresponding change in stellar brightness, as it would happen if the veiling was caused by continuous emission alone.
Indeed this lack of correlation in highly veiled, and thus highly accreting CTTSs, is confirmed by Rei et al. (2018).For a small sample of young active stars, Rei et al. (2018) find that the veiling is strongly line-dependent and is larger in strong photospheric lines and weaker or absent in the weakest ones.This concurs with what we have found in the three metal lines, namely, the stronger lines exhibit stronger reactions to increases in accretion.

Summary and conclusions
We have investigated the accretion variability of RU Lup over time through the H i line luminosities and the veiling present in the photospheric metal lines Li i 6707 Å, Na i 5688 Å, and Mn i 6021 Å.Our main results are the following: -L acc and Ṁacc were observed to vary by up to a factor of ∼2 over the 15 nights.This variability is an intrinsic variation of the accretion rates.In particular, it does not reflect variations related to the stellar rotation.-The average value of the mass accretion rate found is Ṁacc, avg = 4.8 × 10 −8 M yr −1 .This is within the spread of values found in Herczeg & Hillenbrand (2008) (1.8 × 10 −8 M yr −1 ) and Alcalá et al. (2017) (6.7 × 10 −8 M yr −1 ).
-The wind-tracing emission lines ([O i] 6300 Å and [S ii] 6730 Å) are very stable over the time interval of these observations.They do not reflect the short-term changes in the mass accretion rate.-In comparing the models of Romanova et al. (2008), Kulkarni & Romanova (2008, 2009), we establish RU Lup is likely to be in an unstable accretion regime.
-The Li i 6707 Å absorption line is measured to have an average EW Li i = 0.17 ± 0.03 Å and an average r Li i = 3.25 ± 0.20 when the Li depletion in the template star TAP 45 is taken into account.-The temporal variability pattern over the 15 nights for r Li i is strongly correlated to that of L acc (see Fig. 8).-The veiling measured in the photospheric absorption lines Na i 5688 Å (r Na i, avg = 3.27 ± 0.63) and Mn i 6021 Å (r Mn i, avg = 3.29±0.67)demonstrate positive, although different, correlations with the veiling measured in Li i 6707 Å and r Li i, avg = 3.25 ± 0.20.This reveals the relative importance of emission line veiling in each individual line's total veiling.The veiling found in Mn i 6021 Å, for example, increases more rapidly than the other metal lines, given an increase in accretion.
-The correlation we detect between changes in the accretion luminosity and changes in the veiling is strong and provides further evidence that the variability observed is not caused by extinction (the veiling measurements are unaffected by extinction because they are obtained by comparing two normalised spectra).
Because of the strength of the Li i 6707 Å line in CTTSs, its veiling can be an excellent tool in monitoring accretion variability that is independent of extinction and flux calibration.This relationship between variabilities in r Li i and L acc over short timescales can be further explored in other CTTSs to characterise accretion activity regardless of their visual extinction.
A future multi-object study aimed at calibrating this relationship could benefit from the careful selection of template stars that have similar Li abundances as the targets.The correlation between these two quantities found in this paper need not only be applied to the Li i photospheric line; however, its prevalence in TTSs and high S/N make it a better tool than other metal lines that are not always visible.The techniques explained in this paper also provide useful methods to calculate veiling in the absence of a perfect template star, opening the door for more research into the accretion processes in young stars.

Appendix A: Integrated fluxes of accretion lines
Below are the tables containing the integrated fluxes of the accretion-tracing lines.

Fig. 1 .
Fig. 1.ESPaDOnS optical spectrum of RU Lup on the night of 17 June 2011.This spectrum is overlaid with the photometry measured from the ANDICAM observations of the same night.The circles represent the aperture photometry taken in the B, V, R and I bands and the triangles represent the synthetic photometry calculated by convolving the spectrum with the transmission response of ANDICAM in each band.For the B and V bands the triangles and circles are overlapping.
3.4)(Hartigan et al. 1995;Rigliaco et al. 2013;Simon et al. 2016;McGinnis et al. 2018).Taking the flux calibrated [S ii] 6730 Å HVC (the 17 June epoch) as a reference (the [O i] 6300 Å line can also be used, but is less straightforward to analyse because it is in a region of telluric contamination), we re-scaled each epoch by the amount that was necessary to match the flux of its HVC with the flux calibrated one.The lack of temporal variability seen in the resultant[O i] 6300 Å line spectra (Fig.2, left panel) confirms that our re-scaling is a valid calibration, as the [O i] 6300 Å line's HVC

H
i (Hα, Hβ, Hγ and Hδ) and He i (5875 Å, 6678 Å, and 7065 Å).Equivalent widths and integrated fluxes are shown in Tables A.1 and A.2, respectively.Both line profiles and intensities vary with species and (most noticeably in the case of H i) from night to night (see Fig. B.1).In the following we concentrate on the variation of the line luminosity.As seen in the left-hand column of Fig. B.1, the Balmer lines vary significantly over the short timescale of two weeks covered by the ESPaDOnS data.Between the dimmest and brightest nights, there is an increase in integrated flux by a factor of ∼1.7 for Hα, ∼1.9 for Hβ, ∼2.2 for Hγ and ∼2.0 for Hδ.The right-hand column of Fig. B.1 shows the He i emission lines at 5875 Å, 6678 Å and 7065 Å.

Fig. 4 .
Fig. 4. Example of fitting a photospheric (Li i 6707 Å) line in one of the RU Lup spectra (21-06-2011) with the WTTS template from TAP 45 to estimate the veiling present in this region of the spectrum.
. The vast majority of the emission from the [O i] 6300 Å line is blue-shifted and it is completely blue-shifted in the [S ii] 6730 Å.Both lines have a very similar profile with wings reaching very high systematic velocities, namely ∼−300 and ∼−380 km s −1 for the [S ii] 6730 Å and the [O i] 6300 Å lines, respectively.These wing velocity values are slightly higher than those found for RU Lup in the observations of Fang et al. (2018) and Banzatti et al. (2019) in 2008 (∼−270 km s −1 for [O i] 6300 Å) and Whelan et al. (2021) in 2012 (∼−250 and ∼−280 km s −1 , for the [S ii] 6730 Å and the [O i] 6300 Å lines, respectively).All these values for the wings of the HVC in RU Lup are typically higher than those of other TTS with this forbidden emission, often by at least 100 km s −1 , as can be seen in Fig. 17 of Fang et al. (2018) for the case of the [O i] 6300 Å line.These two line profiles show two distinct components: an HVC peaking at ∼−200 km s −1 and a LVC peaking at ∼−17 km s −1 .The LVC of both lines is dominated by a narrow low-velocity component (NLVC), while the HVC contains two components, with one less defined very broad component at ∼−75 and ∼−105 km s −1 for the [S ii] 6730 Å and the [O i] 6300 Å line, respectively (see Fig. D.2).The broad component of the [O i] 6300 Å line has a broader width (roughly 400 km s −1 ) than that of the [S ii] 6730 Å line (roughly 200 km s −1

Fig. 5 .
Fig. 5. L acc variability across the 12 epochs (upper panel) and folded in phase with the rotation of the star (lower panel).The period is 3.71 days (Stempels et al. 2007).The colour of the data points in the lower panel indicates the cycle of rotation of the star.The parameters of the sinusoidal fit are: amplitude = −0.4± 0.1 L and offset = 0.8 ± 0.1 L .The blue dashed line indicates the median value of accretion luminosity of 0.48 L .

Fig. 6 .
Fig. 6.Veiling in the Li i 6707 Å line.r Li i variability across the 12 epochs (upper panel) and folded in phase with the rotation of the star (lower panel).The period is 3.71 days (Stempels et al. 2007).The colour of the data points in the lower panel indicate the cycle of rotation of the star.The parameters of the sinusoidal fit are: amplitude = 1.6±0.1 and offset = 4.7 ± 0.1.The blue dashed line indicates the median value of Li i veiling of 3.23.

Fig. 7 .
Fig. 7. Diagram showing the difference in an observed CTTS (teal) Li i 6707 Å line and that of a WTTS (orange).The deepest absorption feature (dashed light blue) represents a spectrum with no veiling and the same abundance of Li i 6707 Å as the CTTS.The hatched region represents the difference in absorption depth between the CTTS and WTTS.The shaded region represents the difference in absorption depth between the CTTS and the unveiled CTTS spectrum.

Table 1 .
Stellar properties of RU Lup.
Flux calibrated spectra from the 12 epochs taken with ESPaDOnS across the optical regime (4000−7100 Å).Notable emission and absorption lines discussed in this study are identified with their species name.Regions highly contaminated by telluric emission are also indicated.The spectrum in black(17-06-2011)represents the night of observations that we calibrated with simultaneous photometry.The other 11 nights were flux calibrated as described in Sect.2.2.
D.1).The flux of each component varies with a similar trend, indicating that the processes that give A94, page 4 of 17 C. Stock et al.: Accretion variability in RU Lup
Fit representing how the veiling, estimated from the Li i 6707 Å absorption line (including a Li abundance correction -see Sect.4.3), correlates with the L acc logarithmically.The colours of the dots correspond to the same colour scheme used throughout this paper, as shown in Fig.2. by summing ∆ r Li and r Li (obs), which is taking the differing abundance between RU Lup and the template star TAP 45 into account.
Log plots of the veiling measured in the metal lines of Li i 6707 Å, Na i 5688 Å and Mn i 6021 Å in comparison with each other.These plots allow us to visualise how the veiling can vary from one line species to another and that this is indicative of line emission veiling contributing to the total veiling measured in each.

Table A .
1. Integrated fluxes and equivalent widths of the H i lines and their associated errors.Table A.2. Integrated fluxes and equivalent widths of He i 5875 Å and He i 6678 Å and their associated errors.Date F He i 5875 EW He i 5875 F He i 6678 EW He i 6678 yyyy-mm-dd ×10 −13[erg s −1 cm −2 ] [Å] ×10 −13 [erg s −1 cm −2 ] [Å]