© ESO 2018

1. Introduction

Cataclysmic variables (CVs) are binary systems composed of a white dwarf accreting from a low-mass stellar companion filling its Roche lobe (Warner 1995). Cataclysmic variables come in a bewildering variety of types and subtypes corresponding to their variability properties and spectra. The overarching distinction is between the dwarf novae, which show repeated outbursts with an amplitude ≳2 optical magnitudes on timescales of weeks to decades, and the nova-like systems, whose light curves remain roughly steady on the same timescales.

The outbursts of dwarf novae have long been known to originate in the accretion disk surrounding the white dwarf (Smak 1971; Osaki 1974) due to a mechanism identified by Meyer & Meyer-Hofmeister (1981). This instability occurs when the temperature is low enough in the accretion disk that hydrogen recombines. The steep dependence of the opacity with temperature in this regime triggers a thermal and a viscous instability that leads the disk to cycle through two states. In the eruptive state, the disk has a high temperature >10⁴ K, hydrogen is highly ionized, and the mass accretion rate Ṁ is higher than the mass transfer rate Ṁ_t from the companion star. In the quiescent state, the disk has a temperature ≲3000 K, hydrogen is mostly neutral, and Ṁ < Ṁ_t. The disk instability model (DIM) aims at exploring the consequences of this instability on disk accretion and explaining the variety of observed light curves (Osaki 1996; Lasota 2001).

A fundamental test of the DIM is whether it reproduces the distinction between dwarf novae and nova-likes. In nova-likes, the whole disk should be hot enough for hydrogen to be highly ionized. This translates into a minimum mass accretion rate above which a disk of given size is stable. The system is a dwarf novae if this criterion is not met. Hence, the DIM can be tested against observations by deriving the average mass accretion rate and the size of the disk R_out in dwarf novae and novalikes. Importantly, this instability criterion does not depend on the details of how matter is transported in accretion disks, which remains a highly debated issue.

Smak (1982) showed that this criterion roughly separates nova-likes from dwarf novae in the (R_out, Ṁ) plane, using approximate magnitudes and distances for a dozen CVs. Since then, only Schreiber & Lasota (2007) have updated this work, for a sample of ten CVs. The importance of this test was revived by Schreiber & Gänsicke (2002), who pointed out that the HST FGS parallax distance of SS Cyg, the archetypal dwarf nova, led to a revised mass transfer rate that placed this outbursting system in the stable region of the (R_out, Ṁ) plane. This led to much hand-wringing (Schreiber & Lasota 2007; Smak 2010) until a radio VLBI parallax distance (Miller-Jones et al. 2013), confirmed with Gaia (Ramsay et al. 2017), firmly placed the system in the unstable region (as anticipated by Lasota 2008).

Although there are dozens of CVs with relatively well-known binary parameters, a thorough investigation of this test has yet to be performed. Ramsay et al. (2017) proposed to test the DIM by taking the absolute magnitude at a specific stage of a dwarf nova outbursts (at the end of long outbursts or in stand-stills) to compare the derived mass accretion rate to the critical rate for the system, under the assumption that both should be nearly equal. Good agreement was found for the sample of a dozen CVs with distances in the Gaia first data release. Its main advantage is that the required optical magnitude is straightforward to extract from the light curves. However, a major drawback of this approach is that it relies on assumptions on how matter is transported in the disk because one needs to predict when the mass accretion rate happens to be equal to the critical rate during the outburst cycle. Hence, it is a much less fundamental and robust test of the basic premises of the DIM.

Indeed, an important difficulty in using the fundamental test is that it requires extracting the average mass transfer rate from CV optical light curves, which in many cases are unevenly, inhomogeneously, and incompletely sampled. The situation is somewhat easier for X-ray binaries, for which the DIM also applies because X-ray monitors provide continuous coverage and X-rays provide a closer measure of the bolometric luminosity in those systems than optical magnitudes do for CVs, which in outburst emit primarily UV radiation. Coriat et al. (2012) were able to derive average mass accretion rates for ~50 X-ray binaries and showed that the DIM, modified to include X-ray irradiation (van Paradijs 1996), separates consistently transients from steady systems.

We carry out, for the first time on CVs, a similar analysis to that performed on X-ray binaries by Coriat et al. (2012) to test the fundamental prediction of the DIM: the separation of CVs into stable and unstable systems is determined by the orbital period (a proxy for disk size) of the binary and mass transfer rate from the companion. We use a subsample of the CV light curves gathered by Otulakowska-Hypka et al. (2016). We consider all types of dwarf nova outbursts, including so-called super-outbursts that are not supposed to be described by the standard version of the DIM (see, e.g., Lasota 2001, and references therein), but which have the same instability criterion. We first validate the conversion from optical magnitudes to accretion rates against model light curves (Sect. 2). The sample of CVs and our approach to analyzing their light curves is described in Sect. 3. Section 4 discusses the results and possible sources of errors. We conclude on the validity of the DIM.

2. Model

2.1. Converting disk flux to magnitudes

The accretion disk in CVs is geometrically thin and radiatively efficient. We model the disk emission as the sum of local black-body spectra at the local effective temperature T_eff, which is a function of radius R and time. The flux is$$\begin{array}{c}\hfill F_{\nu }=\frac{cosi}{d^2}{\int }_{R_{\mathrm{in}}}^{R_{\mathrm{out}}}{B}_{\nu }2\pi R\mathrm{d}R,\end{array}$$(1)

where$$\begin{array}{c}\hfill B_{\nu }=\frac{2h{\nu }^3}{c^2}{[exp\left(\frac{h\nu }{k{T}_{\mathrm{eff}}}\right)-1]}^{-1},\end{array}$$(2)

i the disk inclination, and d the distance to the binary. The black-body approximation is good when compared to using stellar spectra (Smak 1989). The DIM provides the evolution of T_eff as a function of radius and time. The disk emission is converted to an optical flux in the Johnson V band, i.e.,$$\begin{array}{c}\hfill F_{\mathrm{V}}=\frac{cosi}{d^2}{\int }_{R_{\mathrm{in}}}^{R_{\mathrm{out}}}{w}_{\mathrm{V}}(T_{\mathrm{eff}})\sigma T_{\mathrm{eff}}^{4}2\pi R\mathrm{d}R,\end{array}$$(3)

where w_V(T) is defined as$$\begin{array}{c}\hfill w_{\mathrm{V}}(T)=\frac{1}{\sigma T^4}\frac{\int f_{\mathrm{V}}{B}_{\nu }\mathrm{d}\nu }{\int f_{\mathrm{V}}\mathrm{d}\nu },\end{array}$$(4)

and f_V is the filter response from Mann & von Braun (2015), as implemented in the SVO Filter Profile Service¹. The absolute magnitude is then given by$$\begin{array}{c}\hfill M_{\mathrm{V}}=-2.5log\left[{\left(\frac{d}{10\mathrm{pc}}\right)}^2\left(\frac{F_{\mathrm{V}}}{Z_{\mathrm{V}}}\right)\right],\end{array}$$(5)

where Z_V = 3562.5 Jy is the V band flux zero point in the Vega system. Numerical integration of Eq. (3) for a given radial distribution of T_eff in the disk, using the bolometric correction Eq. (4), gives the absolute magnitude M_V Eq. (5) of the disk. The same procedure can be used for any filter in the database.

2.2. Disk parameters

The ratio of outer disk size to binary separation R_out/a depends only on the mass ratio q. Using Paczyński (1977) or Lin & Papaloizou (1979) gives comparable results. The typical disk size is between (2–7) × 10¹⁰ cm for 80 mn ≤ P_orb ≤ 10 h. The disk circularization radius in units of binary separation a, R_circ/a, also depends only on q.

The minimum inner disk radius is set to the white dwarf radius, parametrized as (Nauenberg 1972)$$\begin{array}{c}\hfill R_{\mathrm{WD}}=7.8\times {10}^8{[{\left(\frac{1.44M_{\odot }}{M}\right)}^{2/3}-{\left(\frac{M}{1.44M_{\odot }}\right)}^{2/3}]}^{1/2}\mathrm{cm}.\end{array}$$(6)

However, the magnetic field of the white dwarf can truncate the inner disk before it reaches the white dwarf surface. Following the standard approach (Frank et al. 2002), the truncation radius is set by$$\begin{array}{c}\hfill R_{\mathrm{mag}}\approx {\left(\frac{{\mu }^4}{GM{\dot{M}}_{\mathrm{in}}^2}\right)}^{1/7},\end{array}$$(7)

where $\mu =B_{\text{WD}}{R}_{\text{WD}}^3$ the magnetic moment of the white dwarf and Ṁ_in the mass accretion rate at the inner radius. The truncation radius is close to or at the white dwarf radius in outburst, when Ṁ_in is high. It can be much larger than the white dwarf radius in quiescence. Aside from changing the outburst properties, this has two main consequences for the present work: first, a cold stable disk can exist if the truncation radius is large and second, the varying inner disk radius changes the relationship between magnitude and mass transfer rate.

Limb darkening changes the inclination dependence compared to the simple cos i dependence Eq. (1). We obtain the inclination dependence by using the magnitude correction derived from more elaborate disk spectral models by Paczynski & Schwarzenberg-Czerny (1980) $$\begin{array}{c}\hfill \mathrm{\Delta }M=-2.5log\frac{F(i)}{\langle F\rangle }=-2.5log[(1+\frac{3}{2}cosi)cosi],\end{array}$$(8)

where 〈F〉 = F(i = 0°)/2 is the flux averaged over all inclinations. This works well for i ≲ 75°.

2.3. Relationship of Ṁ–M_V during a dwarf nova cycle

The goal is to reconstruct the mass transfer rate from the optical light curve. Assuming no mass is lost from the disk except by accretion onto the white dwarf, then the mass transfer rate is equal to the time average of the mass accretion rate Ṁ(R, t) through the disk. Ideally, we would use the bolometric luminosity of the disk to reconstruct Ṁ_in. This is not feasible because of lack of multiwavelength data covering UV (where most of the accretion luminosity is emitted) and/or sufficient time coverage of the outbursts. We must therefore rely on optical magnitudes as a function of time, typically V (Fig. 1).

Fig. 1. Absolute magnitude in U, B, V, R, I bands (blue to red, top panel) and mass accretion rate onto the white dwarf Ṁ(Rin) during a dwarf nova outburst. The model parameters are M = 0.6 M⊙, Ṁt = 1016 g s−1, Rout ≈ 2 × 1010 cm, i = 0°. The dashed line shows the Ṁ(Rin) reconstructed from the V magnitude light curve assuming Ṁ does not depend on radius. The reconstructed mass transfer rate is 9 × 1015 g s−1 after averaging over the outburst cycle.

The radial distribution of T_eff varies in a complex fashion during the dwarf nova cycle as fronts propagate through the disk (see, e.g., Hameury et al. 1998). Schematically, T_eff changes from a steady-state-like R^−3/4 profile in outburst to a nearly flat profile in quiescence. The result is a hysteresis between the optical magnitude and the mass accretion rate (Fig. 2). Ideally, we would build a dwarf nova model for each system so that we can precisely obtain the conversion factor from magnitude to Ṁ_in (Fig. 2). We used a much simpler method, at the price of a modest error in accuracy.

Fig. 2. Relation between the absolute V magnitude and the inner mass accretion rate Ṁin for the model shown in Fig. 1. The disk samples the black curve clockwise during the outburst cycle. The dashed line shows the relationship using the stationary disk approximation.

In the following, we assume that Ṁ does not depend on radius, i.e., that the evolution behaves as a succession of steady states with Ṁ = Ṁ_in. This is not correct, as discussed above, but, as we shall prove, the corresponding error is entirely acceptable given the other uncertainties. For a stationary disk (Frank et al. 2002),$$\begin{array}{c}\hfill \sigma T_{\mathrm{eff}}^4=\frac{3}{8\pi }\frac{GM\dot{M}}{R^3}[1-{\left(\frac{R_{\mathrm{in}}}{R}\right)}^{1/2}],\end{array}$$(9)

such that the absolute magnitude depends only on mass accretion rate Ṁ, white dwarf mass M, inner disk radius R_in, and outer disk radius R_out. The integration to get M_V for given (Ṁ, M, R_in, R_out) is straightforward. We checked that our routine gives Ṁ–M_V relationships that are identical to those shown in Paczynski & Schwarzenberg-Czerny (1980) and Smak (1989), who used the same assumptions, and to the M_V values found by Tylenda (1981), who used a more elaborate spectral model. The relationship can be inverted to give Ṁ as a function of M_V although the rightmost term in the expression of T_eff (related to the notorque inner boundary condition) requires solving an equation by iteration. The Ṁ inferred from the absolute V magnitude is shown as dashed red lines in Figs. 1–2. As we verify below, averaging this Ṁ over an outburst cycle gives a reconstructed mass transfer rate that is close to the true input mass transfer rate Ṁ_t.

2.4. Validation against model light curves

The error made in using the steady-state approach can be quantified with model disk light curves with parameters covering a range of possible values for CVs. Optical light curves (U, B, V, R, I) were obtained by following the evolution of the disk with the code described in Hameury et al. (1998). Table A.1 lists the parameters used for the models. The size of the disk is free to vary during the outburst as the angular momentum from the disk is removed by tidal forces from the companion. The inner disk radius is kept fixed except in the last six models where it is left free to vary according to magnetic truncation with μ ~ 10³⁰ G cm³.

The input mass transfer rate Ṁ_t was compared to the reconstructed Ṁ_filter from the various optical light curves, averaged over an outburst cycle. The reconstruction used the average outer disk radius, the exact value of the white dwarf mass, and set R_in to the value given by Eq. (6) (even when the inner disk radius varied due to magnetic truncation). Inclination effects were not considered. The reconstructed rate nearly always underestimate Ṁ_t, all the more so that the filter is redder, but not by much: for example the average Ṁ_V/Ṁ_t = 0.85 ± 0.10 for the models in Table A.1. This bias, due to the assumptions of the model, is negligible compared to the other sources of error, notably the inclination (Sect. 4.4).

3. Optical light curves

3.1. Source list and system parameters

For the dwarf novae, we analyzed a subset of the light curves from the dataset collected and presented in Otulakowska-Hypka et al. (2016). We required a value for P_orb, a parallax distance d, and at least a five-year timespan of observations to include it in the analysis. The parallax π and its associated uncertainty were taken from the Gaia Data Release 2 (Gaia Collaboration 2016, 2018; Lindegren et al. 2018). The superb quality of the Gaia-based distances remove what was previously the major source of uncertainty in deriving the mass transfer rate. Table A.2 lists the dwarf novae systems and the parameters that we adopted. The same approach was followed for nova-likes (Table A.3), selecting systems with a Gaia parallax and using light curves from the American Association of Variable Star Observers (AAVSO) database (Kafka 2017). AM Her systems (polars) were excluded as there is no disk in those systems, where accretion is entirely controlled by the high magnetic field of the white dwarf. The values for q = M₁/M₂, M₁, M₂, i, and the extinction A_V have been mined from a variety of sources (see the reference list below Table A.2). In some cases these values may vary significantly for the same source; in many others we did not find an estimated value in the literature and we had to supplement those as described below. When several values were available, we generally kept the value from the latest reference.

The italics in Tables A.2–A.3 distinguish the values that we assumed, or derived from other parameters, from those that we took from the literature. If only M₁ or M₂ were available we used the mass ratio q to deduce the other. When neither were available, we estimated M₂ from the evolutionary sequence with P_orb of Knigge et al. (2011), taking an estimated relative error ΔM₂/M₂ = 0.2 and deducing M₁ from q. When q was unavailable, we adopted M₁ = 0.75 M_⊙ for consistency with Knigge et al. (2011). The error on M₁ and M₂ was propagated from the error on q wherever necessary. We only selected systems with i ≤ 80° as the magnitude correction is incorrect for edge-on systems. When we could not find an estimated value for i, we took i = 56.7° ± 20°: this value of i corresponding to zero magnitude correction. We took Δi = ±10° when no error was available on the system inclination i. The median extinction A_V at the location and parallax distance of each object was obtained from the 3D map of Galactic reddening by Green et al. (2018). A rough error was estimated by taking the maximum difference between the median A_V and the 16% or 84% percentile (i.e., the 1σ interval). We used values from the other references listed below Table A.2 for the few binaries whose location is not covered by this map. We did not find values for ST Cha and AT Ara. We took the A_V along the line of sight from Schlafly & Finkbeiner (2011). When no error was available we took ΔA_V = 0.1. We note that the 3D reddening map of Green et al. (2018) gives a range A_V = 0.010 (16% percentile) to 0.063 (84% percentile) in the direction and at the distance of SS Cyg, which is consistent with the A_V = 0.062 ± 0.016 specifically derived for this important source by Ritchey et al. (2013). We used the 3D map median value A_V = 0.017.

3.2. Light curve analysis

Ideally, light curves should be regularly sampled on suborbital timescales and cover several outburst cycles for the dwarf novae. However, most light curves are irregularly sampled with, for example, the frequency of observations increasing during outbursts or during dedicated campaigns. Light curves may also feature orbital modulations or eclipses. The data must be homogenized in time to derive the time-averaged optical magnitude of the binary. To do this, we averaged the data over a one-day sliding window, incrementing by one hour between each bin. This smoothes over orbital modulations without affecting the rise times and decay times to outburst. To avoid large gaps with missing data, we manually selected the portion of the light curve that we deemed long enough to average over multiple outburst cycles while remaining densely sampled. The start and stop days we chose are listed in Table A.2.

We define the filling fraction f of the resulting light curve as the fraction of its bins that contain a measurement. For the dwarf novae, we assumed the missing measurements corresponded to zero flux. We verified that assuming that the flux is equal to an observationally defined nonzero quiescent value makes no difference to our results. However, we may be missing outbursts, which would have a larger impact on the estimated mass accretion rate. With this in mind, we separated our dataset between light curves that have f > 0.5, minimising the impact of missing data, and light curves with f < 0.5, where the uneven sampling causes a larger uncertainty. We discuss the bias this introduces below (Sect. 4.4).

For the nova-likes for which the source variability is more limited, we assumed that the missing measurements were at the average flux seen at other times. Hence, the filling fraction f carries much less weight for these systems. Otherwise, we followed the same analysis as dwarf novae. The time interval over which we analyze the light curves are indicated in Table A.3: for VY Scl systems we took the flux during steady maxima of the light curve. Table A.3 also lists the mean V magnitude of the selected portion of the light curve in order to ease comparison with values in the literature for those steady systems.

The magnitudes were then corrected for extinction, distance, and disk inclination (taken to be in the binary plane). We did not correct for contributions to the V-band flux from the bright spot, companion star, and white dwarf, as we discuss below (Sect. 4.4). The absolute magnitudes were converted to mass accretion rates, following the procedure described in Sect. 2, and time averaged to obtain the mass transfer rate.

We calculated the statistical error on the mass transfer rate through Monte Carlo sampling of the errors on the mass of the binary components, inclination, parallax and A_V. The parallax and extinction probability distributions are taken to be Gaussians whereas the mass distributions and cos i are sampled uniformly. Our statistical error estimate is the smallest 68% confidence interval from the distribution of Monte Carlo sampled mass transfer rates, which is close to a lognormal.

4. Results

4.1. The critical mass transfer rate

Figure 3 compares the calculated mass transfer rates against the critical mass transfer rate (Fig. A.1 includes the system names to help localize them in the diagram). The critical mass transfer rate above which a cataclysmic variable with a disk size R is hot and stable is$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{crit}}^{+}=8.07\times {10}^{15}{R}_{10}^{2.64}{M}_1^{-0.89}\mathrm{g}{\mathrm{s}}^{-1}\end{array}$$(10)

Fig. 3. Mass transfer rates of CVs compared to the stability criterion. Systems above the (red) upper solid line are hot and stable. Systems below the lower (blue) line indicate cold, stable disks if the white dwarf magnetic field B ≥ 105 G. The dashed line represents the expected secular mass transfer rate (Knigge et al. 2011). Square symbols indicate Z Cam type dwarf novae; (red) stars indicate nova-likes. Dwarf novae with f ≥ 0.5 are shown in black and those with f < 0.5 are in gray. Figure A.1 includes system names.

using the fit for a disk with solar composition given in Lasota et al. (2008), where R₁₀ = R/10¹⁰ cm and M₁ the white dwarf mass in solar masses. Similarly, the critical mass transfer rate below which a CV with an inner disk radius R is cold and stable is$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{crit}}^{-}=2.64\times {10}^{15}{R}_{10}^{2.58}{M}_1^{-0.85}\mathrm{g}{\mathrm{s}}^{-1}.\end{array}$$(11)

This is usually too low to be of consequence if the disk extends to the white dwarf surface (low R). However, if the inner disk radius is truncated by the white dwarf magnetic field then, combining with Eq. (7)$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{crit}}^{-}=8.78\times {10}^{13}{\mu }_{30}^{0.85}{M}_1^{-0.70}\mathrm{g}{s}^{-1},\end{array}$$(12)

where μ₃₀ = μ/10³⁰ G cm³ the magnetic moment of the white dwarf. A reasonably strong magnetic field can thus truncate the disk and keep it cold and stable if the mass transfer rate is low. This can be used to place an upper limit on the magnetic field of the white dwarf in a dwarf nova.

Figure 3 shows linear fits to the values of ${\dot{M}}_{\text{crit}}^{+}$ and ${\dot{M}}_{\text{crit}}^{-}$ obtained for each object and the width of the line covers the range of individual values. The plotted ${\dot{M}}_{\text{crit}}^{-}$ corresponds to a white dwarf magnetic field B = 10⁵ G or μ ≈ 3 × 10³² G cm³, which is the typical magnetic field above which CVs are identified as intermediate polars (DQ Her type). For ${\dot{M}}_{\text{crit}}^{+}$, we have$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{crit}}^{+}\approx 3.5\times {10}^{16}{P}_{\mathrm{orb}}^{1.6}\mathrm{g}{\mathrm{s}}^{-1},\end{array}$$(13)

where P _orb is defined in hours. This is roughly recovered by noting that $R_{\text{out}}\propto a\propto P_{\text{orb}}^{2/3}{M}_1^{1/3}$ up to some slowly varying function of q. This also shows that the dependence of ${\dot{M}}_{\text{crit}}^{+}$ on M₁ basically cancels out when R_out is replaced by P_orb in Eq. (10). Except for P_orb, which is usually precisely known, the errors in system parameters have very little influence on ${\dot{M}}_{\text{crit}}^{+}$.

4.2. Dwarf novae

All the dwarf novae are well within or consistent with the unstable region delimited by ${\dot{M}}_{\text{crit}}^{+}$, clearly confirming a key prediction of the DIM. There are eight systems for which the estimated Ṁ_t is above the stability limit. However, for seven of these systems, the estimated Ṁ_t is still compatible with the unstable region within the statistical uncertainty (AT Ara, EM Cyg, ER UMa, SY Snc, V1159 Ori, FO Per, V516 Cyg, and V1316 Cyg). CN Ori is the only system for which the estimated error bar is inconsistent with the critical transfer rate: this system has an estimated Ṁ_t = (5.4–14.8) × 10¹⁷ g s⁻¹ to compare to ${\dot{M}}_{\text{crit}}^{+}=2.7\times {10}^{17}\mathrm{g} {\mathrm{s}}^{-1}$. A likely reason for this discrepancy is the median A_V = 0.48 provided by the 3D reddening map, which is flagged as unreliable for such a nearby source. The extinction could be much lower: indeed the best-fit value from the model of Green et al. (2018) is A_V = 0.083 and Urban & Sion (2006) list A_V ≈ 0. The estimated Ṁ_t becomes compatible with ${\dot{M}}_{\text{crit}}^{+}$ if there is negligible extinction. The error on the inclination, i = 67° ± 3°, may also be underestimated. Given this and the level of the systematic uncertainty on our measurements (see Sect. 4.4), we conclude that none of the systems are incompatible with the stability limit.

SU UMa stars, dwarf novae with superoutbursts, are concentrated below the period gap, as usual, except for VW Vul and ES Dra, which are also classified as Z Cam types in Simonsen et al. (2014). Their classification as SU UMa types in Downes et al. (2001) is likely to be incorrect.

Z Cam types tend to have mass transfer rates that are higher than other dwarf novae. Their light curves have standstills that have long been argued to be due to fluctuations of the mass transfer rate bringing it above the stability limit. One would thus expect Z Cam systems to be close to the stability limit (Meyer & Meyer-Hofmeister 1983; Buat-Ménard et al. 2001). Indeed, fluctuations by a factor 3 would push many of the Z Cam systems very close to or above the stability limit. The lack of Z Cam types below the period gap, even though some systems are very close to the stability limit, suggests that such fluctuations do not occur in these systems, most likely for reasons related to the type of companion star. Above the period gap, some Z Cam systems are somewhat further away than a factor 3 from the stability limit (PY Per, V426 Oph, and WW Cet). The classification of WW Cet as a Z Cam stars is uncertain in Downes et al. (2001) but considered certain in Simonsen et al. (2014). The filling fraction of the PY Per light curve, f ≈ 0.34, may lead us to underestimate Ṁ_t if the flux always remains at a high level. For instance, assuming that the missing flux values are at the average of the measured values increases Ṁ_t by a factor ≈2 (Sect. 4.4). Its disk inclination is also currently unknown and may be high. V426 Oph is listed as a Z Cam, but also as a DQ Her type in Downes et al. (2001) and as a novalike system in the General Catalog of Variable Stars (GCVS). Disk truncation by the white dwarf field may lead us to underestimate the true Ṁ_t (see Sect. 4.3 below). The nova-like classification is more puzzling, although downward fluctuations of the mass transfer rate in V426 Oph might push the system into the cold stable regime with a white dwarf magnetic field B ≈ 10⁶ G. Inversely, not all of the systems that are very close to the stability limit in Fig. 3 are classified as Z Cam types (e.g., FO Per, CY Lyr, and CN Ori), but many U Gem stars could be “unrecognized Z Cam stars” (Warner 1995).

We also find there is considerable variation of Ṁ_t at given P_orb. Some of it is undoubtedly due to inaccurate values of the parameters adopted for some of the systems or biases in the light curve analysis (see Sect. 4.4). Although there is a general trend for higher values of Ṁ_t with increasing P_orb, the values differ significantly from the theoretical secular mass transfer rate derived by Knigge et al. (2011), especially above the period gap. Some of these variations may be due to the mechanism for angular momentum losses from the system, which drives the secular mass overflow rate from the Roche lobe and is much less certain above the gap than below, where it is mainly due to gravitational radiation. In fact, in Fig. 3 the evolutionary tracks correspond to a version of the angular-momentum-loss mechanisms rescaled to correspond to observations, which leads nevertheless to “two glaring inconsistencies” (Knigge et al. 2011) with the observed distribution of dwarf novae and nova-likes. Some of these inconsistencies may be due to variations in the mass transfer rate on timescales longer than the outburst cycle but smaller than the secular timescale, such as cycles in the stellar activity.

Finally, we also find that some binaries are well below ${\dot{M}}_{\text{crit}}^{-}$ (e.g., EI Psc and EZ Lyn), indicating that the white dwarf magnetic field is constrained by the DIM to B < 10⁵ G in those systems or else they would be cold and stable.

4.3. Nova-likes

The nova-likes are clearly at higher mass transfer rates than dwarf novae and are consistently above the DIM stability limit defined by the red line in Fig. 3. Closer inspection reveals six systems in our nova list (Table A.3) with an estimated ${\dot{M}}_{\mathrm{t}}<{\dot{M}}_{\text{crit}}^{+}$ (BK Lyn, BG CMi, V603 Aql, IX Vel, UX UMa, and AE Aqr). BK Lyn is probably misclassified. It is the only novalike below the period gap in Downes et al. (2001) but its recent light curve shows variability typical of a ER UMa dwarf nova (Patterson et al. 2013). We tagged it as a dwarf nova in Fig. 3. The others are statistically consistent with being in the stable region of the diagram, except for AE Aqr.

AE Aqr is a DQ Her type system well known for propelling material outside the system owing to the very fast rotation of the white dwarf (Wynn et al. 1997). In this case, we may be severely underestimating the actual mass transfer rate from the secondary because of non-negligible mass loss. There are other systems classified as DQ Her systems (intermediate polars, e.g., BG CMi). The inner disk is truncated by the white dwarf magnetic field even if it does not rotate fast enough to propel material out of the system as in AE Aqr. We do not take into account this truncation and this affects our estimated Ṁ_t for this type of system. For illustration, consider a steady accretion disk around a 1 M_⊙ white dwarf with Ṁ = 10¹⁷ g s⁻¹ and an outer radius R_d = 4 × 10¹⁰ cm. The absolute magnitude of the disk is M_V ≈ 6.7 if the inner disk is truncated at ≈ 6 × 10⁹ cm by a magnetic field B ≈ 10⁶ G Eq. (7). The estimated Ṁ_t is only 20% lower than the actual value, which is well within our systematic errors, if the disk is assumed to be truncated at ≈ 10⁹ cm, corresponding to B ≈ 10⁵ G instead of the true value B = 10⁶ G. However, the estimated mass transfer rate is ≈ 2 × 10¹⁶ g s⁻¹ if this magnitude is translated to Ṁ assuming the inner disk extends to the white dwarf surface at ≈5 × 10⁸ cm. We thus caution that the presence of a magnetic field B ≳ 10⁵ G may lead us to underestimate the true mass transfer rate in DQ Her systems, and in particular in BG CMi, where B ≈ 4 to 10 × 10⁶ G (Chanmugam et al. 1990). Finally, we note that TV Col, another DQ Her system, shows infrequent and short outbursts but according to Hameury & Lasota (2017) these cannot be attributed to the thermal-viscous instability of the disk.

4.4. Sources of error and bias

The error bars in Fig. 3 take into account statistical errors in the system parameters: M₁, M₂, i, π, A_V. The uncertainty on the distance and disk inclination dominate the statistical error. However, there are other sources of systematic errors and biases.

First, some of these errors and biases result from our theoretical assumptions. The stationary disk assumption in our model biases Ṁ_t by ≈15% to lower values (Sect. 2.4). We also underestimate the true mass transfer rate if the inner disk is truncated by the magnetic field of the white dwarf, as discussed above. We may also underestimate Ṁ_t if a significant percentage of the mass in the disk is lost to a wind, although estimates of this fraction indicate that the typical wind loss rate represents ≲10% of the mass accretion rate in the hot state (e.g., Noebauer et al. 2010). Inversely, the mass transfer rate may be enhanced during the outburst (Hameury et al. 2000; Smak 2008) as a consequence, for example of the irradiation of the donor star, but this is debated (e.g., Viallet & Hameury 2008). We did not compensate for any of these effects.

Second, the light curve analysis also systematically underestimates Ṁ_t. To calculate the average mass transfer rate < Ṁ >, we assume that Ṁ = 0 when there are no measurements and mitigate this source of error by taking light curves with high filling fractions f, i.e., so that outbursts are not missed. Were we instead to assume that Ṁ = < Ṁ > when measurements are missing, the mass transfer rate would be increased to (2 − f) < Ṁ >. This assumption is clearly very unlikely looking at the light curves, but taken at face value the mass transfer rates would be underestimated by 50% with f = 0.5. A more pernicious error is that caused by high flux outliers in a poorly sampled portion of the light curve. A few odd points can thus induce an overestimate of the mass transfer rate. We have minimized this by manually inspecting the light curves to remove glaring outliers (a dozen points altogether). We also investigated the influence of our bin-size (from 0.25 to 12 h) and averaging window size (from 0.25 to 2 days). We find differences ≲50% in Ṁ_t with no obvious trend as to over- or underestimating Ṁ_t.

Third, the extinction may be much more uncertain that we assumed, notably for nearby binaries for which the reddening measurements of Green et al. (2018) are less reliable for lack of stars or for which we used the total extinction along the line of sight. In general, we found good agreement between the values we obtained from the 3D map and those we could collect in the CV literature (e.g., Ak et al. 2007). Yet, very different values of A_V are sometimes found for the same binary. For instance, RW Tri has A_V ≈ 0.3 in Table A.3 and in Benedict et al. (2017) but is listed with A_V ≈ 0.8 in Ak et al. (2007) and with A_V = 0.3−0.7 in Schreiber & Lasota (2007). This uncertainty can affect Ṁ_t for some systems but, in our experience, has no impact on the conclusions.

Finally, we may have overestimated the mass transfer rate by spuriously including contributions from the white dwarf, bright spot, and companion. A heuristic approach to removing these is to subtract some suitably defined minimal flux from the light curve. The assumption is that this minimal flux is dominated by the constant contributions from the stars and bright spot. We chose this flux as the threshold above which 85% of the measurements are contained in the cumulative flux distribution. This choice gave a good estimate of the average quiescent flux. We found that this leads to a decrease of the estimated Ṁ_t of at most 50%, where most systems are barely affected. As a cross-check, we estimated the contribution from the bright spot and the stars theoretically (Fig. 4). The absolute V magnitude of the companion and of the accretion-heated white dwarf are taken from Knigge et al. (2011). For the bright spot, we assumed accretion at the critical rate and computed the absolute V magnitude assuming a temperature of 15 000 K and a luminosity$$\begin{array}{c}\hfill L_{\mathrm{BS}}=\frac{GM{\dot{M}}_{\mathrm{t}}}{2R_{\mathrm{out}}}[1-\frac{R_{\mathrm{out}}}{R_{\mathrm{L}1}}].\end{array}$$(14)

Fig. 4. Absolute V magnitude of the companion star (full line) and of the accretion-heated white dwarf (dashed line) as a function of orbital period (Knigge et al. 2011). The dotted line indicates the absolute magnitude of the bright spot, calculated for Ṁ = Ṁcrit and a white dwarf mass M = 0.75 M⊙.

In nearly all cases the theoretical contribution was fainter than the minimum flux set from the observations. The exceptions are FO Aql, IP Peg, SS Aur, and U Gem, for which the theoretical contribution was ≈2 mag brighter (first two systems) or ≈1 mag brighter (latter two), most likely because we assume ${\dot{M}}_{\mathrm{t}}={\dot{M}}_{\text{crit}}^{+}$ when calculating the contribution from the bright spot, whereas ${\dot{M}}_{\mathrm{t}}\ll {\dot{M}}_{\text{crit}}^{+}$; but it is not clear why other systems do not show the same trend.

We conclude that the underestimation of Ṁ_t due to the model assumption and the incomplete light curve are roughly compensated by the overestimation of Ṁ_t due to contributions from the white dwarf, companion star, and bright spot. Hence, we decided not to include any systematic correction to the estimates presented in this work. The statistical error bars also provide a reasonable estimate of the systematic error on Ṁ_t. In all cases, none of the modifications to our analysis that we explored led to results challenging the DIM. The dwarf novae were always consistently placed within the instability region of the parameter space.

5. Conclusions

We have estimated the mass transfer rate Ṁ_t in ≈130 CVs and found that their separation into stable (nova-likes) and unstable (dwarf novae) systems in the (P_orb, Ṁ_t) plane is fully consistent with the DIM (Fig. 3). Although it was clear that nova-likes were going to have higher average accretion rates than dwarf novae, by virtue of their higher steady optical fluxes, it was not clear that the predicted instability line from the DIM would neatly separate the two subpopulations. We thus confirm a key prediction of the DIM that had not been systematically explored yet. We also confirm that irradiation heating of the disk plays a negligible role in the stability of CVs, unlike X-ray binaries where strong X-ray irradiation modifies the limit between stable and unstable systems (Coriat et al. 2012). One of our motivations was to find systems that would challenge the DIM, much like SS Cyg has in the past before the question of its distance was settled (Sect. 1). We have not uncovered such systems. The only system that clearly does not fit the DIM is AE Aqr, which is well known for its fast-spinning white dwarf propelling matter out of the system.

We focused on analyzing a large sample of sources and very likely traded this with some accuracy in Ṁ_t, especially for the binaries whose light curve contains many gaps. A more careful selection and analysis of the individual light curves, corrected for the various contributions to the optical flux and sources of extrinsic variability such as eclipses², together with a critical appraisal of the system parameters would undoubtedly improve the accuracy of some of the individual Ṁ_t that we derived in Tables A.2–A.3. These (time-consuming) refinements would allow us to quantitatively investigate if or how the location of the binary in the (P_orb, Ṁ_t) plane reflects its variability pattern (sub-class type) and the deviations from the expected secular mass transfer rate. They would also allow an additional test of the DIM, albeit dependent on outburst physics, which is that systems closer to the critical mass transfer rate spend longer in outburst than unstable systems (the transientness diagram of Coriat et al. 2012). At this stage, the uncertainties on system parameters and light curves make it difficult to robustly identify classes of variability with locations in the instability diagram.

The lack of knowledge on system parameters and the biases introduced by the irregular sampling of most of the light curves limit our ability to achieve this goal. The Kepler light curves of the dozen or so CVs that have been observed constitute the ideal dataset. We have not used them for lack of information regarding their system parameters. With those in hand, they would provide a subsample of unprecedented quality. In the future, similar regularly sampled light curves may be expected as a by-product of searches for planetary transits, such as ESA’s Plato and synoptic surveys. For instance, the LSST will provide a huge sample of CV light curves, although the baseline cadence (1 visit every three nights) is not ideal. Even if the Gaia catalog provides distances to these CVs, obtaining orbital periods and constraints on the system parameters from follow-up spectroscopy for this avalanche of data represents a daunting challenge for professional and amateur astronomers.

Acknowledgments

We thank the referee, John Thorstensen, for pointing out to us the 3D reddening map constructed by Green et al. 2018 (http://argonaut.skymaps.info), and Jorge Casares for correcting the mass ratio of SY Cnc. GD and JPL thank the hospitality of the KITP in Santa Barbara, where this work was initiated. GD and JPL acknowledge support from Centre National d’Etudes Spatiales (CNES). This research was supported in part by the National Science Centre, Poland, grant 2015/19/B/ST9/01099 and by the National Science Foundation under Grant No. NSF PHY17-48958. This work has made use of the SVO Filter Profile Service (http://svo2.cab.inta-csic.es/theory/fps/) supported from the Spanish MINECO through grant AyA2014-55216, of the SIMBAD database and of the VizieR catalog access tool, operated at CDS, Strasbourg, France. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. We gratefully acknowledge the variable star observations from the AAVSO International Database contributed by observers worldwide. This research was made possible through their patient monitoring of the activity of CVs over decades.

References

Appendix A: Additional tables and figure

Fig. A.1. Same as Fig. 3 except the name of the system is indicated to the right of the estimated Ṁt.

All Tables

All Figures

Fig. 1. Absolute magnitude in U, B, V, R, I bands (blue to red, top panel) and mass accretion rate onto the white dwarf Ṁ(Rin) during a dwarf nova outburst. The model parameters are M = 0.6 M⊙, Ṁt = 1016 g s−1, Rout ≈ 2 × 1010 cm, i = 0°. The dashed line shows the Ṁ(Rin) reconstructed from the V magnitude light curve assuming Ṁ does not depend on radius. The reconstructed mass transfer rate is 9 × 1015 g s−1 after averaging over the outburst cycle.

In the text

	Fig. 2. Relation between the absolute V magnitude and the inner mass accretion rate Ṁin for the model shown in Fig. 1. The disk samples the black curve clockwise during the outburst cycle. The dashed line shows the relationship using the stationary disk approximation.
In the text

Fig. 3. Mass transfer rates of CVs compared to the stability criterion. Systems above the (red) upper solid line are hot and stable. Systems below the lower (blue) line indicate cold, stable disks if the white dwarf magnetic field B ≥ 105 G. The dashed line represents the expected secular mass transfer rate (Knigge et al. 2011). Square symbols indicate Z Cam type dwarf novae; (red) stars indicate nova-likes. Dwarf novae with f ≥ 0.5 are shown in black and those with f < 0.5 are in gray. Figure A.1 includes system names.

In the text

	Fig. 4. Absolute V magnitude of the companion star (full line) and of the accretion-heated white dwarf (dashed line) as a function of orbital period (Knigge et al. 2011). The dotted line indicates the absolute magnitude of the bright spot, calculated for Ṁ = Ṁcrit and a white dwarf mass M = 0.75 M⊙.
In the text

	Fig. A.1. Same as Fig. 3 except the name of the system is indicated to the right of the estimated Ṁt.
In the text