THE GALACTIC NOVA RATE REVISITED

The Galactic nova rate is important in the study of a variety of astrophysical problems. For example, classical nova explosions are thought to play a significant role in the chemical evolution of the Galaxy (e.g., José et al. 2006, and references therein). In addition to producing a fraction of the ⁷Li and the short-lived isotopes ²²Na and ²⁶Al, novae are believed to be important in the production of the CNO isotopes, particularly ¹⁵N, where novae may account for virtually all of the Galactic abundance of this isotope. Thus, complete models for Galactic chemical evolution necessarily require the nova rate as an input parameter.

Novae may also play an important role as Type Ia supernova (SN Ia) progenitors (Shafter et al. 2015; Soraisam & Gilfanov 2015; Starrfield et al. 2016, and references therein). Indeed, perhaps the most promising SN Ia progenitor extant is the recurrent nova M31N 2008–12a (Tang et al. 2014; Darnley et al. 2015; Henze et al. 2015b). M31N 2008–12a has an extremely short recurrence time of just under a year, possibly as short as 6 months (Henze et al. 2015a), which constrains the accretion rate to $(2\mbox{--}3)\times {10}^{-7}\,{M}_{\odot }$ yr⁻¹ and the mass of the white dwarf to be near the Chandrasekhar limit (Wolf et al. 2013; Kato et al. 2014). These models also suggest that the white dwarf is gaining mass, and that it will reach the Chandrasekhar limit in less than 10⁶ years. The ultimate fate of M31N 2008–12a, as with all similar recurrent nova systems, depends on whether the composition of the white dwarf is CO or ONe. In the former case the system is expected to explode as an SN Ia, while in the latter case electron captures onto ²⁰Ne and ²⁴Mg will result in an accretion-induced collapse and the subsequent formation of a neutron star (Miyaji et al. 1980).

Despite its importance, the Galactic nova rate is not well established. Estimates have varied widely, from as few as 20 to as many as 260 yr⁻¹ (Sharov 1972; della Valle & Livio 1994). Recently, Mróz et al. (2015) have measured a rate of 13.8 ± 2.6 yr⁻¹ for the Galactic bulge alone based on OGLE observations. Global nova rates have been estimated both directly, by extrapolating the observed rate in the vicinity of the Sun to the entire Galaxy (e.g., see Shafter 1997, 2002), and indirectly through comparison with other galaxies (Shafter et al. 2000, 2014; Darnley et al. 2006). Shafter (2002) used the Bahcall & Soneira (1980) model for the stellar density in the Milky Way to extrapolate the rate of novae with $m\lesssim 2$, which was assumed to be complete, to faint magnitudes finding a Galactic rate of 36 ± 13 yr⁻¹. Recently, Schaefer (2014) performed a thorough analysis of the observational selection biases against the discovery of even the brightest novae. After taking these biases into account, Schaefer (2014) makes the surprising suggestion that only 43 ± 6% of Galactic novae with $m\leqslant 2$ are likely recovered.

In this paper, we reconsider the analysis presented in Shafter (2002) and Shafter (1997) by conducting Monte Carlo simulations to estimate the global Galactic nova rate and its uncertainty. The revised analysis includes the increased sample of Galactic novae available since 2000, and makes more plausible assumptions regarding the completeness of the nova sample at bright magnitudes. We conclude by comparing the latest Galactic nova estimates with those recently measured in extragalactic systems such as the nearby spiral M31, and the Virgo elliptical M87.

From 1900 to the end of 2015 there have been a total of ∼250 Galactic novae discovered brighter than m = 10, and we have compiled a list of these novae in Table 1.¹ The spatial positions of the full nova sample from Table 1 are plotted in Figure 1, confirming that the novae are concentrated in the Galactic plane, and in the direction of the Galactic center. The observed asymmetry about Galactic latitude b = 0 confirms the presence of significant and patchy absorption in the disk plane. Clearly, magnitude-limited samples of novae will be incomplete, particularly at fainter magnitudes.

Zoom In Zoom Out Reset image size

Surprisingly, novae at the brightest magnitudes appear to be biased toward the northern hemisphere. In Table 2 we show the number of bright novae ($m\leqslant 5$) discovered north and south of the celestial equator since 1900 as a function of apparent magnitude. At first glance it appears that we may be missing a significant fraction of novae at brighter magnitudes in the southern hemisphere, but the observed sample of novae is small. We can gain some insight into the statistical significance of this apparent asymmetry as follows. For a given magnitude, the probability that N or fewer novae would be discovered in a given hemisphere out of a total of M novae is given by

$$\begin{eqnarray}&&{P}_{\leqslant N,M}=\displaystyle \sum _{N=0}^{N}\displaystyle \frac{M!}{N!(M-N)!}\ {0.5}^{M}.\end{eqnarray} \tag{ 1 }$$

Table 2.  Equitorial Distribution of Bright Nova Discoveries

mlim	${N}_{\mathrm{south}}$	${N}_{\mathrm{north}}$	${N}_{\mathrm{tot}}$	P
2	2	5	7	0.45
3	3	8	11	0.23
4	6	11	17	0.33
5	12	21	33	0.16

Download table as:  ASCIITypeset image

Since the excess could have occurred in either hemisphere, we must multiply the probability given in Equation (1) by two. For m = 2 we have N = 2 and M = 7, which results in ${P}_{\leqslant \mathrm{2,7}}=0.45$. Thus, the fact that only two of the seven novae brighter than m = 2 have erupted in the southern hemisphere is not particularly surprising. Probabilities for m = 3, 4 and 5 are also given in Table 2. It seems clear that given the small number statistics, we cannot claim any statistically significant bias in bright nova detections to any one hemisphere.

Nevertheless, evidence for possible incompleteness at bright magnitudes can be appreciated by considering how the rate of discovery of novae reaching apparent magnitude m, N(m), has varied over time. Figure 2 shows the cumulative distribution of nova discovery magnitudes during the period between 1900 and 2015. Of the seven novae that reached m = 2 or brighter, six were discovered in the first half (58 years) of this period (see Table 3). Only one nova, Nova Cyg 1975 (V1500 Cyg), which reached m = 1.9, has been discovered in the second 58 year interval (actually, in the last 73 years!). We can compute the significance of this result by computing the probability that N or fewer novae with $m\leqslant 2$ would be found within any consecutive 58 year span in the 116 years since 1900. That probability is given by Equation (1), where $M-1$ novae must erupt within the same 58 year window. In this case, where N = 1, we find P = 6 × 0.5⁶ (one nova) + 0.5⁶ (no novae) = 0.11. Assuming that the true nova rate has been constant over time, a Kolmogorov–Smirnov test reveals a similar result, namely only a 7% chance that novae with $m\leqslant 2$ would be distributed as shown in Figure 2. Despite the fact that these probabilities do not rule out 100% completeness for $m\leqslant 2$ at the $2\sigma$ level, the probabilities are small, and suggest that at least one nova was likely missed in recent years. With only seven of a possible eight $m\leqslant 2$ novae being detected since 1900, the completeness becomes ∼88%.

Zoom In Zoom Out Reset image size

Table 3.  Galactic Novae with $m\leqslant 2$

Name	Date (Y, M, D)	l (°)	b (°)	mdisa	d (kpc)	${E}_{B-V}$ (mag)	Refb
V1500 Cyg	1975 8 29	89.8	−0.1	1.9	1.5 ± 0.2	0.43 ± 0.08	1, 2
CP Pup	1942 11 9	252.3	−0.4	0.5	1.5	0.25 ± 0.06	1, 2
DQ Her	1934 12 12	73.2	26.4	1.3	0.40 ± 0.06	0.08	1, 2
RR Pic	1925 5 25	271.0	−25.4	1.0	0.4	0.02	1, 2
V476 Cyg	1920 8 20	87.4	12.4	2.0	1.62 ± 0.12	0.27 ± 0.12	1, 2
V603 Aql	1918 6 8	33.2	0.8	−1.1	0.33	0.08	1, 2
GK Per	1901 2 21	151.0	−10.1	0.2	0.46 ± 0.0.03	0.29	1, 2

Notes.

^aPeak magnitude when available, otherwise discovery magnitude. ^b(1) http://asd.gsfc.nasa.gov/Koji.Mukai/novae/novae.html; (2) Darnley et al. (2012).

Download table as:  ASCIITypeset image

Although it appears counterintuitive, Schaefer (2014) points out that incompleteness at brighter magnitudes may be due in part to the evolution of how amateur astronomers survey the sky, which in recent years has turned to the use of telescopes equipped with CCD detectors rather than memorizing the sky and conducting wide-area visual observations. In addition, seasonal effects (i.e., Sun in Sagittarius) will also have diminished the observed rate of fainter novae concentrated toward the Galactic center (see Figure 3). As mentioned earlier, after considering a variety of selection effects Schaefer (2014) arrived at a completeness of just 43% for novae brighter than m = 2. Because only a summary of this work has been published, it is not possible to critically evaluate the assumptions made in arriving at this value. It does, however, seem quite surprising that more than half of the novae reaching second magnitude since 1900 could have been missed. Whether the completeness is close to 90% as estimated above, or whether Schaefer is correct that we have missed more than half of the second magnitude and brighter novae, one thing seems clear, the assumption of 100% completeness for $m\leqslant 2$ made earlier by Shafter (2002) is likely to be overly optimistic.

Zoom In Zoom Out Reset image size

In the analysis to follow, we constrain the Galactic nova rate by adopting plausible limits on the completeness of bright novae. Given that it seems difficult to understand how the completeness could be lower than the value determined by Schaefer (2014), we have adopted c = 0.43 as a lower limit on the completeness of novae with $m\leqslant 2$. The possibility that all novae with $m\leqslant 2$ have been detected since 1900 provides a hard upper limit of 100% on the completeness. We argue that the best estimate of the completeness lies between these limits, and follows from two considerations. As described earlier, the sharp drop in the number of $m\leqslant 2$ novae observed over the past ∼60 years suggests that at least one out of eight bright novae has likely been missed in recent years. If so, a value of c = 0.88 would seem to offer a reasonable estimate for the completeness of novae with $m\leqslant 2$. This estimate is supported by considering that even the brightest novae will likely be missed if they erupt within 18° (1.2 hr) of the Sun. Based on this correction alone, the completeness drops to ∼90%. Thus, in computing the models described in the following section, we simply take c = 0.9 as our best estimate of the completeness for novae with $m\leqslant 2$. For comparison, we also consider models for c = 0.43, which we take as a lower limit to the completeness of novae reaching second magnitude or brighter.

The annual discovery rate of novae brighter than m = 10 since 1900 as a function of magnitude is shown in Figure 4. The corrected $m\leqslant 2$ nova discovery rates for our estimated completeness of c = 0.9 and the lower completeness advocated by Schaefer (2014) are shown as the blue diamond and the red square, respectively. For comparison, we have also plotted the values of N(m) computed from a sub-sample of novae discovered during the period between 1900 and 1950. We find that the annual discovery rate for novae reaching second magnitude or brighter during this earlier time span is approximately twice that for the full 1900–2015 period, and is consistent with that expected after applying Schaefer's incompleteness estimate.

Zoom In Zoom Out Reset image size

In the analysis to follow, we estimate the Galactic nova rate by extrapolating the local nova rate ($m\leqslant 2$) to the entire Galaxy based on a model consisting of separate bulge and disk components. The bulge component, ${\rho }_{{\rm{b}}}$, is modeled using a standard de Vaucouleurs (1959) luminosity profile, while the disk nova density, ${\rho }_{{\rm{d}}}$, is assumed to have a double exponential dependence on distance from the Galactic center and on the distance from the Galactic plane.

Following Bahcall & Soneira (1980) for the Galactic bulge component we have

$$\begin{eqnarray}&&{\rho }_{{\rm{b}}}(r)={\rho }_{b,0}\,\displaystyle \frac{\exp [-7.669{(r/{re})}^{1/4}]}{{(r/{r}_{e})}^{7/8}},\end{eqnarray} \tag{ 2 }$$

where r is the radial distance from the Galactic center and ${r}_{e}=2.7\,\mathrm{kpc}$ is the scale parameter for the Galactic bulge. The constant ${\rho }_{b,0}$ represents the density of bulge novae at the center of the Galaxy.

For the disk component we can write

$$\begin{eqnarray}&&{\rho }_{{\rm{d}}}(z,x)={\rho }_{d,0}\,\exp [-(z/{h}_{z})-(x-{r}_{0})/{h}_{r}],\end{eqnarray} \tag{ 3 }$$

where ${\rho }_{d,0}$ is the density of novae at the position of the Sun, z is the distance of a nova perpendicular to the Galactic plane, r₀ is the distance from the Sun to the Galactic center, and $x={[{r}^{2}-{z}^{2}]}^{1/2}$ is the distance of a nova from the Galactic center in the plane of the Galaxy.

The parameters h_r and h_z are the scale lengths for the exponential distributions of novae parallel and perpendicular to the Galactic plane, respectively. Bahcall & Soneira (1980) adopted ${h}_{r}=3.5\,\mathrm{kpc}$ for the radial scale length of the disk and ${r}_{0}=8$ kpc for the distance from the Sun to the Galactic center. More recent studies have suggested a slightly shorter scale length for the disk and a somewhat larger distance to the Galactic center. Here we adopt more recent determinations of ${h}_{r}=3.0\pm 0.2\,\mathrm{kpc}$ (McMillan 2011) and ${r}_{0}=8.33\pm 0.11\,\mathrm{kpc}$ (Chatzopoulos et al. 2015).

Several estimates of the scale height of cataclysmic variables (CVs) perpendicular to the plane of the Galaxy have been made over the years. For example, Patterson (1984) found z = 190 pc for CVs in general, while Duerbeck (1984) estimated z = 125 pc for Galactic novae specifically. More recently, Revnivtsev et al. (2008) determined z = 130 pc for X-ray selected CVs, while Ak et al. (2008) find $z=158\pm 14$ pc for a large sample of Galactic CVs. We have made an independent determination of the scale height for Galactic novae using the distance estimates given in Table 4, which have been taken from Darnley et al. (2012) and Downes & Duerbeck (2000). Figure 5 shows the cumulative distribution of scale heights for the 47 novae in this sample. If we assume novae are distributed as $n(z)={n}_{0}\,\exp (-z/{h}_{z})$, then for a cumulative distribution where $N(z)={\int }_{o}^{z}n(z){dz}$, we find that $z={h}_{z}$ at $N={N}_{\max }(1-{e}^{-1})$, where ${N}_{\max }={n}_{0}{h}_{z}$. Taking ${N}_{\max }=47$, we find ${h}_{z}\simeq 250$ pc, which is somewhat larger than previous estimates. In the analysis to follow, we will adopt our estimate of h_z = 250 pc. The adopted value of h_z sets the local nova rate density, ${\rho }_{d,0}$, required to normalize the model to N(2). Since the Galactic nova rate is dependent on this normalization, unlike the value of ${\rho }_{d,0}$, it is not sensitive to the choice of h_z.

Zoom In Zoom Out Reset image size

Taken together, ${\rho }_{b,0}$ and ${\rho }_{d,0}$ determine the relative contribution of the bulge and disk to the overall Galactic nova rate. At the position of the Sun, Bahcall & Soneira (1980) find that the bulge contributes 1/800 of the stellar density. In this case, if we define $\theta (={n}_{\mathrm{disk}}/{n}_{\mathrm{bulge}}$) as the ratio of the specific nova rate of the disk population to that of the bulge population, we find that ${\rho }_{d,0}/{\rho }_{b,0}\simeq \theta /79$. Assuming $\theta =1$ and a nova scale height, h_z = 250 pc, this leads to an integrated disk-to-bulge mass ratio of ∼15. More recent models for the Galaxy suggest that the bulge component makes up a larger fraction of the Milky Way's total mass, and that the disk-to-bulge mass ratio ${M}_{d}/{M}_{b}$ is of order 6 (e.g., see McMillan 2011; Licquia & Newman 2015). Assuming the relative nova rates follow the integrated bulge and disk masses, a value of ${M}_{d}/{M}_{b}=6$ corresponds to a bulge contribution of ∼1/320 to the total stellar density at the position of the Sun. In the calculations to follow we will consider both Bahcall & Soneira (1980) models characterized by $\theta =1.0$ as well as a bulge-enhanced $\theta =0.4$ model. The latter model is equivalent to a $\theta =1.0$ model with a more massive bulge (${M}_{d}/{M}_{b}=6$), where the bulge and disk components produce novae at the same rate per unit mass.

3.1. The Absolute Magnitude Distribution of Novae

Before we can estimate the Galactic nova rate or compute the expected nova rate as function of apparent magnitude, we must specify the absolute magnitude characteristic of Galactic novae. Shafter et al. (2009) have shown that the absolute magnitudes of novae typically range from $M\simeq -5$ to $M\simeq -10$, with the mean peak absolute magnitude for the M31 and the smaller Galaxy samples being given by $\langle {M}_{V}({\rm{M}}31)\rangle \simeq -7.2$ and $\langle {M}_{V}({\rm{G}})\rangle \simeq -7.8$, respectively. To bring the M31 sample up to date, we have redetermined absolute magnitude distribution for the full sample of M31 novae through 2015 May given in the online catalog of W. Pietsch² using the extinction and color corrections given in Shafter et al. (2009). Figure 6 shows the resulting absolute magnitude distribution. The mean of the distribution remains unchanged with $\langle {M}_{V}({\rm{M}}31)\rangle =-7.2$, with a best-fitting Gaussian giving a standard deviation of 0.8 mag.

Zoom In Zoom Out Reset image size

Selection effects likely bias both the Galactic and M31 nova samples. The Galactic nova sample is affected by extinction, and is almost certainly biased toward more luminous (primarily disk) novae. The M31 sample is much larger and contains novae from both M31's bulge and disk components. However, the maximum-light magnitude is based on discovery magnitude rather than confirmed peak magnitude, and thus likely underestimates the absolute magnitude of the novae at maximum light. After taking both of these biases into account we adopt $\langle {M}_{V}\rangle =-7.5\pm 0.8$ as the best estimate currently available for the average absolute magnitude distribution of Galactic novae at maximum light.

A significant advance in our understanding of extragalactic nova rates came when Kasliwal et al. (2011) discovered a significant population of faint, but fast, novae in M31 that did not appear to follow the canonical maximum-magnitude, rate-of-decline (MMRD) relation for classical novae (e.g., see Downes & Duerbeck 2000). A typical example of a faint but fast nova is the recurrent nova M31N 2008–12a mentioned earlier, which only reaches an absolute magnitude ${M}_{V}\,\simeq -6$ at maximum light (${m}_{V}\simeq 18.5$ at the distance of M31), and fades by two magnitudes from peak in less than two days (Darnley et al. 2015). Such novae have likely been missed in the Galaxy and in previous surveys for novae in M31. As a result of their short recurrence times, such novae could make up a significant fraction of the observed nova rate, and, if properly accounted for could shift the peak absolute magnitude distribution to fainter magnitudes. To explore this possibility, we will also take the M31 absolute magnitude distribution at face value, and consider models where we adopt $\langle {M}_{V}\rangle =-7.2\pm 0.8$. Later, we will also consider the possibility that the absolute magnitudes of bulge and disk novae differ.

4.1. Direct Extrapolation

Shafter (2002) estimated the global Galactic nova rate by using the Bahcall & Soneira (1980) model to extrapolate the local nova rate (assumed complete for $m\leqslant 2$) to sufficiently faint magnitudes that the entire Galaxy is covered. If we define R(m) as the distance from the Sun to a nova of apparent magnitude m we have

$$\begin{eqnarray}&&R(m)={10}^{[1+0.2(m-{M}_{V}-{A}_{V}(R))]},\end{eqnarray} \tag{ 4 }$$

where M_V is the absolute visual magnitude at maximum light, and A_V (R) is the visual extinction suffered by a nova at distance R. The value of A_V will be a function of position in the Galaxy, and is approximated here as follows:

$$\begin{eqnarray}&&{A}_{V}(R)={a}_{V}{\int }_{0}^{R}{e}^{-z/{h}_{d}}\,{dR}={a}_{V}R({h}_{d}/z)(1-{e}^{-z/{h}_{d}}),\end{eqnarray} \tag{ 5 }$$

where the constant a_V represents the extinction in the midplane of the disk, and h_d is the scale height of the obscuring dust layer, assumed to drop off exponentially perpendicular to the disk plane. Here, we adopt an average extinction of a_V = 1 mag per kpc in the Galactic midplane, with a scale height perpendicular to the plane, h_d = 100 pc (Spitzer 1978). Equations (4) and (5) are solved iteratively to determine R(m).

The number of novae per year visible to a given magnitude, m, is then given by

$$\begin{eqnarray}&&N(m)=2{\int }_{0}^{2\pi }{\int }_{0}^{\pi /2}{\int }_{0}^{R(m)}({\rho }_{{\rm{d}}}+{\rho }_{{\rm{b}}})\,{R}^{2}\,{dR}\,\cos \,b\,{db}\,{dl},\end{eqnarray} \tag{ 6 }$$

where b and l are the Galactic latitude and longitude, respectively, and ${\rho }_{{\rm{b}}}(r)$ and ${\rho }_{{\rm{d}}}(z,x)$ from Equations (2) and (3) can be cast in terms of R by noting that

$$\begin{eqnarray}&&r={[{r}_{0}^{2}+{R}^{2}-2{{Rr}}_{0}\cos b\cos l]}^{1/2},\end{eqnarray} \tag{ 7 }$$

$z=R\,\sin \,b$, and $x={[{r}^{2}-{z}^{2}]}^{1/2}$.

Figure 7 shows the expected increase in the observed nova rate as a function of apparent magnitude for the full 1900–2015 sample. The extrapolation is accomplished by integrating Equation (4) as a function of apparent magnitude where we assume the average nova has an average absolute magnitude at maximum light of $\langle {M}_{V}\rangle =-7.5$. The solid line shows our model for the expected increase in the nova rate with apparent magnitude. As in Shafter (2002), the model has initially been normalized to the observed value of N(2) assuming 100% completeness for $m\leqslant 2$. The bulge and disk contributions to the overall nova rate are shown as the red dotted and blue dotted–dashed curves, respectively. Given that the bulge contributes so little to the nova density at the position of the Sun, these values essentially reflect the nova rate densities of disk novae. The agreement between the observed rates and the model is quite good for novae brighter than m = 2, with the observed nova rates falling off with increasing apparent magnitude as expected for a disk distribution ($\mathrm{log}\,N\propto 0.4m$). At fainter magnitudes the expected incompleteness in the observed nova rate becomes increasingly apparent. We note that, at magnitudes fainter than $m\simeq 6$, there is a hint that the contribution of bulge novae is beginning to kick in.

Zoom In Zoom Out Reset image size

Despite the excellent fit of the model to novae with $m\leqslant 2$, as we have discussed in Section 2, the observed nova rates for even these bright novae are likely to be incomplete. Figure 8 shows our model fits to the incompleteness-corrected values of N(2), shown by the blue diamond and the red square for c = 0.9 and c = 0.43, respectively. These models assume our best estimate for the mean absolute magnitude of Galactic novae, $\langle {M}_{V}\rangle =-7.5$, and produce global Galactic nova rates of ∼48 and ∼98 per year for the c = 0.9 and c = 0.43 models, respectively. We note that the rate based on an assumed completeness of 43% for novae brighter than second magnitude is almost exactly what we would expect if we normalized the model to the observed rate during the period between 1900 and 1950, as shown in Figure 4, and assumed 100% completeness during that interval. In both Figures 7 and 8, the overall Galactic nova rates have been approximated as in Shafter (2002) by extrapolating the models to sufficiently faint magnitudes that the entire Galaxy is essentially covered. The slightly higher rates compared with those in Shafter (2002) result from our updated values of h_r and r₀.

Zoom In Zoom Out Reset image size

Significant limitations of the direct extrapolation approach adopted by Shafter (2002) and reviewed above are that this method does not provide an uncertainty in the derived nova rates, and that it requires that a specific absolute magnitude for the model novae be specified rather than allowing for a distribution of absolute magnitudes to be considered. A roughly equivalent but superior approach is to extrapolate the local population of novae to the entire Galaxy using a Monte Carlo simulation.

4.2. Monte Carlo Simulations

In our Monte Carlo simulations, we distribute simulated bulge and disk novae following the scaling laws given in Equations (2) and (3) for a range of trial Galactic nova rates, ν. We then record the number of novae brighter than second magnitude, ${N}_{\nu }(m\leqslant 2)$, that are produced over the 116 year period covered by the observations. The magnitude of the ith nova is calculated assuming the extinction given by Equation (5) over a distance R_i from the Sun given by

$$\begin{eqnarray}&&{R}_{i}={[{r}_{i}^{2}+{r}_{0}^{2}-2{r}_{i}{r}_{0}\sin {\theta }_{i}\cos {\phi }_{i}]}^{1/2},\end{eqnarray} \tag{ 8 }$$

where θ and ϕ are the usual spherical polar coordinates with origin at the center of the Galaxy. For the disk component, ${\theta }_{i}=\pi /2-{\tan }^{-1}({z}_{i}/{x}_{i})$, while we assume a uniform distribution in cos θ for the bulge component. In both cases we assume that the Galaxy is axially symmetric (i.e., a uniform distribution in ϕ).

We run the simulation $M={10}^{5}$ times for each trial nova rate and record the number of times ${N}_{\nu ,i}(m\leqslant 2)$ the ith simulation matches the number of novae believed to have reached second magnitude or brighter over the past 116 years. This latter number is simply given by ${N}_{c}(m\leqslant 2)=[{N}_{\mathrm{obs}}/c]$, where c is the assumed observational completeness for $m\leqslant 2$. To account for variation in the number of $m\leqslant 2$ novae observed over the past 116 years, our Monte Carlo simulations sample a Poisson distribution (mean of 7) for N_obs. Similarly, to account for the uncertainty in c, our Monte Carlo simulation samples a distribution c_i which is normally distributed about c with standard deviation ${\sigma }_{c}$. Figure 9 shows the completeness function and the resulting ${N}_{c,i}(m\leqslant 2)$ distribution for our c = 0.9 models, where we have assumed ${\sigma }_{c}=0.2$. The most likely estimate of the global nova rate for a given assumed completeness distribution is then given by the value of ν that produces the largest number of matches between ${N}_{\nu ,i}(m\leqslant 2)$ and ${N}_{c,i}(m\leqslant 2)$. Specifically, the probability of a given nova rate is given by

$$\begin{eqnarray}&&P(\nu )=\displaystyle \frac{{\displaystyle \sum }_{i=1}^{M}{\delta }_{{N}_{\nu ,i},{N}_{c,i}}}{{\displaystyle \sum }_{\nu =0}^{\infty }{\displaystyle \sum }_{i=1}^{M}{\delta }_{{N}_{\nu ,i},{N}_{c,i}}},\end{eqnarray} \tag{ 9 }$$

where δ is the Kronecker delta function.

Zoom In Zoom Out Reset image size

We have run an array of Monte Carlo simulations to explore how the choice of model parameters affects the Galactic nova rate. We initially assumed that the nova rate per unit mass is the same in the disk and bulge (i.e., $\theta =1$), and considered both a disk-to-bulge mass ratio ${M}_{d}/{M}_{b}=15$ as found by Bahcall & Soneira (1980) in their pioneering study of the Milky Way, and a bulge-enhanced model with ${M}_{d}/{M}_{b}=6$. The latter model is equivalent to a Bahcall & Soneira (1980) model with a higher specific bulge nova rate given by $\theta =0.4$, which is the value found by Shafter & Irby (2001) based on the spatial distribution of novae in M31. For comparison, we have also computed a pure disk model ($\theta =\infty$) despite the fact that available observations do not support a heavily disk-dominated nova population.

Figure 10 shows a plot of $P(\nu )$ as a function of the trial global nova rate ν for our most plausible models. Four $c=0.9\pm 0.2$ models are shown representing differences both in the contribution of the Milky Way's bulge and disk components and in the assumed absolute magnitude distribution of novae. The solid red curve shows the nova rate distribution for a Galactic nova population characterized by $\langle {M}_{V}\rangle =-7.5\pm 0.8$ and an assumed disk-to-bulge ratio of 15 ($\theta =1$), while the shaded distribution represents our preferred model with an enhanced bulge component ($\theta =0.4$) that contributes 1/6 that of the Milky Way's disk (e.g., see McMillan 2011; Licquia & Newman 2015). The peak of this distribution corresponds to a most likely nova rate of ${52}_{+32}^{-23}$ novae per year, where the uncertainties are 1σ errors of an assumed bi-Gaussian distribution (Buys & De Clerk 1972). Since the local novae ($m\lesssim 2$) are predominately disk novae, an increase in the bulge-to-disk ratio has the effect of shifting the peak of the probability distribution to somewhat higher nova rates.

Zoom In Zoom Out Reset image size

The broken lines in Figure 10 show nova rate distributions for the two disk-to-bulge ratios assuming the fainter absolute magnitude distribution observed for M31, $\langle {M}_{V}\rangle =-7.2\pm 0.8$ (see Figure 6). These models may be more appropriate for our Galaxy if a significant population of faint and fast novae have escaped detection in previous extragalactic nova surveys (e.g., Kasliwal et al. 2011). If so, $\langle {M}_{V}\rangle =-7.5\pm 0.8$ may overestimate the peak luminosity of the Galactic nova distribution. As expected, a lower assumed mean nova luminosity reduces the volume sampled by $m\leqslant 2$ novae, increases the extrapolation to the entire Galaxy, and results in overall Galactic nova rates as high as ∼65 per year.

In addition to our preferred $c=0.9\pm 0.2$ models, we have also computed model rates based on the low completeness of $c=0.43\pm 0.06$ suggested by Schaefer (2014). As expected, these models result in considerably higher nova rates that reach of order 100 per year in the case of the $\langle {M}_{V}\rangle =-7.5\pm 0.8$ models, or even higher if we adopt the lower mean luminosity of $\langle {M}_{V}\rangle =-7.2\pm 0.8$. As discussed earlier, we consider it extremely unlikely that the completeness of second magnitude and brighter novae could be below 50%. Thus, the nova rates that result from the $c=0.43\pm 0.06$ models likely represent firm upper limits to the Galactic nova rate.

A summary of our full array of models is given in Table 5. Each model is described in columns 1–3, which give the ratio of the specific disk-to-bulge nova rates θ, the disk-to-bulge mass ratio ${M}_{d}/{M}_{b}$, and the corresponding disk-to-bulge nova density ratio at the position of the Sun, (${\rho }_{{\rm{d}}}/{\rho }_{{\rm{b}}}{)}_{\odot }$. The derived nova rates are not sensitive to the adopted model parameters such as the bulge scale parameter, r_e, or the disk scale length, h_r, and height, h_z. The local nova rate density ${\rho }_{d,0}$, however, does depend on the assumed scale height of novae in the disk. Values of ${\rho }_{d,0}$ for three representative values of h_z are given in Table 6.

Table 5.  Galactic Nova Rate Estimates

				$c=0.9\pm 0.2$			$c=0.43\pm 0.06$
			Bulge	Disk	Total	Bulge	Disk	Total
(θ)	(${M}_{d}/{M}_{b}$)	(${\rho }_{{\rm{d}}}/{\rho }_{{\rm{b}}}{)}_{\odot }$	Rate (yr−1)	Rate (yr−1)	Rate (yr−1)	Rate (yr−1)	Rate (yr−1)	Rate (yr−1)
$\langle {M}_{V,\mathrm{disk}}\rangle =-7.5\pm 0.8$, $\langle {M}_{V,\mathrm{bulge}}\rangle =-7.5\pm 0.8$
1.0	15	800	${3.0}_{-1.4}^{+1.9}$	${45}_{-21}^{+28}$	${48}_{-22}^{+30}$	${5.7}_{-1.9}^{+2.5}$	${84}_{-28}^{+37}$	${90}_{-30}^{+40}$
0.4 (1.0)	15 (6)	320	${7.4}_{-3.4}^{+4.6}$	${45}_{-21}^{+27}$	${52}_{-24}^{+32}$	${14}_{-4.6}^{+6.3}$	${84}_{-27}^{+38}$	${98}_{-32}^{+44}$
$\infty$	$\infty$	$\infty$	...	${45}_{-19}^{+27}$	${45}_{-19}^{+27}$	...	${86}_{-30}^{+37}$	${86}_{-30}^{+37}$
$\langle {M}_{V,\mathrm{disk}}\rangle =-7.2\pm 0.8$, $\langle {M}_{V,\mathrm{bulge}}\rangle =-7.2\pm 0.8$
1.0	15	800	${4.0}_{-1.8}^{+2.3}$	${60}_{-27}^{+35}$	${64}_{-29}^{+37}$	${7.6}_{-2.5}^{+3.2}$	${112}_{-37}^{+48}$	${120}_{-40}^{+51}$
0.4 (1.0)	15 (6)	320	${9.7}_{-4.3}^{+6.0}$	${58}_{-26}^{+36}$	${68}_{-30}^{+42}$	${19}_{-6.3}^{+8.0}$	${111}_{-38}^{+48}$	${130}_{-44}^{+56}$
$\infty$	$\infty$	$\infty$	...	${60}_{-27}^{+37}$	${60}_{-27}^{+37}$	...	${112}_{-37}^{+49}$	${112}_{-37}^{+49}$
$\langle {M}_{V,\mathrm{disk}}\rangle =-7.8\pm 0.8$, $\langle {M}_{V,\mathrm{bulge}}\rangle =-7.2\pm 0.8$
1.0	15	800	${2.3}_{-1.2}^{+1.8}$	${34}_{-18}^{+26}$	${36}_{-19}^{+28}$	${4.3}_{-1.4}^{+1.9}$	${65}_{-21}^{+29}$	${69}_{-23}^{+31}$
0.4 (1.0)	15 (6)	320	${5.7}_{-2.6}^{+3.6}$	${34}_{-15}^{+21}$	${40}_{-18}^{+25}$	${11}_{-3.6}^{+4.7}$	${64}_{-21}^{+28}$	${75}_{-25}^{+33}$
$\infty$	$\infty$	$\infty$	...	${35}_{-16}^{+21}$	${35}_{-16}^{+21}$	...	${65}_{-23}^{+30}$	${65}_{-23}^{+30}$

Download table as:  ASCIITypeset image

Table 6.  Local Nova Rate Density Estimates

	c = 1.0	$c=0.9\pm 0.2$	$c=0.43\pm 0.06$
hz	${\rho }_{d,0}$	${\rho }_{d,0}$	${\rho }_{d,0}$
(pc)	(kpc−3 yr−1)	(kpc−3 yr−1)	(kpc−3 yr−1)
$\langle {M}_{V,\mathrm{disk}}\rangle =-7.5$, $\langle {M}_{V,\mathrm{bulge}}\rangle =-7.5$
125	0.18	0.21	0.39
150	0.16	0.19	0.34
250	0.11	0.13	0.24
$\langle {M}_{V,\mathrm{disk}}\rangle =-7.2$, $\langle {M}_{V,\mathrm{bulge}}\rangle =-7.2$
125	0.23	0.26	0.50
150	0.21	0.25	0.49
250	0.15	0.17	0.32
$\langle {M}_{V,\mathrm{disk}}\rangle =-7.8$, $\langle {M}_{V,\mathrm{bulge}}\rangle =-7.2$
125	0.12	0.13	0.28
150	0.11	0.12	0.25
250	0.074	0.082	0.18

Download table as:  ASCIITypeset image

We note that Galactic nova rates determined in our Monte Carlo simulations are slightly lower than those found in the direct extrapolation described in the previous section. This difference is the result of the fact that for a given apparent magnitude, the volume, $V({M}_{V})$, sampled is not directly proportional to M_V, but rather $V({M}_{V})\,\propto \,{10}^{-0.6{M}_{V}}$. Thus, the average volume sampled for a Gaussian distribution of absolute magnitudes is larger than the volume sampled assuming a single average absolute magnitude alone. In other words, $\langle V({M}_{V})\rangle$ > $V(\langle {M}_{V}\rangle )$. The larger volume sampled for novae with $m\leqslant 2$ results in a smaller extrapolation, and therefore a lower overall nova rate.

4.2.1. The Effect of Differing Nova Populations

For some time there has been speculation that that there may be two distinct populations of novae: a bulge population characterized by relatively slow and dim novae, and a disk population characterized by somewhat brighter novae that lie closer to the Galactic plane and generally evolve more quickly (della Valle et al. 1992). The evidence for two populations finds some support from observations showing that novae can be divided into two classes based upon the character of their spectra shortly after eruption. Although all novae display prominent Balmer emission shortly after eruption, the so-called "Fe ii novae," are characterized by relatively narrow (∼1000 km s⁻¹ full-width at half-maximum (FWHM)) Fe iiemission features that exhibit P-Cyg profiles near maximum light, while the "He/N novae" show prominent and broad ($\gtrsim$ 2500 km s⁻¹ FWHM) He i and N emission features, but without the Fe ii emission (Williams 1992). The work of della Valle & Livio (1998) showed that the He/N novae cluster close to the Galactic plane and tend to be fast and bright relative to the Fe ii class.

Given the possible existence of two classes of novae, we have also computed nova rate models that are characterized by different average peak luminosities and specific nova rates in the Galactic bulge and disk. Specifically, we have considered a two-component Galaxy model where disk novae are characterized by $\langle {M}_{V}\rangle =-7.8$ and bulge novae by $\langle {M}_{V}\rangle =-7.2$. Since the $m\leqslant 2$ disk novae used in the normalization are more luminous, they extend to a larger volume of the Galaxy. Thus, both the density of novae in the vicinity of the Sun (${\rho }_{d,0}$) and the extrapolated global Galactic nova rate are significantly reduced. Bulge-dominated ($\theta =0.4$) models produce generally higher rates because the normalization is set by nearby novae that belong overwhelmingly to the disk, whereas disk dominated models ($\theta =\infty$) produce lower rates because there is no contribution from the bulge.

4.2.2. Bulge Rates

In addition to summarizing Galactic nova rates produced by the various models, Table 5 also breaks down the corresponding bulge (and disk) contributions to the overall rate. As with the overall rates, the bulge nova rates also depend on the adopted Galaxy model, absolute magnitude distribution, and completeness at bright magnitudes c. Our $c=0.9\pm 0.2$ models predict Galactic bulge rates ranging between ∼3 yr⁻¹ and ∼10 yr⁻¹ depending on assumed parameters. The lowest bulge rates result from our $\theta =1$, $\langle {M}_{V}\rangle =-7.5$ models, with the $\theta =0.4$, $\langle {M}_{V}\rangle =-7.2$ variations producing significantly higher rates. As expected, the $c=0.43\pm 0.06$ models produce bulge rates that are approximately a factor of two higher. The bulge rate of 13.8 ± 2.6 found by Mróz et al. (2015), would appear most consistent with either the bulge-enhanced ($\theta =0.4$), $c=0.9\pm 0.2$ models or the $\theta =1$, $c=0.43\pm 0.06$ models, which produce bulge nova rates of ∼7–10 yr⁻¹ and 14–19 yr⁻¹, respectively.

4.2.3. The Predicted Nova Rate as a Function of Apparent Magnitude

In addition to estimating the overall Galactic nova rate, our Monte Carlo simulations have been used to predict how the number of Galactic novae should increase with apparent magnitude. Assuming the most probable overall nova rates from Table 5, we show in Table 7 the predicted increase in the number of novae visible as a function of apparent magnitude for our $\theta =1.0$ and $\theta =0.4$ models and assumed completenesses of $c=0.9\pm 0.2$ and $c=0.43\pm 0.06$. These predictions can be compared with the results from several ongoing and planned all-sky surveys such as ASAS-SN³ and ZTF⁴ once they are available, and used to differentiate between the various Galactic nova rate models presented here.

Table 7.  Predicted Visibility of Galactic Novae

	$c=0.9\pm 0.2$		$c=0.43\pm 0.06$
	$\theta =1.0$	$\theta =0.4$	$\theta =1.0$	$\theta =0.4$
${m}_{\mathrm{lim}}$	N (yr${}^{-1})$	N (yr${}^{-1})$	N (yr${}^{-1})$	N (yr${}^{-1})$
$\langle {M}_{V,\mathrm{disk}}\rangle =-7.5$, $\langle {M}_{V,\mathrm{bulge}}\rangle =-7.5$
2.0	0.1	0.1	0.1	0.1
4.0	0.4	0.4	0.7	0.7
6.0	1.5	1.5	2.8	2.8
8.0	5.2	5.4	9.7	10
10.0	14	14	25	27
12.0	23	25	44	48
14.0	32	35	60	66
16.0	38	42	72	79
18.0	42	46	80	87
20.0	45	49	84	92
$\langle {M}_{V,\mathrm{disk}}\rangle =-7.2$, $\langle {M}_{V,\mathrm{bulge}}\rangle =-7.2$
2.0	0.1	0.1	0.1	0.1
4.0	0.4	0.4	0.7	0.7
6.0	1.6	1.6	3.0	3.0
8.0	5.8	5.8	11	11
10.0	16	17	30	32
12.0	29	31	55	59
14.0	41	44	77	84
16.0	50	54	94	102
18.0	56	60	105	114
20.0	59	63	111	121
$\langle {M}_{V,\mathrm{disk}}\rangle =-7.8$, $\langle {M}_{V,\mathrm{bulge}}\rangle =-7.2$
2.0	0.1	0.1	0.1	0.1
4.0	0.3	0.3	0.7	0.7
6.0	1.3	1.4	2.6	2.6
8.0	4.5	4.7	8.7	8.9
10.0	11	12	21	22
12.0	19	20	35	38
14.0	25	27	47	51
16.0	29	33	56	61
18.0	32	36	62	67
20.0	34	38	65	70

Download table as:  ASCIITypeset image

In recent years evidence has been building that extragalactic nova rates, which have been determined primarily through synoptic surveys often with sparse temporal sampling, may have been systematically underestimated. In particular, Curtin et al. (2015) and Shara et al. (2016) have recently argued that the nova rate in the giant Virgo elliptical galaxy M87 is likely 2–3 times larger than previously thought, with the latter authors suggesting that the K-band luminosity-specific nova rates (LSNRs) of all galaxies may be ∼3–4 times higher than the previous average of ∼2 novae per year per 10¹⁰ solar luminosities in K (Shafter et al. 2014).

As mentioned earlier, the work of Kasliwal et al. (2011) has shown that a population of faint and fast novae, which deviate from the canonical MMRD relation, may exist. Owing to their intrinsically low peak luminosities and their rapid declines, such novae have almost certainly been missed in previous magnitude-limited and low cadence synoptic surveys for novae in M31. Because these novae deviate so strongly from the assumed MMRD, they have not been properly accounted for when the surveys were corrected for completeness. In most surveys, the incompleteness, and thus the final nova rates, have likely been underestimated.

The generally accepted nova rate for M31 is ${65}_{-15}^{+16}$ (Darnley et al. 2006). This value has been called into question recently by Chen et al. (2016) who have computed population synthesis models to estimate nova rates in galaxies with differing star formation histories and morphological types. They estimate a global nova rate for M31 of 97 yr⁻¹. In another study, Soraisam et al. (2016) corrected the nova samples of Arp (1956) and Darnley et al. (2006) for bias against the discovery of fast novae in these synoptic surveys. Based on these corrections, they estimate a global nova rate for M31 of order 106 yr⁻¹. Thus, the most recent estimates suggest the nova rate in M31 could be as high as ∼100 yr⁻¹. The stellar mass of M31 relative to the Galaxy is not precisely known, but we can compare estimates for the integrated K-band luminosity, which should reflect the mass difference. For the Galaxy, we adopt ${M}_{V}=-20.6$ (Bahcall & Soneira 1980), while for M31 we have m_V = 3.44 (de Vaucouleurs et al. 1991) and ${(m-M)}_{0}=24.38$ (Freedman et al. 2001), yielding ${M}_{V}({\rm{M}}31)=-20.94$. Since both galaxies have similar morphological types (Sbc), their integrated V − K colors should be similar, and in both cases we adopt $\langle V-{\text{}}K\rangle \simeq 3.25$ (Aaronson 1978). Thus, we estimate ${M}_{K}({\rm{M}}31)=-24.19$ and ${M}_{K}(\mathrm{MW})=-23.85$, which corresponds to a luminosity ratio, ${L}_{K,\mathrm{MW}}/{L}_{K,{\rm{M}}31}\simeq 0.73$. Assuming the relative nova rates follow the relative luminosities, based on a comparison with M31, we expect a nova rate in the Galaxy of between ∼50 yr⁻¹ and ∼70 yr⁻¹, depending on whether we adopt the Darnley et al. (2006) or the recent estimates. Generally speaking, these estimates are in good agreement with an average of the rates given in Table 5. The corresponding LSNRs for both M31 and the Galaxy are $\sim 8.5\pm 1.5$ novae per year per 10¹⁰ ${L}_{\odot ,K}$. This value is in excellent agreement with LSNRs of ${9.2}_{-3.0}^{+2.7}$ and 8.6 ± 2.7 novae per year per 10¹⁰ ${L}_{\odot ,K}$ recently found by Shara et al. (2016) and Mróz (2016) for M87 and the Galactic bulge, respectively.

By considering how the observed nova rate in the Galaxy varies with apparent magnitude, Shafter (2002) estimated the overall Galactic nova rate by extrapolating the observed rate in the vicinity of the Sun ($m\lesssim 2$), which was assumed to be complete, to the entire Galaxy using a two-component disk-plus-bulge model similar to that employed by Bahcall & Soneira (1980) to model the stellar density in the Milky Way. In this paper, we have reconsidered and improved the Shafter (2002) analysis in two important respects. First, we have considered important corrections for the incompleteness in the observed sample of nearby novae with $m\leqslant 2$, and second, we have employed a Monte Carlo analysis to extrapolate the local nova rate to the entire Galaxy. The Monte Carlo analysis has the advantage of allowing us to consider a distribution of nova absolute magnitudes rather than a single mean value as in Shafter (2002) and to better assess the uncertainties in our derived Galactic nova rate estimates. In addition to these two principal improvements, we have also updated the sample of novae that have become available over the past 15 years (albeit with no new novae discovered with $m\leqslant 2$), and updated the values for several Galaxy model parameters.

Galactic nova rate estimates based upon our array of models have been summarized in Table 5. We have considered disk-to-bulge mass ratios, ${M}_{d}/{M}_{b}=15$, as found by Bahcall & Soneira (1980), as well as a more massive bulge model characterized by ${M}_{d}/{M}_{b}=6$. The latter model is equivalent to a Bahcall & Soneira (1980) model where the specific nova rate in the disk is just 0.4 times that of the bulge ($\theta =1.0$ versus $\theta =0.4$). Our models have been normalized assuming correction factors of $c=0.9\pm 0.2$ and $c=0.43\pm 0.06$ for the observed fraction of novae that reached second magnitude or brighter since 1900. The $c=0.9\pm 0.2$ models represent our best estimate of the actual completeness of novae with $m\leqslant 2$ observed since 1900, with the $c=0.43\pm 0.06$ models representing a likely lower limit to the completeness, and thus an upper limit to the nova rate. We have also considered models assuming two different nova absolute magnitude distributions. Based on both Galactic and M31 nova samples, we consider $\langle {M}_{V}\rangle =-7.5\pm 0.8$ to represent the best estimate available for the absolute magnitude distribution of Galactic novae. For comparison, we also considered a somewhat fainter absolute magnitude distribution given by $\langle {M}_{V}\rangle =-7.2\pm 0.8$, which is characteristic of the observed distribution in M31. The latter absolute magnitude distribution may be more appropriate if there is a significant population of faint and "fast" novae that have hitherto largely escaped detection in existing surveys. Our models produce Galactic nova rates that typically range from ∼50 per year to in excess of 100 per year in the case where we are observing only 43% of the novae reaching m = 2 or brighter. We have also considered models where the absolute magnitudes of the disk and bulge nova rates differ. If disk novae are significantly more luminous than bulge novae (e.g., $\langle {M}_{V,\mathrm{disk}}\rangle =-7.8\pm 0.8$ and $\langle {M}_{V,\mathrm{bulge}}\rangle =-7.2\pm 0.8$), then somewhat lower nova rates result, falling in the range of ∼35 to ∼75 per year.

Since we currently know little about how nova properties may vary with population, and the completeness estimates are similarly uncertain, arriving at a definitive estimate for the Galactic nova rate based on existing data is not possible. Taking an average of the most plausible $c=0.9\pm 0.2$ rates from Table 5 (the $\theta =1$ and $\theta =0.4$ models) yields a value of ${50}_{-23}^{+31}$ yr⁻¹, which we adopt as the best estimate currently available for the nova rate in the Galaxy. Rates on the order of 100 per year are possible in the event that we have missed roughly half of the novae that have reached second magnitude or brighter over the last century. Despite the large uncertainties, the rates derived in the present study represent significant increases over recent estimates, and bring the Galactic nova rate into better agreement with that expected based on comparison with the latest results from extragalactic surveys.

I thank an anonymous referee for valuable suggestions that resulted in a significant improvement in the analysis. I have also benefited from fruitful discussions with C. T. Daub, M. Henze, and B. E. Schaefer who first recognized that the historical record of the brightest novae was likely incomplete. C. T. Daub, M. Henze and P. Mróz also provided valuable comments on an earlier version of the manuscript. Financial support through NSF grant AST1009566 is gratefully acknowledged.