Multi-wavelength observations of multiple eruptions of the recurrent nova M31N 2008-12a

We report the optical, UV, and soft X-ray observations of the 2017 − 2022 eruptions of the recurrent nova M31N 2008-12a. We find a “cusp” feature in the r ′ and i ′ band light curves close to the peak, which could be related to jets. The geometry of the nova ejecta based on morpho-kinematic modelling of the H α emission line indicates an extended jet-like bipolar structure. Spectral modelling indicates an ejecta mass of 10 − 7 to 10 − 8 M ⊙ during each eruption and an enhanced Helium abundance. The super-soft source (SSS) phase shows significant variability, which is anti-correlated to the UV emission, indicating a common origin. The variability could be due to the reformation of the accretion disk. We infer a steady decrease in the accretion rate over the years based on the inter-eruption recurrence period. A comparison of the accretion rate with different models on the M WD − ˙ M plane yields the mass of a CO WD, powering the “H-shell flashes” every ∼ 1 year, to be > 1 . 36 M ⊙ and growing with time, making M31N 2008-12a a strong candidate for the single degenerate scenario of Type Ia supernovae progenitor.


INTRODUCTION
Nova eruptions are a consequence of thermonuclear runaway on the surface of a white dwarf (WD) primary in cataclysmic binary systems, resulting in the ejection of material in the range of 10 −7 − 10 −4 M ⊙ (Gehrz et al. 1998;Hernanz &Jose 1998 andStarrfield 1999).Inherently, all novae are supposed to be recurrent, with the primary WD and the secondary red-giant/sub-giant star sustaining all the eruptions.The observed recurrence period of novae can range from 1 year (M31N 2008-12a; Darnley et al. 2014) to 98 years (V2487 Ophiuchi; Schaefer 2010).
M31N 2008-12a is an extraordinary RN whose eruptions have been observed every year in 2008-2023 (Table 1).It was first discovered during its 2008 eruption by Nishiyama & Kabashima (2008), though previous eruptions in 1992, 1993 and 2001 have been retrieved from the archives.Since the 2013 eruption, it has been monitored and studied across different wavelength ranges to understand its short recurrence period (2013 eruption - Darnley et al. 2014;Henze et al. 2014a;Tang et al. 2014;2014-Darnley et al. 2015c;;Henze et al. 2015b;2015-Darnley et al. 2016, 2017b,c;,c;2016. The optical light curve and spectral evolution of this very fast RN were found to be similar during all the eruptions, with Balmer, He, and N lines dominating the spectrum.Light curves showed an extremely rapid rise to maximum (∼ 1 day) followed by a fast linear decline for about 4 days and a plateau with slow decline and jitters from day 4 to 8. The multi-eruption UV light curves were similar, with an initial rapid linear decline followed by slow plateau-like declines.The plateau phase coincided with the supersoft X-ray source (SSS) phase.Like optical and UV light curves, the SSS phase was similar during multiple eruptions.
The 2016 eruption (Henze et al. 2018e), however, deviated from the general trend.It occurred after a longer inter-eruption gap, the optical light curve showed a short-lived cuspy feature, and the UV and X-ray fluxes disappeared relatively early compared to the previous eruptions.The "peculiar" behaviour of the 2016 eruption was suggested to be due to a lower accretion rate prior to the 2016 eruption.
Theoretical models generated to satisfy the short recurrence period and short turn-on time of the SSS phase have allowed constraining the mass of the WD to near the Chandrasekhar limit (Tang et al. 2014;Kato et al. 2014).
Deep Hα and HST imaging revealed the presence of an elliptical (134 × 90 parsecs) super-remnant nebula around M31N 2008M31N -12a (Darnley et al. 2015c).The size and mass of the shell indicate that the system has been undergoing eruptions for ∼ 10 6 years (Darnley et al. 2019c) and would likely do so for ∼ 2 × 10 4 more years before the WD attains M Ch (Darnley et al. 2017c).
This paper discusses the optical photometric and spectroscopic observations during the 2017 − 2022 eruptions and the evolution of UV and soft X-ray emission based on Swift archival data and data from our observations with the AstroSat.We explore spectroscopic modelling to reveal physical processes during outbursts and derive associated physical parameters.The behaviour of the RN M31N 2008-12a during the 2017 − 2022 eruptions is also compared to that of the previous eruptions.We end with a discussion on the recurrence period and its implication on the accretion rate and the mass of the primary WD.

Optical
Photometric and spectroscopic observations were carried out using the following telescopes and instruments.The log of optical observations is given in Table 2.
The GIT images were pre-processed, i.e., bias subtracted, flat-fielded and cosmic-ray corrected by the automated pipeline of the GIT (Kumar et al. 2022).Multiple exposures were obtained every night, and in the case of low SNR, images from the same night and the same filter were stacked using SWarp (Bertin et al. 2002).
For the 2018 data, aperture photometry was performed using an aperture of ∼ 2.5 ′′ , close to the fullwidth half maximum (FWHM) of the stellar profile in the images due to low SNR.For 2020-2022 data, the FWHM was first calculated using SExtractor (Bertin & Arnouts 1996) and subsequently used in Image Reduction and Analysis Facility (IRAF 2 , Tody (1993)) to perform PSF photometry.The magnitudes of the local standard stars given in Darnley et al. (2016) were converted from BV RI to gri using the transformations in Jester et al. (2005) to determine the zero points for photometric calibrations.Aperture photometry was also performed and found to be consistent with PSF photometry, eliminating any systematic differences that could arise from aperture photometry of the 2018 data.

Himalayan Chandra Telescope
The Himalayan Faint Object Spectrograph Camera (HFOSC) 3 mounted on the 2 m Himalayan Chandra Telescope (HCT) located at IAO, Hanle, India was used to obtain images in the V RI bands on 2018 Nov 08 UT.HFOSC is equipped with a 2K × 4K CCD.The pixels correspond to an image scale of 0.296 ′′ pixel −1 , with a FoV of 10 ′ × 10 ′ for the central 2K × 2K region.The images were pre-processed using the standard routines in IRAF.The instrumental magnitudes were obtained using aperture photometry with an aperture set at a radius three times the FWHM.Differential photometry was performed with respect to the local standards (Darnley et al. 2016) to account for the zero points of the images.
Optical spectra obtained using HFOSC (Pavana et al. 2018;Sonith et al. 2021;Basu et al. 2022) are given in Table 2.We used a grism with R ≈ 1200 in the wavelength range of 3500 -7800 Å.Data reduction was performed in the standard manner using IRAF.All the spectra were bias subtracted and cosmic rays were cor- Eruption date (a)  Discovery SSS-on date (b)  Days since Detection wavelength References (UT) mag (Filter) (UT) last eruption (Observatory) (1992 Jan 28)  rected before extraction.Wavelength calibration was carried out using the FeAr arc lamp spectrum.Spectrophotometric standard stars, Feige 110 (2018, 2019 and 2020 eruptions) and Feige 34 (2021 and 2022 eruptions), were used to correct for the instrumental response and bring the spectra to a relative flux scale.Absolute flux calibration was done using zero points obtained from broadband magnitudes based on photometric observations within 3-4 hours of the spectroscopic observation, except in the case of the 2020 observations, which had a gap of around 7 hours.

J.C. Bhattacharyya Telescope
The 2K × 4K UK Astronomy Technology Centre (UKATC) CCD mounted on the 1.3 m Jagadish Chandra Bhattacharya Telescope (JCBT)4 , located at the Vainu Bappu Observatory (VBO), Kavalur, India, was used during the 2018 eruption.It has a 15-micron pixel size corresponding to an image scale of 0.3 ′′ pixel −1 , with a FoV of 10 ′ × 20 ′ .JCBT observed the nova in BV bands 2 days after the eruption.The images were reduced and calibrated following the same steps used for HCT data.

Other data sources
Our observations were combined with publicly available photometric data for analysis, from sources referenced below.

Ultraviolet
Photometric studies were done using images obtained in the ultraviolet (UV) bands from the following two telescopes.

Swift UVOT
High cadence UV imaging data of M31N 2008-12a are available from the Swift (Gehrels et al. 2004) archive 6 .The nova has been monitored by UVOT since 2013 during each eruption.The log of Swift observations between 2017 − 2022 are summarised in Table 3.We have used the uvw2 (1928±657 Å) archival data in this study.The uvot task in HEASOFT (v6.29) was used to extract the magnitudes from a source region of radius 5 ′′ after background subtraction.Since the field is crowded in uvw2 filter, a source-free 10 ′′ radius circle, 80 ′′ southwards of the object, was chosen to estimate the background.The calibration assumes the UVOT photometric (AB) system (Poole et al. 2008 andBreeveld et al. 2011) and is not corrected for extinction.

AstroSat UVIT
6 https://www.swift.ac.uk/index.phpAstroSat (Singh et al. 2014) is a space-based telescope with the Ultraviolet Imaging Telescope (UVIT) as one of its instruments.UVIT observed M31N 2008-12a during its 2019 − 2022 eruptions in F148W (1481±500 Å) filter (see Table 3).The level 1 UVIT data was downloaded from the Indian Space Science Data Center (ISSDC) 7 and reduced using CCDLAB following standard routines presented in Postma & Leahy (2021).The orbit-wise images were registered and merged to obtain a single image with a high SNR on which astrometry was performed.The average PSF size in UVIT images was ∼ 1.5 ′′ across all epochs.We performed PSF photometry with an aperture correction term derived from "good stars" to account for the broad PSF wings in UVIT images.The zero points for photometric calibrations in the AB system were adopted from Tandon et al. (2020) and have not been corrected for extinction.

X-ray
Both Swift and AstroSat observe simultaneously in the UV and X-ray wavelengths.Soft X-ray observations from both facilities were used to study the eruptions during the SSS phase.

Swift XRT
Swift X-Ray Telescope (XRT; Burrows et al. 2005) data were downloaded from the Swift archive, with Observation IDs same as that of the UVOT data (Table 3).For the analysis, HEASOFT (v6.29) with XIMAGE (v4.5.1) and XSELECT (v2.5b) were used following the guidelines summarised by UKSSDC 8 .XRT count rates were determined using both XIMAGE and XSELECT tools provided by HEASARC.Both results followed the same trend and were within 1σ errors of each other.We have presented the results from the XIMAGE sosta/optimize analysis as it corrects the counts for vignetting, dead time loss, background subtraction and the PSF of the instrument.We used XSELECT to extract the spectra for each snapshot.ARF files were generated from the exposure maps, while RMF files were taken from the calibration database.Spectral analysis was performed in XSPEC (v12.12.0) assuming Poisson statistics (cstat) due to low counts.We used ISM abundances given in Wilms et al. (2000) and the Tübingen-Boulder (tbabs) ISM absorption model to account for the intervening medium.

AstroSat SXT
7 https://astrobrowse.issdc.gov.in/astroarchive/archive/Home.jsp 8 https://www.swift.ac.uk/analysis/AstroSat Soft X-ray telescope SXT (Singh et al. 2017), placed in parallel alongside UVIT, is capable of observing in the 0.3 − 8.0 keV range simultaneously with UVIT (Table 3).SXT observed the SSS phase of M31N 2008-12a during the 2020 and 2021 eruptions.Level 2 data was downloaded from ISSDC, and the cleaned event files were merged using SXTTools in Julia9 .The source region was chosen as a circle of radius 7 ′ , smaller than the usual SXT PSF of ∼ 12 ′ so to avoid contamination in the crowded M31 field.We set the bin size to 8 hours and energy range to 0.3 − 2.0 keV in XSELECT (XIMAGE was avoided due to compatibility issues) to attain an adequate SNR for light curve analysis.SXT spectra were extracted using XSPEC from the merged SXT cleaned event files.A new ARF file was generated corresponding to the smaller source extraction region for analysis.

Epoch of eruptions
For a very fast RNe, like M31N 2008-12a, a tight constraint on the eruption time is useful for generating light curve models ( §3.3) and studying its recurrence nature ( §7).Hence, we estimate the epochs of eruption based on available detection and pre-discovery magnitudes and non-detection upper limits for all eruptions since 2017.Even though the exact time of eruption is uncertain, it can be well approximated by the mid-point of first detection and last non-detection in each year.Amateur astronomers' interest in M31 and the increase in survey telescopes over the past decade have made it possible to constrain the eruption date to well within a day.The uncertainty in the eruption dates spans between the first detection and the last non-detection.
The 2017 eruption was discovered just in time to be called the "2017 eruption" on Dec 31.77UT by Boyd et al. (2017).Darnley et al. (2017a) reported spectroscopic confirmations on the same day.The last nondetection was on Dec 31.38 UT at an upper limit of m clear = 19 mags (Naito et al. 2018a).
The 2018 eruption was discovered on Nov 06.80 UT at a magnitude of 19.15±0.05by Darnley et al. (2018b) and confirmed spectroscopically by Darnley et al. (2018c) on the same day.Tan & Gao (2018) reported the last nondetection at > 21.20 mag on Nov 06.54 UT in clear filter.
The 2019 eruption was detected on Nov 06.71 UT by Oksanen et al. (2019) at 19.40 mags.The first spectrum taken on Nov 06.83 UT (Darnley et al. 2019d) confirmed the recurrence of M31N 2008-12a.The last non-detection information was not publicly available for 2019, so we adopted the eruption date provided in Darnley et al. (2019b).Darnley et al. (2020b) discovered the 2020 eruption on Oct 30.89UT.The nova was, however, also detected in images captured 90 minutes before the discovery (Galloway et al. 2020), and we use the pre-discovery detection (m g = 18.74) and non-detection (m g > 19.4) to constrain the eruption time.The discovery was spectroscopically confirmed on the next day (Darnley 2020).
The 2021 eruption was discovered on Nov 14.38 UT by Itagaki et al. (2021) and was spectroscopically confirmed by Wagner et al. (2021).Tan et al. (2021) gave pre-detection upper-limits at m clear > 19.00 mags on Nov 13.96 UT.
The 2022 eruption was discovered on Dec 2.83 UT by Perez-Fournon et al. (2022).It was undetected until Dec 02.61 UT at m L > 19.60 mags (Shafter et al. 2022).Spectra taken on Dec 3.84 UT confirmed the source to be a recurrence of the nova (Darnley & Healy 2022), its 15th successive eruption in as many years.
The estimated eruption dates for 2017-2023 are presented in Table 1 together with those of the previous ones.

Light curve evolution
The optical light curves of the 2017 − 2022 eruptions, based on our observations and publicly available data, are shown in Figure 1.The light curves indicate a rapid rise to the peak magnitude in < 1 day from discovery, followed by a fast decline with t 2 ≈ 2 − 4 days in g ′ r ′ i ′ bands.A brief description of the optical light curve for each eruption during 2017 − 2022 is provided below.
The 2017 light curves show a rapid decline in the first 4 days, followed by a slow decline in all the bands.The decline rate is marginally faster in the r ′ band (0.84 ± 0.12 mag day −1 ) compared to the g ′ band (0.74 ± 0.19 mag day −1 ).However, the decline rates in g ′ and r ′ were measured using only two data points and are within error bars of each other.
The maximum phase during the 2018 eruption appears to be broader.Following the initial rise, the magnitude declined by ∼ 0.4 mag in the BV RI filters, and after a brief halt for about 0.3 days at this level, a marginal increase in the brightness by 0.2−0.3mag lasting ∼ 0.7 days can be seen.We also note that the peak magnitude observed in 2018 in R (18.50 mag) is fainter than the peak R or r magnitudes in the other years.In 2018, the initial decline in the r ′ (0.68±0.03 mag day −1 ) and i ′ (0.53 ± 0.04 mag day −1 ) bands was slower by 0.1 − 0.2 mag per day compared to other years.Further, the decline rates of r ′ and i ′ are slower than g ′ band (0.84 ± 0.02 mag day −1 ) in 2018.
The 2019 data set is sparse and restricted to the initial rise and maximum phases.The nova is brighter in R compared to BV bands during the rise and the peak.It rises about 0.4 mag in all the bands in about 0.5 days from discovery.The peak R magnitude reached in 2019 is m R = 18.12, which is higher than most other eruptions.
The g ′ band traces the rise of the 2020 light curve at ∼ 1.2 mag day −1 while the r ′ band traces the smooth decline from the peak at 0.89 ± 0.05 mag day −1 for 2 days.Limited data only restricts the light curve analysis to the r ′ band.
The 2021 eruption was caught almost a day before it reached its peak.The rise was sharper in the i ′ band compared to the g ′ r ′ bands.It declined rapidly in the g ′ band at 1.25 ± 0.20 mag day −1 but relatively slowly in the r ′ band at 0.88 ± 0.13 mag day −1 for the first 3 − 4 days.The decline rate then slowed with significant enhancement in the r ′ band flux around 6 days after the eruption.The 2022 eruption light curve was similar to previous years.The rise is well captured in the r ′ band with a rate of ∼ 2 mag day −1 , the fastest in the last six years.It then declined speedily in the g ′ band (1.01 ± 0.05 mag day −1 ) but relatively slowly in r ′ band (0.79 ± 0.03 mag day −1 ) and even more slowly in the i ′ band (0.69 ± 0.06 mag day −1 ).
Generally, the nova declines fastest in the g ′ band and comparatively slower in the redder bands.Darnley et al. (2016) combined the 2013, 2014, and 2015 eruptions and found the decline rate to be fastest in the V band at 1.21 mag day −1 during the initial decline phase.The evolution of the optical light curve during the later phases is unavailable as the nova fades beyond the detection limit of the 1−2 m class ground-based telescopes generally used for follow-up observations.
The UV light curves (Figure 2) obtained from spacebased telescopes have a wider time coverage (∼ 20 days) and show a linear decline in the 2017 − 2022 uvw2 magnitudes from day 0 to 3 since the eruption.The plateau phase begins with a re-brightening at ∼ 4 days since the eruption, which is also coincident with the SSS turn-on time (yellow shaded region in Figure 2).The evolution during this phase indicates a gradual de-cline of ∼ 0.15 mag day −1 .However, this decline is not smooth but accompanied by small undulations.Some P-type galactic recurrent novae, such as T Pyx, RS Oph and U Sco, show variability in the V -band during their plateau phase (Strope et al. 2010).Since the optical photometry during this phase is insufficient, we are unable to comment on the presence of a similar variability in the optical bands in the case of M31 2008-12a.We encourage continuous and deep optical monitoring during the plateau phase in future eruptions.
The F148W data indicate a brightness similar to uvw2 during the initial decline phase, but becomes fainter than uvw2 by > 0.5 mag during the SSS phase.
The general trends across all the UV-optical bands in 2017 − 2022 eruptions are more or less similar over the last six years and consistent with the previous eruptions (2013 : Darnley et al. 2014, 2014: Darnley et al. 2015c, 2015: Darnley et al. 2016, 2016: Henze et al. 2018e).However, some deviations of the 2016 light curve were noted and are discussed in §8.

Colour Evolution
The (F148W−uvw2) colour was determined from observations taken on the same day.From Figure 3, it is seen that the (F148W−uvw2) colour becomes bluer at the onset of the SSS phase but is significantly redder during the SSS phase.
In the optical bands, we restrict the colour analysis to only the SDSS primed filters to avoid instrumental and/or filter dependencies of the Bessel filters.Nearsimultaneous observations in g ′ r ′ and r ′ i ′ filters for the same eruption were used to estimate the colours.The colours are plotted together as a function of days since the eruption to bring all the outbursts to the same time scale.The (g ′ − r ′ ) colour linearly increases up to day 3 from the eruption and then decreases.This timeline agrees with the initial rise and linear decline phase of the light curve.The (r ′ − i ′ ) colour shows a steep reddening during the rising phase, which then slows down as the nova follows its initial decline.After 3 days from the eruption, the (g ′ − r ′ ) colour becomes bluer while the (r ′ − i ′ ) colour also tends to be bluer, but due to only one data point between day 3 and day 4, we are unable to confirm this.Beyond day 4, the (r ′ − i ′ ) colour becomes redder when the nova enters its plateau phase.A similar trend was also noted by Darnley et al. (2016) in their colour plots.

Light curve modelling
To understand the temporal evolution of the light curves, we model them by breaking them into three phases corresponding to different decline rates.The  uvw2, g ′ , r ′ , and i ′ band are used.The eruption times presented in Table 1 are used as the reference times, and we measure all other times in days with respect to the reference date of each eruption.The three phases considered are 1.The rise to peak: from eruption to t ≈ 1.5.
2. The initial steep decline: t max ≤ t ≤ 3.5, where t max is the time of maxima 3. The slow decline: t ≥ 3.5 Due to extensive coverage in the uvw2 filter, we could notice that the rate of decline decreased even further beyond 8 days of eruption.The final phase is thus divided into two segments in uvw2.Darnley et al. (2016) employed a similar 4-phase division of all light curves to analyse previous eruptions.
First, we generated models for the combined 2017 − 2022 eruptions and obtained the light curve properties at different phases given in Table 4.Then, we combined the 2013 − 2015 eruptions' data (see §1 for references) with those from 2017 − 2022 and generated overall light curve models spanning from 2013 to 2022.We note here that the data in g ′ is sparse as it was not used in most of the observations before 2016.The 2016 data set has been intentionally excluded as an outlier as it deviated significantly from the general trend of other eruptions (plotted in red points in Figure 4), especially in the UV light curve.The combined light curve models are presented in Figure 4 and the light curve parameters are tabulated in Table 4.
Additionally, Gaussian process (GP) regression techniques were employed to fit the entire light curve for each band.The regression results including a 3σ error range, are shown in blue in Figure 4.

The rise to peak
This phase has been modelled with a quadratic function to trace the rise to the peak and the fall just after.Limited availability of uvw2 data during this phase led to only partial modelling of the rise in this band.On the other hand, the optical bands have dense coverage of the rise and the peak.The uvw2 and g ′ bands show a smooth rise towards the peak and a smooth decline from the peak.In contrast, the r ′ and i ′ bands show a "cusp" just before the peak of the modelled light curve is attained.
On combining the 2013 − 2015 light curves with those from 2017 − 2022, we clearly see the cusp (Figure 4), at least in r ′ and i ′ bands, just before the peak is attained.The cusp-like feature is evident in the 2021 and 2022 light curves (Figure 1) as the data points are dense.The 2018 light curve also indicates the presence of the cusp, although with lesser brightness.The 2017 V band and 2019 R band data also hint at the cusp.Observations post-2016 indicate that the cusp is most likely present during all eruptions.The cusp was first noted by Henze et al. (2018e) in the 2016 eruption in multiple wavebands, who speculated the "cusp" could be an isolated event in 2016 (and 2010?), connected to the short SSS phase and long inter-eruption period of the 2016 event.Alternately, as indicated by Henze et al. (2018e) and Darnley & Henze (2020), the low cadence observations of the 2013-2015 eruptions during the rising phase could have allowed for this feature to be missed.We suggest this feature is a general trend and not connected to the shorter SSS phase of 2016.However, to confirm it, we encourage very early detection and dense observations during the rise phase in all UVOIR bands of future eruptions.
The time of maxima was calculated from the quadratic model fits to the data near the peak.The peak magnitudes and the time of the peak are given in Table 4.

The initial steep decline
The initial decline is very fast in all the bands.The decline rates during this phase are modelled by a straight line fit from t max to t = 3.5 days after the eruption.The uvw2 decline rate of 0.93 mag day −1 is steeper than the uvw1 decline rate (0.78 mag day −1 ) reported by Henze et al. (2018e).In the 2017 − 2022 data, we find that the g ′ band decline rate (0.90 mag day −1 ) is marginally higher than the r ′ band (0.88 mag day −1 ) but when we combine it with the 2013 − 2015 data, the g ′ band decline rate (0.90 mag day −1 ) is less than the r ′ band (0.97 mag day −1 ).Darnley et al. (2016) noted that the decline is fastest in V (1.21 mag day −1 ) whereas the B and r ′ decline rates (in mag day −1 ) are 0.99 and 0.97 respectively.The decline rates in this phase are used to derive the t 2 times given in Table 4.

The slow decline
The slow decline phase is modelled with a linear fit from day 3.5 onward.This phase consists of the plateau in the light curves which is also coincident with the SSS phase in X-rays.Combining all eruptions from 2013 gives a sufficient number of data points in all the filters for modelling except in g ′ .The decline rate during this phase is low in all the bands (see Table 4).The r ′ band also show scatter during this linear decline, but these jitters are more prominent in the uvw2 filter.The t 3 time calculated from the straight line fits are 6.17, 10.83, and 14.15 days in r ′ i ′ and uvw2, respectively.The slowing down of the decline rate can be attributed to the expanding ejecta cooling at t ≥ 4 days from the eruption.It is also reflected in the colour evolution where after day 4, we see the system become redder (see Figure 3 and Figure 2 of Darnley et al. (2016)).Beyond day 8, in uvw2, we see a further decrease in the decline rate and model it with a different slope.Optical photometry is sparse after day 8, but some data points in the i ′ band are available, though not enough for modelling.The i ′ band excess (0.2−0.4 mags) around day 8 is most notable.This bump is traced by GP regression and is shown in blue in Figure 4.
During this phase, we see a secular trend of decreasing flux with undulations on top of it.This scatter from the smooth decline in uvw2 has been discussed in §6.Towards the end of the final decline phase, when the SSS flux drops to zero at t > 18, we see a brief period of UV re-brightening before fading away to quiescent.The line fluxes of the emission features clearly identifiable in the individual spectra are listed in Table 5.Also provided in the table are the FWHM velocities obtained from a Gaussian profile fit to the emission lines (using IRAF).The velocities calculated from the widths of the emission lines have been corrected for the instrumental response by de-convolving with the width of night skylines.The initial velocities within 1 day of eruption are as high as 5000 km s −1 , typical of very fast novae.These observations are consistent with the previous eruption of M31N 2008Darnley et al. 2016Darnley et al. , 2015c)).At around ≳ 1.5 day after the eruption, the emission line widths narrow to 3000 km s −1 .The narrowing of the emission lines could be caused by the expansion and dissipation of the faster moving component (Shore et al. 1996) or by the interaction of the ejecta with the circumbinary material.We find the line velocity decelerating at v exp ∝ t −0.27±0.07 .This is similar to the estimate provided by Darnley et al. (2016) (v exp ∝ t −0.28±0.05 ), who argue that the deceleration is due to the interaction of the ejecta with the circumbinary medium, and that the ejecta is in Phase II of the shocked remnant development (Bode & Kahn 1985).The Hα profile, as seen in Figure 5, indicates that the ejecta geometry is structured and time-dependent.The temporal evolution of the Hα line shows a double-peaked structure, prominent in the 2016, 2019, 2020, and 2021 spectra, taken around 0.9 − 1.5 days after the respective eruptions.Around two days after the eruption, the double-peaked profiles give way to a relatively narrow boxy profile.

Estimation of physical parameters
To understand the physical conditions in the nova ejecta, the spectral synthesis code Cloudy (v17.02;Ferland et al. 2017) was used to obtain a 1D model using the procedure described in Pavana (2020).We generated a 1D model for the best SNR spectra taken on 2018 Nov 8.6 UT and 2019 Nov 7.6 UT.The top left (bottom left) panel of Figure 6 shows the 2018 (2019) synthetic spectrum obtained using a two-component (diffuse+clumps) model.In the 2018 spectrum, the effective temperature and luminosity of the central ionizing source were found to be 1.06×10 5 K and 10 37 erg s −1 respectively.A clump component and a low-density diffuse component of density 10 11 cm −3 and 10 10 cm −3 respectively were used to fit the emission lines in the observed spectrum.The ejected mass and helium abundance from the best-fit modelled spectrum were found to be 7.21 × 10 −8 M ⊙ and 2.47 ± 0.11 He ⊙ respectively using the relations given in Pavana et al. (2019) and references therein.For the 2019 spectrum, the effective temperature and luminosity of the central ionizing source were 7.20 × 10 4 K and 10 37 erg s −1 respectively.A clump component of 2.24 × 10 10 cm −3 and a diffuse component of 1.26×10 8 cm −3 could generate a synthetic spectrum close to the observed one.The ejected mass and helium abundance, in this case, were found to be 1.3 × 10 −8 M ⊙ and 3.09 ± 0.18 He ⊙ respectively.The ejected mass derived from X-rays (see §5.1) and spectral modelling are similar to that reported for the 2015 eruption by Darnley et al. (2016).Overabundance of helium has been estimated in other RNe such as RS Oph, V3890 Sgr, T Pyx (see Anupama & Pavana 2020 and references therein) and V745 Sco (Mondal et al. 2020).Velocities are color coded from 6000 km s −1 (red) to -6000 km s −1 (violet).Note that the scales are different in the ejecta morphology plot.
It was noted that the two-component model was insufficient to generate the synthetic spectrum with a high χ 2 value.The observed spectrum shows N II lines, which were clearly visible once a third, diffuse, component was introduced to the model.This implies that the N II and He lines are clearly originating from different regions with different physical conditions.However, since the optical spectrum of this extragalactic nova has low SNR, modelling with three components is beyond the scope of this work.With these uncertainties, modelling a high SNR spectrum with similar methods in the upcoming eruptions is recommended.
The Hα emission line profile during the 2018 and 2019 eruptions with multiple peaks encouraged us to obtain the morpho-kinematic structure for the ejecta using Shape (Steffen et al. 2011).We carried out the morphokinematic analysis of the Hα (and Hβ for 2018) velocity profile following the procedure described in Pavana (2020).
An asymmetric bipolar structure with bipolar cones and an equatorial ring (Figure 6) with a best-fit inclination angle of 80.75 • ± 1.21 • could generate the synthetic velocity profile of 2018 spectrum.The extended bipolar component stretched up to 4.52 × 10 12 cm along the ejecta axis from the centre while the central bipolar cones (opening angle of ∼91 • ) extended up to 3.62 × 10 11 cm.The inner radius of the equatorial ring and radii of the bipolar cones were 1.27 × 10 11 cm and 5.42 × 10 11 cm, respectively.A similar geometry with a best-fit inclination angle of 79.60 • ± 1.45 • could generate the synthetic Hα velocity profile of the 2019 spectrum shown in the bottom panel of Figure 6.The size of the extended bipolar component and the bipolar cone in the central region (opening angle of ∼40 • ) were 5.58 × 10 13 cm and 6.16 × 10 12 cm along the ejecta axis from the centre respectively.The inner radius of the equatorial ring and the radius of the bipolar cones were 4.12× 10 12 cm and 5.27× 10 12 cm, respectively.It should be noted that the He I (6678 Å) profile is blended with the broad Hα profile, and interestingly, the He I line arises from the inner bipolar cone region.
The outer part of the equatorial ring, the central bipolar cone and the extended bipolar region are discernible in the models shown in Figure 6.The extended bipolar nature suggests a fast-moving polar ejecta along the ejecta axis, i.e. jets, contributing more to the highvelocity hydrogen Balmer emission.

X-ray light curve
Figure 7 shows the light curve of supersoft X-ray emission during the 2017 − 2022 eruptions.The light curves from the previous eruptions are also shown for comparison.The emergence of the SSS phase is marked by the detections at ∼ 8 × 10 −3 counts s −1 , which increases to (3 − 4) × 10 −2 counts s −1 and stays around that level from 8 to 15 days after the eruption.This 'peak' of the SSS phase coincides with the UV light curve plateau region.The mean turn-on and turn-off time of the nova estimated from the mid-points of detections and non-detections of 2014 − 2022 eruptions are 5.06 ± 0.60 and 16.89 ± 0.96 days, respectively from the time of the eruption.The average SSS duration of the nova is 11.83 ± 1.56 days.Through the rise and during the SSS phase, the X-ray emission is variable, while the decline from the SSS phase is relatively smooth.Multiple "dips" are seen in the 2017 − 2022 XRT light curves.One around days 6 − 8, and an even more noticeable one around days 10 − 11.This drop in the count rate is evident in the SXT light curves of 2020 − 2021 eruptions (Figure 7).Variability in the X-ray emission during the SSS phase has also been noted in the previous eruptions by Darnley et al. (2016) and Henze et al. (2018e).The cause of this variability is not yet clear, and further high-cadence observations are required to understand its origin.We also note the unique short and faint nature of the SSS phase in the 2016 eruption compared to other eruptions, which is quite apparent in Figure 7.
To model the rise and decline in soft X-ray flux, we used a simple quadratic function.We included all the observations from 2013, except the peculiar 2016 eruption as its effect was seen most in the SSS phase.We plot all the individual and binned sets in Figure 8. Deviations from the naive quadratic function are evident, especially the peaks at days 9 − 10 and 11 − 12 and the dips at days 7 − 8 and 10 − 11.The prominent features between days 8 and 13 are present in the binned and unbinned data, indicating that these variabilities' causes last for more than half a day.The drop and rise of flux between days 10-11 seem to be a general feature of the SSS phase of M31N 2008-12a.Most of the variability is seen up to day 13, whereafter, the decline is relatively smooth.
X-ray studies of M31 novae (Henze et al. 2010(Henze et al. , 2011(Henze et al. , 2014b) ) have revealed the correlation of ejecta expansion velocity and the SSS t on time.The ejecta mass was calculated from the turn-on times and the ejecta velocities (v exp ) using the relation given in Henze et al. (2014b).A t on time of ∼5 days with an v exp of ≈ 2000 ± 200 km s −1 around this phase gives an ejecta mass range of 4.2 × 10 −8 < M ej,H / M ⊙ < 10.2 × 10 −8 .These are slightly higher than that calculated from optical spectra in §4.2 but less than the average mass accreted in a year.

X-ray spectroscopy
The XRT data was also used to extract spectra by merging two data sets obtained on consecutive days to increase the SNR.We fixed the H column density at 1.4 × 10 21 cm −2 (Darnley et al. 2016) but varied the blackbody temperature and normalization to attain the best-fit values.The time evolution of the SSS temper-ature is shown in Figure 8 for 2017 − 2022 eruptions.Not only does the fluxes peak 10 − 14 days after the eruption, but the temperatures also peak, suggesting a correlation between the SSS flux and temperature.In the 2020 eruption, a temperature fluctuation during the maxima can be seen.Such fluctuations have been reported before by Darnley et al. (2016).This pattern is not seen in other eruptions, possibly due to combining data sets of two consecutive days.In §5.1, it was noted that the rise to the maxima in the SSS phase shows variability, whereas the decline was smooth.The temperature evolution in Figure 8 also shows an asymmetry during the rise and decline of the SSS phase.These could be because of two different underlying causes.The increase in flux and temperature is due to the expansion and thinning of the ejecta, probing the deeper and hotter layers towards the WD surface.Whereas during the later stage, when the obscuring material is already dissipated, the decrease in flux and temperature is because of the residual nuclear burning slowing down and eventually stopping.
Spectra extracted from the merged SXT data is shown in Figure 9. Also shown are the contemporaneous XRT spectra obtained from merged snapshots of two successive days.The data has been restricted to below 2.5 keV for the SXT to avoid background contamination due to its large PSF compared to the XRT.Beyond 1 keV, the flux is too low, owing to the super-soft nature of the source.A faint hard X-ray tail (above 1.5 keV) can be seen in SXT data, but we could not be certain of its origin because of low SNR.The best-fit blackbody temperatures from SXT spectra are slightly higher than the XRT data during similar times (Figure 8).The modelled flux in the 0.3 − 2.0 keV range was similar in the 2020 spectra for both instruments but differed by a factor of 2 (higher in SXT) in the 2021 spectra.As the observations are not continuous and the XRT and SXT epochs do not coincide exactly, the mismatch could be because of the rapid variability seen in flux and temperature during the SSS phase in recurrent novae.

UV -X-RAY CORRELATION?
We noticed the 2016 uvw2 light curve was "shorter and less luminous" compared to the 2017 − 2022 uvw2 measurements during the SSS phase.Henze et al. (2018e) had found the same for soft X-rays and reasoned it to be due to a reduced accretion rate prior to the 2016 eruption.They could not comment on the 2016 uvw2 measurements because of the unavailability of the uvw2 light curve template at that time.This motivated us to find the connection between these two wave bands in the super-soft phase.Since both soft X-ray and UV show different trends during the super-soft phase, it was necessary to detrend them.The uvw2 light curves were detrended with a linear fit as it followed a linear declining trend during SSS phase, whereas the X-rays were detrended with a quadratic function as it followed a rise and a subsequent fall.The detrended light curves of each year from 2017 to 2022 are given in Figure 10.In 2017, we saw the UV and the X-ray fluxes behave inversely between days 7.5 and 13, and a Pearson correlation coefficient (hereafter r) value of -0.76 suggests a strong anticorrelation.In 2018, the anti-correlation lasts shorter from day 8.5 to day 11.5 but is stronger with a r value of −0.86.2019, on the other hand, does not show any strong correlation.In 2020 and 2022, there was a strong anti-correlation from day 9.5 to 14.5 when r value was −0.79 and −0.68, respectively.The 2021 detrended light curves show mild anti-correlation (r = −0.38)which is stronger than 2019 but weaker than the other eruptions during the days 9−14.This anti-correlation seen in most of the eruptions between UV and soft X-ray is strongest from day 8 to day 14 after the eruption, a time corresponding to the maxima of the SSS phase.UV and X-ray flux of most novae have been found to be uncorrelated (Page et al. 2022).Nonetheless, Ness et al. (2009) noted such anti-correlation of UV and X-rays in the detrended light curves for nova V458 Vul, although in the 0.6 − 10 keV X-ray range, just before the start of the SSS phase, and they suggested that it would imply that both the UV and hard X-rays originate from the same region.In HV Cet, the UV and X-ray flux were found to be tied up in phase, which Beardmore et al. (2012) argued to be due to the same cause, the orbital period in their case.V603 Aql also showed correlated UV-X-ray emission, though in this case, it was interpreted as due to X-ray illumination (Borczyk et al. 2003).Since the source of soft X-rays during the SSS phase is the nuclear burning on the surface of the WD, the anti-correlation during the SSS peak, in our case, would hint that the UV radiation origin is also close to the surface of the WD.It is possible that the accretion disk survives each eruption (Darnley et al. 2017b,c).The surviving partial accretion disk would emit UV radiation.The complete reformation of this partial disk could cause the variability seen in both UV and X-ray detrended light curves.The possibility of a wobbly, nascent accretion disk could also cause such a behaviour.

RECURRENCE PERIOD, ACCRETION RATE AND WD MASS
M31N 2008-12a has erupted every year since 2008, making it an exceptional case of the only RN observed 15 times and that too consecutively.This section focuses on the trend of the recurrence period and its relation to the accretion rate and the WD mass.

Increasing recurrence period
The mean recurrence period was reported to be P rec = 351 ± 13 days after the 2015 eruption by Darnley et al. (2016), which was updated to 363 ± 52 days after the outlier 2016 event by Henze et al. (2018e).
Since the 2016 eruption, M31N 2008-12a erupted six more times, and each year, the time gap between two successive eruptions has been more than the mean recurrence period except in 2018 (310.1 days) and 2020 (355.9 days).As of the 2022 eruption, the mean recurrence period is P rec = 363.6 days with a standard deviation of 40.3 days (Figure 11, right panel).On the other hand, the median recurrence period has increased from 347 days in 2016 (Henze et al. 2018e) to 360.5 days in 2022.Figure 11 shows two different sets of histograms, along with their kernel density estimates (KDE), for the recurrence periods between 2008−2015 and 2008−2022.On considering up to the 2015 eruption (grey histogram in Figure 11), the mode of the recurrence period is 340 days and the KDE peaks at 341.28 days (FWHM of 67.27 days).But on incorporating all the eruption information till 2022 (red histogram in Figure 11), the mode of recurrence period shifts to 360 days with the KDE peaking at 358.18 days (FWHM of 96.36 days).Over the last seven years, we see a clear increasing trend in the recurrence period.
To further investigate this matter, we plotted the period as a function of the eruption year in the right panel of Figure 11.We used the GP regression technique to extract the trend in the data set and associate errors with it.The data was modelled using a Matern kernel with a typical length scale of 15 years and an amplitude equal to the median of the recurrence period.We tested our model by applying it to the eruption dates of 2008−2021, and it could predict the 2022 eruption date within 3σ error limits.The 2022 data was subsequently included in the training set to project the upcoming eruptions.Since there is a gap of around 15 years between 1993 and 2008, we did not include the 1993 data point for our modelling, but an extrapolation of our model does seem to incorporate the 1993 eruption within error limits.For comparison, the mean and 1σ of the constant recurrence period model have also been shown in We emphasize here that our model is restricted to only the 'usual' eruptions that follow the trend.Anticipation of 'outlier' events, such as the 2016 eruption, is not feasible.The composition of the accreted material was fixed at H:He:Metals ≡ 70:28:2.We use their results for the massive WD cases (1.20 − 1.35 M ⊙ ), obtain a best-fit power-law relation between the accretion rate and the period and employ the same to infer the accretion rates for each cycle of M31N 2008-12a.These are also shown in the top panel of Figure 12.

Estimating the WD mass
We plot the accretion rates thus obtained corresponding to the periods of M31N 2008-12a in the last 15 years in the M WD − Ṁ parameter space in the bottom panel of Figure 12.For comparison, we over-plot the results from Kato et al. (2014), Prialnik & Kovetz (1995), Wolf et al. (2013), andTang et al. (2014), who predicted the WD mass and accretion rates for a recurrence period of 1 year.
In the M WD − Ṁ plane, when the accretion rate surpasses the critical threshold ( Ṁ > Ṁcr ), the WD exhibits a red giant (RG) like behaviour, undergoing surface mass burning at the critical rate, while any excess material is ejected in the form of optically thick winds (Kato & Hachisu 1994;Hachisu et al. 1996).Conversely, when the accretion rate falls below the critical threshold but remains above the stable H accretion rate, i.e.
Ṁcr > Ṁ > Ṁst , the burning on the WD's surface remains stable, capable of sustaining itself over an extended period as a supersoft X-ray emitter (Kato et al. 2014).Below the stability line, i.e.Ṁst > Ṁ , the accretion rate is insufficient to sustain continuous hydrogen burning.Systems within this parameter range experi-ence "H-shell flashes" or nova eruptions, a category that includes all RNe, including M31N 2008-12a.The limits on Ṁcr and Ṁst taken from Wang (2018), Kato et al. (2014), Wolf et al. (2013), andNomoto et al. (2007) are shown in the bottom panel of Figure 12.The differences in these limits are primarily because of the different techniques employed in modelling.
We infer from Figure 12 (bottom panel) that a WD with mass below 1.30 M ⊙ has a high accretion rate and does not allow the RN phenomenon to occur at the rate of once a year.
We also see that 1.32 − 1.36 M ⊙ WDs do fall into the "H-shell flash" region of some of the models, and any WD with M WD > 1.36 M ⊙ satisfies the necessary criteria of nova eruption ( Ṁ < Ṁst ) for all the models.Thus, the WD mass in M31N 2008-12a is likely to be greater than 1.36 M ⊙ , which would, in turn, allow "H-flash" features at the observed recurrence period of M31N 2008-12a.
However, it should be emphasized that Starrfield (2017) generated models using MESA, which allow for the stability line (and the critical line) to exist for less massive WDs but not for higher mass WDs (∼ 1.35 M ⊙ model shown in Figure 12).It was shown that such massive WDs would show H-flashes initially, followed by He-flashes, and ultimately grow to M Ch .
Strikingly, none of the models could predict the accretion rate derived from accretion disk modelling discussed in Darnley et al. (2017c).The accretion rate was found to vary during eruption, SSS, and quiescence phases.In Figure 12, we show the range of accretion rate onto the WD during quiescence.It may hint that Starrfield's models, which do not predict any upper limit on stable accretion, are better suited for such massive systems.But at the same time, the accretion rates generated by both types of models (with and without an upper limit) are insufficient to match the ones derived from observations of this exceptional RN.The cusp around maxima in the light curves has been suggested to have different origins.It could be due to shock from a secondary ejection (Kato et al. 2009), polar outflow along the line-of-sight, or ejecta-donor interaction (Darnley et al. 2018b).Observational evidence of broad-winged emission features supports the presence of a fast-moving component in the ejecta.Modelling the Hα (and Hβ for 2018) line profiles using Shape could also generate these fast-moving polar ejecta close to the line-of-sight.Darnley et al. (2017b) proposed the presence of these jets using HST data, though they did not confirm it.The photometric and spectroscopic evidence combined with Hα profile modelling presented in this work further strengthens the claim of polar jets emanating close to the line-of-sight, indicating a low-inclination angle for this system.
Optical imaging and spectroscopic modelling of RS Oph revealed the ejecta to be bipolar, possibly associated with jets (Bode et al. 2007;Ribeiro et al. 2009).RS Oph jets were confirmed in radio wavebands (O'Brien et al. 2006;Rupen et al. 2008;Sokoloski et al. 2008;Munari et al. 2022).High-velocity components of emission lines are standard signatures of jets and were observed in V1494 Aql (Iijima & Esenoglu 2003), U Sco (Kato & Hachisu 2003), V6568 Sgr and YZ Ret (McLoughlin et al. 2021a).McLoughlin et al. (2021b) pointed out that jets are usually associated with fast novae, which are essentially linked to massive WDs.M31N 2008-12a falls perfectly into these categories.One possible mechanism for jet formation could include bipolar winds during the nova outburst and subsequent mass ejection into an asymmetric medium.This case is particularly interesting for RNe where the asphericity left behind from previous eruptions triggers material to escape through certain channels.Magnetic field lines near the WD and the accretion disk could also lead to the collimation of wind perpendicular to the disk (Ogilvie & Livio 2001).The accretion disk of M31N 2008-12a is known to be luminous (Darnley et al. 2017c).Such bright disks can also give rise to supersonic winds forming jets (Fukue 2002).

Light curve
The optical light curves from 2017 to 2022 are similar.A sharp linear decline is seen from day 1 since the maximum, followed by an approximately flat but jittery plateau and then the final decline ensues.The evolution is similar to the past eruptions and is close to the light curves of P-class recurrent novae (Strope et al. 2010).
The UV peak is observed before the optical peak in all the eruptions.The UV light curve shows a decline from peak magnitude until the onset of the plateau phase, followed by multiple jitters.The UV plateau phase is consistent with the SSS phase's turn-on time.A flat decline follows it and ultimately ends with a brief period of brightening.The 2016 uvw2 light curve shows considerable deviation from the other eruptions.
The SSS phase turns out to be similar in all the eruptions except for the 2016 eruption, where it ended as early as ∼ 16 days from the eruption.The SSS temperature is strongly correlated to the soft X-ray flux.There is a significant drop in X-ray flux during the 2017 − 2022 eruptions around the same time that was noted in the previous eruptions on day 11 since the eruption.The cause of this drop in flux is yet to be explored.
During the slow decline phase, the 2016 uvw2 light curve (Figure 4) deviated from the general trend.Here, we also point out that for the first time in 2016, detailed uvw2 observations were conducted, and the light curve was found to be similar to the 2015 uvw1 trend.Based on this, Henze et al. (2018e) concluded that the optical and UV evolution in the 2016 eruption were similar to the previous ones, while the peculiarity of the 2016 eruption was reflected only in the X-rays.However, Figures 2 and 4 show that the evolution of uvw2 flux in 2016 differs from all other (subsequent) eruptions.The 2016 uvw2 light curve is fainter, similar to what was also noted for its soft X-ray counterpart ( §5.1).

Decreasing accretion rate
In § 7, we found that the recurrence period shows an increasing trend with time.A modest increase in the recurrence period would imply that either the accretion rate or the WD mass is decreasing over the years (Hillman et al. 2016;Wang 2018).Light curve models provided by Kato et al. (2015) suggested the WD to be as massive as 1.38 M ⊙ , accreting at a rate of 1.3 × 10 −7 M ⊙ yr −1 .Whereas Darnley et al. (2017c) modelled the quiescent phase accretion disk using HST data and showed that the rate of mass accretion could be even higher at (0.6 − 1.4) × 10 −6 M ⊙ yr −1 considering a 50% efficiency.On the other hand, the ejecta masses during each nova cycle (see § 5.1 and § 4.2) are ∼ 10 −8 M ⊙ .Thus, the net mass lost during each eruption is always less than the total mass gained between each eruption.As a result, the WD is growing in mass with time.The absence of Ne lines in the spectra indicates it to be a CO WD.Such WDs can reach > 1.36 M ⊙ only by accreting material.These inferences rule out the possibility of an increasing recurrence period due to decreasing WD mass.We suspect that the accretion rate has been slowly declining over the years (see Figure 11), lengthening the time taken to reach the critical conditions required for thermonuclear runaway reactions to be initiated on the surface of the WD.The following could cause a gradual decrease in the accretion rate: 1.The presence of starspots and increased activity in the secondary (Henze et al. 2018e).
2. The companion slowly running out of gas by supplying material to power the "H flashes" for millions of years (Darnley et al. 2019c).
3. Orbital dynamics can also change accretion rates, especially in violent systems like M31N 2008-12a, where nova eruptions are frequent.
4. Extent of the destruction of the accretion disk during each nova eruption decides the time taken to reform the accretion disk and resumption of accretion.Delayed accretion can also lead to a slowing down of the recurrence period.1.The linear decline post-maximum in the optical light curves is similar to that of the previous eruptions.The evolution of the UV light curve in the 2017 − 2022 eruptions is also similar to the previous eruptions.A rapid decline since the maximum is followed by a plateau phase coincident with the SSS turn-on time.It then follows a secular decline with undulations before dimming beyond the detection limit.A UV rebrightening is also seen towards the end of the SSS phase.
2. The SSS phase features are consistent with previously reported values.The mean SSS turn-on time and turn-off time are 5.1 ± 0.6 days and 16.9 ± 1.0 days since the eruption, respectively.The SSS phase shows X-ray variability, the most prominent being the dip ∼ 11 days after the eruption.
3. The UV and soft X-ray flux are found to be anticorrelated at the peak of the SSS phase, which has not been reported before.This implies that both originate at the surface of the WD and could arise during the reformation of the partially disrupted accretion disk.
4. Balmer, He, and N lines dominate the optical spectra.Hα velocities decelerate from ∼ 5000 km s −1 within 1 day of eruption to ∼ 2000 km s −1 at around 4 days after eruption, consistent with phase II of shock remnant development.
5. The ejecta mass derived from t on and Cloudy modelling is of the order of 10 −8 − 10 −7 M ⊙ , which is consistent with previous estimates derived using different techniques.Compared to the accretion rates derived in this work and previous studies, the ejecta mass is lower than the average mass accreted in a year, implying the WD is potentially increasing its mass.
6.He abundance in the ejecta was found to be high at He/He ⊙ ∼ 2.5 − 3.1, as is the case for most RNe.
7. Hα line morphology indicates an ejecta with an equatorial ring, a slow bipolar conical component, and an extended fast bipolar component along the line-of-sight resembling a jet-like structure.Evidence of a cuspy feature in the light curves near the peak is seen as a general trend after the 2016 eruption in the r ′ and i ′ bands.Together with emission-line modelling, we conjecture the cusp is caused by jets present in the ejecta.The presence of jets in this system was suspected, and we provide strong evidence for its presence here.Such jets could be common in novae systems, especially RNe.
8. We noticed that the recurrence period shows a weak tendency to increase with time, a sign of decreasing accretion rate.9.By comparing the recurrence period with binary evolution models, the mass of the WD is constrained to be > 1.36 M ⊙ .However, we emphasize that none of the models could replicate the observed accretion rate determined in previous studies.Irrespective of that, a CO WD near the M Ch and growing in mass is a good candidate for the single degenerate channel of Type Ia supernova explosions.
Science and Technology, Government of India.

Figure 1 .
Figure 1.Optical light curves of M31N 2008-12a for 2017 − 2022 eruptions.GIT, HCT, and JCBT observations are plotted with publicly available data.Vertical dashed lines in each panel mark the epochs of spectroscopic observations.

Figure 2 .
Figure 2. UVOT light curves of M31N 2008-12a shown for 2013 − 2022 eruptions.UVIT F148W data are overplotted.SSS ton and t off times are highlighted in yellow and red respectively.

Figure 4 .
Figure 4. Light curve properties of 2013-2022 (except 2016) eruptions for uvw2, g ′ , r ′ , and i ′ bands (clockwise from top-left).Outlier data points of 2016 are marked in red.The light curve template, in orange, was generated by stitching the phases in each filter.GP mean and 3σ error functions are shown in blue.
Spectral analysisThe optical spectra with a good SNR are shown in Figure5with important emission features marked.Spectra taken on 2021 Nov 14.8 UT and 2022 Dec 3.73 UT were noisy and have not been used for analysis.All the spectra have been dereddened using E(B − V ) = 0.10 (Darnley et al. 2017b).The spectra were taken within the first three days of the eruption and depict a blue continuum with Hydrogen Balmer and He I (4471 Å, 5876 Å, 6678 Å) emission lines.Some epochs also show the He II (4686 Å) and the N III lines (∼ 4640 Å).Based on a multi-eruption combined spectrum,Darnley et al. (2016) identify several other features in the spectrum.

Figure 6 .
Figure 6.Top (from left): Best-fit Cloudy modelled spectrum (red) overplotted on the observed spectrum (grey) and smoothed spectrum (black dotted lines) of the 2018 eruption.Best-fit Hα and Hβ velocity profiles (red) overplotted on the observed profile (dotted line).Morphology of the ejecta of 2018 eruption obtained from Shape using Hα and Hβ velocity profiles.X-axis is the line-of-sight direction, Y-axis is perpendicular to the plane of sky and line-of-sight.Red represents velocity away from us and violet towards us (2500 to -2500 km s −1 ).Bottom (from left): Same as top panel, but for 2019 eruption.Velocities are color coded from 6000 km s −1 (red) to -6000 km s −1 (violet).Note that the scales are different in the ejecta morphology plot.

Figure 8 .
Figure 8. Top panel: Combined XRT data of the SSS phase during 2013-2022 (except 2016) eruptions.Overplotted are data points binned at 0.5 days, its quadratic fit in black, GP regression and its corresponding 3σ error region in blue.The deviations from the simple quadratic fit are shown below it.Bottom panel: Temperature evolution from XRT data during the SSS phase of 2017 − 2022 eruptions.SXT data points for 2020 and 2021 are over-plotted.Mean turn-on and turn-off times are marked in yellow and red, respectively.

Figure 9 .
Figure 9. XRT and SXT spectra of 2020 (top panel) and 2021 (bottom panel) eruptions.XRT data nearest to the SXT observation dates have been used for better comparison.The exact observation epochs are given in the legends (in days since the eruption).Spectral fitting involved a single black body with ISM absorption models for both SXT and XRT.

Figure 11 .
Figure 11.Left: Distribution of the frequency of eruptions as a function of 'days since the last eruption' binned at 20 days.Black histogram is for 2008 − 2015 eruption data, and black solid line represents the KDE for this data.Red histogram and red solid line represent the same for 2008 − 2022 data.Right: Eruptions marked with stars.Dashed lines and coloured areas represent the mean and error functions of the two models tracking the recurrence period with time.Shaded region is the projection for the 2023 − 2025 eruptions, given we don't see outlier events.Secondary Y-axes represent the accretion rates for 1.38 M⊙ and 1.40 M⊙ WD corresponding to the recurrence period on the primary y-axis.
Figure 11.Both the models could reasonably anticipate the 2023 eruption.However, the GPR model can pick up any underlying data trends and better constrain the change in accretion rate.
In a theoretical study to understand the possible mass growth of a WD accreting matter from a non-degenerate companion,Hillman et al. (2016) have explored a range of accretion rates and derived limits on the accretion rate and on the initial mass that will allow a WD to reach the Chandrasekhar limit.Adopting their relation between D (period) and Ṁ (accretion rate) for hydrogen accretion cases, log Ṁ = − A log D − B, we estimate the accretion rates in the last 15 years for the M31 RN eruptions.Here, the coefficients A and B depend on the WD mass.For each value of M WD , A and B were determined by fitting a linear function to the parameter space of log D and log Ṁ .The accretion rates of WD masses between 1.20 M ⊙ and 1.40 M ⊙ corresponding to the periods of M31N 2008-12a are shown in the top panel of Figure 12.Wang (2018) have used the Modules for Experiments in Stellar Astrophysics (MESA) to model the binary evolution of WD accreting H-rich material from a companion for a range of WD masses and accretion rates.
Figure 12.Top: Observed recurrence period plotted against accretion rate for different WD masses."W", "H" and "S" in the legend indicate accretion rates derived from Wang (2018), Hillman et al. (2016) and Starrfield (2017) respectively.Bottom: MWD− Ṁ parameter space for M31N 2008-12a.Accretion rates from the top panel are shown in red (Wang's), blue (Hillman's) and yellow (Starrfield's).Overplotted in green markers are for the one-year recurrence period taken from the literature.Stable and critical accretion rate limits from four studies are shown as horizontal tracks.The accretion rate range obtained by Darnley et al. (2017c) during quiescence is shown in purple.
More evidence of jets?

5 .
A third body orbiting the M31N 2008-12a CV could perturb the binary motion, changing the accretion rate.Triple systems are known to produce exotic binaries; one such example is T Pyx(Knigge et al. 2022).9. SUMMARY This paper presents the evolution of 2017 −2022 eruptions of M31N 2008-12a in different wavelengths.The main results are summarised as follows.

Table 2 .
Optical spectroscopic and photometric observations of 2018-2022 eruptions of M31N 2008-12a.This is an example table.The full table will be made available in the online format

Table 3 .
AstroSat (UVIT and SXT) and Swift (UVOT and XRT) observations of 2017 − 2022 eruptions of the M31N 2008-12a.The complete table will be made available in online format.

Table 5 .
Flux and FWHM velocities of identified lines in the spectra It is located at the Indian Astronomical Observatory (IAO, Hanle).We acknowledge funding by the IITB alumni batch of 1994, which partially supports the operation of the telescope.Telescope technical details are available at https://sites.google.com/view/growthindia/.This work uses the SXT and UVIT data from the AstroSat mission of the Indian Space Research Organisation (ISRO).We thank the AstroSat TAC for allowing us ToO time to observe this nova from 2019-2022.We thank the SXT and UVIT payload operation centres for verifying and releasing the data via the ISSDC data archive and providing the necessary software tools.We acknowledge the use of public data from the Swift data archive.This work has also used software and/or web tools obtained from NASA's High Energy Astrophysics Science Archive Research Center (HEASARC), a service of the God-dard Space Flight Center and the Smithsonian Astrophysical Observatory.Kulinder Pal Singh (KPS) and G.C. Anupama (GCA) thank the Indian National Science Academy for support under the INSA Senior Scientist Programme.Facilities: HCT:2m, JCBT:1.3m,GIT:0.7m,Swift (UVOT and XRT), AstroSat (UVIT and SXT) Software: CCDLAB (Postma & Leahy 2017), IRAF v2.16.1 (Tody 1993), HEASOFT v6.25, XIMAGE v4.5.1, and XSELECT v2.4e (Nasa High Energy Astrophysics Science Archive Research Center (Heasarc) 2014), XSPEC v12.12.0 (Arnaud 1996), Python v3.6.6 (Van Rossum & Drake 2009), NumPy (Harris et al. 2020), SciPy (Virtanen et al. 2020), Pandas (McKinney et al. 2010), Matplotlib (Hunter 2007), scikit-learn (Pedregosa et al. 2011), astroquery (Ginsburg et al. 2019) APPENDIX