On the Observability of Recurrent Nova Super-Remnants

TL;DR

This study models the growth and observability of nova super-remnants around recurrent novae, revealing how environmental factors and white dwarf properties influence their size, structure, and potential detectability in different astrophysical contexts.

Contribution

It introduces a dynamic simulation of NSRs considering white dwarf evolution, expanding understanding of their formation, structure, and observational signatures around recurrent novae.

Findings

NSRs form as large, low-density cavities with hot ejecta and surrounding cool shells.

Higher environmental density and accretion rates limit NSR size.

Only NSRs around high accretion rate novae are currently observable.

Abstract

The nova super-remnant (NSR) surrounding M31N 2008-12a (12a), the annually erupting recurrent nova (RN), is the only known example of this phenomenon. As this structure has grown as a result of frequent eruptions from 12a, we might expect to see NSRs around other RNe; this would confirm the RN--NSR association and strengthen the connection between novae and type Ia supernovae (SN Ia) as NSRs centered on SN Ia provide a lasting, unequivocal signpost to the single degenerate progenitor type of that explosion. The only previous NSR simulation used identical eruptions from a static white dwarf (WD). In this Paper, we simulate the growth of NSRs alongside the natural growth/erosion of the central WD, within a range of environments, accretion rates, WD temperatures, and initial WD masses. The subsequent evolving eruptions create dynamic NSRs tens of parsecs in radius comprising a low-density cavity, bordered by a hot ejecta pile-up region, and surrounded by a cool high-density, thin, shell. Higher density environments restrict NSR size, as do higher accretion rates, whereas the WD temperature and initial mass have less impact. NSRs form around growing or eroding WDs, indicating that NSRs also exist around old novae with low-mass WDs. Observables such as X-ray and H$\alpha$ emission from the modelled NSRs are derived to aid searches for more examples; only NSRs around high accretion rate novae will currently be observable. The observed properties of the 12a NSR can be reproduced when considering both the dynamically grown NSR and photoionisation by the nova system.

Tables

Table 1. Table 1: Parameters for each run. Columns record the simulation number, initial WD mass, WD temperature, accretion rate, ISM density, spatial resolution, number of eruptions to grow the WD to M Ch subscript M Ch \text{M}_{\mathrm{Ch}} or for the simulation to reach the temporal upper limit of 10 8 superscript 10 8 10^{8} years, the cumulative time of the simulation, and the total kinetic energy released. Run 0 relates to the 10 5 superscript 10 5 10^{5} identical eruptions as modelled by DHO19 . Ejecta characteristics for Run 1 † used a broken exponential/linear interpolation (see Section  4.1 ). Runs 1 ⋆ superscript 1 ⋆ 1^{\star} , 2 ⋆ superscript 2 ⋆ 2^{\star} , 5 ⋆ superscript 5 ⋆ 5^{\star} and 7 ⋆ superscript 7 ⋆ 7^{\star} have the same ejecta characteristics as Runs 1, 2, 5 and 7, respectively, but do not include radiative cooling. Run 22 contains the same nova system as Run 1 but tuned with an ISM density of n = 1.278 𝑛 1.278 n=1.278 to match the ISM predicted in Section  6.3 for the reference simulation WD to grow a NSR to the size (67 pc) of the observed NSR around M 31N 2008-12a.

Run #	${\text{M}}_{\text{WD}}$	$T_{\text{WD}}$	$\dot{M}$	ISM density	Spatial resolution	Number	Cumulative time	Total Kinetic Energy
	(${\text{M}}_{\odot }$)	(K)	(${\text{M}}_{\odot }{\text{yr}}^{-1}$)	($1.67\times {10}^{-24}{\text{ g cm}}^{-3}$)	(AU/cell)	of eruptions	(years)	(erg)
0	n/a	n/a	$1.6\times {10}^{-7}$	1	0004	0,100,000	$1.0\times {10}^5$	$4.5\times {10}^{47}$
1	1	$1\times {10}^7$	$1\times {10}^{-7}$	1	0200	1,900,750	$3.1\times {10}^7$	$2.4\times {10}^{49}$
2	1	$1\times {10}^7$	$1\times {10}^{-7}$	0.1	0200	1,900,750	$3.1\times {10}^7$	$2.4\times {10}^{49}$
3	1	$1\times {10}^7$	$1\times {10}^{-7}$	0.316	0200	1,900,750	$3.1\times {10}^7$	$2.4\times {10}^{49}$
4	1	$1\times {10}^7$	$1\times {10}^{-7}$	3.16	0200	1,900,750	$3.1\times {10}^7$	$2.4\times {10}^{49}$
5	1	$1\times {10}^7$	$1\times {10}^{-7}$	10	0200	1,900,750	$3.1\times {10}^7$	$2.4\times {10}^{49}$
6	1	$1\times {10}^7$	$1\times {10}^{-7}$	31.6	0200	1,900,750	$3.1\times {10}^7$	$2.4\times {10}^{49}$
7	1	$1\times {10}^7$	$1\times {10}^{-7}$	100	0200	1,900,750	$3.1\times {10}^7$	$2.4\times {10}^{49}$
8	1	$1\times {10}^7$	$1\times {10}^{-8}$	1	0200	0,040,343	$1.0\times {10}^8$	$1.3\times {10}^{48}$
9	1	$1\times {10}^7$	$1\times {10}^{-8}$	10	0200	0,040,343	$1.0\times {10}^8$	$1.3\times {10}^{48}$
10	1	$1\times {10}^7$	$1\times {10}^{-8}$	100	0200	0,040,343	$1.0\times {10}^8$	$1.3\times {10}^{48}$
11	1	$1\times {10}^7$	$1\times {10}^{-9}$	1	0400	0,002,094	$1.0\times {10}^8$	$3.1\times {10}^{47}$
12	1	$1\times {10}^7$	$1\times {10}^{-9}$	10	0400	0,002,094	$1.0\times {10}^8$	$3.1\times {10}^{47}$
13	1	$1\times {10}^7$	$1\times {10}^{-9}$	100	4000	0,002,094	$1.0\times {10}^8$	$3.1\times {10}^{47}$
14	1	$3\times {10}^7$	$1\times {10}^{-7}$	1	0200	2,770,545	$4.1\times {10}^7$	$4.6\times {10}^{49}$
15	1	$5\times {10}^7$	$1\times {10}^{-7}$	1	0200	2,029,154	$2.7\times {10}^7$	$5.0\times {10}^{49}$
16	0.65	$1\times {10}^7$	$1\times {10}^{-7}$	1	0200	1,953,955	$3.7\times {10}^7$	$2.5\times {10}^{49}$
17	0.8	$1\times {10}^7$	$1\times {10}^{-7}$	1	0200	1,945,717	$3.6\times {10}^7$	$2.5\times {10}^{49}$
18	0.9	$1\times {10}^7$	$1\times {10}^{-7}$	1	0200	1,933,696	$3.4\times {10}^7$	$2.5\times {10}^{49}$
19	1.1	$1\times {10}^7$	$1\times {10}^{-7}$	1	0200	1,779,622	$2.2\times {10}^7$	$2.4\times {10}^{49}$
20	1.2	$1\times {10}^7$	$1\times {10}^{-7}$	1	0200	1,494,979	$1.0\times {10}^7$	$2.1\times {10}^{49}$
21	1.3	$1\times {10}^7$	$1\times {10}^{-7}$	1	0200	1,149,284	$3.7\times {10}^6$	$1.8\times {10}^{49}$
22	1	$1\times {10}^7$	$1\times {10}^{-7}$	1.278	0200	1,900,750	$3.1\times {10}^7$	$2.4\times {10}^{49}$
1†	1	$1\times {10}^7$	$1\times {10}^{-7}$	1	0200	2,591,344	$2.1\times {10}^7$	$5.3\times {10}^{49}$
1⋆	1	$1\times {10}^7$	$1\times {10}^{-7}$	1	0200	1,900,750	$3.1\times {10}^7$	$2.4\times {10}^{49}$
2⋆	1	$1\times {10}^7$	$1\times {10}^{-7}$	0.1	0200	1,900,750	$3.1\times {10}^7$	$2.4\times {10}^{49}$
5⋆	1	$1\times {10}^7$	$1\times {10}^{-7}$	10	0200	1,900,750	$3.1\times {10}^7$	$2.4\times {10}^{49}$
7⋆	1	$1\times {10}^7$	$1\times {10}^{-7}$	100	0200	1,900,750	$3.1\times {10}^7$	$2.4\times {10}^{49}$

Equations

$$M_{\text{WD},i+1}=M_{\text{WD},i}+\left(m_{\text{ig},i}\times {\eta }_i\right),$$

$$r_{\mathrm{out}}^3=\frac{3S_{\star }}{4\pi n(r{)}^2\beta (r,T)}+r_{\text{in}}^3,$$

Full text

††pubyear: 2023††pagerange: On the Observability of Recurrent Nova Super-Remnants–References

1 Introduction

Recurrent novae (RNe) are a subclass of the cataclysmic variables that experience repeated thermonuclear eruptions on timescales of a human lifetime. Like classical novae (CNe) – systems observed in eruption just once – RNe are interacting binary systems (Walker, 1954; Warner, 1995) containing a white dwarf (WD) and a main-sequence, subgiant, or red giant donor (Darnley et al., 2012). Hydrogen-rich material is expelled from the outer layers of the donor through stellar winds or Roche lobe overflow, following which it accumulates on the surface of the WD usually via an accretion disc. At the base of the accreted layer, compression and heating continually increase until the critical pressure for a thermonuclear runaway (TNR; Starrfield et al., 1972, 1976, 2020) is reached. Once degeneracy is lifted, the accreted envelope is driven upwards by radiation pressure and expands violently, with material travelling faster than the escape velocity of the WD ejected into the surrounding environment as the nova eruption (see, for example, Starrfield et al., 1976, 2020). Mass accretion then continues after (and possibly during; Kato et al., 2017; Henze et al., 2018) the eruption, leading to successive RN eruptions, separated by a recurrence period ($P_{\mathrm{rec}}$) which can vary.

Novae with carbon-oxygen WDs present a compelling single degenerate (SD) pathway to type Ia supernovae (SN Ia; Whelan & Iben, 1973; Hachisu et al., 1999a; Hachisu et al., 1999b; Hillebrandt & Niemeyer, 2000). Multi-cycle nova simulations (Yaron et al. 2005, hereafter Y05; Hachisu et al. 2007; Kato et al. 2015; Hillman et al. 2015, 2016; Starrfield et al. 2021) show that a substantial amount of accreted material is retained on the WD’s surface post-eruption, ultimately growing the WD to the Chandrasekhar (1931) limit (M${}_{\text{Ch}}$) in ${\sim}10^{7-8}$ years (Hillman et al., 2016). The other leading SN Ia pathway is the double degenerate (DD) scenario with two merging WDs (Webbink, 1984; Iben & Tutukov, 1984) yet within both the SD and DD pathways, novae are the brightest proposed progenitor, even at quiescence (Darnley, 2021). Therefore, extragalactic nova population studies can link environmental effects such as star formation and metallicity with SN Ia sub-classes. Alternatively, if the donor evolves such that no donatable material remains in the envelope, then the WD will cease growing and thereby never reach the M${}_{\text{Ch}}$, resulting in an extinct RN (Darnley, 2021).

A CN eruption will eject approximately ${\sim}10^{-4}$ $\text{M}_{\odot}$ of material into its surroundings with typical velocities ranging from a few hundred to several thousand km s*-1* (O’Brien et al., 2001). The interaction of ejecta with different velocities (Aydi et al., 2020b) will shock heat the gas leading to X-ray and radio emission such as that seen in RS Ophiuchi (Bode & Kahn, 1985; O’Brien et al., 1992) and V838 Her (O’Brien et al., 1994). This ejected material then goes on to form a nova shell (see, e.g., Woudt & Ribeiro, 2014; Harvey et al., 2020; Santamaría et al., 2020, 2022). For the ${\sim}10\%$ of Galactic novae with observed shells (see, e.g., Wade, 1990; Slavin et al., 1995; Gill & O’Brien, 1998; Santamaría et al., 2019; Santamaría et al., 2022), their morphologies can inform us of the underlying configuration of the binary. In particular, nova shells are structured with an equatorial waist and polar cones of emission (Hutchings, 1972). This structure forms from the originally near-spherically symmetrical nova ejecta interacting with the material in the orbital plane lost by the donor (see, e.g., Mohamed et al., 2013). Polar blobs, equatorial (and/or tropical) rings as well as knots are common to almost all nova shells; see, for example, DQ Her (Williams et al., 1978), HR Del (Harman & O’Brien, 2003), DO Aql and V4362 Sgr (Harvey et al., 2020) as well as V5668 Sgr (Takeda et al., 2022). In addition, due to the repeating nature of RNe, we have an example of interacting ejecta from successive eruptions producing clumping and shock heating around the RN T Pyxidis (Shara et al., 1997; Toraskar et al., 2013).

Even though the accretion disk surrounding the WD can be altered (Henze et al., 2018) to the point of removal in many cases (Drake & Orlando, 2010; Figueira et al., 2018), it will re-establish after the nova outburst (Worters et al., 2007) in preparation for future eruptions. Consequently, all nova systems are predicted to experience repeated outbursts with substantial variation in recurrence period between systems (Y05). Yet, only the recurrence periods for the known RNe, all contained within the Galaxy (10; Schaefer, 2010; Darnley, 2021), the Large Magellanic Cloud (4) and M31 (19; Darnley & Henze, 2020), have been determined, ranging from 98 years (Pagnotta et al., 2009) down to 1 year (Henze et al., 2015; Henze et al., 2018; Darnley & Henze, 2020). Such short inter-eruption intervals are powered by a combination of a massive WD and a high mass accretion rate (Starrfield et al., 1988).

The most rapidly recurring nova known is M31 N 2008-12a, or simply ‘12a’, (see, e.g., Darnley et al., 2016; Henze et al., 2018; Darnley & Henze, 2020; Darnley, 2021, and references therein). This extreme example erupts annually ($\overline{P}_{\text{rec}}=0.99\pm 0.02$ years; Darnley & Henze 2020) and has the most massive WD known ($\simeq 1.38\,\mathrm{M}_{\odot};$ Kato et al., 2015), likely CO in composition (Darnley et al., 2017a), accreting with a substantial mass accretion rate of ($0.6\leq\dot{M}\leq 1.4)\times 10^{-6}\,\mathrm{M}_{\odot}\,\mathrm{yr}^{-1}$ from a red giant (or clump) companion (Darnley et al., 2014, 2017b).

First associated with 12a by Darnley et al. (2015), this RN is surrounded by a vastly extended nebulosity. Compared to some of the largest Galactic CN shells known such as GK Persei (${\sim}0.5$ pc; Bode et al., 2004; Harvey et al., 2016b), Z Camelopardalis (${\sim}0.7$ pc; Shara et al., 2007) and AT Cancri (0.2 pc; Shara et al., 2012), 12a’s shell has semi-major and -minor axes of 67 and 45 pc, respectively, justifying a nova super-remnant (NSR; Darnley et al. 2019, hereafter DHO19) status. DHO19 ruled out the possibility of the shell being a SN remnant, a superbubble or a fossil H ii region with H$\alpha+$[Nii] imaging and deep low-resolution spectroscopy. Instead, the NSR’s existence was attributed to the cumulative sweeping up of ${\sim}10^{5-6}\text{ M}_{\odot}$ (DHO19) of local interstellar medium (ISM) from many previous nova eruptions.

To test the viability of a RN origin for 12a’s NSR, DHO19 utilised Morpheus (Vaytet et al., 2007) to perform 1D hydrodynamical simulation of $10^{5}$ 12a eruptions. Each of these eruptions ejected $5\times 10^{-8}\ \text{M}_{\odot}$ at a terminal velocity of 3000 km s*-1* over seven days, repeating every 350 days (DHO19). We assign the DHO19 simulation as Run 0 and it will be used as a comparison throughout this work.

Self- and ISM-interaction of the ejecta from each Run 0 eruption formed a huge cavity surrounded by an expanding shell with relative thickness of 22%. An unavoidable consequence of continual eruptions from a central system is the formation of a dynamical structure, be that a nova shell or larger remnant. However, the existence of a dynamical NSR does not necessarily signify a NSR that is observable. Nevertheless, the simulated dynamic remnant of Run 0 was found to be consistent with observations of the 12a NSR (DHO19).

A unique feature of a structure formed from repeatedly interacting eruptions is a continuously shock-heated region located inside the outer shell (DHO19). Extrapolating the growth rate from these simulations to the observed size of the super-remnant, DHO19 suggested an age of $6\times 10^{6}$ yrs. Importantly, the mechanism driving the NSR formation is also growing the 12a CO WD, which Darnley et al. (2017b) predict will surpass the Chandrasekhar limit and explode as a SN Ia in $<$ 20,000 years.

In this paper, we build upon the NSR hydrodynamic modelling presented by DHO19 through consideration of the complete eruption history of a nova system as the WD mass grows from its formation toward the Chandrasekhar mass. We also explore a number of factors, both intrinsic and extrinsic to the nova system that might impact NSR formation, to aid the search for more NSRs. This will be the first attempt to determine if the NSR associated with 12a is unique or whether it is simply the first of the phenomena to be found.

In Section 2 we describe the eruption model used to generate input parameters. We describe the Morpheus hydrodynamic code employed in this paper in Section 3 before outlining each of the separate runs of our main simulations. Various tests conducted after the main simulations are presented in Section 4. We explore the observability of NSRs in Section 5 by modelling emission from the simulations and then compare our simulations to observations of the 12a NSR in Section 6, before concluding our paper in Section 7.

2 Generating Nova Ejecta Properties

The DHO19 simulations of the 12a NSR utilised $10^{5}$ identical eruptions with a fixed recurrence period. While a good approximation for this system during its recent evolution, identical eruptions do not match the expected long term evolution of such a system, whereby the characteristics of the ejecta evolve with the changing WD mass. Therefore, to obtain the properties of a nova system with incrementally changing nova eruptions, we were required to grow a WD (see Section 2.2). We will only describe the model we used to grow the WD for a ‘reference simulation’ as an illustration, however this model was utilised for each of the different WD temperatures and accretion rates. As a reference simulation corresponding to the 12a system, we chose to grow a $10^{7}$ K WD with $\dot{M}=10^{-7}\text{ M}_{\odot}\,\text{yr}^{-1}$ (see Section 2.1 for details), which we then placed within an environment with a hydrogen-only ISM density of $1.67\times 10^{-24}\ \text{g cm}^{-3}$ (1 H atom per cubic centimetre). We refer to this ISM density throughout the paper by the number density $n=1\,\text{cm}^{-3}$ (but drop the units for clarity).

2.1 Parameter space

Y05 provides a parameter space for the characteristics of a nova envelope and the outburst characteristics for an extended grid of nova models with varying WD mass, temperature and accretion rate. This grid runs through all permutations of these parameters and outputs various eruption characteristics such as the mass accreted onto the WD which ignites during the TNR ($m_{\text{acc}}$), the mass ejected from the WD during the nova eruption ($m_{\text{ej}}$) and the duration of the mass-loss phase ($t_{\text{ml}}$) i.e., the timescale of each eruption.

For this study, we use values of $m_{\text{acc}}$, which we equate to the ignition mass ($m_{\text{ig}}$), $m_{\text{ej}}$, and $t_{\text{ml}}$ for WDs with masses 0.65, 1.0, 1.25 and 1.4 $\text{M}_{\odot}$111These WD masses were chosen from the set of WD masses given in Y05, as were the WD temperatures and accretion rates used in our study. We were limited to these accretion rates by the eruption models of Y05 whereby no $t_{\text{ml}}$ is provided for $\dot{M}=10^{-6}\,\mathrm{M}_{\odot}\,\mathrm{yr}^{-1}$., three temperatures of 10 MK, 30 MK and 50,MK, and three accretion rates of $10^{-7}\text{ M}_{\odot}\,\text{yr}^{-1}$, $10^{-8}\text{ M}_{\odot}\,\text{yr}^{-1}$ and $10^{-9}\text{ M}_{\odot}\,\text{yr}^{-1}$ (we consider three temperatures for $\dot{M}=10^{-7}\text{ M}_{\odot}\,\text{yr}^{-1}$ but only 10 MK for the other accretion rate values). To interpolate and extrapolate these points for a continuous set of values for our WD growth model, we required a function that evolved smoothly, behaved as a power law for lower masses, yet which became asymptotic as the Chandrasekhar mass was approached (see Section 4.1 for an alternative approach). The functions we fit to $m_{\text{ig}}$ and $t_{\text{ml}}$ are shown in Figure 1, as well as the continuous function for $P_{\text{rec}}$ (the ratio of $m_{\mathrm{ig}}$ and accretion rate). As we also wish to be consistent with observed characteristics of the nova eruption, we utilised observationally determined relations from Warner (1995) and Henze et al. (2014) to determine a function for the terminal ejecta velocity of the outburst.

2.2 Growing a white dwarf

We grew a $1\,\text{M}_{\odot}$ WD to a $\text{M}_{\text{Ch}}$ WD by accumulating the retained mass from iterated nova eruptions and using the interpolated relationships given in Section 2.1 to obtain properties of each eruption. For this example, a $1\text{ M}_{\odot}$ WD with a temperature of 10 MK experiences approximately 1,900,000 eruptions while growing from $1\text{ M}_{\odot}$ to $1.4\text{ M}_{\odot}$, reaching a recurrence period lower limit of ${\sim}282$ days. This WD mass upper limit of $1.4\text{ M}_{\odot}$ is assumed for all WD scenarios, which we equate to the Chandrasekhar mass ($\mathrm{M}_{\mathrm{Ch}}$) for this study.

A WD is grown (or eroded) according to the amount of accreted material retained (or removed) between eruptions. To model the evolution of the mass accumulation efficiency ($\eta$) over the evolution of a WD, we utilised the values of $m_{\text{ig}}$ and $m_{\text{ej}}$ from Y05 such that $\eta=(m_{\text{ig}}-m_{\text{ej}})/m_{\text{ig}}$ and interpolated between these points for a continuous set of values (see top right panel of Figure 1). The changing mass of the WD can thus be described as:

[TABLE]

where $M_{\text{WD},i}$ is the pre-eruption mass of the WD, $m_{\text{ig},i}$ is the mass accreted by the WD before the eruption, $\eta_{i}$ is the evolving mass accumulation efficiency, and $M_{\text{WD},i+1}$ is the post-eruption mass of the WD. With the initial WD mass being $1\,\text{ M}_{\odot}$, we utilised the relationships found in Section 2.1 to give the associated $m_{\text{ig}}$ value for equation 1. The post-eruption mass was then used as the $M_{\text{WD}}$ value in the next iteration and we continued this until we reached the limiting mass stated previously. We used the output parameters from this iterative model in our simulations. With each iteration, we were also able to use the relationships found in Section 2.1 to illustrate the evolution of a number of parameters including ejecta kinetic energy and momentum in terms of WD mass, recurrence period, elapsed time (from the first eruption), and the number of eruptions.

Utilising the WD growth model, we generated nova ejecta with incrementally changing properties. As the mass of the WD increases, eruptions become more frequent, and ejecta become less massive but with higher velocity in response to the increasing WD surface gravity.

3 Hydrodynamical Simulations

As the net mass loss rate from the WD varies as the WD mass grows, an analytic relation for the growth of the NSR shell cannot be derived. As such, full hydrodynamic simulations are a necessity if we are to understand the evolution of NSRs and their emission characteristics.

As in DHO19, the hydrodynamical simulations in this work were performed with Morpheus (Vaytet et al., 2007) – developed by the Nova Groups from the University of Manchester and Liverpool John Moores University. Morpheus brings together one-dimensional (Asphere; see Vaytet et al., 2007), two-dimensional (Novarot; see Lloyd et al., 1997) and three-dimensional (CubeMPI; see Wareing et al., 2006) codes to form an MPI-OpenMP Eulerian second-order Godunov simulation code that functions with Cartesian, spherical or cylindrical coordinates, and includes radiative cooling and gravity.

The configuration of the nova systems in this work are modelled in the same manner as given in DHO19 such that the mass donor is a red giant exhibiting a continuous wind mass loss rate (after accretion) of $2.6\times 10^{-8}\text{ M}_{\odot}$yr*-1* with a terminal velocity of 20 km s*-1*. These values are assumed to be consistent with the donor in the RS Ophuichi system (Bode & Kahn, 1985), thus are used as representative values with the red giant wind having negligible influence on the NSR evolution. The nova eruption is represented by an instantaneous increase in mass loss and ejecta velocity (the red giant wind’s contribution becomes negligible here) followed by a quiescent period in which only the red giant wind (with decreased mass loss and lower ejecta velocity) is present. Furthermore, unless otherwise stated, each ejection is modelled as a wind with a mass-loss rate and velocity that incrementally increase throughout the simulation as governed by the relationships determined from Y05 models (see Section 2 for details and Figure 1). The eruptions are separated by incrementally decreasing recurrence periods also governed by the aforementioned relationships. True nova ejecta are not spherically symmetric, however largely for computational reasons, we have assumed one-dimensional spherical symmetry for these simulations, effectively modelling the bulk equatorial ejecta (see, e.g., Mohamed et al., 2013). The spatial resolution of the full simulations ($\geq$200 AU/cell) is larger than the expected orbital separation of the WD and the donor (for example, the orbital separation for 12a is ${\sim}1.6$ AU; Henze et al., 2018) so we assume that both are located at the origin. Therefore, interaction between the ejecta and the donor or accretion disk is ignored.

Ideally, we would want to run each complete simulation at a high spatial resolution, however, this is not feasible with temporal and computing constraints. Running the reference simulation (see Section 3.2) several times with varying spatial resolution (and varying number of eruptions), we found that running its full 1,900,750 eruptions at 200 AU/cell would have the same long term structure as a simulation with resolution of 1 AU/cell (the resolution of a test run with 100 eruptions). Consequently, we set a spatial resolution of 200 AU/cell for most of our simulations, while those with lower spatial resolution (as indicted in Table 1) are set in response to the infrequency of eruptions, and therefore lessened impact on resolving the gross NSR structure, within those particular runs.

3.1 Incorporating radiative cooling

Nova ejecta lose energy through radiative cooling, which affects the evolution of any NSR. Therefore, the effects of cooling were tested in DHO19, with a NSR grown from $10^{3}$ eruptions with the inclusion of the radiative cooling module in Morpheus. The cooling model utilised in Morpheus was taken from Raymond et al. (1976, their Figure 1). The cooling rate is given as a function of gas temperature of an optically thin plasma, with no dust or molecules, made up of H, He, C, N, O, Ne, Mg, Si, S, Ca, Fe and Ni. Radiative cooling becomes ineffective below a temperature of $10^{4}$ K. Above $10^{8}$ K, the gas is ionised and only radiates through free-free Bremsstrahlung (Vaytet et al., 2007). Between these limits, cooling is dominated by line-cooling from the metals within the gas (Vaytet, 2009).

DHO19 demonstrated that there was no significant difference between the Run 0 NSR structure with or without cooling (see their Extended Data Figure 4). Cooling was suppressed in the Run 0 NSR as the recurrence period was much shorter than the cooling timescale. Hence, radiative cooling in the full simulation of Run 0 was not included.

In all cases, the NSR evolution presented in this work begins with high mass and low velocity ejecta (see Section 2) leading to less energetic eruptions and, crucially, with long gaps between consecutive eruptions. Therefore, at early times, the recurrence period will be longer than the cooling timescale and, as such, we incorporate radiative cooling in all simulations.

3.2 Reference simulation — Run 1

Our reference simulation, Run 1, models nova eruptions from a growing WD with a temperature $T_{\mathrm{WD}}=10^{7}$ K, with $\dot{M}=10^{-7}\,\text{M}_{\odot}\,\text{yr}^{-1}$, and within a low density ISM ($n=1$). With the varying mass accumulation efficiency, it would take ${\sim}$31 Myr (1,900,750 eruptions) for this WD to grow from $1\text{ M}_{\odot}$ to $\text{M}_{\mathrm{Ch}}$. Run 1 has a spatial resolution of 200 AU. This information, including the total kinetic energy released, is summarised in Table 1 for all simulations in this paper.

Run 1 is presented in Figure 2: the left-hand plot shows the density, pressure, velocity and temperature characteristics of the NSR after all ${\sim}$1,900,000 eruptions; the right-hand plot shows the evolution of the NSR shell outer edge and the inner edge, and the inner edge of the ejecta pile-up boundary (regions of the NSR are outlined in the top left panel of Figure 2).

In the top-left panel of the left-hand plot of Figure 2, we see that the inner and outer edges of the dynamical NSR shell extend to ${\sim}70.5$ and ${\sim}71.3$ pc, respectively – a shell thickness of 1.1%. As can be seen in the right-hand plot of Figure 2, shell thickness varies over the NSR evolution. For example, the shell compresses from $2.72\%$ ($P_{\text{rec}}=50$ years) to $1.14\%$ ($P_{\text{rec}}=1$ year) to $1.10\%$ ($P_{\text{rec}}=282$ days). At all times, this is much thinner than the 12a NSR shell (DHO19), which is $22\%$ from observations and remained at this thickness throughout Run 0 (see Figure 3). The shell thickness evolution during Run 1 is directly related to energy losses via cooling and to the evolution of eruption properties whereby the increasing frequency and kinetic energy of the ejecta drive a compression through the NSR shell.

In Run 1, the higher density found at the NSR shell inner edge ($n\,{\simeq}\,160$) compared to the outer edge ($n\,{\simeq}\,3$), seen in the top panel of Figure 3, is attributed to the contribution from the more recent, more frequent and more energetic eruptions – the rate of change of eruption properties surpasses the dynamic time-scale of the NSR shell at later times. The rate of propagation of the NSR shell into the surrounding ISM, and therefore the outer edge of the shell, remains largely based upon the combined properties of the entire eruption history, whereas the inner edge is shaped by newly arriving material.

As evident in the bottom left panel of the left-hand plot of Figure 2, the velocity of material in the inner cavity is high (${\sim}6.7\times 10^{3}$ km s*-1*) as it is essentially in free expansion. The velocity then drops substantially as the ejecta pile-up region is encountered, with the resultant shock-heating increasing temperatures by five orders of magnitude (see bottom right panel of the left-hand plot). The velocity and temperature in the ejecta pile-up region declines continuously out to the NSR shell as the ejecta encounter previously ejected material and reverse shocks (from the pile-up/inner shell boundary), with the cool outer edge expanding at a relatively low ${\sim}1$ km s*-1*.

Figure 2 provides a comparison between Run 1 and Run 1⋆ (with and without radiative cooling, respectively), to illustrate the significant difference in the NSR size and shell structure. The outer edge of the NSR in Run 1 extends to 71.3 pc yet, without radiative cooling in Run 1⋆, the NSR extends to ${\sim}$90 pc (having swept up around twice as much ISM). This substantial reduction in size can only be attributed to radiative losses within the NSR. Additionally, the radiatively cooled NSR shell from Run 1 is much thinner (${\sim}$1%) than the uncooled equivalent in Run 1⋆ (${\sim}$21%; see Figure 2). This results from the material in the early NSR shell losing energy via radiative cooling and therefore lacking the necessary pressure to maintain its size. This suppresses the early NSR shell formation such that when shell compression takes effect at later times (as increasingly energetic ejecta collide with the inner edge of the shell), the starting point is a thinner shell.

The NSR cavity and ejecta pile up boundary at ${\sim}10$ pc have similar density, pressure, velocity and temperature in Run 1 and Run 1⋆. At later stages, the increased frequency and energy of the eruptions results in the scenario that tends toward the Run 0 regime, whereby there is not enough time for the ejecta or remnant to cool radiatively between consecutive eruptions. Consequently, we see the effects of radiative cooling at the outer edge of the remnant, a relic of the earlier spaced out less energetic eruptions, and the centre of the NSR reflecting the later frequent eruptions. Furthermore, this point can be extended to all of the simulations conducted throughout this paper, whereby the growth and subsequent size of the nova super-remnant is shaped heavily by its early evolution.

So far, we have only considered the final epoch of Run 1, after the full 1,900,750 eruptions (Figure 2). However, to appreciate the changing structure and characteristics of the NSR, we have provided an animation of the Run 1 in Figure 4.

We illustrate in Figure 5 the spatiotemporal analysis of the evolution of the Run 1 NSR in terms of density, pressure, velocity and temperature. The NSR shell in Figure 5 can be identified most clearly in the top left panel as the narrowing light green segment running from bottom left (at ${\sim}0.25$ parsec) to the top right. In addition, the boundary of the ejecta pile-up, separating the cavity and the ejecta pile-up region can be seen as the other apparent line left of the remnant shell, running from the bottom left to the top centre of the panel (this boundary can be seen most clearly in the bottom right panel showing temperature evolution). This radial evolution of the shell and ejecta pile-up boundary directly replicates those seen in the right-hand plot of Figure 2, however here we show how each parameter changes over the full simulation.

The average density of the early NSR shell is approximately $n\,{\simeq}\,6$ for the first $10^{6}$ years of growth (see top left panel of Figure 5). Beyond this epoch, we see the effect of radiative cooling as the NSR shell loses energy and is compressed by the surrounding ISM and incoming eruptions, thereby leading to an increase in the average density within the shell to $n\,{\simeq}\,36$ after ${\sim}3\times 10^{7}$ years. The average density within the ejecta pile-up region is much lower than the surrounding ISM and continuously decreases throughout the evolution, dropping as low as $n=2.4\times 10^{-4}$ by the final epoch. After $10^{6}$ years the mass of the shell is ${\sim}50\ \text{M}_{\odot}$ but then substantially increases to $4\times 10^{3}\ \text{M}_{\odot}$ after $10^{7}$ years and ending with a mass of ${\sim}4\times 10^{4}\ \text{M}_{\odot}$ by the final epoch (${\sim}3\times 10^{7}$ years). This is consistent with the upper limiting shell masses derived from imaging and spectroscopy of the 12a NSR ($7\times 10^{5}\ \text{M}_{\odot}$ and $10^{6}\ \text{M}_{\odot}$ from assuming oblate and prolate geometries, respectively; DHO19, ).

As shown in the top right panel of Figure 5, the average pressure within the NSR shell is initially high as this thin high density region initially forms at high temperature. The pressure within the shell decreases until it matches the average pressure within the pile-up region after ${\sim}2\times 10^{7}$ years. The outer edge of the shell remains at the same pressure for the remainder of the simulation. However, the pressure at the inner edge increases, creating a pressure gradient within the shell. With the average temperature of the ejecta pile-up region increasing monotonically throughout its evolution (see the bottom right panel of Figure 5), the pressure within follows the same trend once that region’s size is established. The average pressure evolution illustrates how the NSR shell compression takes place during an intermediary period. The shell forms initially without compression, is then compressed as it is subjected to pressure gradients and after ${\sim}2\times 10^{7}$ years, the thinner shell remains.

The average temperature of the Run 1 NSR shell falls as a direct result of cooling due to expansion and radiative losses, dropping from an initial $5\times 10^{3}$ K to $40$ K after ${\sim}2.8\times 10^{7}$ years before increasing modestly to $90$ K as later eruptions become more frequent and begin to impact the inner edge of the shell through the pile-up region, leading to compression and re-heating (see bottom right panel of Figure 5). On the other hand, the pile-up region begins with higher temperatures of ${\sim}1\times 10^{6}$ K and continues to experience this temperature throughout before dramatically increasing to ${\sim}2.5\times 10^{8}$ K after the full $3\times 10^{7}$ years, maintaining these extremely high temperatures through shock-heating.

The average velocity of the NSR shell, like the average temperature and average pressure, decreases throughout the evolution before a slight increase for the final $6\times 10^{6}$ years (see bottom left panel of Figure 5). The velocity of the shell’s outer edge at ${\sim}6\times 10^{3}$ years is ${\sim}10\ \text{km s}^{-1}$ and remains below this velocity throughout. However, the velocity of the inner edge does increase due to the more frequent collisions occurring within the pile-up region, leading to a small velocity gradient within the shell. The ejecta pile-up region follows a similar trend but with higher average velocities, a result of increasingly frequent and higher velocity ejecta impacting the ejecta pile-up boundary. As the cavity is essentially a vacuum, the increasing velocities within this region are directly reflecting the increasing velocities of the nova ejecta.

3.3 Varying the ISM density

Here we consider the same nova system as Run 1 ($T_{\text{WD}}=10^{7}$ K; $\dot{M}=10^{-7}\,\text{M}_{\odot}\,\text{yr}^{-1}$), but placed in lower and higher density surroundings. Run 2 is pre-populated by ISM with a lower density of $1.67\times 10^{-25}\text{ g cm}^{-3}$ ($n=0.1$) and the ISM density of Run 5 and Run 7 is $1.67\times 10^{-23}\text{ g cm}^{-3}$ ($n=10$) and $1.67\times 10^{-22}\text{ g cm}^{-3}$ ($n=100$), respectively. We also sampled between these ISM densities with Run 3 ($n=10^{-0.5}\approx 0.316$), Run 4 ($n=10^{0.5}\approx 3.16$) and Run 6 ($n=10^{1.5}\approx 31.6$). As illustrated in Figure 6, the full simulations extend progressively further as the ISM density is decreased (e.g., ${\sim}$116 pc, ${\sim}$43 pc and ${\sim}$26 pc for Run 2 ($n=0.1$), Run 5 ($n=10$) and Run 7 ($n=100$), respectively) and all maintain an exceptionally thin shell due to the suppression of the early shell formation, reminiscent of Run 1. Furthermore, the remnants grown in Run 1, 2, 5 and 7 with radiative cooling are 78.26%, 63.29%, 78.31%, and 77.96%, respectively, of the size of their counterpart without cooling (from Runs 1⋆, 2⋆, 5⋆ and 7⋆) as a direct result of radiative losses from cooling. The relative thickness of the NSR shell varies for each simulation but remains small ($\lesssim 4\%$) for all ISM densities, resulting from the same amount of work done by the same nova system on surroundings that present increasingly higher resistance.

As expected, the density in the NSR cavity and pile-up region increases approximately in-line with ISM density. These regions are not only denser as a result of the ISM environment, but are also more compressed for higher $n$, leading to increased pressure. The velocity of material inside the NSR cavity from Runs 1–7 is identical as in all cases the ejecta are essentially undergoing free expansion. Also, temperatures in this region for each Runs 1–7 all reach the same extreme temperature of ${\sim}1\times 10^{9}$ K, as nova ejecta expanding without resistance collide into earlier ejected matter in the pile-up region, before dropping away to $<10$ K at the nova shell’s inner edge (i.e., the properties in this region don’t strongly depend upon $n$). The growth of the outer edge of the NSR shells within the $n=10$ and $n=100$ ISM follow a similar evolution as that of Run 1 (see the red line on the right plot of Figure 2).

We can summarise our findings for this section as follows: for a given total kinetic energy, an increase in local ISM density results in a smaller nova super-remnant.

3.4 Varying the mass accretion rate

The next six simulations (Runs 8–13) explored NSR evolution while varying accretion rate. We considered a WD with a temperature $10^{7}$ K accreting hydrogen rich material at a rate of $\dot{M}=10^{-8}\ \text{M}_{\odot}\ \text{yr}^{-1}$ as well as a nova with the same WD temperature but with a lower accretion rate of $\dot{M}=10^{-9}\ \text{M}_{\odot}\ \text{yr}^{-1}$, placed within the three ISM densities used in Runs 1, 5 and 7, see Table 1.

Runs 1–7 presumed that accretion was driven by the wind of a giant donor. We include mass loss from the donor between eruptions, although this has no impact upon the results (yet is computationally favourable, see Section 3). As such, we reduce the mass loss rate from the donor in line with any simulated changes to accretion rate for consistency and to ensure that the donor wind does not become important.

The WD growth models for $\dot{M}=10^{-8}\ \text{M}_{\odot}\ \text{yr}^{-1}$ and $\dot{M}=10^{-9}\ \text{M}_{\odot}\ \text{yr}^{-1}$ reveal that the WD loses mass with every eruption; it does not grow towards the Chandrasekhar limit, but is instead eroded. We therefore imposed a temporal upper limit of 100 Myr for the $\dot{M}=10^{-8}\ \text{M}_{\odot}\ \text{yr}^{-1}$ and $\dot{M}=10^{-9}\ \text{M}_{\odot}\ \text{yr}^{-1}$ simulations. The WD growth models indicate that these systems require 40,343 eruptions and 2,094 eruptions, respectively, to reach the temporal upper limit. At which point, these systems would have a recurrence period of ${\sim}$3,000 years and ${\sim}$49,000 years, respectively.

Focusing on Runs 8–10 ($\dot{M}=10^{-8}\ \text{M}_{\odot}\ \text{yr}^{-1}$) presented in the second row of Figure 7, we find that the overall structure of the remnants are similar to those grown with higher accretion rate. The major difference is their much larger size and thicker shells. The shell grown in the lowest density ISM (Run 8; $n=1$) extends to ${\sim}$99 pc, with a shell thickness of ${\sim}$11%, and Run 9 ($n=10$) and Run 10 ($n=100$) grow remnants with radial sizes of ${\sim}$62 pc and ${\sim}$40 pc, and shell thicknesses of ${\sim}$22% and ${\sim}$25%, respectively. These more extended shells are a consequence of the larger amount of kinetic energy ejected by the underlying system and the longer time over which it can act ($1\times 10^{8}$ years compared to ${\sim}3.1\times 10^{7}$ years in Run 1; see Figure 8). The outer edge of the NSR shell follows the same evolutionary trend as seen in Runs 1–7 (in the same manner as the remnant in the right plot of Figure 2).

In Runs 11–13 ($\dot{M}=10^{-9}\ \text{M}_{\odot}\ \text{yr}^{-1}$; $n=1,10,100$, respectively), we see that the NSRs take the familiar shape seen in Runs 1–10 with a very low density cavity preceding a high density shell (see the third row of Figure 7). The remnants grown in the Run 11 ($n=1$), Run 12 ($n=10$) and Run 13 ($n=100$) extend to ${\sim}$75 pc, ${\sim}$48 pc and ${\sim}$26 pc, respectively, and have shell thicknesses of 17%, 34% and 39%, respectively. Yet for each of these runs, the remnant shell is difficult to discern from the surroundings with the peak density within the NSR shell of Run 11, Run 12 and Run 13 reaching only 10.9%, 1.9% and 1.4% beyond that of the prepopulated ISM density, respectively. As expected, the outer shells of the remnants grown in systems with the lower accretion rate ($\dot{M}=10^{-9}\ \text{M}_{\odot}\ \text{yr}^{-1}$) follow the same growth curve over time as previous runs.

The nova eruptions from the systems in Run 11–13 occur infrequently for the vast majority of the evolution, starting with $P_{\mathrm{rec}}\sim$46,600 years when $M_{\mathrm{WD}}=1\,\text{M}_{\odot}$ and increasing to ${\sim}49,000$ years after the full $1\times 10^{8}$ years. Therefore, a combination of low energy eruptions and long recurrence period leads to a very broad, low-contrast shell as the ejecta individually dissipate into the surrounding ISM with minimal pile-up. Dynamically, such a NSR would be difficult to discern from the local environment. However, we would not expect this form of shell to exist around the known RNe as these systems would not (currently) be recognised as recurrent nova with their recurrence periods being $\gg 100$ years (see, for example, Darnley & Henze, 2020).

Equipped with the simulations of NSRs grown from systems with different accretion rates, we find that a lower accretion rate leads to more extended, but less well-defined, NSRs: a direct result of the longer evolutionary timescale.

3.5 Varying the white dwarf temperature

The underlying WD temperature does not have a significant impact on the evolution of most of the various parameters given in Section 2.2. For example, for $\dot{M}=10^{-7}\ \text{M}_{\odot}\ \text{yr}^{-1}$, the evolution of each parameter is very similar throughout, regardless of the WD temperature. Yet, there is a moderate difference in the evolution of the mass accumulation efficiency for the different temperatures. This is also true for the total kinetic energy of the ejecta generated from the entirety of the nova eruptions whereby the 30 MK and 50 MK have approximately twice the kinetic energy output as the cooler 10 MK WD. This is reflected in the set of simulations with the WD temperature varied from 10 MK (Run 1) to 30 MK (Run 14) to 50 MK (Run 15) with $\dot{M}=10^{-7}\ \text{M}_{\odot}\ \text{yr}^{-1}$ and $n=1$. A comparison of the NSR shell, as shown in the fourth row of Figure 7 for the three different WD temperatures, reveals the overall structure of each to be similar, but with the 30 MK WD remnant extending moderately further than the others. The outer edge of the remnant shell for the coolest WD is 71.3 pc and the hottest WD leads to an outer edge of 79.7 pc, whereas the outer edge of the 30 MK WD remnant shell is 97.4 pc. Yet, this informs us that, for the highest accretion rate we have considered, the WD temperature has a small impact on the large scale structure of the NSR in comparison to the effects of ISM density (Section 3.3) and mass accretion rate (Section 3.4).

There are similarities with the evolution of the shell for each WD temperature and at each epoch the density and thickness of the shells are a close match. By analysing how the recurrence period and the total kinetic energy change as the NSR grows in each of these systems, it is apparent that the WD temperature only has a relatively small impact. This may be due to the system having a high accretion rate ($10^{-7}\,\text{M}_{\odot}\,\text{yr}^{-1}$), so being dominated by accretion heating222Y05 accounted for accretion heating within their computations.. Any influence of WD temperature may become more substantial as accretion rate decreases as accretion heating will become less severe and the WD would have more time to cool between eruptions.

A further consideration is that unlike accretion rate and ISM density, which were both varied by factors of 10 and 100, the WD temperatures considered here only vary by factors of 3 and 5. The range we use (10 MK, 30 MK and 50 MK) was initially employed by Prialnik & Kovetz (1995)333Before consequently being adopted by Y05 with the incorporation of lower accretion rates for the cooler WDs. and was chosen to represent two extremes and an intermediate WD core temperature; the lower limit was set as a colder WD delays hydrogen ignition leading to long accretion times (hence more substantial eruptions) and the upper limit accounts for hot WDs being able to quickly reach the conditions for TNR.

We can conclude, for the accretion rate and ISM density ($n$) sampled in Runs 1, 14, 15, that the expected variation in WD temperature has much less impact on NSR evolution than plausible variations in accretion rate or $n$.

3.6 Varying the initial white dwarf mass

So far we have considered nova eruptions generated by a WD growing from $1\ \text{M}_{\odot}$ to $\text{M}_{\mathrm{Ch}}$. Here, we consider a number of different initial WD masses; $0.65\ \text{M}_{\odot}$ in Run 16, $0.8\ \text{M}_{\odot}$ in Run 17, $0.9\ \text{M}_{\odot}$ in Run 18 and $1.1\ \text{M}_{\odot}$ in Run 19 with $\dot{M}=10^{-7}\ \text{M}_{\odot}\ \text{yr}^{-1}$ and $n=1$. This upper initial mass is the upper formation limit for a CO WD (Ritossa et al., 1996). We also sample WDs with masses of $1.2\ \text{M}_{\odot}$ in Run 20 and $1.3\ \text{M}_{\odot}$ in Run 21. The number of eruptions appreciably increases as we lower the initial WD mass, as more eruptions are needed to reach $\text{M}_{\text{Ch}}$ (see Table 1).

A comparison of the NSR shells from these runs, presented in the last row of Figure 7, shows that each remnant becomes marginally larger as the initial WD mass is lowered, as more eruptions lead to more ejecta impacting the surrounding ISM over a longer period of time. The radial size of the NSRs in Runs 16, 17, 18 and 19 ($0.65\ \text{M}_{\odot}$, $0.8\ \text{M}_{\odot}$, $0.9\ \text{M}_{\odot}$ and $1.1\ \text{M}_{\odot}$) almost completely resemble that of the NSR from Run 1 ($1\ \text{M}_{\odot}$) whereas starting with a WD mass $>1.1\ \text{M}_{\odot}$ (in the regime of ONe WDs; Ritossa et al., 1996) such as simulated in Run 20 ($1.2\ \text{M}_{\odot}$) and Run 21 ($1.3\ \text{M}_{\odot}$) does make a difference to the radial size of the NSR. The structure of the shell for each NSR is remarkably similar, with the $0.65\ \text{M}_{\odot}$ WD simulation finishing with a shell thickness of ${\sim}$1.1% compared to ${\sim}$1.2% for the $1.3\ \text{M}_{\odot}$ WD. Each NSR shell also follows a very similar transition, with similar shell widths ratios at the same epochs. The radial growth curves of each simulation follow the same evolution with the $0.65\ \text{M}_{\odot}$ WD taking ten times ($37$ Myr) the time to reach $\text{M}_{\mathrm{Ch}}$ than the $1.3\ \text{M}_{\odot}$ WD ($3.7$ Myr).

We can therefore conclude, that the initial mass of the growing WD has little impact on the final structure of the NSR, much less than the prominent influence of the ISM density (Section 3.3) and accretion rate (Section 3.4).

4 Additional tests

In Section 3, we presented the full set of simulations. Here, we outline several tests of alternative models of the ejecta characteristics.

4.1 Using broken fits to estimate system parameters

For Runs 1–21, we utilised ejecta characteristics determined from our WD growth model. This was based on interpolating between the results of multi-cycle nova evolutionary simulations by Y05 (see Section 2.1). In our work, a smooth function asymptotically approaching $\text{M}_{\mathrm{Ch}}$ was fitted to the Y05 grid.

An alternative way of interpolating between the Y05 grid points is with a ‘knee’ function (e.g., Soraisam & Gilfanov, 2015, their Figure 1), which we replicated by fitting two distinct exponentials (Figure 1). From here, we grew a $1\ \text{M}_{\odot}$ WD with our model as outlined in Section 2.2, but referring in this case to the broken exponential fits.

Eruption parameters evolve in the same way as those from the smooth function fitting, with the main difference being the abrupt ‘knee’ at $1.25\ \text{M}_{\odot}$. The total kinetic energy at the end of the WD growth is ${\sim}5.3\times 10^{-2}$ foe ($10^{51}$ ergs). This is much greater than the total kinetic energy generated from our smooth fitting function in Section 2.2; this ended with ${\sim}2.5\times 10^{-2}$ foe. This reflects the more extreme eruptions later on in this system’s evolution as a direct result of the higher ejecta velocities after the WD has surpassed 1.25 $\text{M}_{\odot}$.

We ran a simulation (Run 1*†) of nova eruptions generated from the two distinct exponential fits, with the same parameters as our reference simulation (Run 1) including $\dot{M}=10^{-7}\ \text{M}_{\odot}\ \text{yr}^{-1}$ and a WD temperature of $10^{7}$ K and $n=1$ ISM (see Table 1). As can be seen in Figure 9, the shell grown from the broken exponential fitting does not grow as large as the shell grown from the smooth fitting. This is a consequence of the much higher mass accumulation efficiency between $1\ \text{M}_{\odot}$ and $1.25\ \text{M}_{\odot}$ (see Figure 1) for the broken exponential fit, resulting in lower levels of ejecta and substantially less kinetic energy during the early stages of NSR growth; this period has a major impact on the proceeding evolution. Beyond $1.25\ \text{M}_{\odot}$, the radial growth curve of the Run 1†* shell deviates from that of Run 1 (at approximately $1.4\times 10^{7}$ years; see the inset of the right panel in Figure 9) as a result of the later eruptions becoming more extreme.

As shown in the inset of the left panel in Figure 9, both the shells in Run 1 and Run 1*†* have a similar structure however the shell in Run 1*†* is thinner, and consequently, comprises a higher density inner edge. This also has a greater impact on the temperature gradient of the shell in Run 1*†*, with the outer edge being much hotter than the inner edge, unlike Run 1.

It is clear that using an alternative interpolation to the values given in Y05 does have an effect on the final simulated NSR. In the case of the Run 1*†*, the shell width is approximately half the shell width of the remnant in Run 1 plus the size of the remnant decreases by a factor of ${\sim}$12%. Whilst a non-negligible difference, we consider the more realistic smooth evolution of system parameters adopted for our study to be a truer representation for NSR simulations. Nevertheless, this does indicate the need for more finely sampled nova model grids.

4.2 Eruption characteristics

Although we are predominantly concerned with the long term evolution (and therefore large scale structure) of a NSR, we explored several eruption characteristics to observe how NSR evolution is affected. Firstly, we know that the timescale of the nova eruption can vary as we see a wide range of SSS periods (see, e.g. Henze et al., 2014). Secondly, shocks play a key role within the nova ejecta, and instead of material being ejected in one event, the eruption contains a number of components with varying masses and velocities (Metzger et al., 2014; Aydi et al., 2020a, b; Murphy-Glaysher et al., 2022).

4.2.1 Eruption duration

As an extension to the Run 0 tests (DHO19), to determine if the duration of a nova eruption affects NSR large scale structure we ran high resolution (${\sim}$4 AU/cell) simulations, each with 1000 eruptions utilising the Run 0 setup with a range of eruption durations: 0.07 d, 0.7 d, 7 d, 70 d and 350 d. For each test, the eruption duration plus the quiescent period match the recurrence period (350 days; e.g., 349.03d + 0.07d or 343d + 7d), with a fixed ejecta velocity of 3000 km s*-1*. We required each test to inject the same total kinetic energy, so the eruption mass-loss rate was decreased to account for the longer timescales. After around 100 eruptions, the inner and outer edges of the NSR shell followed the same evolutionary trend regardless of eruption duration, and even though the NSR pile-up fluctuates more than the shell, they again settle into similar growth rates. This removes eruption duration dependency and indicates that our NSR results are not sensitive to any assumptions made about eruption time-scales.

4.2.2 Intra-eruption shocks

We also wanted to test whether having a non-uniform ejection of material from the nova would affect the large scale structure of the shell. For this, we considered the composition of a classical nova whereby the eruption takes place over a certain timescale and over that time period, the speed of ejection increases (Bode & Evans, 1989; O’Brien et al., 1994; Metzger et al., 2014; Aydi et al., 2020a, b). This implies that the outburst is comprised of a slow wind followed by a faster wind, creating a shock within the ejecta (O’Brien et al., 1994; Metzger et al., 2014; Aydi et al., 2020a, b).

We ran a Run 0-based simulation following 1000 eruptions with a 7 day duration. To incorporate intra-ejecta shocks, we split the ejecta into two separate components. For moderate-speed novae, the ejecta velocities range from 500–2000 km s*-1* but for fast novae, this range is 1000–4000 km s*-1* (O’Brien et al., 1994). As we are considering recurrent nova eruptions and therefore dealing with fast novae, we used the latter range of velocities for this test. We ejected half the mass at 1000 km s*-1* over 3.5 days immediately followed by half of the mass at 4123 km s*-1* over the next 3.5 days, such that the total kinetic energy matched that of a 7 day eruption with an ejecta velocity of 3000 km s*-1*. As the second half of the mass is ejected at a higher velocity than the first, we see intra-ejecta shocks as the later ejecta overtakes and interacts with the earlier ejecta. Again, after around 100 years, the inner and outer edge of the NSR shells created from ejecta with and without intra-eruption shocks follow the same evolutionary trend.

In Sections 4.2.1 and 4.2.2, we have demonstrated that the long term evolution of nova ejecta is not affected by the nova eruption duration nor by the presence of intra-ejecta shocks, and consequently, neither is any NSR. NSR evolution only depends upon the total kinetic energy of the ejecta and the surrounding medium444Here, we are considering a pure adiabatic scenario..

5 Observational predictions

Here, we investigate the evolution of NSR observables, derived from Run 1 (Section 3.2), in part to inform any NSR follow-up observations or searches. The simplest and computationally cheapest way to predict the emission over a full simulation of a NSR is by assuming a pure hydrogen environment. We can thus compute the ionisation fraction ($f$), emission measure (EM), recombination time-scale, X-ray luminosity, and H$\alpha$ emission. In general, an assumption of pure hydrogen provides a good estimate of $f$ throughout the NSR.

5.1 Evolution of emission measure

Assuming pure hydrogen, we employed the Saha (1921) equation to compute $f$ for each NSR cell across all epochs. As the number of free protons in a medium of fully ionised hydrogen is equal to the number of electrons, we define the EM in each NSR cell as the square of electron density ($n_{\mathrm{e}}^{2}$) integrated over the volume of the spherical shell represented by each cell.

The EM from the different Run 1 NSR regions (cavity, ejecta pile-up, shell, and the entire NSR) at each epoch were calculated by integrating over all shells within each region. The mean ionisation fraction ($\bar{f}$) in each region, per epoch, was computed in a similar fashion while also weighting each shell by density.

The evolution of $\bar{f}$ and the total EM for each region is shown in Figure 10. In Figure 11, we show the evolution of $\bar{f}$ and EM for the cavity, ejecta pile-up region and shell alongside the evolution of the mean density and temperature.

As illustrated in Figure 11, the mean temperature of the pile-up region is approximately a few $\times 10^{6}$ K for ${\sim}2.7\times 10^{7}$ years (except during the initial eruption) and begins to increase toward ${\sim}2.8\times 10^{8}$ K during the next ${\sim}3.5\times 10^{6}$ years of the NSR evolution. The density in this region decreases by over a factor of 2 as it grows but the extremely high temperatures maintains $\bar{f}\gtrsim 25\%$. As a result, the EM from this region remains high. Within the cavity, $\bar{f}\gtrsim 1\%$, and so even with density decreasing over time, the emission from this region remains a contributing, albeit fluctuating, factor, until latter stages of NSR evolution.

If we focus on $\bar{f}$ within the NSR shell in Figure 11, we see the effect of recombination as a result of the high densities and cooling. For the first $10^{5}$ years, the NSR shell is fully ionised, here the shell EM is high and the dominant source. After this, $\bar{f}$ in the shell decreases to negligible levels as the material recombines and remains neutral for the majority of the NSR lifetime (from ${\sim}10^{5}$ years to ${\sim}3\times 10^{7}$ years) which, combined with an almost constant mean density during this period, leads to a drop in EM to effectively zero. However, as with the other regions, the late-time frequent highly energetic eruptions begin to re-heat the NSR shell, increasing $\bar{f}$ marginally. The high NSR shell density at this time leads to the NSR shell again contributing to the EM at the end of the simulation.

The evolution of the total NSR EM is shown in the bottom-right panel of Figure 10. The NSR shell initially dominates the EM as this high density region begins to sweep up ISM. After ${\sim}5\times 10^{5}$ years, the average temperature within the shell has decreased enough for the material to recombine, resulting in a dramatic reduction in EM from this region. As a result, the total EM from the NSR becomes dominated by the pile-up region between ${\sim}5\times 10^{5}$ years and ${\sim}3\times 10^{7}$ years, with additional contribution from the fluctuating cavity emission throughout (originating from the eruptions themselves). Once the later stages have been reached (the last ${\sim}5\times 10^{5}$ years), with frequent highly energetic ejecta, the rate of ionisation within the very high density shell (particularly at the inner edge) leads to a substantial increase in EM from this region. However, unlike at early times when EM was dominated by the entire NSR shell, the emission at these later times emanates exclusively from the pile-up region and the inner edge of the shell.

5.2 Evolution of recombination time

The Morpheus code only informs as to the ionisation state of the material based upon the dynamics of the simulation; it does not include radiative transfer. As such, when considering the emission from simulated NSRs, and indeed their observability, we must also take account of recombination timescales ($t_{\mathrm{recomb}}$).

As recombination time depends upon the relative abundances of the gas, from this point on we assume that all material is of Solar composition. While this will be a good approximation for the ISM it will be less so for the ejecta. However, the NSR is predominantly swept up ISM. Abundances from Wilms et al. (2000) were utilised to determine $f$ for H, He, C, N, O, Ne, Na, Mg, Al, Si, P, S, Cl, Ar, Ca, Ti, Cr, Mn, Fe, Co, and Ni within the NSR. We compute the minimum recombination time for all cells of Run 1 by assuming the NSR is fully ionised, thus providing a lower limit on recombination time for each cell.

The recombination time evolution across the entire Run 1 NSR remnant is shown in Figure 12. Here, we see that maximum recombination time within the NSR shell (except for the first epoch considered) is always $\lesssim 3\times 10^{4}$ years and with the peak always being at the inner edge of the shell, the minimum recombination time of the shell approximately corresponds with the peak density. As the evolving WD approaches $\text{M}_{\mathrm{Ch}}$, the amount of ionised mass within the NSR ejecta pile-up region (effectively the entire NSR), reaches ${\sim}10\,\mathrm{M}_{\odot}$, as gas within the pile-up region is heated by the late-time frequent and energetic eruptions. This is once again reflected in the moderate rise of the recombination time at the inner edge of the shell in Figure 12 (the thick black dotted line tracing the NSR shell inner edge). Notably, the mass weighted median recombination time (indicated by the dashed line) remains essentially constant throughout after ${\sim}5\times 10^{6}$ years, hence we adopt $t_{\mathrm{recomb}}=315$ yr throughout the Run 1 NSR shell (the mass weighted median recombination time during the epoch when $P_{\mathrm{rec}}=1$ yr).

5.3 Evolution of X-ray luminosity

Following Vaytet (2009), Vaytet et al. (2011) and DHO19, we compute the EM contribution from each Run 1 NSR spherical shell (as defined by the simulation cells) and then bin the EM contribution into 95 logarithmically divided temperature bins ranging from 149 K to ${\sim}3.9\times 10^{9}$ K (based on the shell/cell temperature). The temperature-binned EMs are used as inputs to XSPEC. Within XSPEC, we utilise the APEC (Smith et al., 2001) model which computes an emission spectrum containing lines for H, He, C, N, O, Ne, Mg, Al, Si, S, Ar, Ca, Fe, and Ni with Solar abundances (He fixed at cosmic) from a collisionally-ionised diffuse gas.

The EM histograms can also be used to broadly explore the evolution of NSR emission as a function of photon energy and hence wavelength. Tracking the emission evolution for the Run 0 NSR (see Extended Data Figure 7 in DHO19, ) reveals that it starts off at high temperatures, emitting mostly in X-rays at ${\sim}$1 keV as in Run 0, as the eruptions are immediately frequent and highly energetic. But, as the NSR shell grows and cools, the EM peak moves toward lower energies, ending in the optical/NIR region (${\sim}2\times 10^{-3}$ keV) after the full $10^{5}$ eruptions. A logarithmic extrapolation of the EM indicates that the present day peak might be in the infrared, around 12–13 $\mu$m, and could be a potential target for JWST (DHO19).

On the other hand, the Run 1 NSR begins with the peak EM at low energies (optical/NIR) due to the long period between the initial low energy eruptions allowing the NSR to cool. The temperature of the NSR as a whole remains low throughout the evolution and the EM peak remains at low energies through all 1,900,750 eruptions.

Separating the EM evolution into the component NSR parts, namely the cavity, pile-up, and shell, provides the contributions from each of these regions. The cavity emission remains relatively low compared to other regions throughout the full evolution. For the first ${\sim}10^{4}$ eruptions, the cavity emits in the optical/NIR regime. However, when the recurrence period approaches one year, the contribution from the cavity, albeit small, branches across to higher energies. This may be attributed to the ejected material colliding with the inner edge of the pile-up region.

Emission levels from the pile-up region are considerably higher than from the cavity and contribute more to the X-ray emission at later times as incoming ejecta continuously shock-heat this region. In fact, after the full 1,900,750 eruptions, a portion of the pile-up region emits in excess of 100 keV. In contrast, the NSR shell emits mostly in the optical at early times before peaking after only $10^{3}$ eruptions, when the majority of the emission lies in the NIR. Beyond this epoch, for the entire evolution of the NSR, the shell contributes a negligible amount of emission and it remains the coolest part of the NSR, largely shielded from the highly energetic material.

We use the EM to predict the evolution of the Run 1 X-ray luminosity. We assumed that our simulated NSR is at a distance of 778 kpc (Stanek & Garnavich, 1998, i.e., within M 31). To remove the impact of single eruptions, we re-bin to a lower temporal resolution. This is illustrated in Figure 13 with comparison to the X-ray luminosity evolution from the NSR created in Run 0.

As shown in the left plot of Figure 13, for the Run 0 NSR, the X-ray luminosity peaks at ${\sim}6\times 10^{31}\,\text{erg}\,\text{s}^{-1}$ after approximately $10^{3}$ years (equivalent to $10^{3}$ eruptions for Run 0). This luminosity then fades to ${\sim}9\times 10^{29}\,\text{erg}\,\text{s}^{-1}$ after $10^{5}$ years/eruptions and with a power-law extrapolation to the latest time, representing present day in DHO19, the total X-ray luminosity drops to ${\sim}3\times 10^{29}\,\text{erg}\,\text{s}^{-1}$. As detailed in DHO19, the X-ray luminosities predicted for the entire NSR evolution lie well below the $3\sigma$ upper limiting luminosity of ${\sim}9\times 10^{34}\,\text{erg}\,\text{s}^{-1}$ constrained by archival X-ray observations (see horizontal dotted line in the left plot of Figure 13).

The Run 1 NSR X-ray luminosity follows an entirely different evolution from Run 0 (see right plot of Figure 13). While X-ray emission is predicted from the onset of Run 0, we predict negligible X-ray emission from the Run 1 NSR until $3\times 10^{5}$ years; $P_{\mathrm{rec}}\lesssim 85$ yr. From that point, the X-ray luminosity rises as the recurrence period falls and the ejecta become more energetic, increasing significantly during the final ${\sim}10^{7}$ years. Starting at ${\sim}2\times 10^{29}\,\text{erg}\,\text{s}^{-1}$ after ${\sim}3\times 10^{5}$ years, the initial X-ray luminosity is dominated by soft emission between 0.3–1 keV.

The influence of the more frequent and energetic eruptions becomes evident over the next 26 Myr as harder emission from shock-heating, with energies between 1–10 keV, reaches ${\sim}1.5\times 10^{30}\,\text{erg}\,\text{s}^{-1}$ after ${\sim}2.7\times 10^{7}$ years (see inset of the right plot in Figure 13), contributing greatly to the total X-ray luminosity of ${\sim}1\times 10^{31}\,\text{erg}\,\text{s}^{-1}$ at this epoch. However, this is still much fainter than typical nova X-ray luminosities such as, for example, M31N 2004-01b, 2005-02a, and 2006-06b with $L_{\text{X}}=(11.1\pm 1.6)\times 10^{36}\,\text{erg}\,\text{s}^{-1}$, $2.6\times 10^{37}\,\text{erg}\,\text{s}^{-1}$ and $(3.6\pm 0.3)\times 10^{36}\,\text{erg}\,\text{s}^{-1}$, respectively555Unabsorbed luminosity between 0.2–10 keV. (see Henze et al., 2010, 2011, for a large sample of M 31 CNe X-ray luminosities). Instead, this X-ray luminosity is more akin to that seen in quiescent novae such as ${\sim}6\times 10^{31}\ \text{erg}\ \text{s}^{-1}$ for RS Ophiuchi (Page et al., 2022).

The NSR X-ray luminosity then continues to increase for the remainder of the evolution, ending with a luminosity of ${\sim}1\times 10^{31}\,\text{erg}\,\text{s}^{-1}$. This is due to hard emission (1–10 keV) becoming increasingly significant, with harder emission between 10–50 keV appearing in the final $4\times 10^{6}$ years. If we consider the $P_{\mathrm{rec}}=1$ yr epoch, the Run 1 NSR X-ray luminosity is ${\sim}9\times 10^{30}\,\text{erg}\,\text{s}^{-1}$ (see inset of the right plot in Figure 13). This is $30\times$ greater than the present day extrapolated luminosity from Run 0.

5.4 Evolution of H$\alpha$ flux

From observations of the 12a NSR, we know such structures should be visible through their H$\alpha$ emission (DHO19). As such, we utilised Run 1 to predict the evolution of H$\alpha$ emission from a NSR in a similar manner to that described in Andersson (2021). The H$\alpha$ luminosity was calculated by convolving the EM histograms with the appropriate temperature-dependent recombination coefficient for the given temperature (from Pequignot et al., 1991). The NSR was placed at the distance of M 31 and we applied extinction of $A_{\text{H}_{\alpha}}=0.253$ to find the H$\alpha$ flux across the simulated NSR. The evolution of the Run 1 NSR H$\alpha$ flux is presented in Figure 14.

The Run 1 NSR H$\alpha$ evolution broadly follows the EM evolution (cf. Figures 10 and 11). Initially, as the early NSR shell sweeps into the ISM, the H$\alpha$ emission (predominantly emanating from the shell) follows a roughly power-law increase, reaching a peak of ${\sim}8\times 10^{-18}\,\text{erg s}^{-1}\,\text{cm}^{-2}$ after $10^{5}$ yr. Beyond this time however, the shell temperature decreases, allowing for recombination and a consequent (power-law-like) drop in H$\alpha$ emission. As described in Section 5.1, between ${\sim}10^{5}$ yrs and ${\sim}3\times 10^{7}$ yrs, the main sources of H$\alpha$ emission are the pile-up region and cavity. The cavity contribution can be seen as the numerous spikes in H$\alpha$ flux, with the later energetic eruptions from the nova colliding with the sparse material within that region. As shown in Figure 14, the last ${\sim}8\times 10^{6}$ years then see a dramatic increase in H$\alpha$ emission, almost exclusively coming from the highly energetic eruptions at this stage impacting the high density inner edge of the formed NSR shell and pile-up region, reaching a maximum of ${\sim}6\times 10^{-18}\,\text{erg s}^{-1}\,\text{cm}^{-2}$ after the full $3.1\times 10^{7}$ years.

As shown in Figure 14, we also modelled the NSR H$\alpha$ flux evolution for Run 8 ($\dot{M}=10^{-8}\,\text{M}_{\odot}\,\text{yr}^{-1}$) and Run 11 ($\dot{M}=10^{-9}\,\text{M}_{\odot}\,\text{yr}^{-1}$) to explore the impact of mass accretion rate on H$\alpha$ observability. Early in their evolution, H$\alpha$ emission from these NSRs follow a similar, but much fainter, evolution to the NSR emission in Run 1 (with $\dot{M}=10^{-7}\,\text{M}_{\odot}\,\text{yr}^{-1}$). However, unlike in Run 1 where the H$\alpha$ flux begins to increase beyond ${\sim}10^{6}$ years, the emission in Run 8 and Run 11 drops away, and continues to do so for the rest of the NSR’s growth. In both Run 8 and Run 11 the H$\alpha$ flux drops to ${\sim}2\times 10^{-20}\,\text{erg s}^{-1}\,\text{cm}^{-2}$ after $1\times 10^{7}$ years (compared to ${\sim}1\times 10^{-19}\,\text{erg s}^{-1}\,\text{cm}^{-2}$ in Run 1) and ends with ${\sim}6\times 10^{-21}\,\text{erg s}^{-1}\,\text{cm}^{-2}$ and ${\sim}4\times 10^{-21}\,\text{erg s}^{-1}\,\text{cm}^{-2}$, respectively, after $1\times 10^{8}$ years.

We can tentatively conclude from these models, that NSR H$\alpha$ emission for systems with high accretion rates is significant early on in NSR growth (younger RNe systems) and again late on in the NSR evolution, from older RNe systems such as the RRNe. Furthermore, the brightest NSRs are the systems containing near-Chandrasekhar mass WDs. However, for systems with lower accretion rates, in which the WD is eroding, the H$\alpha$ emission at latter stages of evolution is orders of magnitude fainter than observed in high accretion systems.

6 Comparing simulations and observations

6.1 Run 1 versus the 12a nova super-remnant: dynamics

To determine how well these simulations recreate properties of the only known NSR, we compare them to observations of the 12a NSR. For this, we will consider the simulated NSR grown from a nova with parameters that most resemble 12a. The 12a mass accretion rate derived from observations is $(6-14)\times 10^{-7}\,\mathrm{M}_{\odot}\,\mathrm{yr}^{-1}$, the closest accretion rate we were able to consider is $10^{-7}\,\mathrm{M}_{\odot}\,\mathrm{yr}^{-1}$, within Runs 1–7. The 12a $P_{\mathrm{rec}}=1$ yr, therefore we compare with simulations at this recurrence period (${\sim}$99.54% through the simulations). At this point, the simulated WD mass is ${\sim}1.396\ \mathrm{M}_{\odot}$.

The most immediate difference we see between observations and the simulations is the NSR radial size and the shell thickness. Within the reference simulation (Run 1; $n=1$), the NSR extends to ${\sim}71.3$ pc compared to the observed 67 pc (DHO19). Furthermore, DHO19 assumed that 12a is located within a high density environment, which leads to a smaller NSR, more closely resembling the Run 7 NSR ($n=100$). The shell thickness of the Run 1 NSR is ${\sim}1\%$, dramatically smaller than the 22% derived from observations of the inner and outer edges of the 12a NSR (DHO19).

As with the first simulation (Run 0) of a NSR, the general shell structure of the NSR in Runs 1 – 7 is reminiscent of the observed shell. They all have a very low density central cavity (not apparent in observations) with freely expanding high velocity ejecta leading up to a very hot pile-up region. Spectroscopy of an inner ‘knot’ in the 12a NSR reveals strong [O iii] emission, indicative of higher temperatures closer to the 12a system (DHO19). In the 12a observations, we see evidence for a high density shell sweeping up the surrounding ISM, which is replicated in Runs 1 – 7. The lack of [O iii] emission in the 12a shell demonstrates that the shell has cooled below the ionisation temperature of O*+* (DHO19).

We can conclude that the simulations that most resemble the 12a NSR, in terms of accretion rate and ISM density (Runs 1 – 7), can replicate the radial size of the NSR that is observed, but not its shell thickness. As a result, we can only conclude that there must be other contributing factors in the evolution and shaping (or geometry) of these structures that we have not yet considered. In particular, we wish to explore the impact of early helium flashes as well as a non-fixed accretion rate on NSR evolution in future work. Additionally, the simulations presented in this work are one-dimensional and so are not susceptible to Rayleigh-Taylor or Richtmyer–Meshkov instabilities. This additional physics is likely to influence the dynamics of the growing shell, through for example, shell fragmentation as seen in Toraskar et al. (2013).

But importantly, our models only simulate the dynamically grown structure, and associated emission of a NSR, they do not (yet) consider additional effects that photoionisation may have on any observed NSR (see Section 6.3).

6.2 Run 1 versus the 12a nova super-remnant: emission

We again explore the epoch of Run 1 that coincides with $P_{\mathrm{rec}}=1$ yr (after $3.04\times 10^{7}$ yr) to predict the X-ray luminosity and H$\alpha$ flux, as in Section 5, to directly compare to the emission from 12a’s NSR.

6.2.1 Emission measure at one year recurrence period

We follow the procedures in Section 5.1 to compute the ionisation fraction ($f$) and emission measure (EM) for the NSR at $3.04\times 10^{7}$ yr (see Figure 15). Here, the entire NSR, up to the inner edge of the shell is fully ionised ($f=1$). The ionisation decreases dramatically, to negligible values, within the shell. This fully ionised state within the cavity (up to ${\sim}10$ pc) can be attributed to the ejecta interaction with the RGW and subsequent free expansion. Within the pile up region (between ${\sim}10-70.2$ pc), gas is continuously impacted by incoming eruptions and shocks resulting in collisional excitation and, consequently, $f=1$. Shocks are also present at the inner edge of the NSR shell (${\sim}70.2$ pc) as gas flows through the pile-up region into the swept up shell. However, further into the shell, toward the outer edge (${\sim}71$ pc), the gas is dynamically shielded from incoming shocks and does not experience a high level of ionisation.

6.2.2 Recombination time at one year recurrence period

In Section 5.2, we computed minimum recombination times throughout the NSR evolution by considering the recombination time for a hypothetical fully ionised NSR. For the epoch of this simulation where $P_{\text{rec}}=1$ yr, we also compute recombination time for the NSR given the $f$ predicted by the dynamic growth.

The recombination times for a NSR dominated by Solar material are illustrated in Figure 16 with the red line. Recombination times throughout the NSR are extremely long, owing to the extremely low density and continuous ejecta–RGW shocks within the cavity (up to ${\sim}10$ pc). Within the pile-up region (${\sim}10-70.2$ pc) the continual shock-heating from colliding ejecta drives the recombination time high. At the inner edge of the NSR shell, where the gas density dramatically increases, we see the recombination time drop to a $2\times 10^{5}$ yrs. Beyond the inner edge (at the front end of the shell), cooler neutral gas forces the recombination time to increase substantially. When considering an already fully-ionised NSR, we still see extremely long recombination times within the cavity and pile-up regions. However, we do see a significant difference within the NSR shell. As before, the recombination time drops dramatically at the inner edge yet now we see $t_{\mathrm{recomb}}{\sim}10$ yr at the inner edge, rising to ${\sim}10^{4}$ yrs at the outer edge.

As a result of the high recombination times within cavity and pile-up regions of the NSR and recombination times in the shell on a par with the travel time for nova ejecta to cross the NSR (${\sim}3.4\times 10^{4}$ yrs for ejecta travelling at ${\sim}2000\,\text{km}\,\text{s}^{-1}$), the NSR shell may exhibit emission induced by photoionisation from the nova eruptions. Furthermore, if the ISM density is low enough (see Section 6.3), then ionising radiation from the central source might traverse the (fully collisionally ionised) inner regions of the NSR with the ability to potentially create an ionised region beyond (or within the shell of) the dynamically grown NSR.

6.2.3 X-ray luminosity at one year recurrence period

The output from the Run 1 NSR at the epoch coinciding with $P_{\mathrm{rec}}=1$ yr was processed and passed to XSPEC. The X-ray luminosity as a function of radius was calculated using the APEC model (without the incorporation of absorption) and is shown in Figure 17.

At the centre of the remnant, there is a high X-ray luminosity from the underlying system due to the nova eruptions, however this is then followed by negligible emission from the cavity as the ejecta are in free expansion. Beyond this cavity, the ejecta begins to impact the higher density pile-up region (up to ${\sim}10$ pc), leading to a significant jump in the X-ray luminosity (${\sim}1\times 10^{22}\,\text{erg}\,\text{s}^{-1}$). As more and more ejecta contribute toward shock-heating the pile-up region further from the centre, we see a continuous increase in X-ray emission up to the inner edge of the NSR shell at ${\sim}70.2$ pc, where $\text{L}_{\text{X-ray}}\simeq 4\ \times 10^{27}\,\text{erg}\,\text{s}^{-1}$. The total predicted X-ray luminosity from the NSR at this epoch is ${\sim}1\times 10^{31}\,\text{erg}\,\text{s}^{-1}$ (see Figure 13). This is consistent with the unabsorbed luminosity upper limit of the NSR associated with 12a derived from archival XMM-Newton observations (${<}9\ \times 10^{34}\ \text{erg}\ \text{s}^{-1}$; DHO19, ).

6.2.4 H$\alpha$ flux at one year recurrence period

We applied the technique set out in Section 5.4 to Run 1 at the epoch corresponding to $P_{\mathrm{rec}}=1$ yr to compare the predicted H$\alpha$ emission to that from the 12a NSR (see Figure 18). Here, we see that there is H$\alpha$ emission from the cavity and increasingly from the ejecta pile-up region, yet this always remains below ${\sim}10^{-27}\,\text{erg s}^{-1}\,\text{cm}^{-2}$. However, as is the case for X-ray emission, the majority of H$\alpha$ flux originates at the inner edge of the NSR shell. Here, the density of hydrogen is extremely high compared to the rest of the NSR and so the large amount of collisional excitation from the impacting ejecta results in high levels of recombination and H$\alpha$ emission of ${\sim}3\times 10^{-18}\,\text{erg s}^{-1}\,\text{cm}^{-2}$, many orders of magnitudes higher than anywhere else across the NSR. The total predicted H$\alpha$ luminosity from the NSR at this epoch is $\text{L}_{\text{H}\alpha}\simeq 3.6\ \times 10^{32}\,\text{erg}\,\text{s}^{-1}$.

6.3 A photoionisation remnant?

As stated in the previous section, the dynamic simulations in this work, with parameters most similar to 12a, replicate the broad observed structure, but not the shell thickness, or potentially observability, of the 12a NSR. But so far, we have only considered the growth and emission of the dynamically formed NSR. However, a proportion of the NSR will be exposed to photoionisation directly from the central system, the accretion disk, the eruptions, as well as any shock emission. As such, we consider here the formation and radial size of the photoionisation remnant, and any dependence upon ISM density. We will assume that material inwards from the NSR shell is fully ionised throughout the evolution as discussed in Section 5.1 and shown in Figure 11.

We show in Figure 19, the dynamical remnant inner (purple) and outer (green) radii for Runs 1–7 with respect to ISM density at the epoch when each of the runs have recurrence periods of one year and assuming that the mass accretion rate is $10^{-7}\ \text{M}_{\odot}\ \text{yr}^{-1}$. We then interpolated these points with a power-law fit.

To estimate the size of any photoionisation region generated by the nova eruptions, we can perform a Strömgren-like analysis as $t_{\mathrm{recomb}}\gg P_{\mathrm{rec}}$ within the NSR shell and the ISM. However, because all material inwards of the NSR shell is always fully shock ionised, instead of a Strömgren sphere, we will have a Strömgren shell. Consequently, the photoionisation region can be estimated thus:

[TABLE]

where $r_{\mathrm{out}}$ is the outer radius of the photoionised region, $S_{\star}$ is the ionising luminosity from the source, $n(r)$ is the number density of the medium, $\beta(r,T)$ is the total recombination rate for Case B recombination (see, e.g., Dyson & Williams, 1980), and $r_{\text{in}}$ is assumed to be the outer edge of the fully ionised region of the NSR. This $r_{\text{in}}$ was determined to be the first point from the center of the NSR in which the ionisation fraction ($\bar{f}$) falls below 100%.

We will take the ionising luminosity from the nova eruptions (or the SSS emission) as the Eddington luminosity of a $1.396\,\text{M}_{\odot}$ WD (the mass of the WD in our models at the time when $P_{\text{rec}}=1$ yr) minus the observed luminosity of the 12a SSS, such that $\text{L}_{\text{Edd}}-\text{L}_{\text{obs}}\approx 41,400\ \text{L}_{\odot}$ for two weeks (the SSS timescale of each eruption; Henze et al., 2018) and assume a spectrum of 15 eV photons, giving a time averaged $S_{\star,\text{SSS}}=6.6\times 10^{48}\,\mathrm{photons}\,\mathrm{s}^{-1}$. Substituting this into the equation 2 along with varying values for $n$ (ISM density) provides us with the width of the ionisation region for Runs 1–7 (see the orange points in Figure 19). We also calculated a similar ionisation region but with the inclusion of the disk luminosity ($5910\ \mathrm{L}_{\odot}$) such that $S_{\star,\text{disk}}=9.4\times 10^{47}\ \mathrm{photons}\ \mathrm{s}^{-1}$ (with this emission present at all times), alongside the SSS emission (see the yellow points in Figure 19). Again, we assumed a $1.396\,\text{M}_{\odot}$ WD, $\dot{M}=10^{-7}\,\text{M}_{\odot}\,\text{yr}^{-1}$, and a spectrum of 15 eV photons to estimate the disk luminosity using $\mathrm{L}_{\mathrm{disk}}=(G\dot{m}\mathrm{M}_{\mathrm{WD}})/\mathrm{R}_{\mathrm{WD}}$. We also considered the contribution of ionising photons from shocks, by computing the shock emission within XSPEC when $P_{\text{rec}}=1$ yr, yielding $S_{\star,\text{shock}}=2.9\times 10^{41}\ \mathrm{photons}\ \mathrm{s}^{-1}$. But this is many orders of magnitude less than $S_{\star,\text{SSS}}$ and $S_{\star,\text{disk}}$ and so was not considered further.

With the two luminosities we do consider, the widths of the ionisation regions (SSS or SSS+disk) produced can be found and are shown in Figure 19 with orange (SSS) and yellow (SSS+disk) points. For all of the NSRs grown in Runs 1–7, the emission from the nova system (eruptions and disk) cannot ionise the NSR shell and so the ionisation regions are fully contained within the remnant shell. This suggests that observations of NSRs should exhibit emission at the inner edge of the NSR shell.

To test this, we created a synthetic sky image (with the inclusion of seeing) to directly compare with observations of the NSR surrounding M 31N 2008-12a. Utilising the outer edge fitting presented in Figure 19, we determined the ISM density required to grow a NSR with the same radial extent as the observed NSR around M 31N 2008-12a (67 pc) to be $n=1.278$. With this, we ran another simulation (Run 22) with the same system parameters as Run 1 but with an ISM density of $n=1.278$. The outer edge of the NSR grown within Run 22 extends to ${\sim}67.4$ pc (as expected) and exhibits a shell thickness of $1.1\%$. The boundary between the cavity and ejecta pile-up region is located at ${\sim}9$ pc and the inner edge of the NSR shell is at ${\sim}66.6$ pc, with a density of $2.6\times 10^{-22}$ ($n\simeq 122$).

Using the same technique as described in Sections 5.4 and 6.2.4, we took the Run 22 NSR at the epoch when $P_{\text{rec}}=1$ yr and predicted its H$\alpha$ emission profile. We then generated a synthetic sky image of H$\alpha$ emission for Run 22 by integrating this H$\alpha$ emission radial profile over the volume of a sphere, collapsing this sphere along one axis into a two-dimensional image before convolving this with a Gaussian with a width of 1 arcsecond to represent the typical seeing at the Liverpool Telescope (LT; Steele et al., 2004). A wedge of this spherical NSR is shown in Figure 20.

As can be seen in Figure 20, the structure does resemble the structure of the observed remnant around M 31N 2008-12a, as seen from the ground with the LT (see Figure 8 in Darnley et al., 2015). Specifically, we can see a negligible measure near the origin of the NSR (light grey) and a very low measure at the transitionary ejecta pile-up region (same light grey section), mimicking the LT observations. Then, at the inner edge of the shell, we see a vastly significant increase in the emission measure (dark grey band) as the ejecta that traversed the pile-up region collides with the extremely high density remnant shell. There is, however, a geometrical difference between the full synthetic sky image which uses a spherically symmetric model and the observed remnant around 12a which is elliptical, likely from an inclined torus or barrel-like structure.

As well as replicating the 12a NSR on the sky, this is the type of structure we would also expect to observe around other novae hosting NSRs, using ground-based facilities. Based on our full suite of simulations of NSRs, we find that detectable remnants can form around novae with very different system parameters and so should actively be searched for around all types of novae, not just those with very short recurrence periods.

7 Conclusions

We have presented a suite of hydrodynamical simulations of recurrent nova eruptions to determine how system parameters such as accretion rate, ISM density, WD temperature and initial WD mass affect the growth of a nova super-remnant. We follow the evolution of the WD from its formation mass up to either the Chandrasekhar mass (for high accretion rate systems) or the mass at a temporal upper limit (for lower accretion rate systems), and evolve the eruption properties as the mass changes. We utilised these simulations to predict the observational signatures associated with NSRs such as X-ray and H$\alpha$ emission, before comparing our simulations with the NSR observed around 12a, including the generation of a synthetic sky image. Here, we summarise the key results:

Dynamic nova super-remnants (NSR) should be found around all RNe, including those with long recurrence periods and lengthy evolutionary times, as the nova eruptions naturally drive their creation. 2. 2.

Unlike the DHO19 study, we find that radiative cooling plays a key part in the formation of dynamic NSRs, and significantly alters the density and thickness of the outer dynamic shell. 3. 3.

The creation of a dynamic NSR occurs whether the WD mass is increasing or decreasing, indicating that NSRs also exist around old novae with low mass WDs. 4. 4.

The evolving eruptions create NSRs many parsecs in radius comprising a very low density cavity, bordered by a very hot pile-up region, and surrounded by a cool, thin, high density shell. 5. 5.

A high density ISM restricts the NSR size, as does a high accretion rate; these parameters have the largest effect on NSR size. 6. 6.

The temperature of the WD and initial WD mass may have much less impact on NSR size, however NSRs grown from ONe WDs (${>}1.1\ \text{M}_{\odot}$) are significantly reduced. 7. 7.

The simulated NSRs can replicate the size of the 12a NSR and can reproduce the associated structure of H$\alpha$ emission. 8. 8.

Only NSRs grown from systems with high accretion rates will currently be observable.

NSR structures may have been overlooked within the Milky Way as they will extend across large regions of sky, far beyond their central RN. Ultimately though, the discovery of a second NSR surrounding another RN would provide strong evidence for an association between RNe and NSRs. NSRs also offer an opportunity to find unknown/unconfirmed RNe, and have the potential to point to ‘extinct’ novae where the donor has been exhausted (Darnley, 2021). Additionally, with the WD in a proportion of these systems being close to $\mathrm{M}_{\mathrm{Ch}}$ with the real possibility to explode as a SN Ia, these phenomena can also provide “a clear and persistent signpost to the progenitor-type of that SN Ia” (Darnley, 2021), and provide a mechanism for the removal of hydrogen from the immediate vicinity of a single-degenerate SN Ia (removing ${\sim}10^{6}\,\mathrm{M}_{\odot}$ of gas tens of parsecs from the central system; Harvey et al., 2016a; Darnley, 2021).

Acknowledgements

The authors would like to gratefully thank our reviewer, Michael Shara, for his insightful suggestions that helped to improve our study and strengthened our conclusions. MWH-K acknowledges a PDRA position funded by the UK Science and Technology Facilities Council (STFC). MWH-K, MJD, EJH and PAJ receive funding from STFC grant number ST/S505559/1. This work made use of the high performance computing facilities at Liverpool John Moores University, partly funded by LJMU’s Faculty of Engineering and Technology and by the Royal Society.

Data Availability

The data in this study can be shared on reasonable request to the corresponding author. This work was conducted with the Morpheus (Vaytet et al., 2007) program and analysed using the Python libraries: Numpy (Harris et al., 2020) and Matplotlib (Hunter, 2007).