The Galactic Nova Rate: Estimates from the ASAS-SN and Gaia Surveys

Though discoveries of Galactic novae date back millennia, the systematic monitoring of the sky for novae by photographic means began around the turn of the twentieth century (Duerbeck 2008). During the first half of the twentieth century, ∼2 novae were discovered per year visually and using photographic plates. With the widespread use of astronomical photography (Pickering 1893, 1895) and the advent of objective prism surveys (Duerbeck 2008), the rate of discoveries increased to ∼3 yr^–1 in the mid-twentieth century. In the 1980s and 1990s, film photography became commonly used by amateur astronomers, and the discovery rate increased to ∼4 yr^–1.

In the 2000s and 2010s, more sensitive large-format CCD and CMOS-based digital cameras became widely available to amateur astronomers, increasing the annual nova discovery rate to ∼8 yr⁻¹ (see Figure 1). It was thought that discoveries were highly incomplete (Liller & Mayer 1987), but the degree to which novae were missed owing to a shallow magnitude limit, low observational cadence, or lack of sky coverage was unknown. Therefore, it was unclear whether professional surveys with systematic observations would produce a large increase in the discovery rate.

Zoom In Zoom Out Reset image size

The All-sky Automated Survey (ASAS-3; Pojmański 2001) was one of the first CCD surveys imaging the entire sky reachable from its observing site and contributed to nova discoveries in the early 2000s. Many surveys at the time were focused on searching for supernovae, asteroids, or transiting exoplanets, so observations of transients in the Galactic plane were lacking, but surveys like ASAS-3 inspired the next generation of wide-field nova searches.

In the 2010s, large sky surveys began to contribute significantly to the discovery of Galactic novae. The Fourth Phase of the Optical Gravitational Lensing Experiment (OGLE-IV; Udalski et al. 2015) began in 2010, with high-cadence I-band observations of the central bulge discovering a large number of candidate bulge novae (Mróz et al. 2015). The New Milky Way Survey (NMW; Sokolovsky et al. 2014) began searching for transients in the northern hemisphere Galactic plane at high cadence down to V ≈ 13.5 mag in 2011. The All-Sky Automated Survey for SuperNovae (ASAS-SN; Shappee et al. 2014) started searching for transients in 2013, and then in 2017 it became the first survey to systematically observe the entire night sky, including the Galactic plane, with nearly daily cadence down to g ≈ 18.5 mag (Kochanek et al. 2017). These high-cadence and systematic observations allowed for fast-declining novae to be discovered anywhere on the sky. The Palomar Gattini-IR survey (PGIR; De et al. 2020f) began surveying the northern hemisphere Galactic plane in 2019 down to J ≈ 15.3 AB mag. As a near-IR survey, PGIR can discover highly extinguished novae not discovered in previous years (De et al. 2021b). These surveys were designed to discover transients, but the astrometrically focused Gaia (Gaia Collaboration et al. 2016) has been a surprising contributor to the discovery of Galactic novae. The broad observing filter, high angular resolution, and all-sky coverage allow Gaia to detect highly reddened novae anywhere on the sky, including the southern Galactic plane (Hodgkin et al. 2021b). These surveys have helped to increase the average discovery rate to 10.5 yr⁻¹ since 2010, with 17 spectroscopically confirmed novae in 2021, the highest number on record (see Figure 1). A higher discovery rate allows for better estimates of the global Galactic nova rate, as there is less sensitivity to model assumptions.

Two methods have been used to estimate the total rate of classical nova production in the Galaxy. The first method extrapolates a sample of discovered Galactic novae based on the estimated completeness (commonly referred to as the Galactic or direct method). The second method estimates the rate in a nearby galaxy and infers the Milky Way rate by scaling on the relative luminosities of the two galaxies (referred to as the extragalactic or indirect method).

The direct method was used to make the first prediction for the total frequency of Galactic nova eruptions. It was estimated that there should be at least, but probably much higher than, R = 50 novae per year (Lundmark 1935). Two decades later, Allen (1954) used a sample of 19 novae to estimate a Galactic rate of R ∼ 100 yr^–1. Kopylov (1955) used 23 novae thought to be within 1500 pc of the Sun observed over a 60 yr period, arriving at a rate of R = 50 Galactic novae per year. In another two decades, Sharov (1972) extrapolated from a sample of eight novae in the solar neighborhood to estimate a nova rate of R = 259 yr^–1 for the entire Galaxy. In a similar fashion, Liller & Mayer (1987) extrapolated from a sample of 17 novae thought to be within a 60° slice of the Galaxy to a total rate of R = 73 ± 24 yr⁻¹. Then, Hatano et al. (1997) assumed that novae in the Galaxy are disk dominated and arrived at an annual rate of R = 41 ± 21 yr⁻¹.

In the 1990s, estimates began to be derived from extragalactic samples of novae using the indirect method. Ciardullo et al. (1990) used the nova rate of NGC 5128 to infer the nova rate in our own Galaxy to be between R = 11 and 46 yr^–1. Van den Bergh (1991a) considered the nova rates in M31 and M33 and the globular cluster population ratios to deduce a Galactic nova rate of R ∼ 16 yr⁻¹. Della Valle & Livio (1994a) derived a Galactic rate of R = 24 yr^–1 by measuring the rates in five other galaxies and assuming that the rates were proportional to the galaxy luminosity. Finally, Darnley et al. (2006) used a survey of M31 to infer a Galactic nova rate of $R={34}_{-12}^{+15}\,{\mathrm{yr}}^{-1}$. Even more recently, Rector et al. (2022) used a 20 yr survey of M31 to indirectly estimate a slightly lower Galactic nova rate of $R={28}_{-4}^{+5}\,{\mathrm{yr}}^{-1}.$

Shafter (1997) is the first in a series of papers that extensively looked at the Galactic nova rate and estimated, by extrapolating samples of novae of varying sector size, limiting magnitude, and time period, that the Galactic nova rate is 35 ± 11 yr⁻¹. Then, Shafter et al. (2000) estimated the nova rates in M51, M87, and M101 to indirectly estimate a Galactic nova rate of $R={27}_{-8}^{+10}\,{\mathrm{yr}}^{-1}$. A couple years later, Shafter (2002) assumed that the discovered sample of m_V < 2 mag novae is complete, extrapolated to a global rate of R =36 ± 13 yr⁻¹, and also derived a rate of R ∼ 25 yr^–1 by comparing the K-band luminosity and nova rate of M31 to the Milky Way. Shafter et al. (2014) indirectly estimated a Galactic nova rate of R ∼ 25 yr^–1 by scaling the rates of 16 galaxies by their K-band luminosities. In the most recent paper of this series, Shafter (2017) assumed that the m_V < 2 mag novae are 90% complete and estimated that the most likely nova rate is $R={50}_{-23}^{+31}\,{\mathrm{yr}}^{-1}$. An even higher rate was predicted a year later when Özdönmez et al. (2018) estimated the local density of nova eruptions to predict an average estimate of the disk nova rate of $R={67}_{-17}^{+21}\,{\mathrm{yr}}^{-1}$, and by combining this with a bulge rate estimate of R = 13.8 ± 2.6 yr⁻¹ (Mróz et al. 2015), they predicted a Galactic nova rate of R ∼ 80 yr⁻¹.

These predictions from the past 40 years, along with the history of the nova discovery rate, are summarized in Figure 1. One noticeable feature is that the extragalactic predictions estimated lower rates on average than the recent predictions from Galactic data. The Hubble Space Telescope (HST) survey of M87 novae (Shara et al. 2016) derived a nova rate for M87 over three times larger than that estimated by Shafter et al. (2000). The authors argue that HST is more sensitive to faint and fast novae and that previous extragalactic nova surveys had underestimated nova rates because they neglected this harder-to-discover class of novae. Assuming this to be the case, the Galactic rate derived from the M31 sample in Darnley et al. (2006) was increased to between ∼50 and ∼70 yr^–1 (Shafter 2017).

Accounting for faint and fast novae could finally make the predictions from both the direct and indirect methods consistent, but it also appears at odds with the only modest increase in the discovery rate after the advent of large sky surveys: the implication would be that we are still discovering fewer than a quarter of the Galaxy's novae. With these new surveys, it is less plausible that time sampling is the main reason that novae are being missed, so another explanation is needed. A reasonable possibility is dust—that many novae are being hidden from optical surveys by foreground extinction. Kawash et al. (2021d) quantified the contribution of interstellar dust extinction to the optical discovery rate of Galactic novae and found that dust can hide ∼50% of the Galactic population from being discovered by observers using V- or g-band filters. This helps explain some—but not all—of the discrepancy between the discovered and predicted rates.

Luckily, the well-defined observing patterns of large time-domain surveys now make it possible to make systematic predictions of the Galactic nova rate by calculating the expected completeness in the survey data. This new era of Galactic nova rate estimates from large sky surveys began with a prediction using data from PGIR, where a sample of 11 highly reddened novae was used to derive a rate of $R={44}_{-9}^{+20}\,{\mathrm{yr}}^{-1}$ (De et al. 2021b). This is consistent with the higher rates that have been published recently, and more predictions from other large sky surveys could help bolster these higher-frequency estimates.

For the first time, in this paper we use data from multiple all-sky surveys to estimate the nova rate of the Milky Way. The use of all-sky surveys reduces the sensitivity of our results to possible differences in nova behavior between the disk and bulge (see, e.g., Della Valle & Izzo 2020).

To date, there are only two all-sky surveys that have reported a nova candidate. First, ASAS-SN became able to scan the entire night sky at a roughly 1-day cadence in 2017. This high cadence was unprecedented, allowing for the discovery of fast novae anywhere on the sky.

Gaia is the only other all-sky survey that reports nova candidates. Gaia is designed to repeatedly scan the whole sky to make astrometric measurements, and its observing pattern allows for the discovery of many transients, including novae. Usually, a pair of observations are taken 106.5 minutes apart and followed up 2–4 weeks later if the field is not Sun constrained. Though the cadence is much lower than ASAS-SN observations, the high angular resolution (006 × 018), the broad G-band observing filter, and the limiting magnitude (G < 19 mag.) make Gaia sensitive to certain novae in the plane that no other survey can detect (Hodgkin et al. 2021b). Hence, these two surveys are sensitive to different types of hard-to-detect novae: ASAS-SN to fast novae and Gaia to highly reddened novae in crowded fields. Making a rate estimation using both surveys allows us to better predict the Galactic nova rate, capitalizing on both surveys' distinct strengths.

In Section 2, we discuss the sample of novae detected by the surveys and the assumptions we make to calculate the discovery rates. In Section 3, we explain how we modeled the population of Galactic novae using a Monte Carlo simulation to estimate the Galactic nova rate. Then, in Section 4 we show the results and compare the simulated detections to the real sample. In Section 5, we compare our results to previous estimates and explore how different model assumptions can change the results. Finally, in Section 6, we summarize our results and the broader implications.

Gaia transients are reported publicly to the Gaia science alerts (GSA) website, ¹⁵ but only 25% of GSA transients have been classified spectroscopically. The vast majority of confirmed candidates are extragalactic supernovae (Hodgkin et al. 2021b), and as seen in Figure 14 of Hodgkin et al. (2021b), there are a large number of Galactic plane transients that are not classified. Therefore, the number of confirmed classical novae reported by GSA is certainly lower than the actual number detected. To quantify this, we have searched the entire archive of GSA transients for missed Galactic nova events.

In 2021 September, we performed a retroactive nova search on the ∼17,500 transients in the GSA website. Those at high Galactic latitude (b > 4°) or those with small measured outburst amplitudes (amp. < 5 mag) are unlikely to be missed Galactic novae (Kawash et al. 2021e). If the quiescent or historic magnitude is too faint to be detected by Gaia, we do not make a cut on the amplitude. We eliminated those at high latitude and those with small amplitudes, and the number of candidates decreased to 435. The bright candidates (G < 14 mag) have all been classified, with 9 out of the 12 transients being classified as classical novae. The faint candidates (G > 14 mag) include six additional confirmed and highly reddened classical novae, but there are also many unclassified transients. In the absence of extinction, an intrinsically faint nova with peak absolute magnitude M_G = −5 mag at a large Galactic distance of 20 kpc would peak at an apparent magnitude of m_G = 11.5, so the faint (G > 14 mag) nova candidates will all be heavily extinguished and appear red. Most reported transients also include a low resolution (R ∼ 100) and uncalibrated (in wavelength and flux) BP/RP spectrum, a tool that has proven to be crucial in identifying highly reddened nova candidates. We examined the spectra of the remaining unclassified candidates to identify any transient with a majority of its flux in the red RP filter or with strong emission lines. This identified five highly reddened Galactic plane transients detected between 2019 and 2021 with no spectroscopic follow-up.

We obtained spectra of these candidates using the Goodman spectrograph (Clemens et al. 2004) on the 4.1 m Southern Astrophysical Research (SOAR) telescope to look for evidence of past nova eruptions. We detected spectral features in Gaia20dfb and Gaia21dwe consistent with past nova eruptions (Kawash et al. 2021a) and list them as confirmed novae in Table 1. Because the follow-up observations were taken months to years after the eruption and the sources are likely highly extinguished, we did not detect flux from Gaia21axf, Gaia20btn, and Gaia20dfc, but we still consider these transients as Galactic nova candidates and list them in Table 2 because their positions, colors, spectral features, and light curves are all consistent with other classified classical novae. Their colors rule out the most common type of classical nova contamination (dwarf novae; Kawash et al. 2021e), but the chance of contamination from microlensing events, young stellar objects, and Be star outbursts remains. Given the rate of confirmation of other candidates, and in lieu of a more complex contamination model, we assign each of these candidates a 50% probability of being a nova for our Monte Carlo simulations described in Section 3.

We assume that this retroactive search has recovered nearly all of the nova candidates reported by Gaia and that any error in the discovery rate caused by missed nova candidates that were reported will be less than Poisson errors. As seen in Table 1, Gaia reported 7, 8, and 12 confirmed novae in 2019, 2020, and 2021, respectively, with an additional 3 unconfirmed candidates listed in Table 2. Assuming Poisson uncertainty for this sample of ∼28.5 novae yields a mean and standard deviation of the discovery rate of 9.5 and 1.8 yr⁻¹, respectively, when running the Monte Carlo simulation over the years 2019−2021. The dates listed as "Gaia t₀" in Table 1 are the epoch of Gaia's first detection, which can lag the start of the eruption by weeks to months because of Gaia's nonuniform scanning law.

The ambitious goal of ASAS-SN to observe the entire night sky daily has provided unprecedented cadence for Galactic observations of transients. ASAS-SN transients are reported publicly to https://www.astronomy.ohio-state.edu/asassn/transients.html. In 2018, ASAS-SN switched from observing in V band to exclusively observing in g band. To allow time for deep g-band reference images to be built, we only calculate the discovery rate between 2019 and 2021. Over this time period, we inspected all the ASAS-SN data for confirmed Galactic novae and found that 26 novae were detected in the transient pipeline. If there is no ASAS-SN value for t₀ or peak brightness listed in Table 1, the nova was too highly reddened and therefore too faint (g ≲ 18.5 mag) to be detected by ASAS-SN. V3730 Oph was detected at g = 16.6 mag but was never flagged in the pipeline as a transient (fainter transients have a smaller chance of being reported relative to brighter ones; see Equation (4)). We therefore do not include it as a nova in the ASAS-SN discovery rate. Kawash et al. (2021e) found that there were four more reported cataclysmic variable (CV) candidates in 2019 that could be luminous enough to be novae if they are distant enough to be significantly obscured by dust extinction. As with the Gaia candidates, we make the simple assumption that these four candidates, listed in Table 2, have a 50% probability of being classical novae and assume that these were the only potentially missed candidates over this time period. Hence, over the 3 yr considered here, there were 26 confirmed and 4 unconfirmed novae, for a total population of 28 novae. Assuming Poisson statistics yields a mean and standard deviation of the ASAS-SN discovery rate of 9.3 and 1.8 yr⁻¹, respectively.

If we combine both surveys and account for overlapping discoveries, there were a total of 39 confirmed classical novae and 7 additional candidates detected between 2019 and 2021. This yields a mean and standard deviation of the discovery rate of 14.2 and 2.3 yr⁻¹, respectively. While this observed nova rate is still far below most predictions, it is still the highest observed Galactic rate ever used to infer the total nova rate.

The positions of the model novae in the Galaxy are derived by assuming that they trace the stellar mass density of the Galaxy. The stellar density model is outlined in the appendix of Kawash et al. (2021d) and includes a thin disk, thick disk, and halo component as described in Robin et al. (2003) and a two-component, elongated, triaxially symmetric bulge from Simion et al. (2017). The ratio of the disk mass to the bulge mass, and therefore the ratio of disk to bulge novae assuming equal production per stellar mass, is 1.7. For each run of the Monte Carlo simulation, we randomly pick 1000 nova positions from a possible 10,000 positions. Because this study requires light curves to be generated in each survey at every position, we limit the total number of possible positions to 10,000. In Kawash et al. (2021d), we explored how sensitive the derived nova rate is to the stellar density model, where it was found to change results by ∼15% for the ASAS-SN rate, subdominant to other uncertainties in our calculation.

The peak absolute magnitude of each nova is assigned by randomly sampling a normal distribution, described with a mean and standard deviation of M_g,G = −7.2 ± 0.8 mag (Shafter 2017). This distribution was measured from M31 novae, and classical novae in the Galaxy could have a different luminosity function (Shafter et al. 2009; Özdönmez et al. 2018). For example, observational biases could suppress the number of bright and faint novae owing to the variable surface brightness of M31. As a test of the effect of changing the luminosity function, we also perform this analysis after doubling the variance to σ = 1.6 mag. This resulted in a <10% increase in the derived rates, leading to the conclusion that our results are not significantly sensitive to the assumed luminosity function.

The extinction along the line of sight to each nova position is calculated from the mwdust package (Bovy et al. 2016). We use the combined19 version of this extinction model, which is built by combining the Marshall et al. (2006) map of the inner Galactic plane, the Green et al. (2019) map of the northern hemisphere, and the Drimmel et al. (2003) map of the southern hemisphere. The various maps take precedence over each other where they overlap in the order they were listed above. The vast majority of the model novae (95%) lie in the Marshall et al. (2006) map region of the sky. For ASAS-SN, we query the extinction in the Sloan Digital Sky Survey (SDSS) g band, as it is the most similar to ASAS-SN's g-band filter. Extinction in the broad Gaia G-band filter is not directly accessible in mwdust, so we estimate this value by first querying the extinction in the SDSS g and i bands. The relationship between these SDSS filters and Gaia's G-band filter was fit to a third-order polynomial with the form

$$\begin{eqnarray}&&G-g=-0.1064-0.4964x-0.09339{x}^{2}+0.004444{x}^{3},\end{eqnarray} \tag{ 1 }$$

where x = g − i (Gaia Collaboration et al. 2021).

An additional concern when estimating the extinction for a model of Galactic novae is that the 3D dust maps that comprise the mwdust package only estimate the extinction out to distances where the colors of stars can be measured, and therefore they do not extend to the largest Galactic distances and the highest extinction regions. If a model nova is beyond the largest distance bin of mwdust, we add extinction by fitting the amplitude of the extinction along a line of sight and assuming that the dust is described with an underlying double exponential distribution that is a function of Galactic radius and height from the plane. We use (r₀, z₀) = (3.0,0.134) kpc for the scale length and height, respectively (Li et al. 2018). The amount of extinction is fit out to the distance where mwdust has information, and then we use the results to extrapolate to the nova distance. For 55% of simulated novae, mwdust extends to the distance of the nova, so no extinction is added. For those remaining novae beyond the 3D dust maps, the mean extinction added in the g and G bands is 0.9 and 0.6 mag, respectively, and the median extinction added in the g and G bands is 0.1 and 0.7 mag, respectively. The distributions of the estimated extinction before and after this correction are shown in Figure 2. The median nova will experience ∼10 mag of extinction in g band and ∼6 mag of extinction in G band.

Zoom In Zoom Out Reset image size

With the peak luminosity, the distance, and the extinction in g and G bands estimated, the peak apparent magnitude of each model nova can be calculated using the distance modulus. The results are shown for 10,000 model novae in Figure 3. The differences between the Gaia and ASAS-SN peak apparent magnitude distributions are entirely due to the different observing filters used; the bluer ASAS-SN g band is more vulnerable to extinction from interstellar dust compared to the broader Gaia G band. This figure shows that extinction is the largest factor in determining the apparent brightness of a Galactic nova. The dashed black distribution ignores extinction, and therefore its variance is caused by variations in the luminosity and distance of a nova; in this case, the standard deviation in the peak apparent brightness is only 1.2 mag. However, when considering extinction, the standard deviation of the peak brightness seen in Gaia's G-band filter is 5.0 mag, and in ASAS-SN's g-band filter the standard deviation is 8.8 mag. Based only on the peak brightness, ∼90% of Galactic novae are bright enough to be detected by Gaia (G < 19 mag), compared to the only ∼60% detectable by ASAS-SN (g < 18.5 mag).

Zoom In Zoom Out Reset image size

The modeling work to derive the apparent magnitude distribution of Galactic novae was largely laid out in Kawash et al. (2021d), but the detection efficiency in the survey data was only broadly estimated. Here we incorporate the model into each survey's data to more accurately estimate the detection efficiency. Each simulated nova is given a random eruption date between 2019 January 1 and 2021 December 31. While in principle novae first discovered at the start of our survey period (2019 January) could have exploded in 2018, and hence outside of our simulation, the predicted and observed number of novae in the first half of January is so low that this potential issue does not have a meaningful (≲1%) affect on our results. Though not entirely symmetric, this effect is also diminished by novae that erupt near the end of the simulation in late 2021 that would be discovered in 2022, after the end of the simulation.

The temporal evolution of the brightness is modeled by sampling a distribution of nova speeds and then constructing the full shape of the light curve from a sample of known nova light curves. We used 93 AAVSO V-band nova light curves from Strope et al. (2010) and 75 Stony Brook/SMARTS V-band light curves from Walter et al. (2012) as light-curve templates. V band was chosen because it was the most well-sampled filter from the two databases and also close to the g and G bands of ASAS-SN and Gaia. Motivated by the typical intrinsic colors of novae (e.g., van den Bergh & Younger 1987) and the fact that extinction is the key factor in determining the apparent brightness (not the luminosity), we assume a flat spectrum (V_peak = g_peak = G_peak) when transforming these templates to g and G bands.

Previous completeness studies have used several methods to estimate how long a nova is detectable, including defining a discrete number of observable days after eruption (Mróz et al. 2015) and assuming a linearly declining light curve (De et al. 2021b). By using real light curves as templates, we do not have to make simplifying assumptions about the shape of the light curves, which can be quite complex (Strope et al. 2010). However, we also carry out the analysis by using linearly declining light curves to explore how sensitive the results are to this property, and we discuss the results in Section 5.2.

We took additional steps to reduce errors in the light-curve templates. A template should only include flux from the nova eruption and not from when the nova has returned to the quiescent state or from nearby background stars. Most of the quiescent magnitudes are listed in Table 1 of Strope et al. (2010) for the novae presented in that work, so we truncate the light curve once it declines to within 3 mag of this quiescent magnitude. Walter et al. (2012) note that there are cases in the SMARTS sample of novae when the fading remnant becomes blended with a background star, and from inspecting these light curves, this does appear to be true in a few cases. We also cut off light curves within a magnitude of any prolonged plateaus many magnitudes below maximum to eliminate this contamination. If the duration of the template light curve is shorter than the survey length, we extend the template by assuming that it will decline linearly at the average pace of the light curve.

Each nova is randomly assigned a decline speed by sampling a lognormal distribution for t₂ (the time it takes a nova to decline 2 mag from peak) with mean and standard deviation equal to 18.7 and 3.2 days, respectively, based on modeling of observed novae and accounting for selection effects (Kawash et al. 2021e).

There has been a large effort to establish novae as standard candles. The relationship between the luminosity and the speed of a nova, commonly referred to as the maximum magnitude versus the rate of decline (MMRD) relation, was once thought to be a tight correlation (Capaccioli et al. 1989; della Valle & Livio 1995), but in recent years many examples of novae that do not follow the published relationships have been found, especially those that are fainter and faster (Kasliwal et al. 2011; Shara et al. 2017a). These faint, fast novae have been suggested to be a reason for a factor of three discrepancy in extragalactic nova rate estimates (Shara et al. 2016). In our primary model, we assume that the luminosity of a nova simply sets an upper bound on the speed or t₂ value. Put another way, the luminous novae are restricted to having small t₂ values, but the fainter novae are allowed all values of t₂. The boundary between the forbidden and allowed values is shown as a blue dashed line in Figure 4, with the allowed values being below this line. Any model nova with a given luminosity that is assigned a forbidden t₂ value has this parameter resampled until an allowed value is found. This forbids luminous, slow novae, as no such example has been found, but allows for any number of faint, fast novae.

Zoom In Zoom Out Reset image size

The maximum magnitude versus rate of decline for 10,000 model novae is shown in Figure 4. The distribution of t₂ is discontinuous because of the discrete t₂ values of the template light curves. Also shown on this plot is the MMRD correlation measured in della Valle & Livio (1995; the red shaded region), real Galactic nova values measured from "Gold" and "Silver" Gaia distances (denoted as the green stars; Schaefer 2018), and various extragalactic novae (Kasliwal et al. 2011; Shafter et al. 2011; Shara et al. 2016). The faint, fast novae are arguably overrepresented in our model compared to observations, so we also rerun the analysis by treating MMRD as a strict correlation and discuss the results in Section 5.2. We explored the detection efficiency of each survey in this parameter space (see Sections 3.4 and 3.5 for details on detection efficiency), and the results are shown in Figure 5. Though faint, fast novae are clearly harder to detect in the model, the difference in the recovery fraction of faint, fast novae compared to novae that adhere to the MMRD relation is no more than a factor of two. This means that it is unlikely that a large population of faint, fast novae exist in the Galaxy that have not been detected by all-sky surveys, unless such novae have a very different spatial distribution and are more embedded in dust than typical novae.

Zoom In Zoom Out Reset image size

Once each model nova is given a decline rate, we assign a light-curve template by finding the Strope et al. (2010) or Walter et al. (2012) light curve closest to the randomly assigned t₂ value. The template is then scaled to the assigned peak apparent magnitude. The model novae now have all of the information needed to inject them into the survey data, and next we discuss how that is performed for each survey.

For each nova position we used the Gaia Observation Forecast Tool ¹⁶ to provide the epochs at which Gaia observed the location, and we assume a fixed detection threshold of G < 19 mag. We then sample the template light curve at the cadence of Gaia, and several examples are shown in the right column of Figure 6. We draw the noise in the light curve from a normal distribution with standard deviation

$$\begin{eqnarray}\begin{array}{rcl}\sigma & = & 3.43779\,\mathrm{mag}-(G/1.13759)\\ & & +{\left(G/3.44123\right)}^{2}-{\left(G/6.51996\right)}^{3}+{\left(G/11.45922\right)}^{4}\end{array}\end{eqnarray} \tag{ 2 }$$

if 13 mag < G < 19 mag and σ = 0.02 mag for G < 13 mag.

Zoom In Zoom Out Reset image size

To determine whether a model nova would be reported by GSA, we start with three basic requirements (outlined in Hodgkin et al. 2021b). First, there need to be detections in both fields of view (FOVs; preceding or trailing) brighter than G < 19 mag within 40 days of each other. Second, the detected brightness needs to exceed that of all stars within a 15 radius in Gaia DR2. Third, there cannot be a G < 12 mag star within a 10'' radius in Gaia DR2.

Even if a transient satisfies these three requirements, there are additional ways for it to be missed in the pipeline. For example, a single source can have multiple source IDs, ambiguous matches in Gaia DR2, multiple sources within the core region, etc. (Kostrzewa-Rutkowska et al. 2018). It is difficult to determine when this will happen in the simulation, so we make a statistical estimate by looking at the reporting efficiency of novae in Table 1. There are 20 novae reported by GSA that were also detected by other surveys (amateurs, ASAS-SN, PGIR, etc.); however, after further inspection of the Gaia database and the Gaia observation forecast tool, there are seven additional novae (V0569 Sct, V1708 Sco, V5693 Sgr, V1710 Sco, V1674 Her, V0606 Vul, and RS Oph) that should have been reported based on the three major criteria listed above. Hence, GSA reported no more than 20/27 ≈ 74% of the novae that passed these criteria, and the true reporting efficiency is likely somewhat lower because of the strong chance of additional candidates that were unreported by other observers. We treat this estimate as a 1σ upper limit with 10% uncertainty, so for each Monte Carlo trial we assign Gaia a reporting probability for novae that pass all of the hard-coded requirements by sampling a normal distribution with mean μ = 0.67 and standard deviation σ = 0.067. This is necessarily a crude model to summarize the complex process that leads to a candidate nova being reported and undoubtedly an important source of uncertainty in our simulations. Future observations of novae and Gaia alerts can help constrain this reporting efficiency.

The cadence and limiting magnitude of ASAS-SN observations were calculated by constructing image subtraction light curves without adding the reference flux from the ASAS-SN data at the 10,000 positions of the model novae. This automatically provides a sampling of the cadence, noise uncertainty (due to the lunar cycle and weather conditions), and contamination from bright, nearby stars. The formal photometric uncertainties reported in ASAS-SN light curves tend to underestimate the true uncertainties (Jayasinghe et al. 2018), and this can particularly be true in the Galactic plane owing to crowding, even with the benefits of image subtraction. We estimated a rescaling of the uncertainties for each light curve by looking at the distribution of the ratio f_i /σ_i , where f_i is the flux and σ_i is the reported error in the flux of each camera at each position. If the error estimates are correct, then the standard deviation of this distribution should be unity. When it is larger, we increase σ_i by the factor needed to make the width of the distribution unity. For many of the random nova light curves, the rescaling is large enough to lead to a ∼1 mag reduction in sensitivity on average. Because this work is looking at Galactic novae that tend to be located at crowded low Galactic latitudes, these affects are more severe than for an extragalactic supernova study.

When injecting the model nova light curves into the ASAS-SN data, we assume that the nova would be detected if it is brighter than these rescaled 5σ upper limits on each particular epoch. Similar to the Gaia analysis, we fit a polynomial to photometric errors in ASAS-SN data, and we add noise to the light-curve templates by sampling a normal distribution with a measured standard deviation of

$$\begin{eqnarray}\begin{array}{rcl}\sigma & = & 0.08\,\mathrm{mag}+0.04(g-13)-0.04{\left(g-13\right)}^{2}/\mathrm{mag}\\ & & +0.02{\left(g-13\right)}^{3}/{\mathrm{mag}}^{2}-0.002{\left(g-13\right)}^{4}/{\mathrm{mag}}^{3}\end{array}\end{eqnarray} \tag{ 3 }$$

if 13 mag < g < 18.5 mag, and σ = 0.02 mag for any g < 13 mag.

When a transient is reported by ASAS-SN, it is cross-checked to catalogs like Gaia and the Panoramic Survey Telescope and Rapid Response System (Pan-STARRS; Chambers et al. 2016) to roughly estimate the outburst amplitude. This helps differentiate classical nova candidates from other CV outbursts, but the large pixel scale of ASAS-SN (8'' pixel⁻¹) can lead transients to have underestimated outburst amplitudes when their position is coincident with another star. To account for this, we also require detections to be 5 mag brighter (Kawash et al. 2021e) than the closest star in Gaia DR2 within half an ASAS-SN pixel, to assure that the transient would be recognized as having a large outburst amplitude.

The goal of the ASAS-SN pipeline is to discover previously unknown transients and variable stars. To avoid continually looking at known variables, the ASAS-SN pipeline does not generate candidate images for flux changes detected within 5 pixels of known Mira, long-period, or semiregular variables. A list of the positions, types, and magnitude range of known variables is acquired from VSX and OGLE and is maintained in the ASAS-SN database. We inspect this same list in the simulation and require that the model novae also satisfy this condition. If a transient is near a different type of variable, a candidate image will be generated; however, if a new transient is found near a known one, it is possible for the new transient to be confused with a known repeater and consequently not reported. To incorporate this into the model, a transient has to be a magnitude brighter than the listed magnitude range of any known variable within half an ASAS-SN pixel of the simulated nova.

The last step for a transient that is detected by ASAS-SN to be reported is that a human needs to flag the source in the data as real and previously unreported. As mentioned above, the probability of this occurring scales with the signal-to-noise ratio (S/N). The images of bright candidates (g < 15 mag) generate alert emails and are filtered into a pipeline dedicated to finding Galactic classical novae. Fainter Galactic transients are commonly reported by ASAS-SN but with lower completeness due to the number of false positives from artifacts increasing at fainter magnitudes. Generally, the closer a candidate is to the detection limit, the less likely the image is to be vetted by a human and reported as a transient. The probability for a single detection of an extragalactic supernova to be reported in ASAS-SN was studied by D. Desai et al. (2022, in preparation) and found to be

$$\begin{eqnarray}\left\{\begin{array}{ll}0.65 & {\rm{S}}/{\rm{N}}\geqslant 12\\ (1+({\rm{S}}/{\rm{N}}-12)/7)\times 0.65 & {\rm{S}}/{\rm{N}}\lt 12,\end{array}\right.\end{eqnarray} \tag{ 4 }$$

where S/N is the signal-to-noise ratio of each detection. In the simulation, each detection that satisfies the above requirements has this probability of being flagged. This detection probability is per individual epoch, so the brighter and slower transients have a much higher likelihood of being reported than the faint and fast ones.

In our primary model, we assume that novae occur in proportion to the stellar mass in both the bulge and disk. Because roughly 40% of the novae are in the bulge region (R < 3 kpc) compared to the 60% of novae that would be considered to be in the disk (R > 3 kpc), this predicts a bulge and disk rate of ≈ 10 ± 2 yr⁻¹ and ≈16 ± 2 yr⁻¹, respectively. There is an independent estimate of the bulge nova rate from the OGLE-IV survey, of 13.8 ± 2.6 yr⁻¹ (Mróz et al. 2015), which is just over 1σ higher than our value.

We find little difference (≲10%) in the detection efficiencies between the bulge and disk in our work, with ASAS-SN predicted to recover ∼37% of disk novae and ∼28% of bulge novae, and Gaia recovering ∼43% of disk novae and ∼40% of bulge novae. This means that our inferred total Galactic nova rate would not meaningfully change for reasonable differing assumptions about the bulge and disk populations of novae (Kawash et al. 2021d). As an example, Shafter & Irby (2001) estimated that the M31 nova rate per unit luminosity of the disk was 0.4 that of the bulge. Assuming a similar ratio for the mass in our model yields bulge and disk rates of ≈15 and ≈11 yr⁻¹, respectively. Hence, even a mild enhancement of the bulge nova rate with respect to the disk nova rate could bring our results into full agreement with the OGLE-IV results, but as stated above, the disagreement is small even with our standard model.

Our inferred rate of disk novae is not consistent even at the 2σ level with the average disk rate estimate of ${67}_{-17}^{+21}$ yr⁻¹ from Özdönmez et al. (2018). This rate was derived from an estimate of the local (within 1–2 kpc) nova population extrapolated to the entire disk. The origin of the disagreement is not immediately clear but could potentially arise if the distances to some of the novae in their sample were underestimated.

The rates derived in this work are perhaps most easily comparable to the PGIR estimate of ${44}_{-9}^{+20}$ yr⁻¹ (De et al. 2021b), the only other direct Galactic estimate made from a systematic large (though not all-sky) time-domain survey. However, the PGIR rate is substantially higher than those derived here, with the overlap in the posterior distribution functions <20%. Here we assess whether the discrepancy can be explained by different assumptions between the two models.

Kawash et al. (2021d) explored the sensitivity of the ASAS-SN-derived nova rates to various model assumptions and found that the rates are most sensitive to the assumed extinction law (${A}_{V}/{A}_{{K}_{s}}=13.44$ vs. ${A}_{V}/{A}_{{K}_{s}}=8.65;$ Nataf et al. 2016) and the stellar density model used to distribute the nova positions (Robin et al. 2003 vs. Cautun et al. 2020). Because a steeper reddening law increases the g-band extinction and the Cautun et al. (2020) stellar density model places more novae closer to the plane (see Figure 1 of Kawash et al. 2021d), the recovery fraction of novae in ASAS-SN decreases, resulting in a 33% increase in the derived rate. Assuming a similar increase for this work would increase the ASAS-SN derived Galactic nova rate to R ≈ 37 yr⁻¹, still consistent with our primary model at the 2.2σ level and now consistent with the PGIR rate at the 1σ level. However, the Gaia rate is much less sensitive to these changes, as it is better at finding extinguished novae in the plane, so a significant discrepancy between the Gaia and PGIR rates would remain.

As seen in Table 1 of Kawash et al. (2021d), the PGIR predicted rate is not sensitive to changes in the distribution of positions, extinction model, or bulge-to-disk ratios of novae. Those rates were derived assuming a detection efficiency of 17%, estimated from the completeness study of PGIR data (De et al. 2021b).

Arguably the largest difference between the PGIR study and this work is the assumed shapes of nova light curves and the adherence to the MMRD relation. We assume that the distribution of speed class is derived from the lognormal distribution of t₂ measured in Kawash et al. (2021e) and then use a real nova light curve with the closest value of t₂ as a template to model the fading from maximum light. De et al. (2021b) use the peak luminosity to find the corresponding t₃ value using the MMRD relationship measured in Özdönmez et al. (2018) and then assume that the apparent magnitude will fade linearly in time. We reran our analysis, this time using the PGIR method to derive the light-curve speed and shape, to see whether this could explain the discrepancy in the results. This led to a small decrease in the derived rates, with the ASAS-SN, Gaia, and joint rates changing to 26 ± 5 yr⁻¹, 21 ± 5 yr⁻¹, and 24 ± 4 yr⁻¹, respectively. This minor decrease is because the MMRD relation measured in Özdönmez et al. (2018) maps the average nova luminosity of M = −7.2 mag to a relatively slow decline t₃ ≈ 73 days, and the strict adherence to the MMRD relation does not allow for faint and fast novae. This makes novae observable for a slightly longer period of time compared to the method used in this work. But overall, the differences in the shape of the simulated nova light curves cannot explain the inconsistencies in the results.

One of the three major requirements for a nova candidate to be identified in the PGIR pipeline is that the transient has to be 3 mag brighter than any 2MASS counterparts within a radius of 10'' (De et al. 2021b). From the positions and peak brightness of PGIR-detected novae listed in Table 1 of De et al. (2021b), it appears that PGIR 20ekz and PGIR 20eig do not meet this requirement, and they could have been discovered retroactively by other means. As seen in Figure 8 from De et al. (2021b), if the sample of PGIR novae was decreased from 11 to 9 novae (by excluding 20ekz and 20eig), it would correspond to a 1σ range for the Galactic nova rate between 27 and 52 yr⁻¹, consistent with both the PGIR rate and the rate derived in this work. Hence, the rates derived in this work are inconsistent at the 1σ level with the PGIR rate, but small changes in either model would appear to bring the results into agreement at a rate of ∼30 yr^–1.

Because of uncertainties in galaxy stellar masses and the potential variation of nova rate with stellar population parameters, there are challenges associated with comparing indirectly derived rates from extragalactic nova surveys with direct rates from Galactic surveys. However, extragalactic surveys have historically been more complete in at least some dimensions, so the exercise has been common practice in the literature. As seen in Figure 1, extragalactic studies typically lead to lower Galactic rate estimates than direct measurements. The microlensing survey of M31 suggested a Milky Way rate of $R={34}_{-12}^{+15}$ yr⁻¹ (Darnley et al. 2006), consistent with the rates derived here. Additionally, the 20 yr survey of M31 recently indirectly derived a Galactic nova rate of $R={28}_{-4}^{+5}$ yr⁻¹ (Rector et al. 2022), also consistent with our rate at the 1σ level. Another common practice is to estimate the linear correlation between the nova rate and the log luminosity of a sample of galaxies. Della Valle & Livio (1994b) studied five galaxies to infer a Galactic nova rate of 24 yr⁻¹, Shafter et al. (2000) studied three galaxies to infer a Galactic nova rate of ${27}_{-8}^{+10}$ yr⁻¹, Shafter et al. (2014) gathered measurements from 16 galaxies to arrive at a Galactic rate of ∼25 yr^–1, and Della Valle & Izzo (2020) compiled measurements from 14 galaxies to infer a Galactic nova rate of ∼22 yr⁻¹. These estimates have large and unreported uncertainties, but they are notably all consistent at the 1σ level with the most likely rate predicted in this work.

Shara et al. (2016) argue that extragalactic nova surveys underestimate nova rates because they neglect to account for faint and fast novae, but the present paper shows that missed faint, fast novae are not an important source of uncertainty in the Galactic nova rate (see Section 3.2 and Figure 5) compared to other factors such as the foreground extinction. This conclusion is consistent with the lack of discovery of a substantial population of faint, fast novae after years of daily monitoring by ASAS-SN. While faint, fast novae undoubtedly exist at some level and are harder to detect, current Galactic data provide no evidence for a large yet-to-be-discovered population of these sources, and it is plausible that extragalactic studies have overestimated their importance.