© The Authors 2023

1 Introduction

Classical T Tauri stars (CTTSs) are young (~1–10 Myr), low-mass (<2 M_⊙) objects surrounded by a circumstellar disk (Hartmann et al. 2016). Their strong magnetic fields truncate the disk at a few stellar radii (typically 5 R_★, Hartmann et al. 1998). The current paradigm for the interaction between the disk and the star is the magnetospheric accretion (Bouvier et al. 2007), in which the material free-falls onto the star following the magnetic field lines. The rich emission line spectrum typical of CTTSs (Joy 1945; Herbig 1962) can be explained in the framework of this model (Hartmann et al. 1994; Muzerolle et al. 1998), as well as the continuum excess flux, which results from the accretion shock at the stellar surface (Calvet & Gullbring 1998).

Young stellar objects (YSOs) are known to be variable, both photometrically and spectroscopically (Joy 1945; Herbst et al. 1994; Hartmann et al. 2016; Fischer et al. 2023). Many different processes can contribute to this variability, such as variable accretion rate, rotational modulation due to stellar spots, circumstellar extinction, and flares (Cody et al. 2014, 2022). The accretion process is often accompanied by outflows (Hartmann et al. 2016; Bally 2016), either in the form of disk winds or jets (e.g., Romanova et al. 2009; Ferreira 2013, and references therein).

Since the spectro-photometric variability of YSOs shows up in a broad range of wavelengths, from the X-rays to the infrared (Appenzeller & Mundt 1989), it is essential to study these objects using a multiwavelength approach, by means of simultaneous observations in different spectral regions. Many works have shown the capabilities of simultaneous spectro-photometry in determining stellar and accretion parameters, unveiling the inner disk structure of CTTSs and studying accretion variability on different timescales (e.g., Bouvier et al. 2007; Alencar et al. 2018; Zsidi et al. 2022a,b; Fiorellino et al. 2022). This requires coordinated monitoring campaigns of a range of instruments, which are notoriously difficult to achieve. The Hubble UV Legacy Library of Young Stars as Essential Standards, (ULLYSES, Roman-Duval et al. 2020; Espaillat et al. 2022) now offers such a possibility. The aim of this program is to obtain low and medium resolution spectra of YSOs covering the wavelength range from the far-UV (~150 nm) to the infrared (~1000 nm) with the Hubble Space Telescope (HST). Together with the accompanying optical program PENELLOPE (Manara et al. 2021) at ESO Very Large Telescope (VLT), it provides an unprecedented spectroscopic dataset to study the accretion variability of YSOs.

PENELLOPE typically provides three high resolution spectra either with the Echelle SPectrograph for Rocky Exoplanets and Stable Spectroscopic Observations (ESPRESSO, Pepe et al. 2021) or the Ultraviolet and Visual Echelle Spectrograph (UVES, Dekker et al. 2000) and one medium resolution X-shooter (Vernet et al. 2011) spectrum per target taken close in time to the ULLYSES observations. The high resolution spectra are needed to study variability in the emission line profiles, while the X-shooter spectrum is used to constrain the stellar and accretion parameters. As much as possible, the ULLYSES and PENELLOPE observations are scheduled during times where the targets are observed with the Transiting Exoplanet Survey Satellite (TESS, Ricker et al. 2014). TESS produces short-cadence, about one-month-long light curves for many YSOs with a spectral response that covers the red/infrared wavelength range (~0.6–1.1 µm). Since most of the excess flux due to accretion is emitted in the near-UV to optical spectral region, at wavelengths shorter than the TESS filter band (Calvet & Gullbring 1998), it is also important to obtain simultaneous multiband photometry (Robinson et al. 2022).

In this work we take advantage of the wealth of data provided by the aforementioned programs to study in detail a single target, chosen because of the simultaneous coverage with all the relevant programs. The target of this study is HM Lup (Sz 72), a CTTS with spectral type M2 and a mass of 0.37 M_⊙ (Alcalá et al. 2014; Manara et al. 2023) located in the Lupus cloud at a distance of 156 pc (Gaia Collaboration 2021). It hosts a disk with a dust mass of 4.09 M_⊕ and a radius of 10.97 au (Manara et al. 2023), inclined by 53 ± 19 deg relative to the line of sight (Ansdell et al. 2016).

The goal of this paper is twofold. First, we aim to study the temperature stratification of the magnetosphere of HM Lup using a range of emission lines. Secondly, we aim to investigate the spectrophotometric variability of the system.

The paper is structured as follows. In Sect. 2 we describe the observations. We report the basic properties of the star in Sect. 3. We present the results based on the optical spectrum of the system in Sect. 4 and the results on the spectrophotometric variability in Sect. 5, while we discuss them in Sect. 6.

2 Observations

For our study, we selected HM Lup as one of the PENELLOPE targets having spectra taken simultaneously with a TESS observation, in this case in sector 38 (2021). HMLup was also observed with TESS in Sector 12 (2019).

2.1 Simultaneous data

We downloaded the TESS light curves of sectors 12 and 38 from the MAST archive¹. The full frame images (FFI) were reduced by the TESS Science Processing Operations Center (SPOC), producing a 10 minute cadence light curve. We define the beginning of TESS Sector 38 observations, MJD₀ = 59 333.363, as the reference time for all the simultaneous spectroscopic and photometric data. Figure 1 shows the TESS light curve in which the pre-search data conditioning simple aperture photometry (PDCSAP) flux (F) was converted into TESS magnitudes (T) using the relation T = −2.5 · log₁₀ F + ZP where ZP = 20.44 is the TESS Zero Point magnitude (Fausnaugh et al. 2021; Vanderspek et al. 2018).

We downloaded BVR_c I_c photometry available for HM Lup from the American Association of Variable Star Observers (AAVSO) International Database². Most of the data were secured from the beginning of TESS observations and span ~140 days. Data from different observers were sometimes obtained during the same day at small temporal distance (Δt ~ 15 min). To have a better view of the variability, we binned the photometry to 0.15 day.

HM Lup was selected as an ULLYSES target and observed on 4 May 2021 with the Space Telescope Imaging Spectrograph (STIS) camera between 1650 and 10200 Å with a resolving power of R ~ 1500. Contemporaneous to the HST data, medium-high resolution optical spectroscopy was obtained in the framework of PENELLOPE (Manara et al. 2021). The journal of the spectroscopic observations is reported in Table 1, while Fig. 1 shows the temporal position of the spectra relative to the 2021 photometric data. High resolution spectra were obtained with ESPRESSO in Pr. Id. 106.20Z8.003 (PI Manara). Since in the first observation the conditions did not fulfill the requirements, namely the seeing was >1.55″, the observation was repeated the day after. However, with its signal-to-noise ratio (S/N) of 12, the first spectrum can still be used. The ESPRESSO spectra cover a wavelength range between 3800 and 7880 Å with a resolution of 140 000. These spectra were flux-calibrated as explained in Appendix A. The X-shooter spectrum taken from Pr. Id. 106.20Z8.004 (PI Manara) is divided into three arms, UVB (3000–5600 Å, R ~ 5400), VIS (5600–10200 Å, R ~ 18 400), and NIR (10 200–24 800 Å, R ~ 11600). In order to achieve this spectral resolution, the observation was carried out with slit widths of 1″, 0.4″, 0.4″ in the UVB, VIS and NIR arms, respectively. The absolute flux calibration correcting for slit losses was achieved using a short exposure taken with a 5″ slit, as explained by Manara et al. (2021). Telluric correction was performed on the ESPRESSO and X-shooter spectra using the molecfit tool (Smette et al. 2015). All the optical spectra were reduced by the PENELLOPE team and made publicly available on Zenodo³.

2.2 Nonsimultaneous data

In addition to the simultaneous spectro-photometric data, we have also included observations of HMLup taken at other epochs. In particular, an X-shooter spectrum obtained on 18 April 2012, which was already presented by Alcalá et al. (2014, 2017) and Frasca et al. (2017) taken from Pr. Id. 089.C-0143(A) (PI Alcalá), and the TESS light curve from sector 12, obtained on 25 May 2019, UT 01:11:43 (MJD₁ = 58 628.050) and shown in Fig. 2.

Fig. 1 Light curves and timing of simultaneous spectroscopy for HM Lup. The topmost panel shows the 2021 10 minute cadence TESS light curve (Sector 38). The inset displays the Lomb–Scargle Periodogram for the TESS light curve. Overlaid in purple is a boxcar smoothed version of the light curve produced by phase-folding the data with the period of the highest peak in the Lomb–Scargle periodogram, 4.79 day. The phase ϕ was computed with this period and ϕ = 0 at the maximum of the TESS light curve, MJD 59 336.97. The vertical dashed lines mark the epochs of the simultaneous spectroscopic observations. Here ES = ESPRESSO and XS = X-shooter. The other panels show the part of the AAVSO BVRcIc photometry that is simultaneous with TESS, with a linear interpolation as a guideline. The large open squares mark the synthetic photometry obtained from the X-shooter spectrum for the four filters. We defined MJD0 = 59 333.363 as the beginning of TESS Sector 38 observation.

3 Stellar parameters

The properties of the system were determined by Alcalá et al. (2017) and Frasca et al. (2017) by fitting the 2012 X-shooter spectrum. While the first paper focused on the spectral type (SpT), the stellar luminosity L_★ and the accretion properties, the second work used the ROTFIT routine (Frasca et al. 2015) to derive the atmospheric parameters of the accreting star, namely the effective temperature T_eff, the gravity log g, the projected rotational velocity v sin i, and the systemic radial velocity RV. Manara et al. (2023) updated the values from Alcalá et al. (2017) using the Gaia DR3 distance and assuming the Herczeg & Hillenbrand (2014) relation between SpT and T_eff and nonmagnetic evolutionary tracks by Feiden (2016). We applied the same procedure to the 2021 X-shooter and ESPRESSO spectra. Except for the RV, the parameters are in agreement with the previous result. The best fit RV value for the 2012 spectrum was 6.9 ± 2.4 km s⁻¹ while we obtained RV = −2.7 ± 1.9 km s⁻¹ from the 2021 X-shooter spectrum, suggesting a possible companion. We measured the RVs in the ESPRESSO spectra, obtaining −3.1 ± 0.7 km s⁻¹, −3.6 ± 0.7 km s⁻¹, −2.5 ± 0.3 km s⁻¹, and −2.7 ± 0.5 km s⁻¹ for the epochs 1, 2, 3, and 4 respectively. Therefore, no RV variation (within 1σ) emerges from the spectra during the 2021 PENELLOPE campaign. We did not find other signs of binarity in the spectra, such as asymmetric absorption lines. In summary from the RV measurements, we cannot exclude that HMLup is actually a long-period (P ≳ 10 d) single-lined spectroscopic binary.

Table 2 reports the stellar parameters of the system. Regarding υ sin i and RV, we adopted the values obtained from an average of the best fit values for the ESPRESSO spectra, given their higher resolution. We computed the stellar radius R_★ inverting the relation $L_{★}=4\pi R_{★}^2\sigma T_{\text{eff}}^4$.

In Fig. 3, we show HMLup in the Ṁ_acc-M_★ diagram together with the whole X-shooter Lupus sample. Both the 2012 measurement of the accretion rate and our anticipated new values (from Sect. 5.2) are included, showing that HM Lup is a strongly accreting CTTS.

Fig. 2 TESS 2019 two-minute cadence light curve (Sector 12) of HM Lup. The inset shows the Lomb-Scargle periodogram. Overlaid in purple is a boxcar smoothed version of the light curve produced by phase-folding the data with the period of the highest peak in the Lomb-Scargle periodogram, 5.41 d. We defined MJD1 = 58 628.050 as the beginning of the TESS Sector 12 observation.

4 The optical spectrum of the system

The optical spectrum of HM Lup is rich in emission lines, the strongest being the Balmer series and the Ca II H & K lines. In addition, permitted emission lines from many other species can be identified, such as He I and singly and doubly ionized metallic elements, for example, Na I, Ca I, Ti I and Ti II, Fe I, and Fe II. We detected more than 150 emission lines from the iron peak elements. This feature is reminiscent of the outburst spectra of EXors (Sicilia-Aguilar et al. 2012, 2017) or the spectra of other strong accretors, such as DR Tau (Beristain et al. 1998). We found outflow signatures in the spectra, such as the [O I] 6300 line in emission and blueshifted absorption in the He I 10830 line. However, these features will not be analyzed in this work. A selection of the observed permitted spectral lines is shown in Fig. 4 for the four ESPRESSO spectra.

The different excitation potentials of the observed transitions and the differences in the emission line profiles indicate the presence of a thermally stratified environment, with multiple regions that contribute to the observed spectrum. In this section, we focus on a single spectrum, the epoch 2 of ESPRESSO (cyan in Fig. 4), to illustrate the stratification of the accretion flow. We chose this spectrum because of its high S/N and the strength of the emission lines from low excitation transitions of neutral and singly ionized species. We discuss the line variability in Sect. 5.3. The atomic parameters for the analyzed emission lines were taken from the NIST Atomic Spectra Database⁴.

Fig. 3 Ṁacc-M★ diagram for X-shooter targets in Lupus, using the values from Manara et al. (2023). The black triangles indicate upper limits on the measured accretion rate. The blue, green and red diamonds mark the accretion rate of HM Lup in 2012 and 2021 (Sect. 5.2). The dashed gray line is a linear best fit to the data.

4.1 Balmerseries, He I, Ca II K

According to the morphological line profile classification by Reipurth et al. (1996), the Balmer series, except for Ha, and the Can K line have a Type IIB profile, that is, a double-peaked emission profile in which the secondary peak exceeds half the strength of the primary peak. The Ha line has instead a flat-topped emission profile. Balmer and Ca II K emission lines have symmetrical wings up to ±400 km s⁻¹. The He I lines have less broad wings, up to ~±200 km s⁻¹, and consist of a narrow component (NC) and a broad component (BC). In the rest of the paper, we focus our analysis on the BC emission.

In the magnetospheric accretion scenario, the BC is expected to originate in the infalling material (Beristain et al. 1998, 2001; Hartmann et al. 2016). The free fall velocity onto a star of mass M_★ starting from rest at the disk truncation radius R_T is $${\upsilon }_{\text{ff}}={\left(2G{M}_{★}\right)}^{1/2}{\left(\frac{1}{r}-\frac{1}{R_{\text{T}}}\right)}^{1/2}.$$(1)

Assuming R_T = 5 R_★ (Gullbring et al. 1998) and the stellar parameters of HM Lup from Table 2, we obtain a free fall velocity of ~175 km s⁻¹ at r = 2 R_★, consistent with the observed line wings of the Ca II K and He I 5876 lines. The wings of the lower Balmer lines, especially Ha and Hβ, exceed these values, possibly due to Stark broadening (Muzerolle et al. 2001; Wilson et al. 2022). Conversely, the higher lines of the series are less affected by Stark broadening and their wings agree well with the wings of He I 5876, as shown in the left panel of Fig. 5 for the H9 line. This suggests a common origin for the BC of the H I and the He I lines. Given the high excitation potentials of the He I lines, which have upper states with energies E_j ≳ 20 eV, a source of ionizing radiation is required (Beristain et al. 2001). In the magnetospheric accretion flow, these conditions are met in the pre-shock region, that is heated by soft X-rays emitted by the accretion shock (Hartmann et al. 2016). Singly ionized calcium, with its ionization potential of 11.87 eV, is almost completely ionized in such conditions. Azevedo et al. (2006) studied the formation of the Ca II infrared triplet lines and showed that the gas departs from LTE conditions. Close to the star, because of the higher dilution factor for the accretion shock radiation, calcium is mostly in the doubly ionized stage (Ca III). However, we observe strong Ca II emission at all ESPRESSO epochs (Fig. 4). This suggests that the Ca II lines are emitted from an outer part of the magnetospheric flow, where calcium is predominantly in the singly ionized stage. The right panel of Fig. 5 shows how the wings of the Ca II K line are less pronounced than the wings of, for example, Hδ, in agreement with this hypothesis.

Fig. 4 Selection of continuum subtracted emission lines in the ESPRESSO spectra of HM Lup. The colors for the ESPRESSO spectra are the same as in Fig. 1. The lines from the 2012 X-shooter spectrum are shown in green for comparison. The vertical dashed lines mark the radial velocity of the system, while the horizontal dashed lines highlight the zero-flux level for each set of lines. The Fe I and Ca I lines were smoothed with a 7 points boxcar filter.

Fig. 5 Continuum subtracted profiles of selected magnetospheric lines in the epoch 2 ESPRESSO spectrum. In the left panel, H9 (3835 Å) and He I 5876. In the right panel, Hδ and Ca II K. The vertical dashed line marks the radial velocity of the system. The H9 line was smoothed with a boxcar filter and multiplied by 0.5 to match the He I 5876 wings. The Ca II K line was rescaled to Hδ in a similar way.

Fig. 6 Reconstruction of the Fe I absorption. All profiles are from the epoch 2 ESPRESSO spectrum. All lines are plotted relative to the rest wavelength of Fe I 3969. The top panel shows the continuum subtracted profile of the Ca II H + Hϵ blend (black) compared to the Ca II K (green) and H8 (blue) lines. The red line is the model for the unabsorbed blend (Eq. (2)). The vertical dashed lines indicate the velocity displacement of Ca II H, Fe I 3969 and Hϵ relative to the rest wavelength of Fe I 3969. The middle panel shows the Fe I 3969 absorption profile, reconstructed as explained in Sect. 4.2, and its best fit. The bottom panel shows the Fe I 4064 and 4132 emission lines.

4.2 Fe I fluorescence

The third plot in the second row of Fig. 4 compares the Fe I 4064, 4132, and 4144 lines. The first two are a doublet, having a common upper level. Although these three lines are from the same multiplet, the integrated flux in the Fe I 4132 line, the weaker line of the doublet, is more than twice the integrated flux of the Fe I 4144 line, as shown in Table 3. This behavior was first recognized by Herbig (1945), who proposed that the unusual strength of this doublet is due to a fluorescence mechanism. The Fe I line at 3969. 26 Å, which shares the upper level with the doublet, is nearly coincident in wavelength with Ca II H (3968.47 Å) and Hϵ (3970.08 Å). Absorption in this line can produce an overpopulation of the y³F⁰ ( J = 3) level, from which the fluorescent emission lines originate. The top panel of Fig. 6 shows that the Ca II H + Hϵ blend is strongly absorbed in the region of Fe I 3969.

We reconstructed the Fe I 3969 absorption profile by creating a model for the unabsorbed emission in the Ca II H + Hϵ blend. To this end, we used the Ca II K and H8 lines as template for Ca II H and Hrespectively. We subtracted the, respectively. We subtracted the continuum in these lines and shifted them in radial velocity to the position of Ca II H and Hϵ relative to the Fe I 3969 rest wavelength. The non-absorbed model profile was computed as a linear combination, $$F_{\text{model}}=0.8F_{\text{CaIIK}}+1.3F_{\text{H8}}.$$(2)

The coefficients were determined by matching the wings of the blended profile (|υ_rad| > 150 km s⁻¹), where absorption from Fe I is not expected, to the wings of the two template lines. The absorption profile was then obtained by dividing the observed Ca II H + Hϵ blend by its model. The result is shown in the middle panel of Fig. 6. We fit this profile with a Gaussian $$F_{\upsilon }=1+C\exp \left[\frac{{\left(\upsilon -{\upsilon }_0\right)}^2}{2{\sigma }^2}\right].$$(3)

The best fit parameters are C = −0.777 ± 0.002, υ₀ = 41.4 ± 0.2 km s⁻¹, σ = 64.8 ± 0.3 km s⁻¹. The line width is in agreement with those observed in the fluorescent lines (bottom panel of Fig. 6).

The Fe I 3969 profile provides constraints on the structure of the circumstellar envelope (Willson 1975). Since the Ca II and H I lines are formed in the pre-shock region (Sect. 4.1), a foreground structure that is opaque in Fe I must be present to produce absorption. From the Gaussian best fit of the profile and assuming a purely extinguishing medium for which F(υ) = exp(−τ_υ), we obtain a minimum optical depth at line center τ₀ ~ 1.5. Therefore, this region is optically thick in Fe I. This is confirmed by the integrated flux ratio of the fluorescent lines, r_FeI, that is equal to 2.6 at epoch 2, much lower than the ratio of the A_ji values, 5.6. The line center is redshifted by ~40 km s⁻¹, indicating that the absorption is produced in infalling material.

Overall, the reconstructed Fe I 3969 absorption profile is similar to the core of the Balmer lines (e.g., Hβ in Fig. 4) in the epoch 1 and 2 ESPRESSO spectra. To compare the two features quantitatively, we fit the Hβ line in the epoch 2 ESPRESSO spectrum with a triple Gaussian model. We forced one Gaussian to fit the broad emission and the other two to reproduce the asymmetric depression. Figure 7 displays the results of the fit. Compared to Fe I 3969, the absorption in Hβ has a blueshifted centroid and extends to lower positive velocities. The difference probably stems from the fact that Hβ is self-absorbed, while Fe I 3969 absorbs against Ca II H and Hϵ. This means that, in the case of the Balmer lines, absorption can take place only if emission at a given radial velocity is intercepted by hydrogen gas moving at the same υ_rad⁵ along the line of sight. On the other hand, the fluorescence phenomenon couples regions having radial velocity differences equal to the shift between the rest wavelengths of Ca II H and Hϵ relative to Fe I 3969, that is, ~−59 km s⁻¹ and ~+62 km s⁻¹ respectively. This indicates the presence of velocity gradients in the accretion flow, in agreement with the magnetospheric accretion scenario. An outer zone of the magnetosphere sees Ca II H emission that is redshifted, leading to the tuned absorption in Fe I 3969 (Gahm 2001). For the same reason, absorption in Fe I can take place due to H I emission that is seen as blueshifted by Fe I. This result also demonstrates that the Ca II emission at moderate blueshifted velocities (−100 km s⁻¹ ≲ v_rad ≲ 0 km s⁻¹) is produced in the infalling gas.

Fig. 7 Comparison between the Hβ and Fe I 3969 absorption components. Top panel: triple Gaussian fit of the continuum normalized Hβ line in the epoch 2 ESPRESSO spectrum. The dashed lines show the three Gaussian components. Bottom panel: COMPARISON between the continuum normalized absorption profiles in Hβ and Fe I 3969. The red line is the Gaussian best fit of the Fe I 3969 from Fig. 6. The blue line is the sum of the two Gaussian absorption components in Hβ, rescaled relative to the continuum to allow the comparison with the Fe I 3969 absorption profile.

4.3 Other metallic lines

Among the numerous metallic lines observed in the spectrum of HM Lup, we selected a set of emission lines that belong to different species and have different ionization potentials and excitation energies, so that they potentially probe different conditions in the accretion structure. Table 3 reports the atomic parameters of the selected transitions. The lines are shown in the second row of Fig. 4. All the lines display a NC+BC structure. The BC has a full width at half maximum (FWHM) of ~ 115 km s⁻¹ and an emission profile that is skewed to the red. The last two columns show the integrated flux for each line in the epoch 2 and 3 ESPRESSO spectra, chosen as representative of the line behavior during and just after the accretion burst. Before integrating, we dereddened the spectra using A_V = 1 mag (Sect. 5.2) and the Cardelli et al. (1989) extinction law with R_V = 3.1. At epoch 2, the emission is strongest in Fe II, followed by the Na I D lines. Fe I lines are weaker than Fe II lines, typically by a factor of ~5–7, except for the Fe I 4064 and 4132 lines (Sect. 4.2). This is different than EX Lup, for which Sicilia-Aguilar et al. (2012) found a ratio of two for the Fe II vs. Fe I during the outburst.

Figure 8 shows the Ca II infrared triplet (IRT) lines from the 2021 X-shooter spectrum, that have profiles similar to those of the metallic lines. The components of the Ca II IRT have relative intensities that differ from the 1:9:5 relation that is expected from the ratios of their gf-values⁶, indicating that these lines are formed in an optically thick environment. This sets a lower limit to the electron density of 10¹¹ cm⁻³ (Hamann & Persson 1992).

Additional information on the optical depth of the medium can be obtained from the flux ratio of lines that share the upper level. We chose the Fe II 4549 and 4352 lines (E_j = 5.55 eV) for this exercise, since they are strong, not blended and near in wavelength. The integrated flux ratio for the epoch 2 ESPRESSO spectrum is r_FeII ≈ 1.16. This value is lower than the ratio of their A_ji values, 2.04, confirming that the region where the metallic lines are formed is optically thick. This complicates the derivation of the properties of the medium in terms of line ratios, since the line emission depends on the escape probability (Sobolev 1960; Kogure & Leung 2007). The comparison between the line ratio and the ratio of the A_ji values for the selected Fe I and Fe II doublets indicate that this region is more optically thick in Fe I than in Fe II. Moreover, the observation of the Fe I fluorescence (Sect. 4.2) suggests that non-LTE (NLTE) effects can be important in the region where the metallic lines are formed. Therefore, the Fe II to Fe I ratio cannot be taken as an unambiguous diagnostic of the actual conditions of this region.

Some constraints on n_e and T can still be placed from line ratios under the assumption of LTE conditions (Appendix B). In the conditions for the coexistence of Mg I, Fe I and Fe II, neutral sodium and calcium are almost completely ionized. In such an environment, the Na I D and Ca I 4227 lines are therefore less optically thick than the Fe I lines. Hence, the Ca II to Ca I ratio is a better temperature indicator than the Fe II to Fe I ratio. The Ca I 4227 line can be compared to one of the IRT lines, for instance the Ca II 8542 line. Another temperature diagnostic for the emitting region is the Mg II to Mg I ratio. The non-detection of Mg II 4481, which has an upper state with E_j = 11.63 eV, provides an upper limit to the temperature. We computed this upper limit by estimating the noise level in the continuum subtracted spectrum at the nominal position of the Mg II line. The atomic parameters of the Mg II 4481 line and the Ca II 8542 line are reported in the last two rows of Table 3. Figure 9 shows the result of the LTE analysis described in Appendix B for Ca and Mg. For 11 ≲ log n_e ≲ 14, the Ca II to Ca I ratio is reproduced with temperatures between ~4000 and 7000 K, in agreement with the density-dependent upper limit on T placed by the nondetection of Mg II 4481 (Fig. 9).

Fig. 8 Continuum subtracted profiles of the Ca II IRT lines in the X-shooter spectrum at MJD – MJD0 = 3.80 days, almost simultaneous to the ESPRESSO epoch 2 spectrum.

Fig. 9 Result of the Saha-Boltzmann analysis for Ca and Mg. Observed line ratios were computed from the epoch 2 ESPRESSO spectrum. Left panel: likelihood for the observed Ca II 8542 vs. Ca I 4227 ratio as a function of log ne and T. Right panel: theoretical Mg II 4481 vs. Mg I 5184 ratios as a function of T for log ne [cm−3] = 11,12,13,14. The red area marks the upper limit on the observed ratio.

5 Spectrophotometric variability

5.1 Photometry

Figure 1 summarizes all available photometric data in 2021, while Fig. 2 shows the TESS Sector 12 light curve from May 2019. We lack AAVSO data in the highest flux peak of the 2021 TESS light curve, that is, in the range 3 days ≲ MJD – MJD₀ ≲ 5 d, which is coincident with the X-shooter spectrum. To complement the multiband light curve at that epoch, we computed synthetic BVR_cI_c photometry from the X-shooter spectrum using the filter transmission curves downloaded from the SVO Filter Profile Service⁷ (Rodrigo et al. 2012; Rodrigo & Solano 2020).

There are differences in the light curve morphology between the two TESS epochs. The 2019 TESS data show no clear periodicity pattern and a tendency for showing dips, while quasi-periodic spikes are observed in the 2021 TESS photometry. To quantify the different morphology, we computed the variability metrics introduced by Cody et al. (2014) for both light curves. The parameter Q describes the degree of periodicity and the parameter M the (a) symmetry around the median value of a given light curve. Q ranges between 0 (= highly periodic) and 1 (= aperiodic), while strongly negative (positive) values of M indicate the tendency of showing bursts (dips) in the light curve. We obtained Q = 0.39 and M = −0.66 for the 2021 TESS light curve and Q = 0.57 and M = −0.18 for the 2019 TESS light curve, suggesting a higher degree of periodicity and a more pronounced bursting behavior in 2021.

We performed a Lomb–Scargle (Lomb 1976; Scargle 1982) analysis on the TESS light curves. The Lomb–Scargle Periodograms (LSPs) for the TESS 2021 and 2019 light curves are shown in the insets of Figs. 1 and 2, respectively. We detected four peaks in both LSPs, all below a False Alarm Probability (FAP) of 10⁻⁵. We chose the periods associated to the highest peak of each periodogram, 4.79 d and 5.41 d for the light curves from 2021 and 2019 respectively, as representative of the periodicity of the light curves. This is shown in Figs. 1 and 2 by visually comparing the light curves to an average profile produced by phase folding the data with those periods and smoothing the phase folded light curve with a boxcar filter with a width of 0.25 in phase, as described by Cody et al. (2014). We estimated the uncertainties on the periods as the standard deviation of a Gaussian function fitted to the peaks in the LSPs. The results are P = 4.79 ± 0.28 days for the year 2021 and P = 5.41 ± 0.35 days for 2019. Although the central value is different, the two periods agree within the uncertainties. The LSP of both TESS light curves present additional, weaker peaks with the same structure. For the case of the 2021 light curve they are at 3. 76, 6. 20 and 9.58 days. The 2019 light curve shows peaks at similar values. The latter is the first harmonic of the 4.79 days period. To understand whether the other two periods are significant, we produced an average profile analogous to the one shown in purple in Fig. 1. The comparison to the observed light curve showed poor match. This suggests that these peaks are aliases.

We studied the periodicity of the AAVSO light curves by means of continuous wavelet analysis, which is a useful tool to determine the frequency content of a signal as a function of time. We chose as wavelet template the Morse wavelet (Lilly & Olhede 2012) with a symmetry parameter γ = 3 and a time-bandwidth product 𝒫 = 90, that produces a time resolution of ~30 days. Since the wavelet analysis works only on evenly spaced data, we linearly interpolated the AAVSO data to 0.25 d centers to remove the gaps in the light curves. Higher cadences did not change the results of the analysis. The continuous wavelet transform (CWT) of the B band photometry is shown in Fig. 10. The light curve and its CWT can be divided into three time segments with different power spectra, as shown in the lower panel of Fig. 10. In the first segment (B1), which is simultaneous to the 2021 TESS light curve, the detected period P_B1 is compatible with the period of the TESS light curve, 4.79 days. In the other two segments the main periods detected are P_B2 = 6.24 days and P_B3 = 8.84 days. In the LSP of the whole light curve, in which the time dependence of the power spectrum is lost, these multiple contributions are averaged out. We considered only P_B3 as statistically relevant, since it has FAP = 2. 5% in the LSP.

The part of the AAVSO data that is simultaneous with TESS data follow the shape of the TESS light curve, as shown in Fig. 1. The relative magnitude variations vary across the bands, being 0.3 mag in R_c and I_c, 0.5 mag in V and 0.7 mag in B. The approximately enhanced variability amplitude in the blue part of the spectrum is compatible with changes in the accretion rate, since the excess flux is expected to peak at shorter wavelengths (Calvet & Gullbring 1998).

5.2 Accretion rate variability

A comparison between the two X-shooter spectra and the STIS spectrum in the 3200–10 000 Å range is shown in Fig. 11. Marked differences in the Balmer continuum and the Balmer jump are evident between the three spectra. In the 2021 X-shooter spectrum the Balmer continuum was ~3 times higher than in 2012. The ratio is reduced to a factor ~1.5 in the Paschen continuum. The flux in the STIS spectrum lies approximately between these two.

We derived the accretion parameters by fitting each spectrum with a combination of a photospheric (nonaccreting) template and a model for continuum emission from a slab of hydrogen representing the accretion shock. This method yields simultaneously with the UV accretion luminosity (L_acc,UV), the spectral type, and extinction (A_v) of the star. For more details on the procedure we refer to Manara et al. (2013). For the 2012 spectrum we used published values obtained with the same method. That spectrum was first analyzed by Alcalá et al. (2014), but in this work we use the updated values from Manara et al. (2023). The resulting values of A_V, L_acc,UV, and Ṁ_acc are listed in Table 4. Typical uncertainties on these parameters are 0.1 mag, 0.25, and 0.45 dex respectively (Manara et al. 2013, 2021; Alcalá et al. 2014). The best fit of the 2021 X-shooter spectrum returned a SpT M2, in agreement with the fit obtained by Alcalá et al. (2014) for the 2012 spectrum, and A_v = 1 mag. The STIS spectrum was fitted in the X-shooter wavelength range. Its low resolution (R ~ 1500) prevents constraining the SpT. Therefore, we fixed it to M2 in the fitting routine and obtained A_v = 1.5 mag and Ṁ_acc = 1.53 × 10⁻⁸ M_⊙ yr⁻¹. Since this spectrum is only ~1.3 d separated in time from the second X-shooter observation, variations in A_V on the order of ΔA_V = 0.5 mag seem unlikely, given that variations in the TESS bandpass are only about 0.2 mag, as shown in Fig. 1. Hence, we fixed A_V = 1 mag in the STIS best fit, obtaining the accretion rate reported in Table 4.

We obtained an independent measure of the accretion rate from the line luminosities, L_line, using their empirical relation with L_acc,UV, calibrated by Alcalá et al. (2017). For the X-shooter spectra, we integrated the available Balmer and Paschen lines, Brγ, Ca II K and infrared triplet (IRT) and He I 5876 and 6678 lines. Since the ESPRESSO spectra cover a narrower wavelength range, we integrated only the Balmer, He I, and Ca II K lines. For the STIS spectrum, we excluded the He I lines because of their low S/N. Before integrating, we dereddened the spectra using the A_V value obtained from the slab model fit of the X-shooter spectrum. Figure 12 shows the comparison between the accretion luminosity derived from the lines, L_acc,lines, and L_acc,UV for the 2021 X-shooter spectrum. Most of the lines are in agreement with the accretion luminosity obtained from the slab model, but the Ca II lines overpredict the accretion luminosity. For this reason, we excluded these lines from the calculation of L_acc,lines. We further discuss this issue in Sect. 6.1. A global value for L_acc,lines was then derived from a weighted average of the values obtained for the individual lines. This mean value was used to compute the accretion rate from the formula Ṁ_acc = 1.25L_accR_⋆/(GM_⋆) (Hartmann et al. 1998). The results are reported in Table 4.

Additional information on the variations in the accretion rate of CTTSs can be derived by measuring the veiling, that is, the excess emission due to the accretion process that makes photo-spheric absorption lines appear less deep (Hartigan et al. 1989). This parameter is defined as the ratio of the excess flux relative to the photospheric flux, r_λ = F_acc(λ)/F_phot(λ). We computed r_λ at 6000 Å for the available spectra using the ROTFIT routine, as described by Frasca et al. (2015) and Manara et al. (2021). The result is shown in the last column of Table 4.

The comparison between L_acc in 2012 and 2021 indicates that HMLup was observed in a more active state in 2021, with an accretion rate ~3–4 times higher. The computed values of Ṁ_acc and r₆₀₀₀ support the hypothesis that the photometric behavior of the system in the main burst of the 2021 TESS light curve was caused by an increase in the accretion rate by a factor of ~1.5.

Fig. 10 Continuous wavelet analysis of the AAVSO B band photometry. The top panel shows the CWT. The inset on the left displays the LSP of the whole dataset, together with a temporal average of the CWT. The orange horizontal dot-dashed lines mark the periods obtained from the LSP analysis, i.e., PB2 = 6.24 days and PB3 = 8.84 day, while the red one indicates the 2021 TESS period, PTESS = 4.79 days. The bottom panel displays the AAVSO B band light curve and its linear interpolation. The insets show the LSP for three different time segments. We excluded the shaded areas from the calculation of the LSPs to highlight the different timescales detected.

Fig. 11 VLT/X-shooter and HST/STIS spectra of HM Lup in the STIS wavelength range. The spectra were smoothed to the HST resolution with a Gaussian filter for clarity.

Fig. 12 Comparison between the accretion luminosity derived from emission lines and that derived from the UV excess in the 2021 X-shooter spectrum. Lacc,UV and its uncertainty are displayed with the red dashed line and the red shaded area. The blue dashed line and the blue shaded area are the value of Lacc,lines and its uncertainty, obtained from a weighted average of the line values. The lines marked in black were excluded from the calculation of the average.

5.3 Line variability

The variability of the Hα, Hβ, Ca II K, and He I 5876 and 6678 lines between the four ESPRESSO spectra is shown in the first row of Fig. 4. The shape of the He I lines is invariant, but the line flux changes. It is highest in the second epoch and lowest in the third and fourth, in agreement with the photometric behavior during the accretion burst. The Hβ and Ca II K lines evolve from a double peaked structure in the first two epochs to a profile with a single red-shifted peak in the third and fourth. Hα always remains flat-topped, suggesting that the double peaks are caused by an optical depth effect. If the temperature and density of the magnetosphere increase, as is likely to happen during the accretion burst, Hβ and the higher lines of the series can develop a self-absorption component (Muzerolle et al. 2001). The same happens for Ca II K which, being a resonance line, is more prone to absorption. On the other hand the Hα line, given its higher opacity, can be strongly optically thick and thermalized in these physical conditions, resulting in the disappearance of the self-absorption component (Hartmann et al. 2016).

The flux emitted in the metallic lines follows the photometric behavior of the system. It is strongest in the first two epochs, while it returns to its quiescent level after the accretion burst, as suggested by the comparison of line fluxes between the third ESPRESSO epoch and 2012 X-shooter spectra in Fig. 4. In the third and fourth epoch, the BC is absent in Fe I and reduced in flux by a factor of ~2–3 in Fe II, Mg I, and Na I. The fact that the emission is stronger in both the Fe I and Fe II lines during the high accretion state suggests that the outburst spectrum is mainly the result of a density enhancement, that is, an overall increase in the number of emitters, rather than a temperature increase. Table 3 shows how the Fe II to Fe I ratio at epoch 2 is lower than at epoch 3. This is compatible with an increase in the electron density during the outburst that favors the LTE population of Fe I relative to Fe II (Appendix B).

To quantify the line variability we fit the Fe II 4924 line in each ESPRESSO epoch with a triple Gaussian model (Appendix C) as shown in Fig. 13. We chose this transition as template for the line variability because it has the highest S/N among the metallic lines. The best fit parameters are reported in Table 5. The line variations can be explained in terms of a broad emission component plus a variable redshifted absorption that reproduces the observed red skewness. The line is symmetric at epoch 1. Then, a redshifted absorption at ~+40 km s⁻¹ appears at epoch 2, it becomes more pronounced, broader and at higher positive velocities at epoch 3, and eventually disappears at epoch 4. The evolution of the redshifted absorption is indicative of optical depth changes in the infalling material, which is likely associated to the rotation of a non-axisymmetric accretion structure, similar to what observed in EX Lup (Sicilia-Aguilar et al. 2012) or CVSO109 (Campbell-White et al. 2021).

5.4 Epoch 4 ESPRESSO spectrum

Although the veiling is constant within the uncertainties between the epoch 3 and 4 ESPRESSO spectra, the accretion rate predicted from the lines is higher in epoch 4 than in epoch 3. This happens because the redshifted emission peak increases relative to the continuum in the higher lines of the Balmer series, as shown in Fig. 4. On the other hand, the He I BC does not change between the two epochs. This further supports the hypothesis that the lines are formed in a stratified environment, both in temperature and in density. The Balmer emission peaks are not linked to the region that produces the veiling, namely, the continuum emission, unlike the wings of the Balmer and He I lines.

Although we lack photometric information for the ESPRESSO epoch 4 spectrum, the phase-folded TESS light curve suggests that during this epoch the system is roughly in the same configuration as in the first two epochs. The phases are indeed ϕ₄ = 1.97, ϕ₁ = −0.11 and ϕ₂ = 0.08. However, the line profiles are remarkably different. The metallic lines are weaker and the double peaked profile is absent in the Balmer lines (Fig. 4). According to the metallic lines and the veiling HMLup has returned to a lower state of accretion in the epoch 4 ESPRESSO spectrum. The behavior of the metallic lines between the epoch 3 and 4 spectra (Fig. 13 and Sect. 5.3) suggests that the enhanced redshifted emission in the Balmer series at epoch 4 is caused by an increase in the emissivity in the outer region of the accretion flow, where these lines are produced. An alternative explanation could be a decrease in the opacity of this region, that has caused redshifted absorption at epoch 3.

Fig. 13 Triple Gaussian fit of the Fe II 4924 line profiles in the four ESPRESSO spectra. See Appendix C for the fit function.

6 Discussion

6.1 Temperature stratification of the accretion flow

The many emission lines in the optical spectrum of HMLup allowed us to probe the temperature stratification of the accretion flow. Figure 14 shows a schematic representation of the structure that emerges from the line analysis. The He I and H I emission wings are produced in the pre-shock region, that we call the inner magnetosphere in Fig. 14, in high temperature conditions (T ≈ 10 000 K). The metallic lines originate in a more external region (near the disk) in which there is coexistence of Na, Mg, Ca, Ca⁺, Fe, and Fe⁺. This region, that we call the outer magnetosphere, has T ≈ 6000 K. The strength of the Ca II lines relative to the H I lines, as well as the observation of the Fe I absorption against Ca II, suggests that the Ca II emission extends further inward, likely between the outer and inner magnetosphere.

According to shock models, the pre-shock region is responsible for the Balmer continuum emission (Calvet & Gullbring 1998); Gullbring et al. 2000). Therefore, the He I lines and the H I line wings are related to the UV excess. Since the Ca II lines are not produced in the pre-shock region, the fact that they overpredict the accretion luminosity might indicate a stronger emissivity from the outer part of the accretion flow in this system, relative to the stars in the Alcalá et al. (2017) sample.

6.2 Origin of the spectrophotometric variability

The analysis of the photometric data indicates that the magnetospheric interaction between the disk and the star is unsteady and highly dynamic in HMLup. Three-dimensional magnetohydrodynamic (MHD) simulations showed that CTTSs may accrete in either a stable or an unstable regime (Romanova et al. 2003, Romanova et al. 2004; Kulkarni & Romanova 2008; Pantolmos et al. 2020). In the stable regime, accretion proceeds in two funnel streams and forms two polar hot spots on the stellar surface. In the unstable regime, the matter accretes in equatorial tongues and forms multiple hot spots on the surface of the star. The light curves are expected to be periodic with the stellar rotation period (P_⋆) in the stable regime and stochastic in the unstable regime. The transition between these two regimes depends on the ratio between the magnetospheric truncation radius R_T and the corotation radius R_co (Blinova et al. 2016). Accretion is unstable if R_T/R_co ≲ 0.71 and stable otherwise. Blinova et al. (2016) showed that when R_T/R_co decreases below ~0.59, unstable accretion becomes ordered. In the ordered regime, the matter accretes in one or two ordered tongues that rotate with the inner disk period.

The comparison between the morphology of the two TESS light curves suggests that the accretion changed from a chaotic regime in 2019 to a more organized regime in 2021, which could be either the stable or the unstable ordered regime. The 2021 TESS light curve was in a quasi-periodic bursting state with P = 4. 79 ± 0. 28 days. In the stable accretion scenario, this period would correspond to P_⋆. The period detected in the third segment of the long term B band AAVSO light curve, P_B3 = 8.84 days, is roughly twice the 2021 TESS period and can be interpreted as its first harmonic, that is easier to detect in the AAVSO light curve given the lower (~ 1 d) cadence. However, using the 4.79 days period as P_⋆ and the υ sin i and R_⋆ values of Table 2 we obtain an inclination i_⋆ = 23 ± 6 deg for the stellar rotation axis, that is, a star that is almost face-on. This is not in agreement with the observed quasi-periodic behavior of the light curve, that is likely the result of rotational modulation and therefore suggests a higher value for i_⋆, as shown for instance by the 1D hydrodynamic simulations of Robinson et al. (2021).

We propose that the quasi-periodic bursting behavior in the 2021 TESS light curve was the result of an increase in the accretion rate that compressed the magnetosphere and drove the system into a regime of unstable ordered accretion. In this scenario, the 4.79 d period detected in the LSP is not the stellar rotation period but it approximately corresponds to the timescale at which matter rotates at the truncation radius. At the boundary between the chaotic and ordered regime, the ratio of the stellar rotation period to the Keplerian rotation period at R_T is equal to (R_t/R_co)^−3/2 = (0.45)⁻¹ = 2.22 (${\omega }_s^{-1}$ in Blinova et al. 2016). Therefore, P_⋆ should be in the 8–10 days range, compatible with P_B3. Assuming this value as P_⋆ we obtain i_⋆ = 47 ± 15 deg, compatible with the observed quasi-periodic modulation. This value is also in agreement with the inclination of the outer disk axis (Table 2) measured by Ansdell et al. (2016). The cycle-to-cycle variability of the bursts in the 2021 TESS light curve further supports this interpretation. According to the simulations by Blinova et al. (2016), the two tongues may carry different amounts of matter and one of them may sometimes disappear.

The simultaneous spectroscopic data made it possible to characterize the photometric variability of the system during the main burst in the 2021 TESS observation in terms of variations in the accretion rate. However, part of the line variability may be due to the rotational modulation of the accretion flow, as suggested by our analysis of the Fe II 4924 line variability and shown in simulations (e.g., Kurosawa et al. 2008; Kurosawa & Romanova 2013). Although the quasi-periodic behavior of the 2021 TESS light curve supports this scenario, our limited spectroscopic coverage does not allow to test the variability of the emission lines as a function of the rotational phase. We found that during the high accretion state the Balmer lines, except for Hα, develop an absorption component and the emission of low-ionization species increases. This is similar to what has been observed in VW Cha by Zsidi et al. (2022b). The dominant effect during the outburst is an increase in the density of the accretion flow. This produces a spectroscopic signature of the colder region near the disk, that is, the enhanced emission in the metallic species.

Fig. 14 Temperature structure of the accretion flow that emerges from the line analysis.

7 Conclusions

We presented a comprehensive spectrophotometric study of the CTTS HMLup. We examined the 2021 TESS light curve and the simultaneous spectroscopy obtained in the framework of the ULLYSES and PENELLOPE programs. The photometric data from 2021 were also compared with the 2019 TESS light curve and the long-term AAVSO monitoring.

The analysis shows that HMLup is a star in which the accretion process is unsteady and rapidly variable. Using different emission lines, we reconstructed the temperature structure of the magnetosphere. The He I lines and the wings of the Balmer lines are formed in the pre-shock region and are related to the UV excess continuum. The emission lines from metallic species (Na I, Ca I, Ca II, Mg I, Fe I, and Fe II) are instead formed in lower temperature conditions, as already shown by Muzerolle et al. (2001) for the Na I D lines. It must be stressed that the line forming region is more complex than our schematic picture. Density and temperature vary continuously along the accretion flow, and NLTE conditions are expected for the gas (Azevedo et al. 2006). A detailed model should take into account the radiative transfer of light through a structure in which the local velocity and optical depth are functions of the distance from the star. The rich emission line spectrum of HMLup makes this system ideal to test the temperature gradient in non-isothermal magnetospheric models.

The photometric behavior of HMLup can be explained as the result of variations in the accretion rate in a system in which the magnetospheric radius is smaller than the corotation radius. At the epoch sampled by the 2019 TESS light curve HMLup was in an unstable and chaotic regime of accretion. During the TESS 2021 observation, the system entered into a regime of ordered accretion due to an increase in the accretion rate. The detected 4.79 d period can be mistaken with the stellar period but it is likely associated to the Keplerian rotation of the matter at the truncation radius. We observed period variations in the long term AAVSO monitoring, in agreement with the simulations by Blinova et al. (2016). In particular, in the third segment of the B band observations, the system has a periodicity of 8.84 d, which could be the actual rotation period of the star. The spectrum of the system during the burst shows the following three main characteristics: (1) emission in the He I lines increases without changes in the shape of the profile, (2) the Balmer lines develop an absorption component, (3) emission in the metallic species increases.

HMLup gives the possibility to monitor a system with an outburst spectrum that is similar to, but less extreme than, that of EXors. Additional observational material, such as an observing campaign capable of covering a timescale of ~10 days with high resolution spectra at nightly cadence and simultaneous photometry, is needed to better understand the nature of the complex variability of this Classical T Tauri Star.

Acknowledgements

This work has been supported by Deutsche Forschungsgemeinschaft (DFG) in the framework of the YTTHACA Project (469334657) under the project codes STE 1068/9-1 and MA 8447/1-1. A.F.R., and J.M.A. acknowledge financial support from the project PRIN-INAF 2019 “Spectroscop-ically Tracing the Disk Dispersal Evolution” (STRADE) and the Large Grant INAF 2022 “YSOs Outflows, Disks and Accretion: towards a global framework for the evolution of planet forming systems” (YODA). C.F.M., J.C.W., and K.M. are funded by the European Union (ERC, WANDA, 101039452). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them. This work benefited from discussions with the ODYSSEUS team (https://sites.bu.edu/odysseus/) (HST AR-16129). The authors acknowledge Eleonora Fiorellino, Suzan Edwards, and René Oudmaijer for suggestions on this work. A.A. acknowledges Steve Shore for valuable discussions on this system. Funding for the TESS mission is provided by NASA’s Science Mission directorate. The authors acknowledge with thanks the variable star observations from the AAVSO International Database contributed by observers worldwide and used in this research, and Elizabeth Waagen for coordinating the AAVSO Alerts. The authors acknowledge the use of the electronic bibliography maintained by the NASA/ADS (https://ui.adsabs.harvard.edu) system.

Appendix A Flux calibration of the ESPRESSO spectra

We flux-calibrated the four ESPRESSO spectra by re-scaling them to match the available BVR_cI_c photometry. Since no photometric data were obtained simultaneously with the spectra, we linearly interpolated the AAVSO light curves to obtain the magnitude values at the ESPRESSO epochs. We treated the flux calibration as a least squares minimization problem, namely, we found the best scaling factor r that matched the BVR_cI_c photometry by minimizing the chi-square function $${\chi }^2\left(r\right)={\displaystyle \sum _{i=B,V,R_{\text{c}},I_{\text{c}}}{\left[\frac{m_i-M_i\left(r\right)}{d{m}_i}\right]}^2}$$(A.1)

where m_i and dm_i are the observed magnitudes and uncertainties and M_i(r) are the r-dependent “model” magnitudes, obtained by multiplying the original ESPRESSO spectra by r and integrating over the passbands i using the respective filter transmission curve. The uncertainty on the scaling factor was estimated as the standard deviation of exp (−χ²). Since the ESPRESSO spectra do not reach 9100 Å, the upper limit of the I_c filter, we appended a rescaled version of the X-shooter spectrum to the ESPRESSO ones long-ward of 7580 Å. The mean ratio of the two spectra between 7200 and 7500 Å was used as scaling factor for X-shooter. We repeated the fit without including the I_c magnitudes, and the scaling parameter did not change within the uncertainties. The obtained scaling factors are reported in Table A.1. To check the reliability of this procedure, we compared the ESPRESSO spectrum at MJD – MJD₀ = 3.83 d to the X-shooter spectrum, which was observed only ~1 h before. We computed the ratio of the two spectra, finding a mean flux ratio of 0.57, in agreement with the scaling factor we obtained from the best fit procedure.

Appendix B Saha and Boltzmann equations

In LTE, the population of the atomic levels of an ion is given by the Saha and Boltzmann equations (e.g., Hubeny & Mihalas 2014). The Saha equation provides the ratio of the populations of two successive ionization stages for a given atomic species $$\frac{N_{\text{j}+1}}{N_{\text{j}}}=\frac{1}{n_{\text{e}}}{\left(\frac{2\pi m_{\text{e}}{k}_{\text{B}}T}{h^2}\right)}^{3/2}\frac{2U_{\text{j}+1}}{U_j\left(T\right)}{e}^{-{\chi }_j/k_{\text{B}}T}$$(B.1)

where T is the temperature, N_j and U_j are the number of atoms and the partition function of the j^th ionization stage, n_e is the electron density, m_e is the electron mass, h is the Planck constant, k_B is the Boltzmann constant and χ_j is the ionization potential of the lower ionization stage. The Boltzmann equation governs instead the population of the atomic levels in a given ion $$\frac{n_j}{N_{\text{ion}}}=\frac{g_j}{U_{\text{ion}}\left(T\right)}{e}^{-E_j/k_{\text{B}}T}$$(B.2)

where n_j is the population of the level with energy E_j and multiplicity g_j relative to the population of the ion N_ion.

To compute the partition functions of the ions we obtained the energy levels and the ionization potentials from the NIST Atomic Spectra Database⁸. If we assume that emission lines are optically thin, we can determine the properties of the emitting region by means of line ratios. The line strength can be defined as $$j_{\text{line}}=\frac{A_{\text{ji}}}{{\lambda }_{\text{line}}}{n}_j\left(n_{\text{e}},T\right)$$(B.3)

where A_ji and λ_line are the Einstein coefficient and the wavelength of the transition and n_j is the population of the upper level of the line, that depends on N_ion through the Boltzmann equation. The Saha equation gives N_ion as a percentage relative to the element population N_elem. The line ratio between two lines from the same species is therefore a function of n_e and T.

In Sect. 4.3, we use this formalism to compute predicted LTE line flux ratios for a grid of n_e and T values, with 10 ≤ log n_e ≤ 14 and 4000 K ≤ T ≤ 8000 K. The observed line ratios were statistically compared to the predicted LTE ratios by maximizing the likelihood function $$ℒ\left(n_{\text{e}},T\right)=\exp \left[-{\left(\frac{r_{\circ }-r_{\text{m}}\left(n_{\text{e}},T\right)}{d{r}_{\circ }}\right)}^2\right]$$(B.4)

where r_o and dr_o are the observed ratio and its uncertainty and r_m is the expected LTE ratio.

Appendix C Triple Gaussian model

In Sects. 4.2 and 5.3 we used a triple Gaussian model to fit the continuum-normalized emission lines in the spectrum of HM Lup. The fit model is $$F_{\upsilon }=1+{\displaystyle \sum _{i=1}^3{C}_i\exp \left[\frac{{\left(\upsilon -{\upsilon }_i\right)}^2}{2{\sigma }_i^2}\right]}.$$(C.1)

For each Gaussian C_i is the amplitude, υ_i is the velocity center and σ_i is the standard deviation.

References

All Tables

All Figures

Fig. 1 Light curves and timing of simultaneous spectroscopy for HM Lup. The topmost panel shows the 2021 10 minute cadence TESS light curve (Sector 38). The inset displays the Lomb–Scargle Periodogram for the TESS light curve. Overlaid in purple is a boxcar smoothed version of the light curve produced by phase-folding the data with the period of the highest peak in the Lomb–Scargle periodogram, 4.79 day. The phase ϕ was computed with this period and ϕ = 0 at the maximum of the TESS light curve, MJD 59 336.97. The vertical dashed lines mark the epochs of the simultaneous spectroscopic observations. Here ES = ESPRESSO and XS = X-shooter. The other panels show the part of the AAVSO BVRcIc photometry that is simultaneous with TESS, with a linear interpolation as a guideline. The large open squares mark the synthetic photometry obtained from the X-shooter spectrum for the four filters. We defined MJD0 = 59 333.363 as the beginning of TESS Sector 38 observation.

In the text

	Fig. 2 TESS 2019 two-minute cadence light curve (Sector 12) of HM Lup. The inset shows the Lomb-Scargle periodogram. Overlaid in purple is a boxcar smoothed version of the light curve produced by phase-folding the data with the period of the highest peak in the Lomb-Scargle periodogram, 5.41 d. We defined MJD1 = 58 628.050 as the beginning of the TESS Sector 12 observation.
In the text

	Fig. 3 Ṁacc-M★ diagram for X-shooter targets in Lupus, using the values from Manara et al. (2023). The black triangles indicate upper limits on the measured accretion rate. The blue, green and red diamonds mark the accretion rate of HM Lup in 2012 and 2021 (Sect. 5.2). The dashed gray line is a linear best fit to the data.
In the text

Fig. 4 Selection of continuum subtracted emission lines in the ESPRESSO spectra of HM Lup. The colors for the ESPRESSO spectra are the same as in Fig. 1. The lines from the 2012 X-shooter spectrum are shown in green for comparison. The vertical dashed lines mark the radial velocity of the system, while the horizontal dashed lines highlight the zero-flux level for each set of lines. The Fe I and Ca I lines were smoothed with a 7 points boxcar filter.

In the text

Fig. 5 Continuum subtracted profiles of selected magnetospheric lines in the epoch 2 ESPRESSO spectrum. In the left panel, H9 (3835 Å) and He I 5876. In the right panel, Hδ and Ca II K. The vertical dashed line marks the radial velocity of the system. The H9 line was smoothed with a boxcar filter and multiplied by 0.5 to match the He I 5876 wings. The Ca II K line was rescaled to Hδ in a similar way.

In the text

Fig. 6 Reconstruction of the Fe I absorption. All profiles are from the epoch 2 ESPRESSO spectrum. All lines are plotted relative to the rest wavelength of Fe I 3969. The top panel shows the continuum subtracted profile of the Ca II H + Hϵ blend (black) compared to the Ca II K (green) and H8 (blue) lines. The red line is the model for the unabsorbed blend (Eq. (2)). The vertical dashed lines indicate the velocity displacement of Ca II H, Fe I 3969 and Hϵ relative to the rest wavelength of Fe I 3969. The middle panel shows the Fe I 3969 absorption profile, reconstructed as explained in Sect. 4.2, and its best fit. The bottom panel shows the Fe I 4064 and 4132 emission lines.

In the text

Fig. 7 Comparison between the Hβ and Fe I 3969 absorption components. Top panel: triple Gaussian fit of the continuum normalized Hβ line in the epoch 2 ESPRESSO spectrum. The dashed lines show the three Gaussian components. Bottom panel: COMPARISON between the continuum normalized absorption profiles in Hβ and Fe I 3969. The red line is the Gaussian best fit of the Fe I 3969 from Fig. 6. The blue line is the sum of the two Gaussian absorption components in Hβ, rescaled relative to the continuum to allow the comparison with the Fe I 3969 absorption profile.

In the text

	Fig. 8 Continuum subtracted profiles of the Ca II IRT lines in the X-shooter spectrum at MJD – MJD0 = 3.80 days, almost simultaneous to the ESPRESSO epoch 2 spectrum.
In the text

	Fig. 9 Result of the Saha-Boltzmann analysis for Ca and Mg. Observed line ratios were computed from the epoch 2 ESPRESSO spectrum. Left panel: likelihood for the observed Ca II 8542 vs. Ca I 4227 ratio as a function of log ne and T. Right panel: theoretical Mg II 4481 vs. Mg I 5184 ratios as a function of T for log ne [cm−3] = 11,12,13,14. The red area marks the upper limit on the observed ratio.
In the text

Fig. 10 Continuous wavelet analysis of the AAVSO B band photometry. The top panel shows the CWT. The inset on the left displays the LSP of the whole dataset, together with a temporal average of the CWT. The orange horizontal dot-dashed lines mark the periods obtained from the LSP analysis, i.e., PB2 = 6.24 days and PB3 = 8.84 day, while the red one indicates the 2021 TESS period, PTESS = 4.79 days. The bottom panel displays the AAVSO B band light curve and its linear interpolation. The insets show the LSP for three different time segments. We excluded the shaded areas from the calculation of the LSPs to highlight the different timescales detected.

In the text

	Fig. 11 VLT/X-shooter and HST/STIS spectra of HM Lup in the STIS wavelength range. The spectra were smoothed to the HST resolution with a Gaussian filter for clarity.
In the text

Fig. 12 Comparison between the accretion luminosity derived from emission lines and that derived from the UV excess in the 2021 X-shooter spectrum. Lacc,UV and its uncertainty are displayed with the red dashed line and the red shaded area. The blue dashed line and the blue shaded area are the value of Lacc,lines and its uncertainty, obtained from a weighted average of the line values. The lines marked in black were excluded from the calculation of the average.

In the text

	Fig. 13 Triple Gaussian fit of the Fe II 4924 line profiles in the four ESPRESSO spectra. See Appendix C for the fit function.
In the text

	Fig. 14 Temperature structure of the accretion flow that emerges from the line analysis.
In the text