Multiwavelength Observations of Multiple Eruptions of the Recurrent Nova M31N 2008-12a

Photometric and spectroscopic observations were carried out using the following telescopes and instruments. The log of optical observations is given in Table 2.

Table 2. Optical Spectroscopic and Photometric Observations of 2018–2022 Eruptions of M31N 2008-12a

Date	Telescope	Instrument	Grating	Exp
(UT)				(s)
2018 Nov 7.8	HCT	HFOSC	Gr7	2700
2018 Nov 8.6	HCT	HFOSC	Gr7	3600
2019 Nov 7.5	HCT	HFOSC	Gr7	3000
2020 Oct 31.6	HCT	HFOSC	Gr7	3600
2021 Nov 14.8	HCT	HFOSC	Gr7	2100
2021 Nov 15.6	HCT	HFOSC	Gr7	3600
2022 Dec 3.7	HCT	HFOSC	Gr7	3600
2022 Dec 4.5	HCT	HFOSC	Gr7	3600

Date	Telescope	Filter	Magnitude	Exp
(UT)				(s)
2018 Nov 7.58	JCBT	B	19.29 ± 0.28	1800
2018 Nov 7.57	JCBT	V	18.78 ± 0.10	1200
2018 Nov 8.57	HCT	V	19.27 ± 0.13	300 × 3
2018 Nov 8.56	HCT	R	19.11 ± 0.11	180 × 3
2018 Nov 8.55	HCT	I	18.83 ± 0.20	150 × 3
2018 Nov 8.76	GIT	$g^{\prime}$	19.60 ± 0.10	300 × 3
2018 Nov 8.77	GIT	$r^{\prime}$	19.15 ± 0.09	300 × 3
2018 Nov 8.74	GIT	$i^{\prime}$	19.28 ± 0.15	300 × 3

Note. The full table of the photometric data will be made available in machine-readable format. The spectroscopic data are not included in the machine-readable table.

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

2.1.1. GROWTH-India Telescope

2.1.2. Himalayan Chandra Telescope

2.1.3. J.C. Bhattacharyya Telescope

2.1.4. Other Data Sources

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

2.2.1. Swift UVOT

2.2.2. AstroSat UVIT

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.

2.3.1. Swift XRT

2.3.2. AstroSat SXT

For a very fast RN, like M31N 2008-12a, a tight constraint on the eruption time is useful for generating light-curve models (Section 3.3) and studying its recurrence nature (Section 7). Hence, we estimate the epochs of eruption based on available detection and prediscovery magnitudes, and the nondetection 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 nondetection 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 nondetection.

The 2017 eruption was discovered just in time to be called the "2017 eruption" on December 31.77 UT by Boyd et al. (2017). Darnley et al. (2017a) reported spectroscopic confirmations on the same day. The last nondetection was on December 31.38 UT at an upper limit of m_clear = 19 mag (Naito et al. 2018).

The 2018 eruption was discovered on November 6.80 UT at a magnitude of 19.15 ± 0.05 by 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 November 6.54 UT in a clear filter.

The 2019 eruption was detected on November 6.71 UT by Oksanen et al. (2019) at 19.40 mag. The first spectrum taken on November 6.83 UT (Darnley et al. 2019d) confirmed the recurrence of M31N 2008-12a. The last nondetection 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 October 30.89 UT. The nova was, however, also detected in images captured 90 minutes before the discovery (Galloway et al. 2020), and we use the prediscovery detection (m_g = 18.74) and nondetection (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 November 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 mag on November 13.96 UT.

The 2022 eruption was discovered on December 2.83 UT by Perez-Fournon et al. (2022). It was undetected until December 2.61 UT at m_L > 19.60 mag (Shafter et al. 2022). Spectra taken on December 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.

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–4 days in $g^{\prime} r^{\prime} i^{\prime}$ bands. A brief description of the optical light curve for each eruption during 2017–2022 is provided below.

Zoom In Zoom Out Reset image size

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^{\prime}$ band (0.84 ± 0.12 mag day⁻¹) compared to the $g^{\prime}$ band (0.74 ± 0.19 mag day⁻¹). However, the decline rates in $g^{\prime}$ and $r^{\prime}$ 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 BVRI filters, and after a brief halt for about 0.3 day at this level, a marginal increase in the brightness by 0.2–0.3 mag lasting ∼0.7 day 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^{\prime}$ (0.68 ± 0.03 mag day⁻¹) and $i^{\prime}$ (0.53 ± 0.04 mag day⁻¹) bands was slower by 0.1–0.2 mag day⁻¹ compared to other years. Further, the decline rates of $r^{\prime}$ and $i^{\prime}$ are slower than $g^{\prime}$ band (0.84 ± 0.02 mag day⁻¹) 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 day from discovery. The peak R magnitude reached in 2019 is m_R = 18.12, which is higher than most other eruptions.

The $g^{\prime}$ band traces the rise of the 2020 light curve at ∼1.2 mag day⁻¹ while the $r^{\prime}$ band traces the smooth decline from the peak at 0.89 ± 0.05 mag day⁻¹ for 2 days. Limited data only restricts the light-curve analysis to the $r^{\prime}$ band.

The 2021 eruption was caught almost a day before it reached its peak. The rise was sharper in the $i^{\prime}$ band compared to the $g^{\prime} r^{\prime}$ bands. It declined rapidly in the $g^{\prime}$ band at 1.25 ±0.20 mag day⁻¹ but relatively slowly in the $r^{\prime}$ band at 0.88 ± 0.13 mag day⁻¹ for the first 3–4 days. The decline rate then slowed with significant enhancement in the $r^{\prime}$-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^{\prime}$ band with a rate of ∼2 mag day⁻¹, the fastest in the last 6 yr. It then declined speedily in the $g^{\prime}$ band (1.01 ± 0.05 mag day⁻¹) but relatively slowly in $r^{\prime}$ band (0.79 ± 0.03 mag day⁻¹) and even more slowly in the $i^{\prime}$ band (0.69 ± 0.06 mag day⁻¹).

Generally, the nova declines fastest in the $g^{\prime}$ 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⁻¹ 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 space-based 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 rebrightening 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 decline of ∼0.15 mag day⁻¹. 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.

Zoom In Zoom Out Reset image size

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 6 yr 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 Section 8.

The (F148W − uvw2) color was determined from observations taken on the same day. From Figure 3, it is seen that the (F148W − uvw2) color becomes bluer at the onset of the SSS phase but is significantly redder during the SSS phase.

Zoom In Zoom Out Reset image size

In the optical bands, we restrict the color analysis to only the SDSS primed filters to avoid instrumental and/or filter dependencies of the Bessel filters. Near-simultaneous observations in $g^{\prime} r^{\prime}$ and $r^{\prime} i^{\prime}$ filters for the same eruption were used to estimate the colors. The colors are plotted together as a function of days since the eruption to bring all the outbursts to the same timescale. The ($g^{\prime} -r^{\prime}$) color 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^{\prime} -i^{\prime}$) color 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^{\prime} -r^{\prime}$) color becomes bluer while the ($r^{\prime} -i^{\prime}$) color 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^{\prime} -i^{\prime}$) color becomes redder when the nova enters its plateau phase. A similar trend was also noted by Darnley et al. (2016) in their color plots.

To understand the temporal evolution of the light curves, we model them by breaking them into three phases corresponding to different decline rates. The ${uvw}2,g^{\prime} ,r^{\prime} ,\ \mathrm{and}\ i^{\prime}$ bands 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 }\leqslant t\leqslant 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 four-phase division of all light curves to analyze 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 Section 1 for references) with those from 2017 to 2022 and generated overall light-curve models spanning from 2013 to 2022. We note here that the data in $g^{\prime}$ are 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.

Zoom In Zoom Out Reset image size

Table 4. Multi-eruption Light-curve Model Parameters

Filter	${t}_{\max }$ a	${m}_{\max }$	t2	t3	Decline Rates (​​​​mag day−1)
	(days)	(AB)	(days)	(days)	${t}_{\max }\,\mathrm{to}\,3.5$	3.5 to ∼ 8	>8
2017–2022
uvw2	0.66 ± 0.26	18.82 ± 0.30	4.41 ± 0.51	14.15 ± 0.41	0.93 ± 0.03	0.17 ± 0.02	0.09 ± 0.01
g'	0.86 ± 0.04	18.27 ± 0.16	2.13 ± 0.04	⋯	0.90 ± 0.02	⋯	⋯
r'	0.86 ± 0.09	18.21 ± 0.27	2.37 ± 0.05	5.87 ± 5.61	0.88 ± 0.02	0.19 ± 0.18	⋯
i'	0.88 ± 0.12	18.30 ± 0.62	3.26 ± 0.20	⋯	0.56 ± 0.04	⋯	⋯
2013–2022 (Excluding 2016)
uvw2	0.66 ± 0.26	18.82 ± 0.30	4.41 ± 0.51	14.15 ± 0.41	0.93 ± 0.03	0.17 ± 0.02	0.09 ± 0.01
g'	0.86 ± 0.04	18.27 ± 0.16	2.13 ± 0.04	⋯	0.90 ± 0.02	⋯	⋯
r'	0.97 ± 0.02	18.44 ± 0.03	2.29 ± 0.01	6.44 ± 0.08	0.97 ± 0.01	0.21 ± 0.02	⋯
i'	1.14 ± 0.06	18.63 ± 0.10	2.42 ± 0.02	10.83 ± 0.32	0.88 ± 0.01	0.10 ± 0.01	⋯

Note.

^a Days since eruption.

Download table as:  ASCIITypeset image

Additionally, Gaussian process regression (GPR) 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.

3.3.1. The Rise to Peak

3.3.2. The Initial Steep Decline

3.3.3. The Slow Decline

The optical spectra with good S/N are shown in Figure 5 with important emission features marked. Spectra taken on 2021 November 14.8 UT and 2022 December 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 3 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 multieruption combined spectrum, Darnley et al. (2016) identify several other features in the spectrum. While these features are not seen in the individual spectra presented here, a few faint features could be identified in the merged spectra taken at similar epochs after outbursts. He i 4922 Å and N ii 5679 Å were detected in the combined spectra of 2018 (1.14 days) and 2020 (1.11 days), and He ii 4686 Å, He i 5016 Å, N ii 5679 and 6346 Å, and Raman O vi 6830 Å could be identified in the late phase merged spectra of 2018 (1.94 days), 2022 (2.03 days), and 2016 (2.23 days).

Zoom In Zoom Out Reset image size

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 deconvolving with the width of night skylines. The initial velocities within 1 day of eruption are as high as 5000 km s⁻¹, typical of very fast novae. These observations are consistent with the previous eruptions of M31N 2008-12a (Darnley et al. 2015c, 2016; Henze et al. 2018e). At around ≳1.5 days after the eruption, the emission-line widths narrow to 3000 km s⁻¹. 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 }\propto {t}^{-0.27\pm 0.07}$. This is similar to the estimate provided by Darnley et al. (2016) (${v}_{\exp }\propto {t}^{-0.28\pm 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 2 days after the eruption, the double-peaked profiles give way to a relatively narrow boxy profile.

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 S/N spectra taken on 2018 November 8.6 UT and 2019 November 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⁵ K and 10³⁷ erg s⁻¹ respectively. A clump component and a low-density diffuse component of density 10¹¹ and 10¹⁰ cm⁻³ respectively were used to fit the emission lines in the observed spectrum. The ejected mass and helium abundance from the best-fit modeled spectrum were found to be 7.21 × 10⁻⁸ 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⁴ K and 10³⁷ erg s⁻¹ respectively. A clump component of 2.24 × 10¹⁰ cm⁻³ and a diffuse component of 1.26 × 10⁸ cm⁻³ 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⁻⁸ M_⊙ and 3.09 ± 0.18 He_⊙ respectively. The ejected mass derived from X-rays (see Section 5.1) and spectral modeling are similar to that reported for the 2015 eruption by Darnley et al. (2016). An 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).

Zoom In Zoom Out Reset image size

It was noted that the two-component model was insufficient to generate the synthetic spectrum with a high χ² 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 S/N, modeling with three components is beyond the scope of this work. With these uncertainties, modeling a high S/N 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 morpho-kinematic 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 8075 ± 121 could generate the synthetic velocity profile of 2018 spectrum. The extended bipolar component stretched up to 4.52 × 10¹² cm along the ejecta axis from the center while the central bipolar cones (opening angle of ∼91°) extended up to 3.62 × 10¹¹ cm. The inner radius of the equatorial ring and radii of the bipolar cones were 1.27 × 10¹¹ cm and 5.42 × 10¹¹ cm, respectively. A similar geometry with a best-fit inclination angle of 7960 ± 145 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¹³ cm and 6.16 × 10¹² cm along the ejecta axis from the center respectively. The inner radius of the equatorial ring and the radius of the bipolar cones were 4.12 × 10¹² cm and 5.27 × 10¹² 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 high-velocity hydrogen Balmer emission.

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⁻³ counts s⁻¹, which increases to (3–4) ×10⁻² counts s⁻¹ 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 nondetections 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), 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.

Zoom In Zoom Out Reset image size

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 general features of the SSS phase of M31N 2008-12a. Most of the variability is seen up to day 13; whereafter, the decline is relatively smooth.

Zoom In Zoom Out Reset image size

X-ray studies of M31 novae (Henze et al. 2010, 2011, 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⁻¹ around this phase gives an ejecta mass range of 4.2 × 10⁻⁸ < M_ej,H/M_⊙ < 10.2 × 10⁻⁸. These are slightly higher than that calculated from optical spectra in Section 4.2 but less than the average mass accreted in a year.

The XRT data were also used to extract spectra by merging two data sets obtained on consecutive days to increase the S/N. We fixed the H column density at 1.4 × 10²¹ cm⁻² (Darnley et al. 2016) but varied the blackbody temperature and normalization to attain the best-fit values. The time evolution of the SSS temperature is shown in Figure 8 for 2017–2022 eruptions. Not only do 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 Section 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 toward 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 are shown in Figure 9. Also shown are the contemporaneous XRT spectra obtained from merged snapshots of 2 successive days. The data have 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 supersoft 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 S/N. The best-fit blackbody temperatures from SXT spectra are slightly higher than the XRT data during similar times (Figure 8). The modeled 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.

Zoom In Zoom Out Reset image size

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 (gray 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 until 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 7 yr, we see a clear increasing trend in the recurrence period.

Zoom In Zoom Out Reset image size

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 were modeled using a Matern kernel with a typical length scale of 15 yr 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 were subsequently included in the training set to project the upcoming eruptions. Since there is a gap of around 15 yr between 1993 and 2008, we did not include the 1993 data point for our modeling, 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 Figure 11. Both of 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.

We emphasize here that our model is restricted to only the usual eruptions that follow the trend. An anticipation of outlier events, such as the 2016 eruption, is not feasible.

In a theoretical study to understand the possible mass growth of a WD accreting matter from a nondegenerate 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 $\dot{M}$ (accretion rate) for hydrogen accretion cases, $\mathrm{log}\,\dot{M}=-A\mathrm{log}\,D-B$, we estimate the accretion rates in the last 15 yr 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 $\mathrm{log}\,D$ and $\mathrm{log}\,\dot{M}$. 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.

Zoom In Zoom Out Reset image size

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.

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.

We plot the accretion rates thus obtained corresponding to the periods of M31N 2008-12a in the last 15 yr in the ${M}_{\mathrm{WD}}\mbox{--}\dot{M}$ parameter space in the bottom panel of Figure 12. For comparison, we overplot the results from Kato et al. (2014), Prialnik & Kovetz (1995), Wolf et al. (2013), and Tang et al. (2014), who predicted the WD mass and accretion rates for a recurrence period of 1 yr.

In the ${M}_{\mathrm{WD}}\mbox{--}\dot{M}$ plane, when the accretion rate surpasses the critical threshold ($\dot{M}\gt {\dot{M}}_{\mathrm{cr}}$), the WD exhibits an RG-like behavior, 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., ${\dot{M}}_{\mathrm{cr}}\gt \dot{M}\gt {\dot{M}}_{\mathrm{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., ${\dot{M}}_{\mathrm{st}}\gt \dot{M}$, the accretion rate is insufficient to sustain continuous hydrogen burning. Systems within this parameter range experience H-shell flashes or nova eruptions, a category that includes all RNe, including M31N 2008-12a.

The limits on ${\dot{M}}_{\mathrm{cr}}$ and ${\dot{M}}_{\mathrm{st}}$ taken from Wang (2018), Kato et al. (2014), Wolf et al. (2013), and Nomoto 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 modeling.

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 ($\dot{M}\lt {\dot{M}}_{\mathrm{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 modeling 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. Modeling 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, although they did not confirm it. The photometric and spectroscopic evidence combined with Hα profile modeling 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 modeling 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 wave bands (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).

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 (Section 5.1).

In Section 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⁻⁷ M_⊙ yr⁻¹. Whereas Darnley et al. (2017c) modeled 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⁻⁶ M_⊙ yr⁻¹ considering a 50% efficiency. On the other hand, the ejecta masses during each nova cycle (see Sections 5.1 and 4.2) are ∼10⁻⁸ 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.  
  the orbital dynamics can also change accretion rates, especially in violent systems like M31N 2008-12a, where nova eruptions are frequent;
• 4.  
  the 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;
• 5.  
  a third body orbiting the M31N 2008-12a cataclysmic variable (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).