Accretion Funnel Reconfiguration during an Outburst in a Young Stellar Object: EX Lupi

EX Lupi, a low-mass young stellar object, went into an accretion-driven outburst in 2022 March. The outburst caused a sudden phase change of ∼112° ± 5° in periodically oscillating multiband lightcurves. Our high-resolution spectra obtained with the High Resolution Spectrograph (HRS) on board the Southern African Large Telescope also revealed a consistent phase change in the periodically varying radial velocities (RVs), along with an increase in the RV amplitude of various emission lines. The phase change and increase in RV amplitude morphologically translates to a change in the azimuthal and latitudinal location of the accretion hotspot over the stellar surface, which indicates a reconfiguration of the accretion funnel geometry. Our three-dimensional magnetohydrodynamic simulations reproduce the phase change for EX Lupi. To explain the observations, we explored the possibility of forward shifting of the dipolar accretion funnel as well as the possibility of the emergence of a new accretion funnel. During the outburst, we also found evidence of the hotspot’s morphology extending azimuthally asymmetrically with a leading hot edge and cold tail along the stellar rotation. Further, our high-cadence photometry showed that the accretion flow has clumps. We also detected possible clumpy accretion events in the HRS spectra that showed episodically highly blueshifted wings in the Ca ii IR triplet and Balmer H lines.


INTRODUCTION
In the last three decades or so, magnetospheric accretion has been well-established in the low-mass young stellar objects (hereafter, YSOs; see reviews in Bouvier et al. 2007;Hartmann et al. 2016).The phenomenon of magnetospheric accretion links the accreting matter from the circumstellar disk to the stellar surface.Originally developed for the highly magnetized accreting stars like white dwarfs, neutron stars, etc., the theory of magnetospheric accretion was later applied to the YSOs (Ghosh & Lamb 1978;Koenigl 1991).The observation of Zeeman splitting in the photospheric lines of the YSOs indicated a surface magnetic field strength of the order of 1-3 kG in the YSOs (Johns-Krull et al. 1999, 2001).The magnetic field truncates the circumstellar disk at inner disk radius equal to a few stellar radii.The material from the inner disk follows the stellar magnetic field lines and gets accreted on the stellar surface near the magnetic poles.The pathway of magnetic field lines channeling the accretion matter is called the magnetic accretion funnel (or accretion funnel).This freely falling material impacts the star near the photosphere at supersonic speed, creating a shocked and hot region (∼ 10 4 K) called a hotspot (Lamzin 1998;Romanova & Owocki 2015;Hartmann et al. 2016).The success of the magnetospheric accretion mechanism lies in its ability to explain the periodically varying lightcurves by the presence of the hotspot on the rotating star and the explanation of high-velocity absorption in the spectral lines, like Hα, in terms of absorption by cool material falling through the accretion channel (Herbst et al. 1994).The interaction of the stellar magnetosphere with the disk is so complex that it gives rise to all sorts of accretion phenomena with varying thermal and morphological structures of the accretion funnel and the hotspot.
The three-dimensional (3D) and 2D magnetohydrodynamic (MHD) simulations have revealed different properties of the magnetospheric accretion.Simulations have shown that the disk is stopped by the magnetosphere at the magnetospheric radius r m where the magnetic pressure of the magnetosphere matches the matter pressure in the disk (Romanova et al. 2002;Zanni & Ferreira 2009;Romanova & Owocki 2015;Takasao et al. 2022).3D simulations of accretion onto stars with the tilted dipole magnetosphere have shown that in many instances, matter accretes onto a star in two funnel streams (as predicted by theoretical work) and forms the banana-shaped spots which have the inhomogeneous density and flux distribution (Romanova et al. 2003(Romanova et al. , 2004;;Espaillat et al. 2021), which can be azimuthally asymmetric (Espaillat et al. 2021).The spots may shift along the surface of the star if the inner disk rotates more rapidly than the closed magnetosphere (Bachetti et al. 2010;Kulkarni & Romanova 2013).Simulations also have shown a possibility of unstable accretion where matter penetrates through the magnetosphere in the form of temporary equatorial filaments, or "tongues" which create hotspots at lower latitudes (Kulkarni & Romanova 2008, 2009;Kurosawa & Romanova 2013;Blinova et al. 2016;Burdonov et al. 2022).This regime is favorable if the magnetosphere rotates slower than the star and conditions for the magnetic Rayleigh-Taylor instability are satisfied (e.g., Spruit et al. 1995).
The magnetic field of stars may be more complex than the tilted dipole (e.g., Donati et al. 1999).3D MHD models of accretion onto stars with complex magnetic fields have shown that the dipole component typically dominates in the disk-magnetosphere interaction, while the higher multipoles determine the shape of the spots on stellar surface (e.g., Long et al. 2011Long et al. , 2012)).
More recently, 3D simulations were performed in a general case where both the magnetic and rotational axes of the star are tilted about the rotational axis of the disk.Simulations have shown that the inner disk tends to be shifted towards the equatorial (rotational) plane of the star and other features associated with the tilt of the rotational axis.(Romanova et al. 2021).
For a single star, only the mass accretion rate or magnetic field strength and topology can vary (Johns-Krull et al. 2022;Finociety et al. 2023), which would then change the accretion funnel geometry.Any change in the accretion funnel geometry would alter the thermal and morphological structure of the hotspot.Long et al. (2012); Kulkarni & Romanova (2013) predicted that there could be changes in size, shape, azimuthal and polar positions and thermal structure of the hotspot, during an increase in accretion rate.The magnitude of these changes depends upon the magnetic field structure, rotation period, magnetic obliquity, etc. (Kulkarni & Romanova 2013).Thus, the stars with varying accretion rates become an ideal laboratory for comparisons of simulations with observations.
A class of YSOs, FUors and EXors, undergo episodic accretion-driven outbursts where the accretion rate increases by a factor from a few to 100.FU Orionis is the prototype of the outbursting YSOs, FUors.FUors undergo outbursts on the timescale of years to decades followed by a quiescence for decades to century (Herbig 1977).The EXors undergo outbursts on the timescale of months followed by a quiescence for years to decades.EXors are more frequent repeaters of the outburst, with a lesser brightness increase than FUors (Audard et al. 2014).EX Lupi, the prototype of EXors, is a touchstone star as it has undergone multiple outbursts.It underwent three major outbursts (∆V ∼ 5 mag) along with several small outbursts since its first reported outburst observation in the 1940s (McLaughlin 1946;Herbig 1977;Bateson et al. 1991;Herbig et al. 2001;Herbig 2007;As-pin et al. 2010;Zhou & Herczeg 2022;Kóspál et al. 2022;Cruz-Sáenz de Miera et al. 2023;Wang et al. 2023).Recently, EX Lupi went into an accretion-driven outburst in 2022 March (Zhou & Herczeg 2022;Kóspál et al. 2022;Cruz-Sáenz de Miera et al. 2023;Wang et al. 2023).The state of high accretion rate remained for about 100 days before returning to its pre-outburst photometric state.Photometrically, the star brightened by ∆g ∼ 2 mag.The wealth of high-cadence photometry and highresolution spectroscopy during the outburst made EX Lupi an ideal laboratory to enhance our understanding of the magnetospheric accretion phenomenon.
Generally, the star's disk is inclined with respect to the line of sight by the inclination angle (i ).Sipos et al. (2009) performed a spectral energy distribution (SED) fit onto the quiescent state fluxes of EX Lupi.The authors found that the inclination angle between i=20 • to i=40 • produces similar results.Goto et al. (2011) modelled CO fundamental bands at 4.6-5.0µm and found that i=40 • -50 • produces good fits to the observation.Further Kóspál et al. (2014) and Sicilia-Aguilar et al. (2015) found vsin(i) = 4.4 km/s which corresponds to an inclination angle of 25 • .5.Hales et al. (2018) determined an inclination angle of 32 • .4± 0 • .9 by fitting a 2D elliptical Gaussian over the 1.29 mm continuum image of EX Lupi star-disk system taken with the Atacama Large Millimeter/submillimeter Array (ALMA).Hales et al. (2018) also calculated the disk inclination angle of 38 • ± 4 • from the modelling of CO isotopologues.This CO isotopologues modelling assumed a smaller value of stellar mass of 0.5M ⊙ for EX Lupi.Recent work by Sicilia-Aguilar et al. (2023) used an intermediate inclination angle of 20 • -40 • to explain color modulation in lightcurves with a single hotspot model.
The authors further modelled the width of Ca II 8662 Å lines and found that an inclination angle of around 45 • best matches the observation.Though dust continuum traces the outer disk, we adopted an inclination angle of 32 • .4± 0 • .9 as it complies with Sipos et al. (2009); Sicilia-Aguilar et al. (2015).Wang et al. (2023) also used a value of i=30 • to illustrate the two-component model of the hotspot.The inclination angle and the magnetic obliquity are such that (90 • -i > Θ = 13 • ) only one hotspot is visible.
Since EX Lupi is located inside a star-forming region, it could be reddened by the dust extinction.Alcalá et al. (2017) found A V = 1.1 mag for EX Lupi by doing a chi-square minimization of the X-shooter spectra with a sum of the stellar photosphere and slab model for the hotspot.Other studies have suggested A V = 0 mag (e.g.Sipos et al. 2009).Recently, Wang et al. (2023) estimated an extinction of A V =0.1 mag towards EX Lupi.We used the conservative value of A V = 1.1 mag in our paper.
We elaborate on the observation and reduction of the spectroscopic and photometric data in Section 2. Analysis of the observed data is done in Section 3 and we conclude this section by listing the key observational results.In Section 4, we discuss the observational key results in detail and make further inferences about EX Lupi.We conclude our work in Section 5 by enumerating the key understandings about EX Lupi and the accretion process in general.Appendices describe the periodicity of EX Lupi, spectral line fittings, an alternate hypothesis of phase change in lightcurves and radial velocity and estimation of clumps' masses.

OBSERVATIONS AND DATA REDUCTION
The spectroscopic and photometric observations we obtained from various telescopes are elaborated below.

HRS-SALT
Beginning almost at the peak of the outburst (∼ 2022 March 15), we observed EX Lupi spectroscopically on 16 epochs, over the next two and half months, with the High Resolution Spectrograph (HRS) in Medium Resolution (MR) mode (R = 37000; Bramall et al. (2010); Crause et al. (2014)) on the 11m Southern African Large Telescope (SALT; Buckley et al. 2006).The spectra were collected under 'The SALT Transient Programme' (2021-2-LSP-001) (Buckley 2017).HRS covers a wavelength range of 3740 Å -8780 Å using blue and red cameras.The wavelength coverage of the blue camera is 3740 Å -5560 Å while the red camera covers 5450 Å -8780 Å.The 2D to 1D extraction and wavelength calibration were done by the MIDAS pipeline (Kniazev et al. 2016a,b).Barycentric correction was applied to the spectra using a Python-based module: barycorrpy (Kanodia & Wright 2018).The instrumental response was estimated from the observation of the stars HD074000, HD38949 and HD205905 on 2022 March 02, 2022 April 12 and 2022 May 15 respectively.The HRS-SALT spectra of these standard stars were divided by the fluxcalibrated spectra of the respective stars taken from the CALSPEC database (Bohlin et al. 2014).The average of the three observations after respective median normalization was taken as the instrument response curve.EX Lupi spectra were divided by this instrument response curve to get the response corrected spectra.The spectra were flux-calibrated against the corresponding epoch's photometric flux.Since ASAS-SN g-band observations were more densely sampled, we used it for interpolating the g-band flux to the epochs of SALT HRS observations.Interpolation was done using a Gaussian Process with a combination of the Matern-5/2 and exponentialsine-squared kernel in the Python module, George (Ambikasaran et al. 2015).B, V and I -band magnitudes from our TMMT multiband observations (see Section 2.2) were found to correlate with ASAS-SN g band observations.Hence, we used this correlation to convert the interpolated g-band magnitudes to interpolated, B, V and I magnitudes.We found that flux calibration of spectra using the TMMT I-band flux produces spectral flux consistent with TMMT B-band in the blue part of the HRS spectra.Thus, the TMMT I-band interpolated fluxes were then used to scale the HRS spectra of EX Lupi.The spectra were further corrected for extinction (A V =1.1 mag, Alcalá et al. 2017).The extinction curve was generated from Cardelli et al. (1989) via the Python module named dust-extinction for R V = 3.1.The error for each pixel value in spectra was taken to be the square root of the effective number of photoelectrons detected.A typical signal-to-noise ratio (SNR) of 75 was attained at around 5000 Å.

TMMT
We followed up EX Lupi during the 2022 March outburst and the post-outburst period using the Threehundred MilliMeter (300 mm) Telescope (TMMT, Monson et al. 2017) at Las Campanas Observatory in Chile.TMMT uses an Apogee Alta U42-D09 CCD Camera and we observed EX Lupi at a very high cadence for 52 nights with TMMT, collecting 1812, 1817 and 1817 frames in B, V and I bands, respectively.During the outburst, we observed for 43 nights starting from 2022 March 13 (JD = 2459651.7106)to 2022 June 20 (JD = 2459750.7426)and 1776, 1781 and 1782 frames were taken in B, V and I bands, respectively.The rest of the data were observed in 2023 January, covering 9 nights.The observation was conducted nearly at a one-day cadence, and inside each night, we conducted three-minute cadence exposures for about 4-5 hours.The photometric observation logs are tabulated in Table 1.Only a portion of the logs are provided here, while a complete log is available online in machine-readable format.

LCRO
Along with TMMT, we also followed up EX Lupi during the 2022 March outburst and the post-outburst period using the 305 mm Las Campanas Remote Observatory (LCRO) telescope at the Las Campanas Observatory in Chile.LCRO uses an FLI Proline 16803 CCD Camera with g ′ , r ′ , and i ′ photometric bands.We observed EX Lupi over 85 nights with a total of 933, 928 and 922 frames in g ′ , r ′ , and i ′ bands respectively.EX Lupi was observed on 11 nights during the outburst (2022 March 13 (JD = 2459651.7093)to 2022 March 25 (JD = 2459663.9104))and for 74 nights in the postoutburst state starting from 2023 January 24 (JD = 2459972.8347)to 2023 May 14 (JD = 2460078.9175).

CTIO
We conducted multiband photometric observations of EX Lupi with an FLI camera on the Prompt6 telescope at Cerro Tololo Inter-American Observatory (CTIO), Chile.The telescope operates under the Skynet Robotic Telescope Network (Martin et al. 2019).Prompt6 is a 0.4m aperture telescope with a field of view of 15.1 ′ × 15.1 ′ .It has U, B, V, R, Hα, OIII and SII filters.We observed EX Lupi in B and V bands from almost the peak of the outburst, 2022 March 12 (JD=2459650.768)to the subsequent 60 days, 2022 May 11 (JD=2459710.5854),collecting 51 epochs of nearly one-day cadence photometric data.The images were corrected for dark currents and were also flat-corrected.

ASAS-SN
High-cadence g and V -band lightcurves are used from All-Sky Automated Survey for Supernovae (ASAS-SN, Shappee et al. 2014).The lightcurve we used in our study covers from 2017 May 27 (HJD = 2457900.7206 2 ) to 2023 September 6 (HJD = 2460194.3491).For reference purposes, the time = 0 of all the calculations related to lightcurves is taken at HJD = 2457900.7206(see Section 3). 2 The HJD (Heliocentric Julian Date), BJD (Barycentric Julian Date), and JD (Julian Date) are not corrected to the same scale.
The maximum difference of 8 min among them corresponds to an insignificant ∼ 0 • .27 of EX Lupi's rotation.Our phase calculation on lightcurves has a gross error of ∼ 3-10 • .

AAVSO
We have used the high cadence multiband U, B, V, R and I lightcurves from the American Association of Variable Star Observers (AAVSO) data archive.The data used for our study covers the whole outburst starting from 2022 March 13 (JD = 2459651.9002)to the subsequent quiescence up to 2022 October 17 (JD = 2459870.4878).A significant amount of AAVSO data was contributed by the Remote Observatory Atacama Desert telescope (ROAD; Hambsch 2012) and Perth Observatory's R-COP telescope.The ROAD telescope is a fully automated telescope located in Chile that obtains nightly photometry in Astrodon B, V, I and Clear bands for a wide range of astronomical projects.It consists of a 40-cm f/6.8 Optimized Dall Kirkham telescope and uses an FLI ML16803 CCD camera with a 4k×4k array with pixels of 9 µm in size.Data was reduced using a custom pipeline and then published in the AAVSO database.The R-COP is a telescope partnership that stands for Remote Telescope Partnership: Clarion University -Science in Motion, Oil Region Astronomical Society, and Perth Observatory.The R-COP telescope is a part of the Skynet Robotic Telescope Network.

Photometry and Radial Velocities from previous literatures
We have also used the RV estimates for EX Lupi from the literature.Kóspál et al. (2014); Sicilia-Aguilar et al. (2015); Campbell- White et al. (2021) presented the RVs of EX Lupi by cross-correlating the spectra with template spectrum and by fitting the emission spectral lines.These RVs cover an interval from the year 2007 to 2019.The time evolution of the RV phase is taken from Campbell- White et al. (2021).We have used the RV amplitudes from Sicilia-Aguilar et al. (2015); Campbell- White et al. (2021).We have also used multiband high cadence photometry at the onset of the outburst from Kóspál et al. (2022).

ANALYSIS AND RESULTS
The central part of this work is to understand the phase evolution of the various waveband lightcurves and the periodically modulated RVs.The techniques used to measure the phases of the lightcurves and the RV are discussed below in Sections 3.1 and 3.2.

Phases of the lightcurves
A start time (t = 0) is needed to define a reference point for the calculation of the phase of the periodically varying lightcurves.We took t = 0 as the first epoch of our ASAS-SN lightcurve.It is JD = 2457900.7206(hereafter startJD).All the phase estimates are done on the light curves after converting from magnitude to flux units.Lightcurve cutouts are made at an interval of 40 days.A straight line is fitted to each individual lightcurve cutout to remove any linear trend.Then the cutout is cross-correlated with a sine wave of the form sin( 2π(t−t0) P + ϕ) where P is the rotation period of EX Lupi and the t is the epoch of observation in JD, and t 0 is startJD.The phase ϕ is varied over [0, 2π) at an interval of 1 • , and the corresponding Pearson correlation coefficient is calculated.The value of ϕ which corresponds to the maximum estimate of the Pearson correlation coefficient is the phase of the lightcurve cutout.This phase ϕ is associated with the first epoch of the cutout.
We calculated the phase along with the phase error estimates for every lightcurve cutout by performing the Monte Carlo method described below.Random light curves are generated by sampling from a Gaussian noise distribution centered at the original light curve values of each epoch.The standard deviation of this Gaussian distribution is taken as a quadrature sum of the lightcurve flux error and stellar activity error.We assigned the stellar activity error to the robust standard deviation of the residual of a sine curve fit onto the original lightcurve cutout.The standard deviation of the Gaussian distribution is dominated by the stellar activity error in almost all the cases.Then, phase is calculated on this simulated lightcurve cutout as described in the previous paragraph.The simulation is run 1000 times.The mean and standard deviations of the phase estimation from 1000 runs are taken as phase and phase error, respectively.
To confirm our phase estimates are not driven by any flares, for our phase analysis of the ASAS-SN lightcurve, we also manually selected clean durations of the lightcurve devoid of any flares and other variability.Three such epochs were identified in the quiescent state.Three more epochs of phases are estimated on manually selected regions during the state of outburst.Similarly, the phase estimation is done for other lightcurves of EX Lupi.
The AAVSO U, B, V, R, and I-band lightcurves are first converted into flux units (erg/s/cm 2 / Å) following zero-point flux in Bessell et al. (1998) and then the phase calculation is done.The phase calculations on the CTIO B and V -band lightcurves are also done after converting them into the flux unit.The CTIO lightcurves cover 60 days, and so a single phase is calculated onto the whole lightcurve.TESS lightcurves in sector-12 and sector-39 have also been used to estimate the phases.We calculated the phases of the g ′ , r ′ , i ′ and z ′ band lightcurves, taken from Kóspál et al. (2022), at the onset of the 2022 March outburst.The phases are calculated one year after the outburst using LCRO g ′ , r ′ and i ′ band lightcurves.The Sloan, g ′ , r ′ , i ′ and z ′ band, magnitudes are converted into flux unit by converting them to AB magnitudes and then using AB zeropoint flux density of 3631 Jy.
The ASAS-SN g-band lightcurve phases during the quiescence, prior to the 2022 March outburst, are similar at around 100 • -110 • , as shown in the upper panel of Figure 1.The ASAS-SN g-band phase (green circle) around 2019 is 102 • .9± 2 • .7 and the close-by sector-12 TESS lightcurve (inverted orange triangle) has a phase value of 103 • .5 ± 0 • .4 while the phases of ASAS-SN g-band lightcurve and TESS sector-39 lightcurves, around 2021, are 111 • .3± 3 • .2 and 113 • .0 ± 0 • .7,respectively.In Figure 1, the grey circles are the phases from the ASAS-SN g-band lightcurve while the green circles are phases estimated on manually selected cleaner flare-free ASAS-SN g-band lightcurve cutouts.Thus, the phases from the TESS lightcurves are consistent with the nearby epoch ASAS-SN phases.The multiband phases from Kóspál et al. (2022) lightcurves (diamonds) at the onset of the outburst have values consistent with that of the TESS sector-39 lightcurve.g ′ -band phase from Kóspál et al. (2022) is 115 • .4± 4 • .2 while maximum phase difference among the g ′ , r ′ , i ′ and z ′ bands is 10 • .0.
At the peak of the outburst, the phase of the ASAS-SN lightcurve increased to 223 • .9± 3 • .8.The B (236 • .0 ± 7 • .6)and V -band (225 • .3± 8 • .8)phases derived from CTIO (squares) during the peak of the outburst are consistent with the ASAS-SN g-band phase, supporting the phase increase.AAVSO multiband phases, U (241 2), also agree well with the corresponding band phases from CTIO and ASAS-SN, during the peak of the outburst.We estimated the phase change during the outburst by comparing the ASAS-SN phase values (green circle, in Figure 1) from close to the peak of the outburst (223 • .9± 3 • .8)and around 2021 (111 • .3± 3 • .2).The ASAS-SN lightcurve yields an increase in phase of 112 • .6 ± 5 • .0.
During the decay of the outburst (see Section 4.5), the phases decayed differentially in multiband AAVSO (dots) and ASAS-SN lightcurves.Among AAVSO multibands, the U -band phases decreased smallest while I-band phases decreased by largest phase values.From the peak of the outburst to the subsequent post-outburst state, the U -band phases remained largest followed by B, V , R and I band phases, sequentially (see Figure 1).The U and B band phases stayed close to their peak outburst phase values, after the outburst was over.Further, the lightcurve phases could not regain their preoutburst value even after one year of the outburst.The first LCRO phases (pentagon) estimated on 2023 February 7 showed a larger than pre-outburst phase value in the g ′ -band.Further, the LCRO g ′ , r ′ and i ′ bands showed a phase difference among the bands: g ′ (217 ) and i ′ (151 • .9± 4 • .4),The ASAS-SN g-band phase values are within ∼10 • to the close epoch LCRO g ′ -band phases.However, the latest ASAS-SN phase value from 2023 July to 2023 September has returned to the pre-outburst quiescence phase value (see grey circles in Figure 1).

Phase of the radial velocity
We took 16 epochs of high-cadence spectra with HRS on SALT, during the 2022 March outburst.Some of the high energetic emission spectral lines of He I, He II, Fe I, Fe II, Mg I, and Si II, which are visibly distinguishable, are modelled with multi-Gaussians + a constant continuum (see Table 2 for energy levels of the lines).The modelled spectral lines have an SNR as low as ∼15 on a few epochs of quiescence.These high-energy lines are expected to originate in the hotspot and hence carry information about the hotspot dynamics.As inferred from the periodically varying lightcurve, the hotspot comes and goes out of our view as the star rotates3 , and thus a periodicity in the RVs is expected for the spectral lines originating in the hotspot.From the fitted lines, the RV of the Gaussian representing the narrow component (NC) of the spectral lines are phase-folded and fitted with a sine wave of the form A sin( 2π(t−t0) where A, P , ϕ and V 0 are the amplitude of the RV, rotation period of the star (7.417 days), phase of the RV and an arbitrary constant respectively.Here, t is the epoch of the spectra (in JD) and t 0 is startJD.The ϕ of this sine fit gives the RV phase at our startJD.To avoid outliers affecting the sine fit, we did a leave-oneout analysis.One RV point is removed at a time, and sine is fitted.This is repeated for each of the 16 RV points.For each spectral line, the median of each of the fit parameters (A, ϕ and V 0 ) and best-fit errors are taken as the value of the parameter and respective error.The phases cover a range of ∼ 110 • -155 • except for the He II 4686 Å, which has a phase of 176 • .Except He II 4686 Å, all the spectral line RVs are clumped around in a small phase band as shown in Figure 2. We also found that the RV amplitudes, inferred from the aforementioned sine fit onto the RVs, have increased consistently for all the spectral lines with respect to the quiescent state RV amplitudes (see Table 2).RV amplitude is also largest for He II 4686 Å.We do not understand the cause behind the large phase and large RV amplitude of He II 4686 Å.
Since we do not have immediate pre-outburst spectra from HRS-SALT, we compared our RV phases with RV estimates from the literature to calculate any preto post-outburst phase changes in periodically varying RVs, similar to a phase change observed in the photometry.2012) in the 2008 outburst spectra.There is no apparent phase difference among the RVs from pre-to post-2008 outburst (see Figure 2 in Kóspál et al. 2014).We estimated the phases of the absorption and emission line RVs by fitting a sine curve to this dataset as done in Section 3.2 except that the zero point (t=0) is taken as JD = 2454309.6147(in order to be consistent with Campbell- White et al. (2021), see paragraph below).For the emission line RVs, the phase is 327 • .3± 13 • .3 while for the absorption line, the RV phase is 114 • .0 ± 3 • .1.Considering EX Lupi's rotation period of 7.417 days, these RV phases are translated from JD = 2454309.6147to our startJD.The re- spective RV phases for the emission and absorption lines are 30 • .5 ± 13 • .3 and 176 • .0 ± 3 • .1 respectively (at our startJD).This emission line phase is plotted in Figure 2 as a blue pentagon.For comparing the phase shift inferred from the emissions lines with that of the photometric lightcurves, we have overlaid a black curve representing the 112 • phase shift inferred from the photometric lightcurves on the same plot.Sicilia-Aguilar et al. ( 2012) estimated the RVs for metallic lines during the state of outburst in the year 2008.The outburst RVs were estimated on 5 epochs spread over 56 days in 2008 April -2008 June.We could not fit a sine curve onto these RVs due to insufficient data points and thus we do not have any direct measure of phase during the major outburst of the year 2008.
Furthermore, Campbell-White et al. ( 2021) did multi-Gaussian fittings to various narrow metallic and He I spectral lines from the spectra observed for EX Lupi during its quiescence.The authors also computed the phase of the RVs by fitting a sine function (the t=0 is taken as JD = 2454309.6147)to the RVs of various time intervals to understand the time-evolution of the RV phases and they found that the phases had some scatter during the time interval of the recorded data (JD = 2454874.8-2457563.9)but remained confined to a 90 • quadrant (see Figure 8 in Campbell- White et al. 2021).To compare the quiescent state RV phase with that of the 2022 March outburst state RV phase, we take the latest RV phases (JD = 2457548.9-2457563.9)from Campbell-White et al. ( 2021) and re-evaluate them at our startJD.A typical quiescence phase value for He I 5016 Å is 42 • ± 7 • .RV phases for other lines are tabulated in Table 2.These values are consistent with the emission line RV phase of Kóspál et al. (2014).For a few spectral lines, Campbell-White et al. ( 2021) also gave RV phases for the time period JD = 2458634.9-2458644.0,which translates to FeII 4924 (152 • ± 60 • ), Fe II 5018 (181 • ± 28 • ) and HeII 4646 (-29 • ± 110 • ), at our startJD.These phase values differ significantly among various spectral lines in contrast to the consistency among vari-ous spectral line phases for JD = 2457548.9-2457563.9and hence we considered the RV phase values for timeinterval JD = 2457548.9-2457563.9as pre-outburst RV phases.All these quiescent state RV phases are plotted by grey squares in Figure 2.This suggests that RV also showed a phase change during the 2022 March outburst, consistent with the photometric phase shift.

Mass accretion rate
Our multi-epoch high-resolution HRS-SALT spectra of EX Lupi are used to estimate the mass accretion rate evolution during-and post-outburst.The evolution of the spectral line over the outburst is shown by plotting the Ca II infrared triplet (hereafter, Ca II IRT) line over the ASAS-SN g-band lightcurve in figure 3. The Ca II IRT lines are plotted on the epochs of observed spectra.The broad component (BC) of Ca II IRT lines faded in tandem with the decay of the outburst (see Figure 3).The mass accretion rate is calculated by estimating the line luminosities of good accretion estimator lines (as marked in Alcalá et al. 2017).They are Ca II (K) (3933.66Å), H8 (3889.049Å), H β (4861.325Å), H γ (4340.464Å), H δ (4101.734Å), He I (4026.191Å), He I (4471.48Å), He I (5875.621Å), He I (6678.151Å) and He I (7065.19Å).A spectral cutout is prepared by manually selecting a wavelength window around every spectral line to avoid any nearby lines.The local continuum of the spectral line is estimated by fitting a straight line on the 2 Å window at both edges of the spectral cutout.Then, the local continuum subtracted spectral line is integrated.The integrated flux is scaled by 4πd 2 * to get the line luminosity at the stellar surface, where d * is the distance to EX Lupi.The line luminosity is converted to accretion luminosity following the relations developed in Alcalá et al. (2017).Assuming that the energy of the infalling matter is completely converted into radiation, the accretion luminosity is converted into a mass accretion rate following the Equation 1.
where Ṁ and L acc are mass accretion rate and accretion luminosity.R * and M * are the radius and mass of EX Lupi and G is the gravitational constant.R in is inner disk radius which is taken to be 5R * .The estimated mass accretion rates from various spectral lines are shown in Table 3.The mass accretion rates for all the aforementioned good accretion tracer lines are similar except for Ca II (K).So, for calculations concerning accretion rates, the average accretion rate for all lines except Ca II (K) will be used.The mass accretion rate estimated from various lines has a spread of about a factor of 2-3 (see Table 3).The average accretion rates on JD = 2459654.4899(around the peak of the outburst) and 2459736.539(almost around quiescence) are 3.2±0.3× 10 −8 M ⊙ /yr and 1.0±0.1 × 10 −8 M ⊙ /yr, respectively.An average from the last three epochs gives a quiescence mass accretion rate of 7.7±0.4× 10 −9 M ⊙ /yr that indicates a factor of 4 increase in mass accretion rate during the outburst.Wang et al. (2023) estimated a mean mass accretion rate of 1.74 × 10 −8 M ⊙ /yr during the peak of the outburst and 3.33 × 10 −9 M ⊙ /yr after the completion of the outburst.The authors used U -band excess flux over the photospheric level to estimate the mass accretion rate.The authors used a different extinction of A V =0.1 mag, and this caused a small difference in mass accretion rates estimates from that of ours.Cruz-Sáenz de Miera et al. ( 2023) performed slab modelling and estimated a factor of seven change in the mass accretion rate from quiescence to the peak of the outburst which is consistent with that of our result and from that of Wang et al. (2023).
From the detailed analysis of the long-term multiband photometry and multi-epoch high resolution spectra, the main observational conclusions are enumerated below.

O1:
The RV phase change was not observed among the spectra from pre-and post-outburst of the year 2008.
O2: Lightcurve and RV phases changed during the 2022 March outburst in multiband data.
O4: While the multi-band lightcurves had similar phase values during pre-outburst, they all started showing different phase values among themselves during the outburst.
O5: U, B, and V band phases didn't revert after the completion of the outburst, while R and I did revert immediately.
O6: Phase decreased differentially in optical and infrared lightcurves but not in U and B bands.2021) for the time-interval JD = 2457548.9-2457563.9.Blue pentagon is from Kóspál et al. (2014), for emission lines.The azimuthal position of the markers corresponds to the RV phase (at our startJD), and their radial position corresponds to the latitude inferred from the respective RV amplitude (see Section 4.9).The much-deviated orange-yellow color circle is He II 4686 Å.A grey curve of 112 • is plotted anti-clockwise from the phase of the blue pentagon, corresponding to the phase change estimated from the lightcurves.The same region is shaded with a light red color.
The following sections shall scrutinize these observational results.

DISCUSSIONS AND INTERPRETATION
To aid the discussion, we start by defining coordinate systems to characterize the phenomenon of phase change.

Coordinate System
We define two coordinate systems, one fixed on the rotating star ((X,Y,Z) ≡ (r,θ,φ), called fixed frame) and the other one fixed with respect to the hotspot , variable frame).Z and Z ′ axes are aligned to the rotation axis of the star, and the X-Y and X ′ -Y ′ planes coincide.The angle between X and X ′ axes is the phase of the lightcurve at any given epoch, which is ∼70 • at our startJD (in V -band, see the first diamond grey phase point in Figure 1).The physical interpretation of the RV phases and lightcurve phases are different.The same value of RV and lightcurve phases do not indicate the same location of the hotspot on the stellar surface.The photometric phase is maximum (90 • of the sine fit) when the hotspot is facing us (90 • in Figure 4) while the RV phase is maximum (90 • of the sine fit to RVs) when the hotspot is at the tangential right side point (180 • in Figure 4).Redshifted RV is positive in our convention.Thus, the hotspot location is given by the photometric lightcurve phase value, or equivalently, the RV phase + 90 • .

Detection of the azimuthal movement of hotspot via phase measurement
The phase of a lightcurve depends upon the location of the hotspot with respect to the observer's line of sight.An increase in phase implies that the hotspot is observed at an azimuthal position larger than 90 • , at a time when it is expected to be at 90 • .Since the rotation period of EX Lupi is found to be consistent with 7.417 days (see Appendix A), the location of the hotspot should move forward over the stellar surface to explain the phase change in the lightcurve.In terms of the coordinate system, (X ′ ,Y ′ ,Z ′ ) frame rotated anticlockwise with respect to (X,Y,Z) frame when viewed from the top (i.e., from above in Figure 4, this rotation direction is a convention used in this paper).
The hotspot is an observable component of matter carried by the funnel streams.The freely falling accreting matter along the accretion funnel defines the location of the hotspot over the stellar surface.A change in the hotspot location over the stellar surface indicates a change in the accretion funnel geometry, which depends on the relative rotation of the magnetosphere and the inner disk and also on the stellar magnetic field configuration.For EX Lupi, Sicilia-Aguilar et al. ( 2015) favour a dipolar magnetic field configuration, since it is a low mass (Gregory et al. 2014), fully convective star.The magnetic field configurations of typical YSOs are observationally found to be dipolar and also multipolar (e.g., Johnstone et al. 2014).Owing to the unavailability of the magnetic field configuration of EX Lupi, the discussions concerning the accretion funnel geometry shall assume a dipolar magnetic field configuration.
In the subsection below, we propose a hypothesis to interpret the azimuthal movement of the hotspot over the stellar surface.

Hotspot's movement due to shifting of the accretion funnel
The accretion onto EX Lupi, with dipolar magnetic field topology, is expected to be dominated by two accretion funnels which dump the infalling matter near magnetic pole, at high latitudes.The hotspot's azimuthal position on the stellar surface depends upon the inner disk radius (which is r m ) and r c , i.e., effectively on the relative angular velocity between the inner disk and the star.It can be understood as follows: the Keplerian angular velocity falls as ∝ r −3/2 .The inner disk radius decreases when the accretion rate increases (see Equation 3), and the infalling matter will have higher angular velocity than the prior accretion state.Since the star and its magnetosphere rotate slower, they will start lagging behind the faster rotating infalling matter.Thus, the matter falls at a larger longitude, along the disk rotation, when the accretion rate increases.This creates the hotspot at a larger longitude along the disk rotation.Observationally, it would appear as a phase change (increase) in the periodic lightcurve and RV.
As the hotspot moves to a larger longitude over the stellar surface, it would result in an azimuthal forward shifting of the dipolar accretion funnel.Hence, we name this phenomenon 'shifting of the accretion funnel'.Kulkarni & Romanova (2013) performed 3D MHD simulations of magnetospheric accretion and combined a range of parameters to derive relations governing the shape and size of the hotspot.We used Equation 6from Kulkarni & Romanova (2013) to estimate the azimuthal location of the hotspot at various instances of the outburst in EX Lupi: where r m , the magnetospheric radius is defined by the balance of disk ram pressure and magnetic field pressure.Since the azimuthal position of the hotspot depends upon the r m (Equation 2), we estimate r m for EX Lupi at the epochs we have the accretion rates from our HRS-SALT by equating the disk ram pressure to the magnetic field pressure as follows where r m and R * are pre-defined parameters.B 3 is the magnetic field strength in units of kG.Ṁ−8 is mass accretion rate in units of 10 −8 M ⊙ /yr.M 0.5 and R 2 are stellar mass and radius in units of 0.5M ⊙ and 2R ⊙ , respectively (Bouvier et al. 2007).k is a factor for correcting the spherical to disk-mediated accretion, where k ≈ 0.5 (e.g., Long et al. 2005) 4 .The values of the stellar parameters are described in Section 1.The mass accretion rates used in Equation 3 are estimated from the spectral lines, for 16 epochs (see Section 3.4 and Table 3).Substituting the value of r m /R * in Equation 2, we calculated the model predicted azimuthal positions of the hotspot on 16 epochs.

Hotspot location during the rise of the outburst
We estimated the predicted relative change in the hotspot's azimuthal position by subtracting the model predicted quiescent state azimuthal position from that of the model predicted outburst state position.Using the accretion rate from our first epoch spectra (JD = 2459654.4899),we estimated the hotspot's azimuthal position predicted by the model in the outburst state as 74 • .3± 6 • .3, as this epoch is close to the outburst peak (see Figure 3).To estimate the quiescent state azimuthal position of the hotspot, we averaged the hotspot's azimuthal position inferred from the last three epoch spectra (-43 • .2± 5 • .5)as they correspond close to the quiescence (JD = 2459732.5452,2459735.5317,2459736.539).Thus, the calculated change in the hotspot's azimuthal position from quiescence to the peak of the outburst is 117 • .5 ± 8 • .3.This is consistent with observations of phase increase by 112 • (see Section 3.1 and Figures 1 and 2).
We would like to note that, due to the periodic variability in accretion rate towards the tail end of the outburst decay, the choice of how we calculate the accretion rate of the quiescent state slightly alters this calculation.For instance, if we choose the azimuthal position estimate from our last spectroscopic observation as the quiescence hotspot location, then the phase change is 88 • .9± 10 • .6.Due to this uncertainty of the true quiescence accretion rate, the above-estimated phase change (using equation 2) would have a difference of 15 • -30 • , depending upon the choice of the quiescence accretion rate.

Hotspot location during the decay of the outburst
In order to understand the evolution of the hotspot location during the outburst decay, a straight line is fitted onto the model predicted azimuthal locations (ϕ) estimated from Equation 2. This model predicted phase decay curve is plotted in Figure 1 as a black dashed line.For visual clarity, the first epoch of the predicted phase decay line is shifted to match the closest AAVSO-U band phase.One can see in Figure 1 that the AAVSO-U band phases did not revert as per the dashed line while R and I bands followed the model predicted decay curve.
Our hypothesis of 'shifting of the accretion funnel' predicts a phase variation in tandem with the accretion rate variation following Equation 2. Hence, this hypothesis could explain O1: The increased phase could have reverted back after the year 2008 outburst was over.The first quiescent state spectrum was taken nearly five months after the completion of the 2008 outburst (Kóspál et al. 2014), and thus, the RVs estimated among the pre-and post-outburst states would have missed any phase change.This hypothesis also explains the O2 and O3: hotspot's longitude position changed during the outburst in the forward direction of disk rotation.However, this hypothesis fails to explain O5 and O6 completely: the U, B, and V band phases did not revert back after the completion of the outburst.
To explain O5, we hypothesize a thicker and warmer disk post-outburst.Outbursts are expected to heat up inner parts of the disk.The detections of 2.29 micron CO band heads in emission in some YSOs are believed to be an indication of this phenomenon (Lorenzetti et al. 2009).There has also been evidence of outbursts heating up outer disks and annealing silicates (in EX Lupi) ( Ábrahám et al. 2009), as well as moving the water snowline outward (in V883 Ori, albeit this is much more powerful FUor) (Cieza et al. 2016).The heated inner disk will be puffed up, and the disk will also start receding from the star with a decrease in the accretion rate.However, in the geometrically thicker disk, the funnel matter flows from higher altitudes of the disk, where the lifting force is smaller.This may help to support a high azimuthal shift of the hotspot for a longer time, compared with the case of a thinner disk.The funnel would carry most of the accreted mass and hence the corresponding hotspot will emit mostly in shorter wavelengths.This could explain a delayed phase decay in the AAVSO-U, B, and V bands.In this scenario, the phases will decay on the cooling timescale of the heated inner disk.
The observed phenomenon of persistent phase shift after the outburst (O5) is one of the most difficult observations to explain.The above scenario of the disk heating and thickening is one of the possibilities.In Appendix C we discuss a few other hypotheses that may explain the persistent phase shift in U, B and V bands after the outburst.

3D MHD simulation of the accretion onto EX Lupi
We simulated the accretion onto a star with parameters of EX Lupi to understand the phase shifts of the hot spots.We used our earlier developed 3D MHD Godunov-type code with a cubed sphere grid, which has been developed earlier and used to study the magnetospheric accretion (e.g., Kulkarni & Romanova 2013;Blinova et al. 2016).Here, we used this code to model funnel streams and spots in EX Lupi.We developed a model of a star with parameters of EX Lupi as discussed in Section 1: mass and radius of the star, 0.6M ⊙ , 1.6R ⊙ , its period: 7.417 days and corresponding corotation radius: 8.44R * .The tilt of the magnetospheric axis Θ = 13 • and the stellar magnetic field B = 3 kG.
The outburst in accretion rate is caused by some processes that are not well understood and which we cannot presently model numerically due to the unknown physics of the process and its possible complexity.Instead, we modelled parts of the process corresponding to different moments of time in the outburst.The result of simulations (the location and the shape of spots) depends on several dimensionless parameters: the truncation (magnetospheric) radius r m /R * , the corotation radius r c /R * (or, ω s ), and the tilt of the magnetosphere Θ. Earlier, Kulkarni & Romanova (2013) investigated shapes and phases of spots in models at a wide range of these parameters and derived different dependencies and correlations.Typically, they combined results for a wide range of magnetospheric tilts, taking, e.g., values in the range from Θ = 2 • − 5 • up to Θ = 20 • − 30 • .Here, we recalculated the model using a tilt of the magnetosphere Θ = 13 • corresponding to EX Lupi and at several values of r m /R * which may correspond to different stages of EX Lupi outburst.We also used twice as fine a numerical grid compared with Kulkarni & Romanova (2013) to ensure even more accurate simulations.To speed up the simulation, the star is placed inside the inner boundary (see Figure 5) such that the inner boundary of the simulation region is 2R * .
Initially, we placed the inner disk away from the star at the corotation radius.This provided more precise initial equilibrium (e.g., Romanova et al. 2003).Subsequently, the disk matter moved towards the star and started accreting.This led to the gradual rise of the accretion rate shown in the bottom left panel of Fig. 5.
Simulations show that when the disk matter starts approaching the star, it is stopped by the magnetosphere at the distance of r m ≈ 6 − 7R * and formed two funnel streams (see 3D views of matter flow in Figure 5; second row from the top, leftmost panel at t = 20.At this time, the magnetosphere is relatively large, and this moment may correspond to the state which is close to the pre-outburst state.The top left panel shows that the hotspot shifted only slightly from the µ − Ω plane (which is the plane of the tilt of the magnetic axis).Subsequently, at t = 33, the accretion rate increased, the inner disk moved closer to the star, and the spot was shifted by a larger angle, ∼ 35 • (see 3D plot and spot at t = 33 ).Later, at t = 40, the accretion rate increased, and the inner disk moved the base of the funnel stream to a larger angle such that the spot was shifted by ∼ 85 • − 100 • .
At moments of enhanced accretion rate at t > 50), we observed the transition from stable to the mildly unstable regime where both funnels and tongues were observed (see the rightmost panels in the same figure at t = 58.).Instability started to appear when the magnetospheric radius decreased to values of r m ≈ 5.5-5.8 and the fastness parameter ω s = (r m /r c ) 3/2 ≈ 0.53-0.57which corresponds to mildly unstable regime (Blinova et al. 2016) 5 .In this regime, matter flows through different channels however most energetic part of the spot with the highest kinetic energy of the flow was determined by the modified funnel stream and was shifted by ∼ 100 • .The 3rd and 4th rows from the top show the slices of density distribution in the µ−Ω plane and in the perpendicular plane (where the flow is expected in case of strong phase shift).One can see that at t = 20 (small shift), the funnel streams are stronger in the µ−Ω plane, while at t = 58 (strong shift), streams are stronger in the perpendicular plane.
We calculated the lightcurve for observer located at the inclination angle 32 • (see bottom right planel of Figure 5).Here, we suggested that all kinetic energy of matter falling to the surface of the star is re-radiated isotropically.The rotation of the spot leads to variation of the amount of radiation towards the observer, and the lightcurve shows maxima and minima.We chose one period of stellar rotation (shown between dashed vertical lines in Figure 5 for lightcurve) and show the location and shape of spots at different moments in time and different phases of stellar rotation (see Figure 6).One can see that at the inclination angle of 32 • the observer will always see the spot.The maximum in the lightcurve corresponds to the time t = 36.8when the spot is located closer to the center of the star, facing us, and the two minima correspond to the times t = 33.1 and t = 40.5 when the spot is closer to the edge of the star.

Chromatic phase shift during the outburst period
The analysis of Section 3.1 showed that the TESS and nearby epoch ASAS-SN phases were consistent during the quiescence of pre 2022 March outburst.Further at the onset of the outburst, g ′ , r ′ , i ′ and z ′ band phases are all consistent within 10 • (see Section 3.1 and Figures 1 and 7).This implies that the azimuthal location is the same for the regions emitting across the wavebands.Therefore, the hotspot was azimuthally confined and nearly single-temperature on the EX Lupi during the pre-outburst quiescence.However, the AAVSO U, B, V, R and I-band lightcurves observed during the peak of the outburst (see Section 3.1) showed finite differences among the hotspot longitude position inferred for U, B, V, R and I-bands (see O4 and Figure 7).The inferred hotspot longitude positions for wavebands among AAVSO U, B, V, R and I have a decreasing trend.U -band lightcurve peaks first followed by B, V, R and I-band lightcurves, sequentially.The maximum phase difference between U to I band is 36 • .The longitude inferred from the CTIO B-band lightcurve is also larger than that of the CTIO V -band lightcurve by 10 • .7.Similarly, Espaillat et al. ( 2021) also found phase differences among UV and optical lightcurves of the star GM Aurigae where the peak of the UV lightcurve precedes that of the optical lightcurve.The authors suggested it to be a consequence of azimuthal asymmetry in the hotspot.For EX Lupi, this indicates a temperature gradient in the hotspot where the front end, along the rotation, is hottest, subsequently followed by cooler and cooler regions.Spread in multiband phases during the outburst could also refer to a chromatic azimuthal enlargement of the hotspot.This could further be interpreted as a broadening of the accretion funnel to channel a large amount of mass during the outburst.
This chromatic phase difference could also be a consequence of the appearance of many azimuthally nearby  funnels during the outburst.Blinova et al. (2016) showed that the accretion becomes unstable when ω s ≤ 0.6 for magnetic obliquity of 5 • (close to that of EX Lupi).Under this condition, many equatorial mattercarrying tongues penetrate the magnetosphere and form chaotic spots near the magnetic equator.But under this condition, the lightcurves are expected to be irregular.The lightcurves become periodic with smaller periods at other extremes when w s ≤ 0.45.We did a periodogram with the declining lightcurves after the outburst, and we didn't find any statistically significant period value smaller than the 7.417 days.We also calculated the w s at the epochs of our spectra and found that the accretion lies in the mildly unstable regime for the initial epochs and then turns to the stable regime of accretion in the later epochs.Upon varying the value of 'k' in between 0.5 and 0.7 in Equation 3, EX Lupi stayed around mild unstable (appearance of a few equatorial tongues with profound dipolar funnel) to the stable regime (accretion is channelled by dipolar funnel).
Not only did the inferred longitude from multiband lightcurves become different during the outburst, but in the subsequent post-outburst state, the inferred longitude of the wavebands with longer effective wavelengths decreased much faster than that of the wavebands with shorter effective wavelengths.

Evolution of chromatic phase shift post-outburst
During the decay of the outburst, the phases decreased differentially among U, B, V, R, and I-bands (see Figure 7).The slopes of the multiband phase decays are-U :-0.153,B :-0.263, V :-0.504,R :-0.758,I :-0.824 (degree/day).As we discussed in Section 4.3.2,R and I band lightcurves followed the model predicted phase decay but the warmer and thicker inner disk could have delayed/slowed the phase decay among the shorter wavelengths.Interestingly, EX Lupi maintained chromatic phase differences even after the outburst was over.We found chromatic phase differences in our LCRO lightcurves taken from 2023 February to 2023 May.We suspect that the occasional small bursts and variability seen in the ASAS-SN lightcurve during the post-2022 March outburst kept the inner disk hot and thick.A hot disk could maintain the larger amount of matter flow along the 'shifted funnel' and support the hotspot at a larger azimuth that corresponds to emission at a shorter wavelength.However, a smaller amount of matter might have been flown by the azimuthally lagging funnels showing the phase decay at the longer wavelengths.These chromatic phase differences showed longterm sustainability of a temperature gradient in the hotspot as it was observed among LCRO g ′ , r ′ and i ′ band lightcurves, one year after the outburst.It explains the persistent chromatic phase differences in the post-outburst state.If this hypothesis is true, it suggests that the chromatic phases should return to quiescent values in the relaxation timescale of the puffed-up hot inner disk.To verify this, we recommend a further high-cadence multiband observation of EX Lupi during its quiescence, akin to that of 2019.Further, to visualize the phase difference among g ′ , r ′ and i ′ bands, we introduce Lissajous figures for g ′ -r ′ , r ′ -i ′ and g ′ -i ′ band lightcurves.The lightcurves are normalized before plotting the respective Lissajous figures.The data were averaged for one-day observations before fitting the sine curve for the plots in Figure 8. Two sine curves with a δ phase difference, when plotted onto the X-and Y-axes, produce an ellipse.The ratio of the semi-minor to semi-major axis varies between 0 (δ = 0 • , the ellipse becomes a straight line) and 1 (δ = 90 • , the ellipse becomes a circle).The ratio of the semi-minor to semi-major axis is maximum for g ′ -i ′ than the other two combinations, indicating a larger phase difference between g ′ and i ′ bands (see Figure 1) than the phase differences between g ′ -r ′ and r ′ -i ′ .
During the final stages of this manuscript preparation, an independent study by Wang et al. (2023) also reported the chromatic phase shift in EX Lupi postoutburst.

Larger phases after one year of the outburst: LCRO and ASAS-SN lightcurves
The phase values from the multiband LCRO and ASAS-SN lightcurves after many months of the outburst are essential to understanding the long-timescale evolution of the system.Our LCRO observations in g ′ , r ′ and i ′ taken one year after the 2022 March outburst show the light curve phases were still around 200 • , 180 • and 150 • in g ′ , r ′ and i ′ band lightcurves respectively (see Figures 1 and 7).These LCRO multiband phases are comparable to the respective band phase values of the AAVSO observations in September 2022 (g ′ ∼ V , r ′ ∼ R and i ′ ∼ I; see Figure 1).
These elevated phases for more than a year after the March 2022 outburst are in sharp contrast to the major outburst in 2008.After the 2008 outburst, there was no indication of an elevated phase after about 5 months of the outburst (Kóspál et al. 2014).If a similar phase change had happened in the 2008 outburst, it relaxed to the quiescent phase much faster.ASAS-SN g-band lightcurve shows the flux from the EX Lupi did not fully return to pre-outburst level after the March 2022 outburst (See lower panel of Figure 9).Between March 2022 and August 2023, there have been numerous low amplitude brightening.Sometimes brightening as high as ∼1.2 mag in g-band during February and March 2023.If our hypothesis is true, these higher levels of accretion activity might have kept the inner disk warm and thicker to support the elevated azimuth of the accretion funnel.

Phase decay after one and half year of the outburst: ASAS-SN view
The ASAS-SN g-band lightcurve of EX Lupi covering from 2023 January 12 (JD = 2459956.8644)to 2023 September 06 (JD = 2460194.3491)showed the larger phases.The phase of the lightcurve remained similar to the value calculated just after the outburst (see Figures 1, 7 and 9).But, the latest 50 days of the ASAS-SN gband lightcurve from 2023 July 16 (JD = 2460140.7206)to 2023 September 06 (JD = 2460194.3491)shows the phase to have decreased to the pre-outburst value (see Figure 9).This latest phase value is 103 • ± 6 • while the pre-outburst phase values estimated for the manually selected clean lightcurve cutouts (see Section 3.1), epochs close to the TESS sector 12 and sector 39 observations, are 102 • ± 3 • and 111 • ± 3 • , respectively.The bottom panel of the Figure 9 shows the brightness level variation of the EX Lupi in the pre outburst, during the outburst and post-outburst state.EX Lupi did not completely return to the pre-outburst brightness level for more than a year of the post-outburst state.The ASAS-SN g-band lightcurve shows the phase of EX Lupi's lightcurve returned to the pre-outburst level in tandem with the flux (see Figure 9) during September 2023.

Change in the radial velocity amplitudes: Possible latitudinal movement of the hotspot
A sine curve model fit the phase folded RVs estimates the phase and amplitude of the RV for the specific spectral line (see Section 3.2).The fitted parameters are tabulated in Table 2.The RV amplitudes from previous existing literature are also tabulated in the same table.In Table 2, columns 8 and 10 contain the ratio of the RV amplitudes from our estimation to that of Sicilia-Aguilar et al. (2015) and Campbell- White et al. (2021) respectively.The ratios are larger than 2 in column 8 and mostly larger than 1.5 in column 10.Thus, there is a significant RV amplitude enhancement during the outburst.For a hotspot on a rotating sphere with a period P and radius R, the RV amplitude can be modelled as: where i and Λ are the inclination angle of the rotation axis with respect to the line of sight and latitude of the hotspot, respectively.Considering Equation 4, an increase in the RV amplitude can be a consequence of four possible causes: 1) the hot region where lines form is at a higher altitude from the stellar surface, along the accretion funnel.2) The accreted matter is channeled through funnels to the magnetic pole, along with the equatorial tongues dumping matter near the equator.3) The hotspot is more longitudinally confined during the outburst than during the quiescence.4) The hotspot's latitude decreased.We scrutinize each of the possibilities below.
Can 'hot region' at a higher altitude explain the increased RV amplitude?
Taking He I 5875 Å as a reference, the RV amplitude increased by a factor of 2.38 (see column 8 in Table 2).
Considering other parameters are the same in Equation 4, the hot region where spectral lines form would have to be at a radial distance of 2.38R * along the accretion funnel (if the quiescent state spectral lines formed at the stellar surface (R * )).In contradiction to this hypothesis, the increased accretion rate would decrease the size of the pre-and post-shock region by the increase of infalling matter's ram pressure during the outburst state.We calculated the size of the pre-shock of the accretion-shocked region to be 0.0025R * << R * (Equation 10 in Hartmann et al. 2016); where we used mass accretion rate of 2 × 10 −8 M ⊙ /yr and the hotspot's fractional area of 5 percent of the total stellar surface.Thus, the explanation of the RV amplitude increase by the spectral line formation at higher altitudes along the funnel is highly unlikely.

Can an unstable accretion regime explain the increased RV amplitude?
The Rayleigh-Taylor instability could develop at the magnetosphere-inner disk boundary during the states of high accretion rate, which produces equatorial mattercarrying tongues and dumps matter near the stellar equator (Blinova et al. 2016).These equatorial tongues can co-exist with the dipolar accretion funnel.If equatorial tongues carry a substantial amount of accreted mass, it could result in hotspots near the equator.Thus, the effective RV amplitude could be an amalgamation of the lines forming at the near-polar and near-equator hotspots.It will effectively increase the RV amplitude.Our analysis in Section 4.4 found that the equatorial tongues could have developed on the initial few epochs of spectra observation around the outburst peaks.However, our RV amplitude is an aggregate estimate of 16 epochs over 80 days, and those one or two epochs would not alter the estimates much.Also, the most energetic flow of matter is determined by dipolar funnel on these epochs (see top right panel in Figure 5) but we found consistent RV amplitude increase in all the high excitation energy lines, like He I, as well as in low excitation energy lines, like Fe I (see Table 2).Thus, an unstable accretion regime could not explain the increased RV amplitude during the small outburst of 2022 March.However, this phenomenon could produce apparent RV amplitude variations in stars under a highly unstable accretion regime.
Can a varying size of the hotspot explain the increased RV amplitude?A change in hotspot size will change the average RV amplitude.For instance, if the hotspot is spread uniformly across all azimuth, due to symmetry, any redshift from one side of the star will be cancelled by the blueshift from the other side of the star, effectively showing no RV changes.But, if the hotspot is confined within a small azimuthal region, the spectral lines produced in it will show redshift and blueshift as it rotates with the star.To understand the relation between the size of the hotspot and the RV amplitude, we simulated hotspots of various sizes on the stellar surface with three parameters: azimuthal size (∆φ), latitudinal size (∆Λ) and latitude location (Λ 0 ).The hotspot is sampled onto a grid (latitude and azimuth angle) with a grid size of 0.5 • .Each grid point of the hotspot is assigned an RV (A RV sin(φ)), where φ is the longitude of the grid.We approximated the RV of any spectral line forming in this hotspot as a uniformly weighted mean of RVs calculated over the hotspot grid.To simulate the peak RV amplitude of the RV time series, we placed the hotspot at the tangent of the stellar surface to our line of sight (at 180 • in Figure 4).By placing the hotspot at the tangential point, only half of the hotspot's azimuthal size contributed to the RV.The azimuthal and latitude sizes of the hotspot were varied between 10 • -100 • and 10 • -30 • respectively.And we further varied the hotspot's latitude location from 55 • to 80 • (see Figure 10).Hotspot's latitudinal size is limited by the latitude location.Each panel in Figure 10 shows RV amplitude variation as a function of one parameter (X-axis label) while keeping the other two parameters free.
We found that the increase in the azimuthal size decreases the RV amplitude (right panel in Figure 10).The RV amplitude of around 1.0 km/s at ∆φ=10 • (see right panel of Figure 10) decreases as azimuthal size is increased to 100 • .Similarly, an increase in latitudinal size decreases the RV amplitude but is not apparent in the middle panel of Figure 10 because of a smaller range of latitude sizes.But, one can see in the left panel that decreasing the latitude position of the hotspot, increases the RV amplitude.Taking He I 5875 Å as a reference (RV amplitude = 2.45 km/s from our HRS spectra is plotted as a dotted horizontal line in Figure 10), we found that the hotspot latitude location below ∼65 • is able to reproduce the observed RV amplitudes in our analysis.The RV amplitude of ∼1 km/s, as shown in Sicilia-Aguilar et al. (2015), can be reproduced with the hotspot's latitude location of ∼80 • along with the latitude and longitude sizes of 10 • -15 • and ∼20 • , respectively (see Figure 10).A larger inclination angle of 45 • (Goto et al. 2011) does not change our qualitative understanding of hotspot's latitude change by 10 • -15 • .
During an outburst, an accretion funnel is expected to become broader and enlarged to transport more mass.Hence, the outburst should have increased the hotspot size, effectively decreasing the RV amplitude.Thus, the scenario of increasing RV amplitude due to tighter con-finement of the accretion funnel is also highly unlikely.
Can a hotspot at a lower latitude explain the increased RV amplitude?
The simulation in the previous section showed that the decrease in the hotspot's latitude could explain the increment of the RV amplitude.Considering a decrease in the hotspot's latitude while keeping the other two parameters invariant during the outburst, we calculated the hotspot's latitude using equation 4. The HRS-SALT spectra covered ∼ 80 days, and hence our calculated hotspot latitude in Table 2 is an aggregate of all epochs.The calculated hotspot latitude from HRS-SALT spectra and the prior available RV amplitudes are tabulated in Table 2.The average latitude differences from that of Sicilia-Aguilar et al. (2015) and Campbell- White et al. (2021) to our HRS-SALT data are 16 • and 6 • , respectively.Latitude estimates from our Mg I 5167 Å and Si II 6347 Å lines showed a larger values than that in Campbell- White et al. (2021), removing these two points gives an average latitude decrease of 9 • during the 2022 March outburst.
The dipolar magnetic field configuration also predicts a decrease in the hotspot's latitude during the state of increased accretion rate.The dipolar field lines follow a parametric equation r/ sin 2 (θ) = const, where θ is measured from the magnetic axis and r is radial distance (Kulkarni & Romanova 2013).The dipolar magnetic field structure is such that the inner closed field lines converge on the stellar surface at lower latitudes.The increased accretion would cause the inner disk to come closer to the stellar surface (see Equation 3).The decreased inner disk radius would let the infalling disk matter to follow the inner magnetic field lines, eventually dumping matter at the lower latitudes.Thus, a decrease in latitude during an accretion outburst is consistent with the dipolar accretion funnel.Assuming a dipolar magnetic field topology, EX Lupi's hotspot latitude is calculated using (Equation 3 in Kulkarni & Romanova 2013) where Θ = 13 • for EX Lupi.R * and r m are already defined.Since the disk would also compress the magnetic field lines, forcing them to converge at some higher latitude, the spot latitude inferred from Equation 5 is a lower limit.Using Equation 5, we calculated the lower limit of the spot latitude among the spectral lines listed in Table 3 to be 49  2021).Thus, a decrease in the hotspot's latitude explains the increase in the RV amplitude.However, a detailed simulation study is required to understand the variation of hotspot latitude during the state of variable accretion rate.
Thus, we conclude that the RV amplitude increment is a consequence of an effective decrease in the stellar hotspot's latitude.

A possible azimuthal asymmetry in the accretion
funnel: Blue shifted wing in the Ca II IRT lines During the outburst, the Ca II IRT lines are seemingly composed of three Gaussian components: a Narrow Component (NC), and two Broad Component (BC) Gaussians; one on the blue side and the other on the red side of the NC (see Figure 3).The chromospheric origin of the NC is suggested by Hamann & Persson (1992); Batalha et al. (1996).The broad component is believed to be formed in a highly turbulent, physically extended region outside the stellar surface (Hamann & Persson 1992;Batalha et al. 1996;Sicilia-Aguilar et al. 2015).There is a high correlation between the existence of broad components in Ca II IRT and the high accretion rate/optical continuum, as illustrated in Figure 3.We modelled the Ca II IRT lines with a sum of three Gaussians: two broad Gaussians each leftward and rightward of a central narrow and small Gaussian.On three epochs of JD = 2459682.4127, 2459683.4179 and 2459712.604 (30.2 and 29.2 days apart, ≈ an integer multiple of EX Lupi's rotation period, 7.417 days), the Ca II IRT lines showed the presence of a prominent blue shifted wing (see 3rd, 4th and 8th epoch spectra in Figures 3 and  11).We see similar structures in H-Balmer lines also  2).The colored and grey circles represent RV amplitude larger and smaller than 2.45 km/s, respectively.The X-axis is in degree, and Y-axis is in km/s.during these epochs.The velocity of the blue-shifted wings are -116.6km/s, -71.4 km/s and -81.4 km/s on JD = 2459682.4127,2459683.4179and 2459712.604,respectively.To visualize rotational phase-based dynamics, the Ca II IRT line profiles are plotted in a polar diagram (see Figure 11).The polar position of the Ca II IRT in the diagram is the hotspot's azimuthal location at the respective epochs of observation.Since the photometric phase from the AAVSO U -band lightcurve is nearly stable even after 5 months of the outburst, EX Lupi's surface is mapped onto the phase values using the AAVSO U -band phase during the outburst.The 3rd, 4th, and 8th epoch spectra (epochs, where the intense blue shifted wings are seen) are observed when the hotspot is approaching us and is almost tangential to our line of sight (0 • in the Figures 4 and 11, pictorially, the star is shown to be rotating anti-clockwise,).As we discussed in Section 4.4, the accretion is dominated by a dipolar funnel, and it approaches us with maximum velocity when the hotspot approaches us.Thus, the blue shifted emission components in spectral lines could imply a hot blob/material coming towards us with maximum projected RV at specific phases (around 0 • in Figure 11, when the funnel approaches us with maximum speed).The freefall timescale of material falling onto the stellar surface from the inner disk is only a few hours.However, the signature of a hot blob appeared even after 4 rotation cycles of the star and hence the blob cannot be expected to be fixed in the accretion funnel.One possibility could be that the mass accretion itself is clumpy.And, we see the extreme blue wings in Ca II IRT lines when the clump is falling through the accretion funnel and the funnel is approaching us at the highest projected velocity.The clumpiness being a random event could explain why not every Ca II IRT line has the same level of extreme blue-shifted wings.
Moreover, the blob has to be hot enough to produce an emission comparable to the rest of the Ca II IRT line emission regions.A detailed analysis is needed to understand the thermo-physical dynamics of the funnel, which is beyond the scope of this paper.We plan to understand the dynamics by doing spectral line analysis in our future work.However, we can estimate the distance of the hot blob from the stellar surface, assuming the measured velocity of the blob to be a sum of radial infall velocity and Keplerian rotation velocity.
Since the infalling matter doesn't exactly follow the Keplerian velocity but slows down, we calculated the lower limit of the radial distance by assuming the blob to be co-rotating with the star: where r is the radial distance of the hot blob.r m and i are pre-defined parameters.φ is the phase of  4.11.Accretion is not streamlined but clumpy: Photometric perspective During the outburst, on multiple nights, we conducted high-cadence multiband photometric observations spanning a few hours.These observations revealed the presence of sudden enhancement in brightness structures on the time scales of a few hours (see Figure 12 for a few examples).These few hour-timescale brightness enhancements follow the same path in color-color and color-magnitude diagrams as followed by the EX Lupi in its quiescent state to outburst state transition (see Figure 13).It suggests that these hour timescale events are related to the accretion process.Gullbring et al. (1996) also reported similar hourly time scale events for T Tauri star BP Tauri.It was understood to be a consequence of accretion of the mass clumps (Gullbring et al. 1996).Figure 12 shows that the time scales are not the same for each event, which may suggest the difference in the clumps' masses.We estimated the duration of the rise and fall of these events by visual inspection of the B-band lightcurve.In some of the cases, the estimated duration is a lower limit, as the whole event couldn't be observed.To estimate the asymmetry of the profiles, we measured the rise and fall times for each individual event.The rise and fall timescales show a linear relation with the slope and intercept of 0.75±0.13and 0.27±0.18,respectively (see Figure 14).The relation shows that the fall timescale of the events is longer than that of the rise timescale by a factor of ∼1.3.
We further estimated the masses of these clumps by converting the brightness increase to accretion luminosity and then to the total mass (see Appendix D).The estimated masses of the clumps are given in Table 4.The clumps' masses range from 0.1 to 19.2 × 10 −5 lunar mass (an equivalent of ∼ 300 to 57000 Halley's comet mass) Assuming a random distribution, 12 such events observed out of 42 high cadence photometric observations using TMMT implies that the clumps were ∼ 29% of the time.In contrast, out of 74 epochs of observations with LCRO one year after the outburst from January 2023 to May 2023, only a very few observations hint at such events (one such event is shown in the bottom right plot in Figure 12).The LCRO observations in 2023 are relatively sparse, which makes it hard to find and characterize the clumpy events.The other possibility could be that the clumpiness of the accretion is dependent upon the mass accretion rate.During the quiescent state mass accretion rate, the accretion could be less clumpy.The clumps could start appearing as random events as the accretion rate increases.In extremely large accretion rates, the falling clumps could become too frequent to be termed continuous mass-infall.We calculated the rotational phase of EX Lupi during the epochs of these clumpy events to search for any rotational-phase correlation.The clumpy accretion is observed at all rotational phases of EX Lupi.

SUMMARY AND CONCLUSIONS
The 2022 March outburst of EX Lupi was observed at an unprecedented high-cadence multiband photome-  shows the order number of the events.The solid black line is reddening vectors plotted at an interval of AV = 0.5 mag (Cardelli et al. 1989) 2. The accretion-driven 2022 March outburst in the EX Lupi caused a forward movement of the hotspot by 112 • ± 5 • .This forward movement is understood as a consequence of the high accretion rate.The increased accretion rate brings the disk closer to the stellar surface, thus increasing the Keplerian angular velocity of the inner disk material to be accreted onto the star.The slowly rotating star and its magnetosphere start to lag behind the faster-rotating infalling matter.The matter falls on the stellar surface at a larger longitude.Thus, the hotspot moves forward during the outburst.The forward movement of the hotspot is consistently observed in our high cadence, high resolution spectroscopic, and multiband photometric monitoring with HRS-SALT, LCRO, TMMT, and CTIO, along with the available archival high cadence multiband photometry from AAVSO and ASAS-SN.
3. The predicted accretion rate-dependent hotspot's azimuthal position is not observed in the outburst decay state.The red lightcurves (R and I) show the phase decaying to the pre-outburst state value, as per the prediction.But, the shorter wavelength waveband's (U, B, and V ) phases stay at higher values.One particular promising hypothesis is that the disk remained heated and thick for a longer period after the outburst, resulting in sustaining the higher longitude of the accretion hotspot.
4. Our analysis shows a predominantly confined single hotspot in the pre-outburst state gets spread azimuthally and covers a large azimuthal area during the outburst.The hotspot in the post-outburst state has an azimuthal temperature gradient.The front end of the hotspot is hotter, while the lagging end has a cooler temperature.
5. Our SALT-HRS spectra show an increase in the RV amplitudes during the outburst.We discuss in detail multiple hypotheses to explain this observation in Section 4.9.The decrease in the hotspot's latitude from ∼ 75 We recommend further follow-up of this touchstone star EX Lupi to understand the long-term evolution of EX Lupi's magnetosphere.Further simulations of the inner accretion disk and the magnetosphere need to be developed to understand the complex behaviours seen in these systems.data with respect to our HRS-SALT data, respectively.The detected spectral lines are cross-checked against the lines found during the earlier outburst by Sicilia-Aguilar et al. ( 2012) and the respective energies are taken from the NIST Atomic Spectra Database (Kramida et al. 2022).Fe II 4923.9216 5.823-3.9033.57±0.66150.2±11.1 3.37±0.4852.4±9.0 Table 3. Mass accretion rate estimated from the HRS-SALT spectra of EX Lupi for 16 epochs.The mass accretion rates are given in units of 10 −9 M⊙/yr.using the Lmfit-module (Newville et al. 2016).The phase-folded RV showed a periodic variation.The phase of the periodically varying RV is plotted in Figure 2. Though the RV inferred from the multi-Gaussian fit revealed a periodic nature, the line width of the NC Gaussian did not show any periodic modulation.The Ca II IRT lines (8498 Å, 8542 Å, and 8662 Å) are hard to model.They seem to be composed of three Gaussians: blue-shifted broad emission line, central NC, and red-shifted BC.We modelled the Ca II IRT lines by a sum of three Gaussians to estimate the velocity of the blue-shifted wings on the three epochs (see Section 4.10).The three epochs where a highly blue-shifted feature is seen in the Ca II IRT line are shown in Figure 11.Detailed line profile modelling and analysis are beyond the scope of this paper and are planned in our upcoming work.
In Figure 17, we show a sample of multi-Gaussian fit to the spectral lines.

C. ALTERNATE HYPOTHESES FOR PHASE CHANGE: ACTIVATION OF A NEW ACCRETION FUNNEL
The phases of the wavebands with shorter effective wavelengths (U and B) being almost constant and at larger longitudes could also be a consequence of the permanent reconstruction of some new accretion funnel at the observed longitude.During the rise of the outburst, it is possible that some global changes happened in the star-disk system, which led to the flow of infalling matter streams along new paths of lowest energy.The new funnels that form at new azimuthal positions can stay there for a significant time, depending upon the decay timescale of the changes.A completely new funnel can be activated at a different azimuthal location under two possible scenarios: 1) permanent reconstruction of the magnetic field and/or 2) a tilt in the inner disk caused by some cloudlet capture.
Permanent reconstruction of the magnetic field: If the magnetic field gets reconstructed during the outburst, then the phase of the lightcurve can change.This phase change would be completely random and will remain permanent owing to the new magnetic field configuration.The observations have previously shown a varying magnetic field configuration (Johns-Krull et al. 2022;Finociety et al. 2023), with suggestions of varying convection processes inside the star.No such varying magnetic field reconstruction has been inferred in EX Lupi for almost one and a half decades otherwise it would have appeared in terms of lightcurve phase change.It is highly unlikely that the intrinsic parameters governing the magnetic field structure could change during the 2022 March outburst but remain intact during the gigantic outburst of the year 2008.This hypothesis opens too many questions than solving the existing ones.

Tilt in the inner disk:
Another probable reason for the permanent change in the hotspot's azimuth could be the tilt of the inner disk rotation axis with respect to the pre-outburst inner disk rotation axis.A stream of clouds falling onto the circumstellar disk could cause a tilt in the disk.This tilt could be completely random depending upon the angular momentum, energy, impact parameter, etc. of the infalling stream.Thus, a phase change of the 112 • could be a random event rather than a systematic forward motion of the hotspot along the rotation as suggested by Kulkarni & Romanova (2013).One of the solutions to the misaligned short-period extrasolar planets is suggested to be a consequence of the tilt in the circumstellar disk caused by the capture of passing by a stream of cloud/gas (Thies et al. 2011). Dullemond et al. (2019) suggested the possibility of gravitational instability in the disk caused by the cloudlet capture, eventually leading to an FU Orionis-like event.The authors further suggested that the timescale of EX Lupi-like event is too small to be explained by a long-time event of cloudlet capture.Romanova et al. (2021) did 3D simulations on the accreting stars with a tilted inner disk with respect to the stellar rotation axis.The new tilted inner disk evolves on its viscous timescale after the initial tilt settlements.In this scenario, the phase will not return for a very long time.This hypothesis also fails to explain the differential phase evolution in AAVSO multiband data in the post-outburst state.

D. ACCRETION RATE ESTIMATION FROM PHOTOMETRY
We calculated the mass of the accretion clumps (see Section 4.11) by estimating the mass accretion rate from the photometry under simplistic assumptions.The AAVSO U, B, V, R and I-band lightcurves are converted into the flux unit using Bessell et al. (1998).Then the sets of AAVSO multiband photometry are chosen such that each band is observed apart not more than 0.05 a day.The fluxes are multiplied by 4πd 2 * .The extinction correction is done for A V =1.1 mag.Then the upper bound on the stellar flux is estimated as the minimum of the fluxes in each photometric band during the quiescence of the post-outburst.This estimated stellar flux is subtracted from the respective band lightcurves to get an estimate of the hotspot flux.Then, a Blackbody is fitted onto the hotspot flux using the Astropy module in Python.The temperature of the Blackbody fit is kept fixed at 10000 K following Hartmann et al. (2016), while the scale is kept variable with an initial value of R 2 * .The integrated luminosity of the best-fit blackbody is taken as the accretion luminosity, which is then converted to the mass accretion rate using Equation 1. Figure 18 shows the accretion rate estimated with the photometry along with the average accretion rate inferred from the spectroscopic lines.The accretion rate estimated from photometry is larger than that from spectroscopic emission lines by a factor of approximately 5 to 6 (see Figure 18).Figure 18 also shows the EX Lupi lightcurve during the outburst for a clear picture of the estimated increased accretion with the brightness enhancements during the outburst.
Similarly, we estimated the mass accretion rate during the 'clumpy events' observed in the TMMT B, V and I lightcurves.TMMT B, V , and I band lightcurves are converted into flux units and corrected for extinctions of A V =1.1 mag.The estimates of the stellar fluxes in the B, V and I bands are the same as that in the AAVSO lightcurves.Then, the mass accretion rates are integrated for the timescales of the clump to get an estimate of the clump's mass (see Table 4).The masses of the clumps range from 1 to 192 × 10 −6 Lunar mass.Murty Foundation, is a not-for-profit organization dedicated to the preservation and celebration of culture, science, and knowledge systems born out of India.The Murty Trust is headed by Mrs. Sudha Murty and Mr. Rohan Murty.

3. 3 .
RV phase change in between the 2008 outburst and the 2022 outburst Kóspál et al. (2014); Sicilia-Aguilar et al. (2015); Campbell-White et al. (2021) have presented RVs of EX lupi for around a decade starting from the year 2007.Kóspál et al. (2014) presented the absorption and narrow emission line RVs, covering 5 years across July 2007 and July 2012.This time period covers a major outburst of the year 2008.The authors estimated emission line RVs by the multi-Gaussian fittings of 133 lines identified by Sicilia-Aguilar et al. (

Figure 1 .
Figure 1.Top panel : The Longitude of the Hotspot estimated with respect to the fixed frame is plotted.Various colors and markers represent estimated hotspot longitudes from the multiband lightcurves and the RVs.Longitudes inferred from AAVSO lightcurves are shown for U and I-bands only.Longitude inferred from RV is plotted only for He I 5875 to make the figure less crowded; however, all the RVs are plotted in Figure 2. A straight black dashed line is the model-predicted decay rate of the phase, starting with the AAVSO-U band phase at the outburst.Bottom panel : ASAS-SN V and g-band lightcurves are plotted by blue and green dots respectively.The time period of 2022 March outburst is shaded by light red color.A more detailed post-outburst phases of multiple wavelengths are shown in Figure 7

O7:
The ASAS-SN g-band phase decreased to preoutburst value during the time interval of 2023 July 16 (JD = 2460140.7206)to 2023 September 06 (JD = 2460190.7206)O8: The RV amplitude of the emission lines increased during the outburst.

Figure 2 .
Figure 2. RV phases estimated for various spectral lines are plotted.The circles correspond to the RV phases from HRS-SALT spectra of EX Lupi during the 2022 Outburst.The grey squares are from Campbell-White et al. (2021) for the time-interval JD = 2457548.9-2457563.9.Blue pentagon is fromKóspál et al. (2014), for emission lines.The azimuthal position of the markers corresponds to the RV phase (at our startJD), and their radial position corresponds to the latitude inferred from the respective RV amplitude (see Section 4.9).The much-deviated orange-yellow color circle is He II 4686 Å.A grey curve of 112 • is plotted anti-clockwise from the phase of the blue pentagon, corresponding to the phase change estimated from the lightcurves.The same region is shaded with a light red color.

Figure 3 .
Figure 3.The Ca II IRT spectra are plotted over ASAS-SN g-band lightcurve during the outburst.Each spectrum is plotted on the respective epoch of observation, shown by a vertical dashed grey line.Ca II lines at 8498.02 Å, 8542.09Å and 8662.14Å are shown by red, green and blue colors, respectively.The broad component disappeared from the Ca II IRT lines with the decay of the excess continuum from the outburst, around after 8th epoch.The red dotted vertical lines are the epochs of observation of hour timescale photometric enhancement in the TMMT lightcurves (see Section 4.11) .

Figure 4 .
Figure 4. Azimuthal angle sign convention to be followed in this work.0 • and 180 • are the tangential left and right points, respectively, while 90 • is facing us.The two arrows corresponding to Ω and µ are rotational and magnetic axes, respectively.The star is shown to be rotating anticlockwise when viewed from above the Ω axis.

Figure 5 .
Figure 5. Results of 3D MHD simulations of the model with parameters of EX Lupi.1st row from the top: The distribution of the kinetic energy of the infalling matter at the surface of the star at different moments of time t counted from the beginning of simulations.Time is measured in days.Red and dark blue colors show areas of the highest and lowest kinetic energy, respectively.Vectors µ and Ω show the directions of the magnetic and rotational axes of the star.The white arrow shows the phase shift.Spots are shown at the inclination angle of the observer i = 32 • and at the zero phase of stellar rotation.2nd row: 3D views of matter flow at different moments of evolution.One of the density levels is shown for clarity.Lines are selected magnetic field lines.3rd row: Density distribution in the µ − Ω plane.4th row: Density distribution in the plane perpendicular to µ − Ω plane.Red and dark blue colors show the highest and lowest densities, respectively.Bottom left panel: the accretion rate onto the surface of the star.Bottom right panel: The lightcurve calculated for inclination angle i = 32 • .

Figure 6 .
Figure 6.Hot spots at different phases of stellar rotation during one rotational period of the star.The color background shows the flux of the kinetic energy density on the surface of the star.µ and Ω are vectors of the magnetic and rotational axes of the star.

Figure 7 .
Figure 7.The evolution of multiband phases from the onset of the outburst to the subsequent post-outburst are shown.A straight line fit (dotted line) is plotted on the first three phase points of AAVSO multiband phases (during outburst decay) as they follow a linear trend.The phase values of the CTIO lightcurves are shown near the peak of the outburst.Other phases are plotted from Kóspál et al. (2022), ASAS-SN, and LCRO lightcurves.The grey dash-dot lines are the model predicted phase decay curve estimated from the accretion rate during the outburst, as discussed in Section 4.3.The model-predicted phase decay curve is shifted to match the AAVSO phases at the epoch of the first spectrum observation, shown for U and R bands.The color-code for the various bands are: cyan (U ), blue (B), green (V, g ′ ), orange (R, r ′ ), brown (I, i ′ ) and pink (z ′ ).

Figure 8 .
Figure 8. Lissajous figures for the LCRO g ′ , r ′ and i ′ band lightcurves during JD = 2460127 -2460145.Left column: Normalized lightcurve cutouts along with their best sine fit are plotted for g ′ , r ′ and i ′ band by green, orange and brown colors, respectively.The photometry is 1-day averaged.Middle column: The yellow dots are the Lissajous figure for the sine fit onto the respective lightcurves, with green arrows pointing to the nearest next epoch.Right column: The normalized flux-flux plots (Lissajous figures) for combinations among g ′ , r ′ and i ′ bands (blue dots).Black arrows connect to the nearest next epoch observation data.The Lissajous figures from the middle panel are shown in this panel for comparison.

Figure 9 .
Figure 9. Top panel : ASAS-SN g-band phase evolution is shown during the pre-and post-2022 March outburst.The Y-axis is the phase in degrees, and the X-axis is time.The grey circles are phases from the ASAS-SN g-band lightcurve while the green circle is a manually estimated phase from ASAS-SN g-band lightcurve cutout (see Section 3.1 and Figure 1).The orange-colored inverted triangle is the phase from the TESS sector-12 lightcurve.The horizontal line attached to the grey circles shows the range of epochs used for phase estimation.Bottom panel : ASAS-SN g-band lightcurve is plotted with sections of colors showing the data used for phase estimation for the respective phase values in the top panel.The Y-axis is flux in mJy, while the X-axis is time in years.Only a portion of the outburst is shown.The outburst flux increased to ∼ 90 mJy.

Figure 10 .
Figure 10.RV amplitudes are plotted against the hotspot's latitude location (left panel ), latitude size (middle panel ), and longitude size (right panel ).Each plot of RV amplitudes vs a parameter has other two parameters as variables.The scatter in RV amplitude, say around 80 • is because of various latitudinal (10 • -30 • ) and azimuthal (10 • -100 • ) sizes acquired by the hotspot.A horizontal dotted line represents RV amplitude of 2.45 km/s, for He I 5875 Å (see Table2).The colored and grey circles represent RV amplitude larger and smaller than 2.45 km/s, respectively.The X-axis is in degree, and Y-axis is in km/s.
the hotspot.Any velocity component perpendicular to the disk is neglected.r m and φ are estimated for the given epoch.Solving the equations 6 and 7, we get a bound on the radial distance of the hot blob.The radial distances are 0.80-1.52R* , 0.50-4.03R* and 0.54-3.08R* for the epochs JD = 2459682.4127,2459683.4179and 2459712.604,respectively.A larger inclination angle would reduce these radial distances of the hot blobs.The estimated radial distances for the hot blobs are smaller than the inner disk radius and thus it further supports the idea that the blue-shifted wings in Ca II IRT lines could be representing the blobs of matter falling through the accretion funnel.

Figure 11 .
Figure 11.Polar plot showing the Ca II IRT spectra at the azimuthal location of the hotspot on EX Lupi's surface.Grey sparsely dashed radial spikes show the hotspot location at the epoch of spectra.Black-colored numerals show the order of the spectra observation.A blue dashed vertical line shows the mean position of the blue wing, as estimated from the 3component Gaussian fit to the spectral lines.The velocities of the blue shifted wing on spectra marked 3, 4 and 8 are -116.6km/s, -71.4 km/s and -81.4 km/s, respectively.Ca II IRT lines at 8498.02 Å, 8542.09Å and 8662.14Å are shown by blue, green, and red colors respectively.Grey dash-dot lines on spectra 3, 4, and 8 mark the zero velocity position.The continuous radial spikes are at an interval of 30 • .The Ca II IRT lines are plotted at different radial positions for visual clarity.

Figure 12 .
Figure 12.Top panel: TMMT lightcurve cutouts showing hourly photometric flux enhancement.The X-axis is in Hours for a particular date and the Y-axis is scaled and shifted by mean of respective band fluxes.B, V, and I band lightcurve cutouts are plotted by blue, orange and red dots respectively.The Savitzky-Golay filter is applied over B, V , and I cutouts and is plotted with respective color lines.The last column shows the scaled and shifted fluxes from the rest of the observations where no clumps were detected.The colors of the data points represent the epoch of observation, with violet being the first and red being the last epoch of TMMT observation.The top panel of the last column shows B-band fluxes while the middle and bottom panels show V and I band observations respectively.Bottom panels: LCRO g ′ , r ′ and i ′ band lightcurves showed similar hour timescale photometric enhancement on 2022 March (left, during outburst) and 2023 March (right, quiescence of post-outburst).tryas well as high-resolution multi-epoch spectroscopy.We searched for a change in the hotspot's azimuthal and latitudinal position during the 2022 March outburst and compared it with the hotspot's position during the outburst of the year 2008.The conclusions of our work are enumerated below:1.The hotspot's azimuthal position has been stablefor about one and a half decades, starting from the epochs preceding the year 2008 outburst.The

Figure 13 .
Figure 13.Color-color (top) and color-magnitude (bottom) plot during the outburst are represented by grey dots.Colored dots represent the hour timescale events.All 12 visually identified events are shown by color spectrum.The color barshows the order number of the events.The solid black line is reddening vectors plotted at an interval of AV = 0.5 mag(Cardelli et al. 1989)

Figure 14 .
Figure 14.The rise timescale (hours) is plotted against the fall timescale (hours) and a linear fit is shown by the dashed line.The various colors mark the various epochs of the event.The circular and square markers denote the events as observed by TMMT and LCRO telescopes, respectively.

Figure 15 .
Figure 15.EX Lupi lightcurve over a baseline of 6 years.Insets of periodic oscillations in the photometry are shown.The Y-axis is in mJy while the X-axis is in HJD.The two middle insets (green and blue) during the rise and fall of the outburst show photometric oscillations overlaid on the increase and decrease of the photometric continuum.

Figure 16 .
Figure 16.Periodogram for post-outburst multiband photometric data.The dotted horizontal lines correspond to a level of 10 % False-Alarm-Probability for corresponding data.The legend shows the estimated period from frequency corresponding to maximum power.The black dashed vertical line shows the location of the period = 7.417 days.

Figure 17 .
Figure 17.Multi Gaussians fit to the EX Lupi spectral lines are plotted.The pink color curve is the combined fit of the spectral line, while the blue color represents the NC of the line and green represents the BC (or the redshifted, for He II).All the spectral lines correspond to the peak of the outburst, JD = 2459654.4899.The spectra were calibrated before fitting as discussed in Section 2.1.The last profile plot is modelled as a sum of four Gaussians: blue represents Mg I 5167, green represents Fe I 5169 and cyan represents Mg I 5172 and red represents the broad component.

Figure 18 .
Figure 18.Top panel : The accretion rate inferred from the emission lines are shown by colored dots while the average accretion rates are plotted by black circles.Middle panel : The mass accretion rate inferred from the photometry is shown by black dots.The blue dotted vertical lines are the epochs of our HRS-SALT spectra observation, shown in the upper panel.Bottom panel : The scaled AAVSO U -band and the scaled ASAS-SN g-band lightcurves are plotted to show the correspondence of the epochs of the outburst.

Table 1 .
The photometric and spectroscopic observation logs from TMMT, LCRO, CTIO and SALT-HRS observations.This table contains a portion of the whole data.The JDs of TMMT, LCRO, and CTIO data are rounded off to two decimal points to show all the multiband data points simultaneously for close epoched observations.The online table contains a complete data set.
Sicilia-Aguilar et al. (2015)g the post-outburst quiescence of the year 2008.Cruz-Sáenz deMiera et al. (2023);Wang et al. (2023)showed that the prior mass accretion rate estimates using 10% of Hα line-width are underestimated.We used a new mass accretion rate estimate of 1.34 × 10 −9 M ⊙ /yr fromWang et al. (2023)to calculate the hotspot latitude on 2010 May 4. This accretion rate is calculated from the line luminosity of Br γ in X-Shooter spectra, using relations fromAlcalá et al. (2017), and hence it is best suited for comparison with our mass accretion estimates.The hotspot latitude during the quiescence on 2010 May 4 is 59 • .According to this theoretical dipolar magnetic field model, the expected latitude change between the quiescence of the post-2008 outburst and during the 2022 March outburst is ∼ 10 • .This is consistent with the observed difference in the hotspot latitude between the 2022 March outburst and that of the quiescence fromSicilia-Aguilar et al. (2015)and Campbell-White et al. ( Sicilia-Aguilar et al. 2012ated spot latitudes from prior estimates (e.g.,Sicilia-Aguilar et al. 2012 • to ∼ 65 • explains the RV amplitude increase.Outburst in EX Lupi caused the hotspot to move azimuthally forward and latitudinally down.The phase of the EX Lupi's lightcurves (g ′ ,r ′ and i ′ bands from LCRO and g-band from ASAS-SN) remained almost constant to the outburst state phase value up to almost one and half years after the outburst.But, the ASAS-SN g-band phase value for the lightcurve from 2023 July 16 (JD

Table 2 .
Sine fit parameters to radial velocities derived from the HRS-SALT spectra of EX Lupi during its 2022 March Outburst.a indicates that the RV phases are re-estimated at our startJD.b and c indicate the RV amplitude ratio from Sicilia-Aguilar et al. (2015) and Campbell-White et al. (2021) .0