Abstract

The prototype dwarf nova SS Cyg unexpectedly exhibited an anomalous event in its light curve in the early few months of 2021 in which regular dwarf nova-type outbursts stopped, and only small-amplitude fluctuations occurred. Inspired by this event, we have performed numerical simulations of light curves of SS Cyg by varying mass transfer rates and varying viscosity parameters in the cool disk. We have also studied the effect of gas-stream overflows beyond the outer disk edge in the light curve simulations. We have confirmed that the enhanced mass transfer is not likely to be responsible for the 2021 anomalous event or its forerunner. We have found that the enhancement of the viscosity in the disk may reproduce the forerunner of that event but may not be sufficient to explain the 2021 anomalous event, although the latter result might be particular to the thermal equilibrium curve we used. Within our simulations, a model of the gas-stream overflow with a slightly higher mass transfer rate than that of our standard model reproduces light curves similar to the 2021 anomalous event. We suggest that the gas-stream overflow is necessary to reproduce that event. The gas-stream overflow may also be responsible for the abnormally high X-ray flux during the normal quiescent state in SS Cyg.

1 Introduction

Dwarf novae (DNe), one subclass of non-magnetic cataclysmic variables (CVs), are eruptive variables that show sudden brightening called “outbursts” of amplitudes of 2–6 mag and with a typical repetition cycle of a few weeks to several months. CVs are close binary systems composed of the primary white dwarf (WD) and the secondary low-mass star. An accretion disk is formed around the WD by the transferred mass from the secondary star [see Warner (1995) for a review]. The DN outbursts are caused by sudden brightening in the accretion disk, which is, in turn, now thought to be caused by instabilities in the accretion disk. This model is called the disk instability model (DIM), which was first proposed by Osaki (1974). In this model, it is considered that the disk stores mass provided by the secondary star in the quiescent state, and that the accumulated mass accretes on to the WD in the outburst state. The physical mechanism of this instability was found to be the thermal–viscous instability triggered by partial ionization of hydrogen (Hōshi 1979). The disk experiences relaxation oscillations between the hot state with high accretion rates and the cool state with low accretion rates (i.e., the thermal limit cycle). The numerical simulations have been successfully performed along the DIM by several different groups to reproduce the observed light curves of DNe (see, e.g., Cannizzo 1993a; Osaki 1996; Hameury 2020 for reviews).

Non-magnetic CVs above the period gap are classified into three subclasses: SS Cyg-type stars, Z Cam-type stars, and nova-like stars (NLs). SS Cyg stars repeat dwarf-nova outbursts, and NLs experience no outbursts. In the DIM, the critical mass transfer rate (⁠|$\dot{M}_{\rm crit}$|⁠) distinguishes these subclasses. If the mass transfer rate is lower than |$\dot{M}_{\rm crit}$|⁠, the disk becomes thermally unstable, and the system repeats outbursts, which is recognized as a SS Cyg star. If the mass transfer rate is higher than |$\dot{M}_{\rm crit}$|⁠, the disk is in a thermally stable hot state, and the system shows constant high luminosity, which is observed as an NL. Z Cam stars are the intermediate class between SS Cyg stars and NLs, having mass transfer rates close to |$\dot{M}_{\rm crit}$|⁠. They repeat DN outbursts at some time and from time to time show a phenomenon called “standstill”, which is a constant luminosity state, for some indefinite intervals (months to more than 1 yr), as if they are a hybrid class between SS Cyg stars and NLs. There exists another class of CVs, called IW And-type stars (or “anomalous Z Cam-type stars”). The standstill in the ordinary Z Cam stars is terminated by fading to the quiescent state, while the standstill in IW And stars (or, more precisely speaking “quasi-standstill”) is rather terminated by brightening to small-amplitude outbursts. IW And stars exhibit a characteristic light variation; a repetition of a quasi-standstill (i.e., a mid-brightness interval with small-amplitude oscillatory light variations) terminated by brightening (Kato 2019).

The brightest dwarf nova, SS Cyg, has been monitored for more than 100 years by visual observations and optical photometry. This system has repeated DN outbursts with quiescence intervals of ∼1 month for a long time (see the AAVSO historical light curve),1 and has been recognized as the prototype of SS Cyg-type stars. However, SS Cyg unexpectedly began to show an anomalous light curve, i.e., something like “standstill”, from the end of 2021 January as if it was imitating Z Cam stars (see figures 1 and 2). Figure 1 exhibits the recent light curve of SS Cyg over five years beginning from BJD 2458000 while figure 2 shows the expanded one around the anomalous event. When this happened in the early few months of 2021, we suspected it could be something like a standstill in Z Cam stars because the brightness of the star stayed more or less constant and the mean light level increased to ∼10 mag (vsnet-alert 25453). However, it soon turned out that the anomalous event was not a genuine standstill, but the system showed rapid oscillatory light variations with small amplitude in which the maxima and minima seemed to approach each other. We call this event “the 2021 anomalous event in SS Cyg” in this paper, a phrase already used in Kimura et al. (2021) (Event B there). The X-ray luminosity during the anomalous event was ∼3 × 1033 erg s−1, which is about five times higher than that in the normal quiescence, and something unusual happened over wide wavelengths (see Kimura et al. 2021). Such an event had never been observed in SS Cyg in its long history of observations.

Optical light curve of SS Cyg after BJD 2458000. The data are taken from the AAVSO archive 〈http://www.aavso.org/data/download/〉. We use the V-band photometry and the visual observations.
Fig. 1.

Optical light curve of SS Cyg after BJD 2458000. The data are taken from the AAVSO archive 〈http://www.aavso.org/data/download/〉. We use the V-band photometry and the visual observations.

Expanded optical light curve of SS Cyg around the 2021 anomalous event, adopted from figure 2 in Kimura et al. (2021). The data are taken from the AAVSO archive 〈http://www.aavso.org/data/download/〉 and the Variable Star Network (VSNET), and the same as that used in Kimura et al. (2021). The filled circles and open squares denote the V-band and RC-band photometric data, respectively.
Fig. 2.

Expanded optical light curve of SS Cyg around the 2021 anomalous event, adopted from figure 2 in Kimura et al. (2021). The data are taken from the AAVSO archive 〈http://www.aavso.org/data/download/〉 and the Variable Star Network (VSNET), and the same as that used in Kimura et al. (2021). The filled circles and open squares denote the V-band and RC-band photometric data, respectively.

Let us scrutinize the overall light curve of SS Cyg extending nearly five years, shown in figure 1. In its earliest phase, the light curve of the year 2018, in particular, was that of a typical one in SS Cyg in which eight outbursts were observed in one year with the minimum around 12 mag and the maximum reaching 8.5 mag, so that the outburst amplitude was ∼3.5 mag and the long and short outbursts alternated. However, something unusual began to occur in 2019–2020. The optical and X-ray quiescent levels simultaneously and gradually rose, beginning around 2019 August (see also Kimura et al. 2021). One of the authors (M. Kimura) noticed this and requested intensive optical photometry of SS Cyg by the AAVSO team (Waagen 2020) in 2020 September. Although nobody knew what was going on in SS Cyg at that time, this phenomenon in 2019–2020 could be a forerunner of the anomalous event in 2021. As mentioned above, the 2021 anomalous event then started to occur from the end of 2021 January.

However, we can see from figure 1 that the 2021 anomalous event was short-lived (for a couple of months) and the system seemed to return gradually to the normal state. After about 60 days from the beginning of the anomalous event, the outburst amplitude became gradually larger and the outburst frequency became gradually lower. The X-ray luminosity was also decreasing. The main part of the anomalous event seems to be the first 60 d after the long outburst in 2021 (see the horizontal solid-line bar in figure 2) with saw-tooth-like light variations (rapid-rise, slow-decline type) of amplitude 1 mag, and the anomaly of this event was gradually fading after this part, though we cannot determine exactly when this event ended. Currently (in 2022 November), SS Cyg seems to have returned back to the normal state (see the upper panel of figure 1). We can see that the quiescent flux showed wavy modulations on long timescales of ∼100 d (see figure 1). A general tendency is that the brighter the quiescent level is, the lower the outburst amplitude is, and the higher the outburst frequency is. The 2021 anomalous event occurred when the quiescent level was maximized in these undulations.

Motivated by the 2021 anomalous event, we started performing numerical simulations of the light curve of SS Cyg to reproduce this anomalous event. Since we initially thought that the 2021 anomalous event could be a genuine standstill of Z Cam stars, we tried, for instance, to increase the mass transfer rate up to the critical mass transfer rate. However, this anomalous event was not a genuine standstill. Nevertheless, this event in itself is definitely unusual and it will undoubtedly provide a great opportunity to reconsider the nature of standstill and/or standstill-like phenomena in DNe. In this paper, we have extensively performed the light curve simulations of SS Cyg by taking this anomalous event in mind. To do so, we examine the response of the light curve with varying mass transfer rates and the varying viscosity parameter of α. In the previous observational paper by Kimura et al. (2021), they suggested that enhancement of viscosity in the disk could cause the anomalous event and its forerunner. They suggested possible reasons for the enhancement of viscosity in subsection 4.2 of that paper. Another possibility may be the effect of gas-stream overflow because the gas-stream overflow can keep the inner disk in the hot state as discussed by Kimura, Osaki, and Kato (2020a). We, therefore, study the effect of the possible gas stream overflow beyond the outer disk rim. Although the reason for the sudden occurrence of the stream overflow is unknown, the change in the thermal state of the outer disk and vertical oscillations of the disk induced by the collision of the stream with the disk might cause the overflow (Hessman 1999).

Cannizzo (1993b) has already studied the light curve response of SS Cyg to varying input parameters, in particular, on the mass transfer rate and the viscosity parameter, extensively. Since such an anomalous event was not expected at that time, the parameter ranges in his simulations were not wide enough to examine the present anomalous event. Furthermore, the treatments of the outer disk radius, tidal torque, and tidal dissipation are different between ours and those of Cannizzo (1993b); we treat variable outer disk radius in our simulations, while he used the fixed disk radius. The variable disk radius and the treatment of tidal torques and of the tidal dissipation play the essential role in our case.

This paper is structured as follows. In section 2, we describe the method of our numerical simulations and set the standard model parameters for SS Cyg. In section 3, we present the simulation results for varying mass transfer rates, varying viscosity parameters, and varying overflow rates in relation to the 2021 anomalous event in SS Cyg. Finally, we discuss our results in section 4 and summarize our findings in section 5.

2 Method of numerical simulations and standard model parameters

In this paper, we perform numerical simulations for the time evolution of an accretion disk to study outburst light curves in CVs. The method of our numerical simulations is basically the same as that described in Ichikawa and Osaki (1992) and Kimura et al. (2020b). We calculate the time evolution of a geometrically thin and axisymmetric disk (but “non-tilted” in this paper). We adopt the one-zone approximation in the vertical direction, by assuming a geometrically thin disk in hydrostatic equilibrium. We formulate mass, angular momentum, and energy conservation laws for the accretion disk in the radial direction, and solve them by a hybrid method of explicit and implicit integration. We make some minor modifications to the numerical code of Kimura et al. (2020b), the details of which are described in appendix 1. The disk radius is time-varying in our formulation to conserve the total angular momentum of the disk.

In this study, we assume that the tidal truncation radius works as a solid brick wall. When the disk reaches the tidal truncation radius and tries to expand further, the expansion is stopped at the tidal truncation radius by the strong tidal torques and the disk cannot expand further. In our previous paper (Kimura et al. 2020b), the removal of the extra angular momentum was taken into account, but the extra tidal heating was not considered. In this paper, we take into account the extra tidal heating at the outer disk edge due to the tidal truncation as one of the heating functions. Its formulation and effect on the simulations are presented in appendix 2.

The binary parameters used in this paper for SS Cyg are as follows: the orbital period (Porb) is 0.27512973 d, the WD mass (M1) is |$0.94\, M_{\odot }$|⁠, the mass of the secondary star (M2) is |$0.59\, M_{\odot }$|⁠, the binary separation (a) is 1.427 × 1011 cm, the inclination angle (i) is 45°, the tidal truncation radius (rtidal) is 0.356a = 5.08 × 1010 cm, and the the Lubow–Shu radius (or the circularization radius, rLS) is 0.105a = 1.50 × 1010 cm (Paczynski 1977; Hessman et al. 1984; Hill et al. 2017). In calculating the V-band magnitude, we use the distance of SS Cyg, which is 114.25 pc measured by the Gaia parallax (Bailer-Jones et al. 2018). We also add the brightness of the Roche-lobe filling secondary star with the temperature of 4200 K, which corresponds to 12.8 mag, and the brightness of the WD of 50000 K, which corresponds to 15.3 mag (see Hameury et al. 2020). We regard the flux from these two objects as constant for simplicity.

One of the most important factors for simulations of dwarf-nova outbursts is the selection of the viscosity parameter of Shakura and Sunyaev (1973), α:αhot for the hot disk in outburst, and αcool for the cool disk in quiescence, which affects the “S”-shaped thermal equilibrium curve characterizing the outburst behavior. If the ratio αhotcool is larger, the outburst amplitude becomes larger (Smak 1984; Mineshige & Osaki 1985). The viscous timescale is longer for smaller α, so that the outburst duration becomes longer as αhot decreases, and the recurrence time of outbursts increases as αcool decreases. We have performed a parameter search to reproduce the characteristics of a typical observed light curve in SS Cyg by our code. The parameters adopted here for our standard model for SS Cyg are |$\dot{M}_{\rm tr} = 10^{16.9}\:$|g s−1, αhot = 0.15, and αcoolhot = 0.1, in order to reproduce an outburst interval (the quiescence duration) of ∼20–40 d, an outburst amplitude of ∼4 mag, and a duration of outbursts of ∼10–20 d. The radial dependence of αcool is not implemented in this paper because it is not necessary in the case of reproducing normal outbursts in SS Cyg, while we did implement it in Kimura et al. (2020b).

The “S”-shaped thermal equilibrium curve that we use in this paper is calculated on the basis of figure 3 of Kimura et al. (2020b), which is a simple approximate expression for figure 2 of Mineshige and Osaki (1983). The left-hand panel of figure 3 in this paper presents examples of thermal equilibrium curves for four representative annuli of the disk with αhot = 0.15 and αcoolhot = 1.0. We exhibit in the right-hand panel of the same figure the thermal equilibrium curves at r = 2.0 × 1010 cm with various αcoolhot values. As αcoolhot is increased, the cool branch rises upwards in the vertical direction, the critical surface density Σmax above which there is no stable solution in the cool branch decreases, and the unstable intermediate branch becomes shorter.

(Left) Thermal equilibrium curves at r = 109, 109.5, 1010, and 1010.5 cm for our basic parameter sets in section 2. Here, αhot = 0.15 and αcool/αhot = 1.0. (Right) Comparison of thermal equilibrium curves at r = 2.0 × 1010 cm for αhot = 0.15 with different αcool/αhot. These equilibrium curves are calculated by $2F = Q_1^{+}$ [see Kimura et al. (2020b) for the expressions of F and $Q_1^{+}$].
Fig. 3.

(Left) Thermal equilibrium curves at r = 109, 109.5, 1010, and 1010.5 cm for our basic parameter sets in section 2. Here, αhot = 0.15 and αcoolhot = 1.0. (Right) Comparison of thermal equilibrium curves at r = 2.0 × 1010 cm for αhot = 0.15 with different αcoolhot. These equilibrium curves are calculated by |$2F = Q_1^{+}$| [see Kimura et al. (2020b) for the expressions of F and |$Q_1^{+}$|].

We note here that different authors used different thermal equilibrium curves in past studies since different authors calculated the vertical structure of the disk in different ways, such as in the treatment of convective energy transport and opacities used, and so on. Our thermal equilibrium curve is based on the computations in Mineshige and Osaki (1983) in which the vertical structure was solved, and the convective energy transport was taken into account by the mixing length theory, but molecular opacities were not considered. As seen in figures 3 and 14 of Mineshige and Osaki (1983), the thermal equilibrium curve takes a ξ-shaped curve and has two unstable branches (i.e., two local maxima in Σ in the cold branch). As discussed by Cannizzo (1993c), the lower maximum is determined by the opacity, while the upper maximum is determined by convection [see figure 2 of Pojmanski (1986)]. The end of the cool branch in our thermal equilibrium curve is the lower maximum of the two. As seen in figure 14 of Mineshige and Osaki (1983), the lower maximum in Σ is dominant when the convective energy transport is inefficient, while the upper maximum is dominant when the convective energy transport is very efficient. For example, the thermal equilibrium curve used in Hameury et al. (1998) corresponds to the case of very effective convection (see figure 1 of that paper). In their thermal equilibrium curve, the maximum effective temperature in the cold branch is independent of αcoolhot, which is higher than ours. The main source of the opacity in the cool branch of our thermal equilibrium curve is negative hydrogen, and Σmax in the cool branch of our thermal equilibrium curve corresponds to the point at which the disk changes from the optically thin state to the optically thick state, and there is still controversy whether the disk in the cold state is optically thin or thick. Molecular opacities were not considered in our thermal equilibrium curve, a certain drawback in our calculations. However, molecular opacity is important if the temperature is lower than ∼2500 K [see figure 2 of Cannizzo & Wheeler (1984)]. Our thermal equilibrium curve thus tends to have the maximum effective temperature of the cool disk lower than those in other works such as Hameury et al. (1998).

The inner edge of the disk, rin, is also one of the important parameters. In the quiescent state, it is thought that the inner disk is truncated via evaporation of the disk mass by the coronal siphon flow from theoretical and observational studies (Meyer & Meyer-Hofmeister 1994; Balman & Revnivtsev 2012). In this case, if rin becomes larger, inside-out outbursts are suppressed, and the X-ray flux in quiescence increases (Cannizzo 1993b). On the other hand, when the outburst occurs, the inner hole is thought to be refilled by the accretion flow and the inner edge of the disk is thought to extend down to the WD surface during outburst. In our standard model, we adopt rin to be 5 × 108 cm, which is approximately the WD radius for a WD mass of |$0.94\, M_{\odot }$| (Nauenberg 1972; Provencal et al. 1998).

We divide the accretion disk between rin and the outer disk edge into N concentric annuli for our finite difference scheme, where N is variable in time because the disk radius is variable in our formulations. The radial distribution of meshes is equally spaced in |$\sqrt{r}$| and the maximum number of concentric annuli, denoted here by N0, is 200, which occurs when the disk outer edge reaches the tidal truncation radius.

The results of our numerical simulations for normal outbursts in SS Cyg by using the parameter sets as described above are displayed in figure 4. The V-band magnitude is estimated by the method in Dubus, Otulakowska-Hypka, and Lasota (Dubus et al. 2018) under the assumption that the radiation from the disk is multi-temperature blackbody emission. The flux from the bright spot is included in the same manner as described in Kimura et al. (2020b). We see that the recurrence time, the outburst interval, the outburst duration, and the repetition of a short outburst and a long outburst are consistent with the observations of SS Cyg (Cannizzo & Mattei 1992). The cycle length of long outbursts is ∼100 d in observations, and this is sensitive to the number of short outbursts between two long outbursts. In the case of |$\dot{M}_{\rm tr} = 10^{16.9}\:$|g s−1, the average cycle length of long outbursts in figure 4 is ∼100 d, which is consistent with observations. In quiescence, radiation from the secondary star and that from the hot spot dominates over that of the disk, in our case, just in the same way as the other simulations by Cannizzo (1993b) and Hameury et al. (2020). We note here that the shoulder-like feature (or the precursor) in the long outburst is barely visible in the V-band light curve. This feature is produced when the disk’s outer radius reaches the tidal truncation radius, and the tidal torques and tidal dissipation are greatly increased. We regard this model as the standard model in our simulations. We perform numerical simulations by changing various parameters from the standard model in the following sections.

Results of our simulations for $\dot{M}_{\rm tr} = 10^{16.9}\:$g s−1, αhot = 0.15, αcool/αhot = 0.1, and rin = 5 × 108 cm. From top to bottom: the light curve for the bolometric luminosity of the system, that for the V-band magnitude, the disk radius in units of the binary separation, the total disk mass, and the total disk angular momentum. The blue line in the top panel approximately represents the expected X-ray flux in quiescence where rin = rWD. The green and orange lines in the second panel represent the V-band magnitude of the disk and the bright spot, respectively. The black line in the same panel stands for a sum of the flux from these two components plus the WD and the secondary star.
Fig. 4.

Results of our simulations for |$\dot{M}_{\rm tr} = 10^{16.9}\:$|g s−1, αhot = 0.15, αcoolhot = 0.1, and rin = 5 × 108 cm. From top to bottom: the light curve for the bolometric luminosity of the system, that for the V-band magnitude, the disk radius in units of the binary separation, the total disk mass, and the total disk angular momentum. The blue line in the top panel approximately represents the expected X-ray flux in quiescence where rin = rWD. The green and orange lines in the second panel represent the V-band magnitude of the disk and the bright spot, respectively. The black line in the same panel stands for a sum of the flux from these two components plus the WD and the secondary star.

We have checked if the numerical resolution is enough in the case of N0 = 200 by performing simulations with coarser or finer grid numbers. We display in figure 5 the V-band light curves and the time evolution of the disk mass of our simulations with N0 = 100, 200, 400, and 600. We see that the number of short outbursts between two long outbursts is 1 in the case of a coarse grid with N0 = 100 (the top part of each panel), but it becomes 2 in the three simulations for finer grids with N0 ≥ 200. The light curve and the time evolution of the disk mass are almost the same and indistinguishable among these three simulations (see the lower three parts of each panel). This result signifies that a resolution with N0 ≥ 200 is required but the grid number N0 = 200 is enough to investigate simulation results.

V-band light curves of the system (left) and time evolution of the disk mass (right) at various resolutions. The green line in the left-hand panel represents the disk brightness.
Fig. 5.

V-band light curves of the system (left) and time evolution of the disk mass (right) at various resolutions. The green line in the left-hand panel represents the disk brightness.

The X-ray luminosity in quiescence is estimated based on the evaporation model. In this model, the geometrically-thin disk is truncated at some inner radius denoted by revap, below which the disk expands to the hot spherical corona and the central WD is surrounded by the optically thin hot corona with a temperature of ∼107 K. The X-ray luminosity is then given by
(1)
where |$\dot{M}_{\rm acc}$| is the mass supply rate and rWD is the WD radius. The mass supply rate to the hot corona is estimated by the mass accretion rate at revap. In this paper, we assume that revap is equal to rin. Here, we assume that the WD spin is much slower than the Keplerian velocity and that the accreted gas does not convey any energy inside the WD by dissipating all energy in the hot corona. Since the local accretion rate at a given radius, r, is expressed as |$\dot{M}_{\rm acc}(r) = 2 \pi r \Sigma (-v_r)$|⁠, where Σ is the surface density, vr is the radial velocity by the viscous diffusion, and (−vr) ∼ αcool(h/r)cs, the local accretion rate in quiescence at a given radius, |$\dot{M}_{\rm acc}(r)$|⁠, is estimated roughly proportional to αcoolTcr2.5. The X-ray luminosity is highly dependent on the inner truncation radius, revap = rin, and αcool.

We can estimate X-ray luminosity by using simulated mass accretion rates at the inner disk edge [|$\dot{M}_{\rm acc}(r_{\rm in})$|]. The blue line in the top panel of figure 4 approximately represents the X-ray luminosity in quiescence where rin = rWD. The estimated value is ∼4 × 1029 erg s−1 on average, which is about three orders of magnitude lower than the observed X-ray luminosity ∼6 × 1032 erg s−1 (Wheatley et al. 2003; Ishida et al. 2009). This discrepancy between the observed and theoretically predicted X-ray luminosity was already pointed out by Wheatley, Mauche, and Mattei (2003). To increase the X-ray luminosity in the optical quiescence in model calculations, one might increase either αcool or rin, or both. Although rin is constant in our simulations, the inner disk would be truncated during quiescence in SS Cyg. For instance, Balman and Revnivtsev (2012) suggested that rin is ∼5 × 109 cm by their observations. We have tried simulations with rin = 5 × 109 cm and estimated the X-ray luminosity to be ∼5 × 1031 erg s−1 during quiescence. This is, however, still less than a tenth of the observational value.

3 Simulations for the 2021 anomalous event in SS Cyg

3.1 Enhanced mass transfer model

We here explore the possibility of enhanced mass transfer as a cause for the 2021 anomalous event in SS Cyg. As mentioned in the introduction, we originally thought that this anomalous event could be a kind of the standstill phenomenon seen in Z Cam stars. It is widely believed that the Z Cam-type standstill is produced by fluctuations in mass transfer rates, where the mass transfer rate is enhanced above the critical one above which the accretion disk is in the hot stable state (Meyer & Meyer-Hofmeister 1983; Lin et al. 1985; Buat-Ménard et al. 2001).

We first perform simulations by gradually increasing the mass transfer rate from that of the standard model in order to see what happens in such a case. Figure 6 illustrates our result of such simulations. In figure 6, the mass transfer rate is increased after 200 d with a function of |$\dot{M}_{\rm tr} = 10^{16.9} [1 + 0.0025 (t - 200)]$|⁠, where t is the time in units of days. We see from figure 6 that the cycle length of outbursts decreases with the increase in mass transfer rates until ∼800 d, but it begins to increase after that, mainly because the duration of the long outburst becomes longer. Eventually, the system brightness does not drop anymore and the disk enters the hot steady state. It is interesting to note that the system approaches to the steady state by increasing outburst interval but not decreasing outburst amplitude. In fact, the outburst amplitude for the disk component remains large even in an outburst just before entering the steady state (see the green line in figure 6). We see that |$\dot{M}_{\rm crit}$| is around 1017.4 g s−1. The V-band magnitude in quiescence gradually increases with the increase in mass transfer rates, which is mainly due to the increasing luminosity of the bright spot.

Results of our simulations in which the mass transfer rate is gradually increased with time. The upper and lower panels represent the mass-transfer rate and the V-band magnitude light curve, respectively. See the caption of figure 4 for descriptions of the black and green lines in the lower panel.
Fig. 6.

Results of our simulations in which the mass transfer rate is gradually increased with time. The upper and lower panels represent the mass-transfer rate and the V-band magnitude light curve, respectively. See the caption of figure 4 for descriptions of the black and green lines in the lower panel.

Figure 7 shows the relation between the cycle length of long outbursts and |$\dot{M}_{\rm tr}$|⁠. Here the cycle length is defined by the interval between two consecutive long outbursts because the system repeats one cycle with this timescale for a given |$\dot{M}_{\rm tr}$| (see figure 4). We see that the system approaches a steady state as the cycle length of long outbursts increases in our simulations. The duration of long outbursts becomes infinity if the mass transfer rate approaches 1017.4 g s−1. We note here that the cycle length abruptly increases and approaches very rapidly to infinity in figure 7. The transition from the dwarf nova outbursting stage to the nova-like stage with a steady disk occurs within a very narrow range in the mass transfer rate: a situation for a favorable condition for the Z Cam phenomenon.

Variations of the cycle length of long outbursts against $\dot{M}_{\rm tr}$.
Fig. 7.

Variations of the cycle length of long outbursts against |$\dot{M}_{\rm tr}$|⁠.

We may note that similar behavior was also found between the supercycle length and the mass transfer rate in the thermal–tidal instability model for the SU UMa-type stars by Osaki (1995), who tried to explain ER UMa stars: DNe undergoing frequent superoutbursts (Kato & Kunjaya 1995; Robertson et al. 1995). In this model, when the mass transfer rate is increased in numerical simulations, the recurrence time decreases at first and reaches a minimum value at a certain point but it begins to increase after that because the duration of superoutbursts increases. The cycle length eventually becomes infinity, and the system becomes a nova-like star.

By taking into account the above-mentioned results, we next try to simulate the 2021 anomalous event in SS Cyg based on the enhanced mass transfer burst model. We perform two simulations by switching the mass transfer rate from the standard value of |$\dot{M}_{\rm tr}=10^{16.9}\:$|g s−1 to some enhanced ones after 240 d. The two enhanced mass transfer rates adopted here are 1017.4 g s−1 and 1017.3 g s−1. We mark these mass transfer rates by the numbers 1 and 2 in figure 6.

The results are shown in figure 8 for the V-band light curves of these two simulations. We see that the system enters into an (almost) constant-luminosity state, such as a standstill in the case of the highest mass transfer shown in the upper panel in figure 8. However, the simulated light curve did not look like the observed 2021 anomalous event of SS Cyg, as the simulated constant light level was as high as the maximum of the short outburst before the switching in mass transfer rates and it was much higher than the observed level averaged in the 2021 anomalous event of SS Cyg. In the second case, shown in the lower panel of figure 8, the outbursts still occur after enhanced mass transfer. The outburst amplitude in simulations was, however, ∼3 mag, which is much higher than the observed one (∼1 mag) for the 2021 anomalous event in SS Cyg. Furthermore, long outbursts inevitably occur in such a case because of enhanced mass transfer. The simulated light curves do not match well with the observed one for the 2021 anomalous event in SS Cyg. We thus conclude that the enhanced mass transfer is very unlikely to be a cause for this event.

Time evolution of the V-band magnitude of the system for the enhanced mass transfer burst model. In the upper and lower panels, the transfer rates are raised up to 1017.4 g s−1 and 1017.3 g s−1 after 240 d, respectively.
Fig. 8.

Time evolution of the V-band magnitude of the system for the enhanced mass transfer burst model. In the upper and lower panels, the transfer rates are raised up to 1017.4 g s−1 and 1017.3 g s−1 after 240 d, respectively.

3.2 Enhanced viscosity model

In order to reproduce the clear-cut outburst and quiescence in DNe by the thermal viscous instability model, it is known that we need to choose αhot in the hot state higher by a factor of 3–10 than αcool for the cool state (Meyer 1984; Osaki 1989). It is widely accepted that the viscosity in the hot ionized disk is produced by the magneto-rotational instability (MRI; Balbus & Hawley 1991). On the other hand, the origin of viscosity in the cold disk is not known yet. In the cold disk, the gas is mostly neutral, and the magnetic field may decouple from the disk matter, and the MRI turbulence could not be maintained in such a case. Under such a circumstance, several different possible sources for the viscosity in the cold disk are discussed, but so far, no model is widely accepted. It is a common practice that light curve simulations are performed by choosing the viscosity parameter αcool for the cold disk in such a way to reproduce observed light curves of DNe.

We here investigate the possibility of enhanced viscosity in the cold disk for the 2021 anomalous event in SS Cyg. We first perform simulations by gradually increasing the ratio of viscosity parameter αcoolhot from the standard value of 0.1 to 1.0 to see what happens. Figure 9 shows the result of our simulations, where αcool is gradually increased during 900 d from the standard value from date 200 d to αhot. We see in figure 9 that slow-rise outbursts (inside-out outbursts) occur more frequently as αcool increases. Small reflares emerge on the fading tail of outbursts. We also find that the flat-bottomed quiescence disappears. In our simulations, we have chosen the inner edge of the disk near the WD surface. However, the inner disk in quiescence is, in reality, most likely truncated by the evaporation of the disk gas via coronal siphon flow (Meyer & Meyer-Hofmeister 1994) to produce the hot optically thin corona near the WD, which produces X-ray radiation. The small-amplitude inside-out outbursts are mostly likely suppressed, and the quiescent duration may be longer in such a situation. Also, small reflares on the fading tail of outbursts may disappear by taking into account the variation of rin (Menou et al. 2000).

Results of our simulations in which αcool is gradually increased with time and finally becomes comparable with αhot. (Upper) Ratio of viscosity parameters αcool/αhot. (Lower) Light curve for the V-band magnitude. See the caption of figure 4 for a description of the black and green lines in the lower panel.
Fig. 9.

Results of our simulations in which αcool is gradually increased with time and finally becomes comparable with αhot. (Upper) Ratio of viscosity parameters αcoolhot. (Lower) Light curve for the V-band magnitude. See the caption of figure 4 for a description of the black and green lines in the lower panel.

As the viscosity in the cold disk is increased further in figure 9, the outburst cycle becomes shorter, and the quiescent brightness level increases. Also, the difference in the amplitude and duration between long and short outbursts become smaller. This behavior looks similar in some sense to the observed light curve before the 2021 anomalous event in SS Cyg. As discussed in Kimura et al. (2021), low-amplitude and slow-rise outbursts seem to have frequently occurred after around BJD 2458700, the long outbursts were suppressed, and the cycle length of outbursts became shorter (see figure 1).

We understand the effect on the light curve due to enhanced viscosity from the change in the “S”-shaped thermal equilibrium curve. If the viscosity in the cool state is enhanced, the unstable branch in the thermal equilibrium curve becomes shorter (see the right-hand panel of figure 3). Then the upper limit of the effective temperature in the quiescent state increases. This change also prevents cooling and heating waves from propagating over the entire disk and causes the thermal instability frequently in local regions. The innermost region of the disk stays for a longer time in the hot state, and the thermal instability frequently occurs in other regions. Besides, more mass is conveyed inwards, and the thermal instability is easily triggered in regions other than the outermost disk. This is why we confirm small-amplitude and slow-rise outbursts and small reflares in the fading tail in the simulated light curves.

Let us now try to simulate the 2021 anomalous event in SS Cyg by artificially increasing the value of αcoolhot. Figure 10 illustrates the results of our simulations in which αcoolhot is increased from 0.1 to 1.0 instantaneously on 240 d. We mark this switching ratio as the number 1 in figure 9. We see in figure 10 that the minimum level in the simulated light curve increases around 0.5 mag after the switching in this ratio. The saw-tooth-like oscillatory light curves (or outbursts with rapid-rise and slow-decay ones) result with an amplitude of 2 mag. Qualitatively speaking, the model light curve after the switching event in figure 9 looks like the 2021 long outburst and its subsequent anomalous phenomenon in SS Cyg: three saw-tooth-like outbursts have occurred within about 60 days of the local minimum at 275 d. However, quantitatively speaking, the oscillatory amplitude of about 2 mag in our simulations is larger than the observed one with about 1 mag in the 2021 anomalous event, and the increment in the minimum brightness in our simulations is smaller than the observed one of about 1.5 mag.

Results of our simulations for the enhanced viscosity model. Here, αcool/αhot jumps from 0.1 to 1.0 on 240 d. From top to bottom: the ratio of viscosity parameters αcool/αhot, the light curve for the V-band magnitude, and the disk radius in units of the binary separation.
Fig. 10.

Results of our simulations for the enhanced viscosity model. Here, αcoolhot jumps from 0.1 to 1.0 on 240 d. From top to bottom: the ratio of viscosity parameters αcoolhot, the light curve for the V-band magnitude, and the disk radius in units of the binary separation.

We have tried two more simulations basically in the same way as that mentioned above but with different increments in αcool, and we compare them in figure 11. In the middle panel, αcoolhot is increased to 0.8, while it is increased to 0.6 in the bottom panel. The final αcoolhot values in these two simulations are denoted by the numbers 2 and 3 in figure 9. We see that the larger the increment is, the smaller the amplitude of outbursts after the long outburst starting from ∼240 d is, and the shorter the recurrence time of outbursts is. We have tried the same simulations with different rin and α values and confirm that the main features, as mentioned above, are always reproduced.

Time evolution of the V-band magnitude of the system for the enhanced viscosity model. Here, αcool is enhanced on 240 d. The top panel is the same as the middle panel of figure 10. In the middle panel, αcool/αhot is raised up to 0.8. In the bottom panel, αcool/αhot is raised up to 0.6.
Fig. 11.

Time evolution of the V-band magnitude of the system for the enhanced viscosity model. Here, αcool is enhanced on 240 d. The top panel is the same as the middle panel of figure 10. In the middle panel, αcoolhot is raised up to 0.8. In the bottom panel, αcoolhot is raised up to 0.6.

Before we tried our simulations, we had anticipated that the simulated light curves would result in small-amplitude fluctuations in the case of constant α for the hot and cool states. However, our results shown in figure 10 demonstrate that fairly clear outbursts occur rather than small-amplitude fluctuations even in the case of αcoolhot = 1.0 (in particular, see the green line in figure 10 for the disk component only) . For instance, Mineshige and Osaki (1985) showed that the simulations with αcoolhot = 1.0 result in small-amplitude fluctuations. Also, Meyer (1984) and Mineshige and Shields (1990) explored the conditions for clear-cut outbursts in which the heating and cooling fronts propagate all the way to the disk from the outer edge to the inner edge. One such condition is that Σmaxmin is larger than 2, where Σmin is the critical surface density at the downward transition from the hot state. The reason is that if Σmaxmin is less than 2, the heating and cooling fronts cannot travel a long distance as they are reflected within a short distance. As seen in the S-shaped thermal-equilibrium curves shown in the left-hand panel of figure 3, Σmaxmin is slightly less than 2 in our case of αcoolhot = 1.0, and thus our cases do not satisfy this condition.

To see what happens in our simulations in the case of constant α shown in figure 10, we show figures for the time evolution of the transition front. Figure 12 illustrates three-dimensional plots of the propagation of a cooling wave in a part of simulation results with αcoolhot = 1.0, which are displayed in figure 10. The three-dimensional plots for the temperature distribution in the lower panels show a part of the light curves indicated by the blue line in the upper panels, and the right-hand figures correspond to a locally expanded version of the left-hand panel with further details. We see from the left-hand panel of figure 12 that the cooling front propagates inward and reaches to the deep inner disk, to r < 1010 cm. However, its propagation is not smooth, but it goes locally back and forth, which is seen as a small ripple-like variation in the decay light curve. The right-hand panel of figure 12 is its locally expanded version, and it shows this feature very clearly. However, the cooling front never reaches the innermost region of the disk, but it is reflected in the middle of the disk as a heating front. The innermost part of the disk thus remains always in the hot state in our simulation.

Three-dimensional figures of cooling front propagation (lower panels) and corresponding light curves (upper panels) of the simulation displayed in figure 10. (Left) Time evolution of the disk temperature from 278.5 to 291 d per 0.5 d. The corresponding time period in the lower panel is represented by the blue line in the upper light curve. (Right) Detailed version of a part of cooling front propagation. The time evolution of the disk temperature from 284.4 to 289 d per 0.1 d and the corresponding light curve.
Fig. 12.

Three-dimensional figures of cooling front propagation (lower panels) and corresponding light curves (upper panels) of the simulation displayed in figure 10. (Left) Time evolution of the disk temperature from 278.5 to 291 d per 0.5 d. The corresponding time period in the lower panel is represented by the blue line in the upper light curve. (Right) Detailed version of a part of cooling front propagation. The time evolution of the disk temperature from 284.4 to 289 d per 0.1 d and the corresponding light curve.

The reason why our simulation results in a clear outburst with amplitudes as large as 2 mag in the case of constant α is not certain, but the results in the case of constant α may be very sensitive to the detailed profile of the S-shaped thermal-equilibrium curve used in the simulations. The value of Σmaxmin in our case is less than 2 but very close to 2. Besides that, the maximum effective temperature of the cool branch is as low as 4000 K in the outer disk in our case (see the left-hand panel of figure 3). The separation of the effective temperature between the hot branch and the cool branch in the thermal equilibrium curve in our case is thus wider than those in other works, as described in section 2. These two effects combined could be the reason for this. We here note that a clear outburst found in our case is not a genuine full-scale outburst because the innermost part of the disk always remains in the hot state. The problem is then how far (how deeply) the cooling front propagates in the disk, and this may depend on the thermal equilibrium curve used, so the amplitude of oscillatory light variations may be different with a different thermal-equilibrium curve. This cannot, however, be easily tested because different authors use different thermal-equilibrium curves.

Figure 13 illustrates three-dimensional plots of the propagation of a heating wave, a corresponding one to figure 12 for a cooling wave. As in figure 12, the three-dimensional plots for the temperature distribution in the lower panels correspond to a part of the light curves indicated by the red line in the upper panels, and the right-hand figures are a locally expanded version of the left-hand figures . We see from figure 13 that the heating wave, which is generated by the reflection of the cooling wave in the middle of the disk, propagates outward (i.e., the inside-out type) and its propagation is rather smooth in this case. Its right-hand panel is a locally expanded version of the left-hand panel and it clearly shows no ripple-like structure, in contrast with the case of the cooling front. It was found that the heating wave propagated faster than the cooling wave, in this case, resulting in a rapid-rise and a slow-decay light curve. It is generally believed that the rise and decay of light curves are symmetric in the inside-out type outbursts. The reason why we have obtained a rapid-rise and slow-decay type light curve in our case may probably be because the heating front propagated smoothly while the cooling front did not.

Three-dimensional figures of heating front propagation (lower panels) and corresponding light curves (upper panels) of the simulation displayed in figure 10. (Left) Time evolution of the disk temperature from 290.9 to 293.5 >d per 0.1 d. The corresponding time period in the lower panel is represented by the red line in the upper light curve. (Right) Detailed version of a part of heating front propagation. The time evolution of the disk temperature from 290.7 to 292.2 d per 0.025 d and the corresponding light curve.
Fig. 13.

Three-dimensional figures of heating front propagation (lower panels) and corresponding light curves (upper panels) of the simulation displayed in figure 10. (Left) Time evolution of the disk temperature from 290.9 to 293.5 >d per 0.1 d. The corresponding time period in the lower panel is represented by the red line in the upper light curve. (Right) Detailed version of a part of heating front propagation. The time evolution of the disk temperature from 290.7 to 292.2 d per 0.025 d and the corresponding light curve.

In the previous subsection, we have demonstrated in the case of αcoolhot = 0.1 that the system approaches a hot steady state by increasing the duration of the long outburst (see figures 6 and 7), but not by decreasing the outburst amplitude, when we increase the mass transfer rate from that of the standard case up to the critical one. As described above, the outburst amplitude becomes smaller when the viscosity in the cool state is enhanced. We wonder whether or not there is another possibility (i.e., the second possibility) that the system might approach to a steady state by decreasing the amplitude of outburst and finally with vanishing outburst amplitude in the case of αcoolhot ≃ 1.0, when we increase mass transfer rate with time. Figure 14 illustrates our simulations in the case of αcoolhot = 1.0 in which the mass transfer rate is increased with time. We see from this figure that the amplitude of outbursts becomes smaller with the increase in mass transfer rates but never vanishingly small when the system approaches to a hot steady state. By comparing this case with the case where αcoolhot = 0.1, the outburst amplitude is smaller by ∼1.0 mag just before the system enters into a constant luminosity state in this case. Besides, |$\dot{M}_{\rm crit}$| is lower in this case since the disk easily becomes thermally stable with higher ratios of αcoolhot. However, the approach to a hot steady state is basically the same between the two cases, that is, the duration of long outburst becomes infinity. Other simulations we have performed give basically the same result. The second possibility seems to be ruled out in our formulation, where the tidal truncation and variable disk radius are taken into account.

Results of our simulations in which the mass transfer rate gradually increases with time in the case of αcool/αhot = 1.0. See the caption of figure 4 for a description of the black and green lines in the lower panel.
Fig. 14.

Results of our simulations in which the mass transfer rate gradually increases with time in the case of αcoolhot = 1.0. See the caption of figure 4 for a description of the black and green lines in the lower panel.

3.3 Possibility of gas-stream overflow

The most important observational feature of the 2021 anomalous event in SS Cyg is an elevation of brightness level in quiescence, that is, from 12 mag in the ordinary quiescence to ∼10–11 mag during the anomalous event. If the entire disk stays in the cool state with low viscosity in quiescence, such an anomalous event is obviously not possible. We need to have a hot state at least in a part of the disk, in particular in the inner part of the disk. In this respect, gas-stream overflow could be another possible solution to this problem. In fact, Kimura et al. (2020b) performed light curve simulations in which the effect of gas-stream overflow plays an important role in the case of tilted disks in relation with the IW And-type phenomenon. Even in the case of a non-tilted disk, a part of the gas stream from the secondary star can overflow above and below the disk edge to penetrate into the inner part of the disk, as discussed first by Lubow and Shu (1976). Hessman (1999) showed that the gas stream could overflow the disk edge, in particular in the case of a cold disk in quiescence. Kunze, Speith, and Hessmann (2001) performed the smoothed particle hydrodynamic simulations on the stream overflow, showing that a significant fraction of the gas stream can overflow the disk edge and penetrate into the inner disk.

As for the light curve simulations for DN outbursts based on the thermal–viscous instability model, Schreiber and Hessman (1998) studied the effects of the stream overflow upon the outburst light curves. However, their study was rather limited for a special case. In this study, we perform light curve simulations to explore the effects of the stream overflow for much wider ranges in parameter space.

To treat the gas-stream overflow, we change the mass supply pattern to the disk. If the gas-stream overflows the outer disk rim, a part of the stream enters the outer disk edge, and the rest of it flows into the inner disk. In this study, we assume that the overflowing gas entering into the inner disk basically goes to a region around the Lubow–Shu radius (i.e., the circularization radius). There are two reasons why we adopt this assumption. First, the hydrodynamic simulations by Kunze, Speith, and Hessmann (2001) demonstrate that the overflowing gas stream settles at a region around the circularization radius, thus justifying our assumption. The second reason concerns a problem of numerical difficulty which Kimura et al. (2020b) encountered in the case when the gas stream reaches the region of r < rLS in the cold disk. Its details were discussed in subsection 6.3 of Kimura et al. (2020b). In short, the disk matter inside of rLS, to which the overflowing gas is added, tends to move towards rLS in the cold disk, since the transferred gas has the specific angular momentum of |$\sqrt{G M_1 r_{\rm LS}}$|⁠. This causes thinning of the disk matter there, which results in a stoppage of calculations. To avoid this numerical difficulty, we assume in our formulation that the mass of overflowing gas is exclusively supplied to a region around rLS as a delta-function-like form. We therefore set the mass supply rate s(r) as displayed in figure 15. Here, |$\dot{M}_{\rm edge}$| and |$\dot{M}_{\rm in}$| are the mass transferred to the outer disk rim and that transferred to the region near rLS, respectively. We define the ratio of |$\dot{M}_{\rm in}$| to the total mass transfer rate as foverflow, and |$\dot{M}_{\rm edge} = (1 - f_{\rm overflow}) \dot{M}_{\rm tr}$| and |$\dot{M}_{\rm in} = f_{\rm overflow} \dot{M}_{\rm tr}$|⁠. In this formulation, s(r) is expressed as follows:
(2)
 
(3)
where ra is the border between the two meshes which is closest to rLS − 0.01a, rb is that closest to rLS + 0.01a, and rc is that closest to rN − 0.02a.
Mass supply pattern for the gas-stream overflow. The mass transfer rate is divided in $\dot{M}_{\rm in}$ and $\dot{M}_{\rm edge}$ in this case.
Fig. 15.

Mass supply pattern for the gas-stream overflow. The mass transfer rate is divided in |$\dot{M}_{\rm in}$| and |$\dot{M}_{\rm edge}$| in this case.

We have carried out simulations in the case of the gas-stream overflow in which we gradually increase foverflow with time while the other parameters are kept the same as those of the standard model. Figure 16 illustrates our simulation results. Here, we calculate the V-band magnitude of the bright spot by assuming that its luminosity is a sum of |$0.25\,\, G M_1 \dot{M}_{\rm edge} / r_N$| and |$0.25\,\, G M_1 \dot{M}_{\rm in} / r_{\rm LS}$|⁠. The size of each spot is postulated to be |$2\%$| of the disk size as in Kimura et al. (2020b). The temperature of the bright spot is determined by TBS = [LBS/(σδS)]1/4, where LBS is the luminosity of the bright spot, σ is the Stefan–Boltzmann constant, and δS is the size of each bright spot. Although it depends on |$\dot{M}_{\rm edge}$| and |$\dot{M}_{\rm in}$|⁠, it is typically ∼10000 K. As foverflow increases, the cycle length of long outbursts becomes longer and a new type of outburst with a small amplitude of about 1 mag begins to appear beside the ordinary short and long outbursts, its number increasing with the increase in foverflow. All of the outbursts become the inside-out type, and the quiescent level increases with the increase in foverflow because of an increase in brightness of the inner bright spot. When foverflow becomes greater than 0.8, the inner part of the disk stays in the hot state most of the time. In the special case of foverflow = 1.0, all of the mass transferred from the secondary is supplied to the inner disk and no mass is supplied to the outer disk edge, i.e., |$\dot{M}_{\rm edge}$| is zero. In such a case, the outer disk edge expands and reaches the tidal truncation radius even in the cool state and the mass slowly accumulates there. Eventually, the heating wave produced by the inside-out outburst propagates outward and reaches the outer edge of the disk, i.e., the tidal truncation radius, producing a long and large outburst, as seen as the last outburst in figure 16. The outer edge of the disk remains at the tidal truncation radius even when the outburst has ended and the outer part of the disk returns to the cool state in this case.

Our results of simulations in the case of gas-stream overflow in which foverflow is gradually increased with time. The lower panel shows the V-band light curve and the upper panel shows the corresponding foverflow value. See the caption of figure 4 for a description of the black and green lines in the lower panel.
Fig. 16.

Our results of simulations in the case of gas-stream overflow in which foverflow is gradually increased with time. The lower panel shows the V-band light curve and the upper panel shows the corresponding foverflow value. See the caption of figure 4 for a description of the black and green lines in the lower panel.

Let us now examine if the effect of gas-stream overflow could explain the 2021 anomalous event of SS Cyg. We have performed simulations by artificially increasing foverflow from 0.0 to 0.5 (i.e., half of the mass is supplied to the inner disk and the other half to the outer edge) on 240 d, and the results are displayed in figure 17. We mark this foverflow value as the number 2 in figure 16. We see from figure 17 that the cycle length of outbursts becomes longer after switching foverflow because the mass transferred to the outer disk edge becomes lower. On the other hand, inside-out-type low-amplitude outbursts frequently occur since a large amount of mass is directly transferred to the inner disk, and the total duration of the quiescent state thus becomes shorter. We have done two more simulations by increasing foverflow to other values on 240 d. Figure 18 illustrates three light curves and the middle panel is the same as that of figure 17. We see from the top panel of figure 18 that inside-out outbursts do not occur in this case with low foverflow but the decaying time from the outburst maximum becomes longer as compared to that before the switching of the foverflow value. On the other hand, if foverflow is high, the quiescent duration becomes shorter because inside-out outbursts frequently occur, so that luminosity dips appear in the light curve (see the bottom panel of figure 18). None of the three light curves shown in figure 18) look like the 2021 anomalous event of SS Cyg, because the minimum level of light curves in all three simulations remained near 12 mag; we must conclude that the effect of gas-stream overflow is not enough to explain this particular event.

Results of our simulations in the case where foverflow is amplified on 240 d. From top to bottom: variations in foverflow, the light curve for the V-band magnitude, and the disk radius in units of the binary separation. Here, foverflow is 0.5 after 240 d.
Fig. 17.

Results of our simulations in the case where foverflow is amplified on 240 d. From top to bottom: variations in foverflow, the light curve for the V-band magnitude, and the disk radius in units of the binary separation. Here, foverflow is 0.5 after 240 d.

Time evolution of the V-band magnitude of the system for the overflow model. The middle panel is the same as that in the second panel of figure 17. In the top and bottom panels, foverflow is raised up to 0.2 and 0.8 after 240 d, respectively.
Fig. 18.

Time evolution of the V-band magnitude of the system for the overflow model. The middle panel is the same as that in the second panel of figure 17. In the top and bottom panels, foverflow is raised up to 0.2 and 0.8 after 240 d, respectively.

When we performed our simulations with the effect of gas-stream overflow, we came to realize the similarity between our calculations and those by Kimura et al. (2020b) in which the gas-stream overflow occurred because of a tilted disk [see, e.g., figure 17 of Kimura et al. (2020b)]. Kimura et al. (2020b) simulated light curves produced by the thermal instability in the accretion disk tilted out of the binary orbital plane by implementing the mass supply pattern to the tilted disk under an assumption that the stream trajectory was calculated as the ballistic one of a gas particle. Although details in the mass supply pattern are different between these two cases, the resultant variation in light curves is similar to each other where a higher foverflow in our case corresponds to higher tilt angles there. By comparing two figures, figures 16 and 17 of Kimura et al. (2020b), we find that the inner part of the disk is persistently in the hot state (possibly a necessary condition for explaining the 2021 anomalous event of SS Cyg) if the mass transfer rate is a little higher than that of the standard model in this paper. To confirm this, we have performed the same type of simulations as figure 16 but with a higher mass transfer rate than that of our standard model, and the results are shown in figure 19. We see in figure 19 that the minimum light level is elevated to around 10 mag if foverflow is sufficiently high because the inner part of the disk is persistently in the hot state in such a case. This behavior is close to that of corresponding cases for tilted disks [see figure 16 of Kimura et al. (2020b)].

Same as figure 16 but $\dot{M}_{\rm tr} = 1.5 \times 10^{17}\:$g s−1. See the caption of figure 4 for a description of the black, green, and orange lines in the lower panel.
Fig. 19.

Same as figure 16 but |$\dot{M}_{\rm tr} = 1.5 \times 10^{17}\:$|g s−1. See the caption of figure 4 for a description of the black, green, and orange lines in the lower panel.

The mass-transfer rate in SS Cyg may be a little higher than |$\dot{M}_{\rm tr} = 10^{16.9}\:$|g s−1, a value that we have adopted in the standard model (see section 2). We have thus tried one more simulation in which foverflow is amplified to 0.8 on 240 d in the case of |$\dot{M}_{\rm tr} = 1.5 \times 10^{17}\:$|g s−1. The results are displayed in figure 20. We see from the middle panel of figure 20 that the system shows low-amplitude fluctuations and occasionally large brightening after 240 d. Also, the average flux level is as high as 10 mag during oscillations because the inner disk always stays in the hot state and does not enter the cool state. This behavior is similar to the first ∼60 d phenomenon of the 2021 anomalous event in SS Cyg (see also figure 2), though the amplitude of oscillations is a little smaller, and the duration of oscillations is a little shorter than observations.

Results of our simulations in the case where $\dot{M}_{\rm tr} = 1.5 \times 10^{17}\:$g s−1 and foverflow changes on 240 d. Here, foverflow is increased to 0.8 after 240 d. The corresponding foverflow is denoted as a dashed line in figure 19.
Fig. 20.

Results of our simulations in the case where |$\dot{M}_{\rm tr} = 1.5 \times 10^{17}\:$|g s−1 and foverflow changes on 240 d. Here, foverflow is increased to 0.8 after 240 d. The corresponding foverflow is denoted as a dashed line in figure 19.

To reproduce the increment of the V-band magnitude in the first ∼60 d light curve of the 2021 anomalous event in SS Cyg against the normal quiescence, the inner disk should stay in the hot state during this event. It is unlikely that the WD or the secondary star or the bright spot becomes more than an order of magnitude brighter. The maximum brightness of the cool disk is estimated from the maximum effective temperature of the cool state to be less than 11.5 mag in the V band, which is much fainter than the minimum V-band magnitude (=∼10.5 mag) of the observed light curve during the first ∼60 d of the anomalous event (see figure 2). Our simulation displayed in figure 20 at least satisfies the condition that the innermost part of the disk never returns to the cool state.

4 Discussion

4.1 Does the 2021 anomalous event originate from the enhanced viscosity?

Kimura et al. (2021) suggested that the enhancement of viscosity in the cool state would reproduce the 2021 anomalous event and its forerunner in SS Cyg. We find that some of the light curves in our simulations of the enhanced viscosity model are similar to the observed light curve during the forerunner of the anomalous event. It is, however, difficult to reproduce the first ∼60 d behavior of the 2021 anomalous event because our simulated light curve still shows an outbursting behavior with amplitude as large as ∼2 mag even in the case of αcoolhot = 1.0, contrary to the observed light fluctuations with amplitude ∼1 mag. As discussed in subsection 3.2, the light curve behavior in the case of the same α for the hot and cool states may be sensitive to the S-shaped thermal-equilibrium curve used in the simulation. If one uses other S-shaped curves with a shorter intermediate branch, the outburst amplitude may become ∼1 mag, but this is another problem.

The gas-stream overflow may play an important role in reproducing the 2021 anomalous event. When the stream overflow is taken into account in our simulations, small-amplitude inside-out outbursts begin to occur (see also figures 16, 17, and 18). If the mass transfer rate is a little higher than that of the standard model in our simulations, the stream overflow generates a light curve similar to the anomalous event (see figure 20). The overflow is likely to have occurred at least during the first ∼60 d of the anomalous event.

4.2 Can the gas-stream overflow explain the observed X-ray luminosity?

As pointed out in section 2, there is a large discrepancy between the observed X-ray luminosity in quiescence of SS Cyg and that predicted from the standard model in which the inner edge of the cold disk extends down near the WD surface. As discussed in section 2, we adopt the evaporation model, in which the geometrically thin viscous disk is truncated at the inner edge, rin, below which the disk expands to the hot spherical corona. In such cases, if rin or αcool is increased, |$\dot{M}_{\rm acc} (r_{\rm in})$| increases, which makes the predicted X-ray luminosity higher. However, there must be some limit to doing so. We have tried simulations where the inner edge of the disk is chosen to be rin = 1.0 × 1010 cm by assuming that a large part of the inner disk is evaporated, and calculated the X-ray luminosity in quiescence to be ∼2 × 1032 erg s−1. However, this X-ray luminosity is still lower than the observed one. It is noted here that the X-ray luminosity in quiescence estimated here depends on the thermal equilibrium curve used. The thermal equilibrium curve in this paper is based on the computations in Mineshige and Osaki (1983). If another thermal equilibrium curve such as that of Hameury et al. (1998) is used in which the convective energy transport is much more efficient in the vertical structure calculations, the X-ray luminosity of the cold disk could become larger, as large as a factor 3 of our estimate (Lasota et al. 2008).

The main difficulty in explaining the fairly high X-ray flux observed in quiescence in SS Cyg lies in the difficulty in the first place of transferring a large amount of mass from the outer disk to the inner disk by viscous diffusion because of low viscosity in quiescence. To solve this difficulty, one possible solution is to consider a situation in which part of the overflowing gas is directly transferred to the inner disk by the gas stream and to the inner edge of the disk from where the disk matter evaporates into the hot corona (i.e., a kind of shortcut in mass supply). However, we do not here discuss any specific mechanism for how the overflowing gas stream is converted there to the hot coronal matter. We simply assume that the disk matter at the inner disk edge evaporates into the hot corona. For example, the observed X-ray luminosity during normal quiescence is reproduced if |$\sim\!\! 5\%$| of the transferred mass is deposited directly to the inner edge and then to the hot corona via the gas-stream overflow. Even if |$\sim\!\! 5\%$| of the transferred mass overflows the disk surface, the outburst behavior does not change so much in comparison with that of the standard model. The gas-stream overflow may alleviate the discrepancy between the observed and simulated X-ray luminosity.

We next consider the X-ray luminosity during the 2021 anomalous event. The light curve displayed in figure 20 is closest to the 2021 anomalous event in SS Cyg among our simulations. In this case, the inner disk always stays in the hot state, and the accretion rate at rin during the mid-brightness interval after 240 d is as high as ∼4 × 1016 g s−1 on average. The X-ray luminosity is calculated to be ∼5 × 1033 erg s−1 by equation (1), which is a little higher but consistent with the observed X-ray luminosity. Although we have used the result in figure 20 in this calculation, the above discussion is not limited to this specific case in this simulation. The X-ray luminosity in the 2021 anomalous event would be reproduced by other simulations where the inner disk keeps the hot state and the averaged optical flux is as high as ∼10 mag because a large amount of mass accreted from the inner disk edge to the X-ray emitting hot corona. In this case, it is not required that the overflowing gas is directly injected into the hot corona to explain the observed X-ray luminosity. Nevertheless, the stream overflow is needed in our simulations to reproduce the light-curve behavior similar to the first ∼60 d light curve of the 2021 anomalous event in SS Cyg and to increase the mass accretion rate at the inner disk edge.

5 Summary

We have performed numerical simulations of the light curves for the 2021 anomalous event in SS Cyg by varying mass-transfer rates, varying viscosity in the cool state, and varying overflow rates of the gas stream. Our findings are listed below.

  • If the mass-transfer rate increases, the cycle length of long outbursts becomes shorter initially. If the mass-transfer rate is enhanced further, the duration of long outbursts gradually lengthens, and the cycle length becomes longer. Finally, the long outburst persists all the way, and the system enters the hot steady state.

  • In the enhanced mass-transfer model, the outburst amplitude remains large even when the system approaches to the hot steady state, and we do not find any cases in which the outburst amplitude diminishes and finally vanishes with the increase in the mass-transfer rate.

  • It is found that when the mass transfer rate is gradually increased, the transition in light curves from the DN-type outburst state to the nova-like state occurs rather suddenly within a very narrow range in mass transfer rate; a condition favorable to explain the Z Cam phenomenon.

  • If the viscosity in the cool state is enhanced, inside-out outbursts frequently occur, and the quiescent level increases.

  • If αcool is chosen to be equal to αhot, it is found that the inner disk never drops to the cool state for parameters corresponding to our standard model for SS Cyg. Even in this case, fairly large amplitude oscillatory light variations are found to occur with an amplitude of 2 mag, contrary to a naive expectation that small-amplitude fluctuations may result in such a case.

  • We have studied the effect of gas-stream overflow beyond the outer disk rim in our light curve simulations where the mass by overflowing gas is assumed to be supplied to the inner disk around the circularization radius. If the overflow rate is increased in our simulations, a new type of inside-out type outbursts with an amplitude of 1 mag occurs, and its number increases. It is found that no models with gas-stream overflow can reproduce the light curve of the 2021 anomalous event of SS Cyg with parameters of our standard model.

  • Within our simulations, only a model in which the gas-stream overflow is considered, and the mass transfer rate is a little higher than that of our standard model, may reproduce a light curve with fluctuations having amplitudes of less than 1 mag as observed in the 2021 anomalous event of SS Cyg.

  • The X-ray flux calculated from our simulations with the standard parameters is much lower than the observed X-ray flux during normal quiescence in SS Cyg. If a few percent of the overflowing gas is provided to the inner edge of the disk, the observed flux could be explained. On the other hand, the direct injection of the overflowing gas to the inner edge of the disk is not necessary to reproduce the observed X-ray flux in the case of the 2021 anomalous event since the accretion rate at the inner disk edge is high enough.

We conclude that the enhanced mass transfer cannot reproduce the 2021 anomalous event and its forerunner in SS Cyg. We have confirmed that the enhancement of viscosity in the cool state may reproduce the observational features of the forerunner of SS Cyg, as suggested in Kimura et al. (2021). However, it is not enough to reproduce the light-curve behavior of the first 60 days of the 2021 anomalous event. The gas-stream overflow may be necessary to explain this event. Also, the stream overflow may play an important role in reproducing the observed X-ray flux during normal quiescence.

Acknowledgements

M. Kimura acknowledges support by the Special Postdoctoral Researchers Program at RIKEN. We are thankful to many amateur observers for providing a lot of data used in this research. This work was financially supported by Japan Society for the Promotion of Science Grants-in-Aid for Scientific Research (KAKENHI) Grant Numbers JP20K22374 (MK) and JP21K13970 (MK).

Appendix 1 Details of the modifications of our numerical code

In Kimura et al. (2020b), we treated the tilted disk and mixed logarithmic and linear meshes. In this study, we set N0 concentric annuli from the innermost disk radius, rin, to the tidal truncation radius, rtidal, equally spaced in |$\sqrt{r}$| as used in Cannizzo (1993b). We first set N0 to be 200 and the number of meshes (N) in the time-dependent calculation is variable. The center of each annulus is given by
(A1)
 
(A2)
The mass supply pattern to the disk is also modified. We basically input the mass transferred from the secondary star to the outer disk edge with the width of drs. The source term s(r) is given by
(A3)
where |$\dot{M}_{\rm tr}$| represents the mass transfer rate. We input the heating by the energy dissipation of the gas stream from the secondary star to the outermost NS meshes between rNdrs and rN. In this paper, drs is 0.02a in the case of N0 = 200, and a typical value of NS is 7 in this paper.

Appendix 2 Extra tidal heating due to the tidal truncation

In our formulation, we assume that the tidal truncation radius works as a solid brick wall, that is, when the disk reaches the tidal truncation radius and tries to expand further, it is assumed that the expansion stops at the tidal truncation radius due to the strong tidal torques. As for the tidal torque, our formulation is basically the same as that in our previous paper (Kimura et al. 2020b). The tidal torque is given by the following equation [see equation (7) of Kimura et al. (2020b)] if the disk radius is less than the tidal truncation radius:
(A4)
where ν is the kinematic viscosity (see Kimura et al. 2020b). In this paper, we adopt the constant cω to be 0.4 s−1, which is much smaller than that (cω ≃ 9.5 s−1) expected in the case when the disk radius in the hot steady state just reaches the tidal truncation radius in this formulation. The tidal torque by equation (A4) is thus weak and we have to remove the extra angular momentum if the disk radius tries to exceed the tidal truncation radius. In our previous paper (Kimura et al. 2020b), the removal of the extra angular momentum at the tidal truncation radius was automatically taken into account, but the extra tidal heating due to this effect was not considered there.
Here we take into account this heating effect, and we implement in our code the extra tidal heating at the tidal truncation radius. To do so, we need to know the extra tidal torque by the tidal truncation (expressed here as Dext) in its explicit form. Once we know this extra tidal torque, we can easily calculate the extra tidal heating. Let us now consider a fictitious situation in which the tidal truncation were not exerted and so the disk were allowed to expand beyond rtidal and to reach rN. The total angular momentum of the disk is given by
(A5)
where the quantities with and without the subscript “new” denote before and after a time step. Here, |$h_{N-1/2}^{\rm new}$| is defined as |$\sqrt{G M_1 (r_{N-1} + r_N)/2}$|⁠. The other symbols appearing in equation (A5) are explained in appendix 1 in Kimura et al. (2020b).
In our formulation, it is assumed that the disk radius is limited to the tidal truncation radius by the extra tidal torque when the disk tries to exceed the tidal truncation radius. In such a case, the conservation of the disk angular momentum given after a time step (dt) by equation (A5) is written as follows :
(A6)
where |$J_{\rm disk}^{\rm new^{\prime }}$| and |$h_{N-1/2}^{\rm new^{\prime }}$| are the total angular momentum of the disk and the specific angular momentum of the largest annulus, respectively, in the case when the disk stops at the tidal truncation radius due to Dext. The quantity |$h_{N-1/2}^{\rm new^{\prime }}$| is given by |$\sqrt{G M_1 (r_{N-1} + r_{\rm tidal})/2}$|⁠.
By comparing these two cases, we then obtain Dext that is given by
(A7)
That is, the extra tidal torque, Dext, is obtained from a fictitious calculation in which the disk is allowed to expand beyond the tidal truncation radius.
We now calculate the extra tidal heating by assuming that it is exerted in the outermost region of the disk with the width of drd. We set drd to be 0.02a in this study. If we introduce a quantity, dext, which is given by
(A8)
we then add the dissipation rate by the extra tidal torque to the outermost Nd meshes between rNdrd and rN by
(A9)
where Ωi − 1/2 is the Keplerian angular velocity at ri − 1/2 written by |$\sqrt{G M_1 {r_{i-1/2}^3}}$| and ωorb is the orbital angular velocity, respectively.

We may note here that drd affects the number of short outbursts sandwiched by two long outbursts. As drd becomes smaller, the number of short outbursts between two long outbursts becomes smaller. This is because an outburst easily enters a long plateau stage with smaller drd because of strong tidal torques at the outer disk edge when the disk radius reaches the tidal truncation radius. We exhibit an example of light curves in the case of drd = 0.005a in figure 21.

Same as figure 4 but in the case of drd = 0.005a.
Fig. 21.

Same as figure 4 but in the case of drd = 0.005a.

Footnotes

References

Bailer-Jones
 
C. A. L.
,
Rybizki
 
J.
,
Fouesneau
 
M.
,
Mantelet
 
G.
,
Andrae
 
R.
 
2018
,
AJ
,
156
,
58

Balbus
 
S. A.
,
Hawley
 
J. F.
 
1991
,
ApJ
,
376
,
214

Balman
 
Ş.
,
Revnivtsev
 
M.
 
2012
,
A&A
,
546
,
A112

Buat-Ménard
 
V.
,
Hameury
 
J.-M.
,
Lasota
 
J.-P.
 
2001
,
A&A
,
369
,
925

Cannizzo
 
J. K.
 
1993a
,
Adv. Ser. Astrophys. Cosmology
,
9
,
6

Cannizzo
 
J. K.
 
1993b
,
ApJ
,
419
,
318

Cannizzo
 
J. K.
 
1993c
, in
Accretion Disks in Compact Stellar Systems
, ed.
Wheeler
 
J. C.
(
Singapore
:
World Scientific Publishing
),
6

Cannizzo
 
J. K.
,
Mattei
 
J. A.
 
1992
,
ApJ
,
401
,
642

Cannizzo
 
J. K.
,
Wheeler
 
J. C.
 
1984
,
ApJS
,
55
,
367

Dubus
 
G.
,
Otulakowska-Hypka
 
M.
,
Lasota
 
J.-P.
 
2018
,
A&A
,
617
,
A26

Hameury
 
J. M.
 
2020
,
Adv. Space Res.
,
66
,
1004

Hameury
 
J.-M.
,
Knigge
 
C.
,
Lasota
 
J.-P.
,
Hambsch
 
F. J.
,
James
 
R.
 
2020
,
A&A
,
636
,
A1

Hameury
 
J.-M.
,
Menou
 
K.
,
Dubus
 
G.
,
Lasota
 
J.-P.
,
Hure
 
J.-M.
 
1998
,
MNRAS
,
298
,
1048

Hessman
 
F. V.
 
1999
,
ApJ
,
510
,
867

Hessman
 
F. V.
,
Robinson
 
E. L.
,
Nather
 
R. E.
,
Zhang
 
E.
 
1984
,
ApJ
,
286
,
747

Hill
 
C. A.
,
Smith
 
R. C.
,
Hebb
 
L.
,
Szkody
 
P.
 
2017
,
MNRAS
,
472
,
2937

Hōshi
 
R.
 
1979
,
Prog. Theor. Phys.
,
61
,
1307

Ichikawa
 
S.
,
Osaki
 
Y.
 
1992
,
PASJ
,
44
,
15

Ishida
 
M.
,
Okada
 
S.
,
Hayashi
 
T.
,
Nakamura
 
R.
,
Terada
 
Y.
,
Mukai
 
K.
,
Hamaguchi
 
K.
 
2009
,
PASJ
,
61
,
S77

Kato
 
T.
 
2019
,
PASJ
,
71
,
20

Kato
 
T.
,
Kunjaya
 
C.
 
1995
,
PASJ
,
47
,
163

Kimura
 
M.
 et al.  
2021
,
PASJ
,
73
,
1262

Kimura
 
M.
,
Osaki
 
Y.
,
Kato
 
T.
 
2020a
,
PASJ
,
72
,
94

Kimura
 
M.
,
Osaki
 
Y.
,
Kato
 
T.
,
Mineshige
 
S.
 
2020b
,
PASJ
,
72
,
22

Kunze
 
S.
,
Speith
 
R.
,
Hessman
 
F. V.
 
2001
,
MNRAS
,
322
,
499

Lasota
 
J.-P.
,
Dubus
 
G.
,
Kruk
 
K.
 
2008
,
A&A
,
486
,
523

Lin
 
D. N. C.
,
Faulkner
 
J.
,
Papaloizou
 
J.
 
1985
,
MNRAS
,
212
,
105

Lubow
 
S. H.
,
Shu
 
F. H.
 
1976
,
ApJ
,
207
,
L53

Menou
 
K.
,
Hameury
 
J.-M.
,
Lasota
 
J.-P.
,
Narayan
 
R.
 
2000
,
MNRAS
,
314
,
498

Meyer
 
F.
 
1984
,
A&A
,
131
,
303

Meyer
 
F.
,
Meyer-Hofmeister
 
E.
 
1983
,
A&A
,
121
,
29

Meyer
 
F.
,
Meyer-Hofmeister
 
E.
 
1994
,
A&A
,
288
,
175

Mineshige
 
S.
,
Osaki
 
Y.
 
1983
,
PASJ
,
35
,
377

Mineshige
 
S.
,
Osaki
 
Y.
 
1985
,
PASJ
,
37
,
1

Mineshige
 
S.
,
Shields
 
G. A.
 
1990
,
ApJ
,
351
,
47

Nauenberg
 
M.
 
1972
,
ApJ
,
175
,
417

Osaki
 
Y.
 
1974
,
PASJ
,
26
,
429

Osaki
 
Y.
 
1989
,
PASJ
,
41
,
1005

Osaki
 
Y.
 
1995
,
PASJ
,
47
,
L11

Osaki
 
Y.
 
1996
,
PASP
,
108
,
39

Paczynski
 
B.
 
1977
,
ApJ
,
216
,
822

Pojmanski
 
G.
 
1986
,
Acta Astron.
,
36
,
69

Provencal
 
J. L.
,
Shipman
 
H. L.
,
Høg
 
E.
,
Thejll
 
P.
 
1998
,
ApJ
,
494
,
759

Robertson
 
J. W.
,
Honeycutt
 
R. K.
,
Turner
 
G. W.
 
1995
,
PASP
,
107
,
443

Schreiber
 
M. R.
,
Hessman
 
F. V.
 
1998
,
MNRAS
,
301
,
626

Shakura
 
N. I.
,
Sunyaev
 
R. A.
 
1973
,
A&A
,
24
,
337

Smak
 
J.
 
1984
,
Acta Astron.
,
34
,
161

Waagen
 
E. O.
 
2020
,
AAVSO Alert Notice
,
720
,
1

Warner
 
B.
 
1995
,
Cataclysmic Variable Stars
(
Cambridge
:
Cambridge University Press
)

Wheatley
 
P. J.
,
Mauche
 
C. W.
,
Mattei
 
J. A.
 
2003
,
MNRAS
,
345
,
49

This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://dbpia.nl.go.kr/journals/pages/open_access/funder_policies/chorus/standard_publication_model)