A Self-Consistent Model for a Full Cycle of Recurrent Novae -- Wind Mass Loss Rate and X-Ray Luminosity

TL;DR

This paper develops a self-consistent model for recurrent novae, accurately calculating wind mass loss rates and X-ray luminosity, improving understanding of nova light curves and evolution.

Contribution

It introduces a combined method using wind solutions and stellar evolution codes for precise modeling of nova outbursts.

Findings

Accurate wind mass loss rates determined during nova outbursts.

Calculated light curves for WDs of 1.2 and 1.38 solar masses.

Identified phases of X-ray emission during nova evolution.

Abstract

An unexpectedly slow evolution in the pre-optical-maximum phase was suggested in the very short recurrence period nova M31N 2008-12a. To obtain reasonable nova light curves we have improved our calculation method by consistently combining optically thick wind solutions of hydrogen-rich envelopes with white dwarf (WD) structures calculated by a Henyey-type evolution code. The wind mass loss rate is properly determined with high accuracy. We have calculated light curve models for 1.2 and 1.38 M_sun WDs with mass accretion rates corresponding to recurrence periods of 10 and 1 yr, respectively. The outburst lasts 590/29 days in which the pre-optical-maximum phase is 82/16 days, for 1.2/1.38 M_sun, respectively. Optically thick winds start at the end of X-ray flash and cease at the beginning of supersoft X-ray phase. We also present supersoft X-ray light curves including a prompt X-ray flash and later supersoft X-ray phase.

Tables

Table 1. Table 1 Summary of Recurrent Nova Models

SubjectaaThe WD mass and mass accretion rate are given parameters and recurrence period and other quantities are the results of calculation.	units	1.2  $M_{\odot }$	1.38  $M_{\odot }$
${\dot{M}}_{\mathrm{acc}}$	${10}^{-7}{M}_{\odot }$ yr-1	2.0	1.6
$P_{\mathrm{rec}}$	yr	9.9	0.95
$\log T_{\mathrm{core}}$	K	7.93	8.03
$M_{\mathrm{acc}}$	${10}^{-7}{M}_{\odot }$	17	1.4
$M_{\mathrm{ig}}$	${10}^{-7}{M}_{\odot }$	26	2.0
$L_{\mathrm{nuc}}^{\max }$	${10}^6{L}_{☉}$	2.2	3.9
$\log T_{\mathrm{nuc}}^{\max }$	K	8.11	8.23
$t_{\mathrm{nova}}$bbduration of the flash with $L_{\mathrm{ph}}>{10}^4{L}_{☉}$.	day	590	29
$t_{\mathrm{wind}}$	day	130	19
$t_{X-\mathrm{flash}}$ccduration of the X-ray flash and SSS phase: $L_X(0.3-1.0$ keV) $>{10}^3{L}_{☉}$ for 1.2$M_{☉}$, but $>{10}^4{L}_{☉}$ for $1.38M_{☉}$.	day	110	0.77
$t_{\mathrm{SSS}}$ccduration of the X-ray flash and SSS phase: $L_X(0.3-1.0$ keV) $>{10}^3{L}_{☉}$ for 1.2$M_{☉}$, but $>{10}^4{L}_{☉}$ for $1.38M_{☉}$.	day	440	7.8
$t_{\mathrm{rise}}$dd the time elapsed between the time of the peak nuclear luminosity and the time of the peak mass loss rate.	day	82	16
$t_{\mathrm{decay}}$eeoptical decay time calculated as $t_{\mathrm{decay}}=t_{\mathrm{nova}}-t_{\mathrm{rise}}$.	day	510	13

Equations

$${\dot{M}}_{\mathrm{evol}}^{n+1}={\dot{M}}_{\mathrm{evol}}^n\exp [10(\log R_{\mathrm{ph}}^n-\log R_{\mathrm{ph}}^{n-1})].$$

$$U\equiv \frac{d\ln M_r}{d\ln r}=\frac{4\pi r^3\rho }{M_r},$$

$$V\equiv -\frac{d\ln P}{d\ln r}=\frac{G{M}_r\rho }{rP},$$

$$L_{\mathrm{Edd}}=\frac{4\pi cG{M}_{\mathrm{WD}}}{\kappa },$$

Full text

††slugcomment: to appear in the Astrophysical Journal

1. INTRODUCTION

A nova is a thermonuclear runaway event on a mass-accreting white dwarf (WD) (Nariai et al., 1980; Iben, 1982; Prialnik, 1986; Sion et al., 1979; Sparks et al., 1978). Theoretical light curves of nova outbursts have been calculated mainly based on a steady-state approximation (i.e., optically thick wind theory: Kato & Hachisu, 1994) with successful reproduction of characteristic properties, such as multiwavelength light curves including supersoft X-ray phase, of a number of classical novae (Hachisu & Kato, 2006, 2007, 2010, 2014, 2015, 2016a, 2016b; Hachisu et al., 2008; Kato et al., 2009). In spite of these successful results, this method is inapplicable when a steady-state is not a good approximation. Such failures occur for very fast evolutions in recurrent novae like U Sco and RS Oph and pre-maximum-optical phases of all kinds of novae. Thus, we need to take a different approach to follow the entire evolution of short period recurrent novae, which we aim in the present work.

Multicycle nova outbursts have been calculated with Henyey-type evolution codes. These codes, however, meet numerical difficulties when the nova envelope expands to a giant size. To continue numerical calculation beyond the extended stage, various authors have adopted various mass loss schemes and approximations (Prialnik & Kovetz, 1995; Kovetz, 1998; Denissenkov et al., 2013; Idan et al., 2013; Wolf et al., 2013a, b; Kato et al., 2014, 2015; Tang et al., 2014). These works, however, paid little attention in reproducing reliable nova light curves, so there is still room for improvement.

Recently, X-ray flash of a nova, a brief X-ray brightening before the optical maximum, came into attention with realistic observational plans (Morii et al., 2016; Kato et al., 2015, 2016). Kato et al. (2016) showed that outbursts of recurrent novae are so weak that their pre-optical-maximum evolution could be very slow. They reported the nondetection of X-ray flux during 8 days prior to the 2015 outburst of M31N 2008-12a. They estimated that the X-ray flash had already occurred long ($9-16$ days) before the optical maximum. Such unexpected slow evolution of the very early phase, in spite of the very fast optical decline, indicates the fact that we do not fully understand nova outbursts yet.

Recurrent novae of very short recurrence periods ($P_{\rm rec}<10$ yr) were recently reported: 0.5 or 1 yr for M31N 2008-12a (Darnley et al., 2014, 2015, 2016; Henze et al., 2014, 2015; Tang et al., 2014; Darnley et al., 2016), 5 yr for M31N 1963-09c (Williams et al., 2015), $\sim 6$ yr for LMCN 1968-12a (Darnley et al., 2016), and 10 yr for U Sco (e.g., Schaefer, 2010). Recurrent novae are binaries harboring a massive white dwarf (WD). Short recurrence period novae occur in very massive WDs of $M_{\rm WD}\gtrsim 1.2~{}M_{\odot}$ with very high mass accretion rates of $\dot{M}_{\rm acc}\gtrsim 1.5\times 10^{-7}M_{\odot}$ yr*-1* (Prialnik & Kovetz, 1995; Wolf et al., 2013a, b; Tang et al., 2014; Kato et al., 2014). Such massive WDs are considered to be candidates for Type Ia supernova (SN Ia) progenitors (Hachisu et al., 1999a, b; Hachisu & Kato, 2001; Hachisu et al., 2010; Han & Podsiadlowski, 2004; Li & van den Heuvel, 1997; Kato & Hachisu, 2012). SNe Ia play very important roles in astrophysics as a standard candle for measuring cosmological distances (Perlmutter et al., 1999; Riess, 1998) and as main producers of iron group elements in the chemical evolution of galaxies (Kobayashi et al., 1998). However, their immediate pre-explosion progenitors remain unclear (e.g., Maoz et al., 2014; Pagnotta & Schaefer, 2014). Thus, studies of very short recurrence period novae are very important.

Our aim of this work is to establish a self-consistent numerical method for modeling a full cycle of nova outbursts. We combine the structure of optically thick winds with a model obtained by our Henyey-type code at each phase of a nova, which gives a mass loss rate and surface values of temperature and luminosity that are required to calculate multiwavelength light curves. In order to obtain successful fittings in our first attempt we have chosen a $1.2~{}M_{\sun}$ WD with a very high mass-accretion rate of $\dot{M}_{\rm acc}=2.0\times 10^{-7}~{}M_{\odot}$ yr*-1* because the nova explosion is weak and expansion is slow, so the envelope structure is close to that of steady-state envelope. We have also calculated wind structures fitted to the evolution models of a $1.38~{}M_{\sun}$ WD with 1 yr recurrence period published in Kato et al. (2016). Section 2 describes our numerical method, and Section 3 presents our results on the $1.2~{}M_{\sun}$ WD in detail. The $1.38~{}M_{\sun}$ WD model is briefly described in Section 4. Discussion and conclusions follow in Sections 5 and 6, respectively.

2. NUMERICAL METHOD

We calculated WD structures from the center of the WD up to the photosphere. We adopt the optically thick winds (e.g., Kato, 1983, 1985; Kato & Hachisu, 1994) as the main mechanism of mass loss during the nova outburst. The optically thick winds have been applied to a number of novae and successfully reproduced light curves and characteristic properties (e.g., Hachisu & Kato, 2006, 2007; Hachisu et al., 2008; Hachisu & Kato, 2010, 2014, 2015, 2016a, 2016b; Kato et al., 2009). The interior structures of WDs are calculated with a Henyey code (Section 2.1). When the optically thick winds occur the structure of mass-losing envelope is calculated as described in Section 2.2. We use the $UV$ plane for our fitting process of outer and inner structures (Section 2.3) in the way as described in Section 2.4.

2.1. Time-dependent Calculation

Evolution models of mass-accreting WDs are calculated by using the same Henyey-type code as in Kato et al. (2014, 2016) and Hachisu et al. (2016). Accretion energy outside the photosphere is not included. We have calculated multicycle nova outbursts until the characteristic properties of flashes become unchanged. The evolution is started with the equilibrium model for a given mass accretion rate (Kato et al., 2014), in which the WD core is thermally in a steady state. Thus, the WD core temperature is not a free parameter.

In a later phase of a shell flash, the photospheric radius (as well as the luminosity) increases to the point where the outer boundary condition of our evolution code fails because the density becomes too small at $T\sim 2\times 10^{5}$ K (at the Fe opacity peak, Iglesias & Rogers, 1996). To avoid the numerical difficulty, we start mass loss, which hinders the growth of the radius, when the photospheric radius $R_{\rm ph}^{n}$ at time-step $n$ exceeds a radius $R_{0}$. The mass loss rate for the next step $\dot{M}_{\rm evol}^{n+1}$ is calculated as

[TABLE]

At the beginning of mass loss $\dot{M}_{\rm evol}^{n}$ is replaced with $\dot{M}_{0}$ ($<0$), which is an arbitrary parameter to give an initial mass loss rate at $R_{\rm ph}=R_{0}$. This equation gives a rapid increase (decrease) in $|\dot{M}_{\rm evol}|$ when the photospheric radius is increasing (decreasing). Equation (1) is used until $R_{\rm ph}$ becomes smaller than a radius $R_{1}$, then $\dot{M}_{\rm evol}$ is switched to the accretion rate given at the start of the calculation. In the $1.38~{}M_{\sun}$ model, we restarted the accretion immediately after the mass loss stopped. This is because we adopt internal models of Kato et al. (2016). In the $1.2~{}M_{\sun}$ model, we restarted accretion when the outburst is almost finished, to compare the duration of the supersoft X-ray phase with that obtained from the static sequence which does not include mass-accretion (see Section 3.3).

The values of the parameters, $\dot{M}_{0}$, $R_{0}$, and $R_{1}$ are chosen case by case. The details of these parameters are explained in Section 5.2.

2.2. Wind Mass Loss

Shortly after the onset of unstable nuclear burning, the envelope expands and optically thick winds are accelerated. The occurrence of winds is detected at the photosphere using the BC1 surface boundary condition in Kato & Hachisu (1994). To obtain envelope structures with winds, we solve the equations of motion, mass continuity, radiative energy transfer by diffusion, and energy conservation assuming steady-state conditions and spherical symmetry (Kato & Hachisu, 1994). The boundary conditions of the wind solution are set at the photosphere and critical point (=sonic point) (Kato & Hachisu, 1994), which determines two values, i.e., the temperature and radius of the critical point, ($T_{\rm cr}$, $R_{\rm cr}$). We calculated the optically thick wind solutions down to a chosen inner boundary, i.e., where the temperature rises to a specified value (e.g. $T=T_{\rm f}$) and obtain the envelope mass of the wind ($M_{\rm env}^{\rm wind}$) between the photosphere and the point of $T=T_{\rm f}$. In short, for a chosen temperature ($T=T_{\rm f}$) a wind solution is characterized by $T_{\rm cr}$ and $R_{\rm cr}$. The other quantities are automatically calculated, such as the mass loss rate, wind velocity, and local luminosity ($L_{r}$). In general $L_{r}$ is the sum of radiative and convective luminosities. In our calculation convection does not occur at the matching point and $L_{r}$ is always equal to the radiative luminosity.

2.3. Fitting in the $UV$ Plane

In the fitting procedure we use the $UV$ plane (Chandrasekhar, 1939; Schwarzschild, 1958; Hayashi et al., 1962), where $U$ and $V$ are the homology invariants defined as

[TABLE]

and

[TABLE]

where $M_{r}$ is the mass within the radius $r$.

In general, the fitting point of the wind solution ($U_{\rm w}$, $V_{\rm w}$) moves downward in the $UV$ plane if we increase the wind mass loss rate (i.e., increase $T_{\rm cr}$ or $R_{\rm cr}$). On the other hand, the fitting point of Henyey-code solutions ($U_{\rm f}$, $V_{\rm f}$) moves upward for a larger mass loss rate. Thus, we find a unique matching point for the same wind mass loss rate and envelope mass in the $UV$ plane. Matching $U$ and $V$ at $T_{\bf f}$ guarantees the continuation of the mechanical structure at the fitting point so that the radius and density are automatically matched there (see Equations (2) and (3)). In addition, as there is no energy source/sink around the fitting place the local luminosity $L_{r}$ is automatically fitted.

2.4. Fitting Procedure

In the wind phase we need several steps of iteration until we obtain a final result in which the wind mass loss is consistently taken into account. Our fitting procedure is as follows.

(1) First, we calculate nova outbursts with the Henyey code assuming a time-variable mass loss rate until the shell flashes reach a limit cycle. We use the last one cycle in the following procedure.

(2) When the mass loss occurs on the Henyey code, we calculate optically thick wind solutions with a specified parameters, ($T_{\rm cr}$, $R_{\rm cr}$), down to the fitting place $T=T_{\rm f}$. If the fitting condition ($U_{\rm w}$, $V_{\rm w}$) = ($U_{\rm f}$, $V_{\rm f}$) is not satisfied, we change $T_{\rm cr}$ and $R_{\rm cr}$ until the condition is satisfied. We need typically 20 iterations done by hand.

(3) If the envelope mass $M_{\rm env}^{\rm wind}$ of the wind solution does not match to that of the interior solution we change the fitting place $T_{\rm f}$ and repeat the same process until the two envelope masses are sufficiently close to each other. A few to several iterations are needed for the process.

(4) We replaced the outer part ($T<T_{\rm f}$) of the evolution model with a well matched wind solution. Then, we obtain a continuous structure from the WD core up to the photosphere that consistently includes the optically thick wind at a specified time during the wind phase. The photospheric temperature, velocity, and wind mass loss rate etc. are uniquely determined.

(5) We need to do such a procedure (2) - (4) for each timestep. However, our wind phase contains $\sim 10,000$ timesteps, which is far beyond the capacity of manual handlings. For this reason, we select about 30 epochs and apply the iteration procedure (2)-(4) by hand. Thus, we obtain stellar structure models for one cycle of nova outburst.

(6) The mass loss rates of the optically thick winds, which varies time to time, is generally different from those assumed in step (1). We repeated the process (1)-(5) assuming different mass loss rates $\dot{M}_{0}$ in Equation (1) until the calculated wind mass loss rates well match to the assumed values. It takes several iterations.

Among the above process, the most difficult process is step (1). A Henyey-type code often meets numerical difficulties without assuming large mass loss rates when the envelope expands. If we assume too large mass loss rates, however, the procedure (1)-(6) does not converge.

Note that the fitting point between the wind and the interior should be carefully chosen especially in the first trial. The fitting point should be deep enough to be insensitive to treatments in the surface region, especially when we assume inadequate mass loss rates in the trial process. This point also should be shallow enough to ensure the time-dependent term $\epsilon_{\rm g}$ (gravitational energy release rate per unit mass: see equation (1) in Hachisu et al., 2016) is negligible.

Our procedure (1)-(6) currently needs many iterations by hand, and thus, takes much human-time. For this reason, we have calculated only 1.2 and 1.38 $M_{\sun}$ models. We certainly need to improve our fitting procedure.

3. 1.2 $M_{\odot}$ Model

We have calculated two models of 1.2 and 1.38 $M_{\odot}$. The 1.38 $M_{\odot}$ model are almost the same as those partly published in Kato et al. (2016). Thus, we first describe the 1.2 $M_{\odot}$ model and then summarize for the 1.38 $M_{\odot}$ results focusing on the difference. For the mass loss parameters of the $1.2~{}M_{\odot}$ model discussed in Section 2.1 we adopt $\log R_{0}/R_{\sun}=-1.25$ (corresponding to the case that we start mass loss at $\log T_{\rm ph}$ (K)= 5.54) and $\log R_{1}/R_{\sun}=-1.36$ (corresponding to the case that we stop mass loss at $\log T_{\rm ph}$ (K)=5.60). The resultant mass loss rates are shown later in Figure 7.

Table 1 summarizes the results of our calculations. The WD mass and mass-accretion rate are given parameters. The others are the results of calculation: the recurrence period $P_{\rm rec}$, temperature at the WD center $T_{\rm core}$, mass accreted after the previous outburst until the beginning of the current outburst (=mass transferred from the companion) $M_{\rm acc}$, ignition mass $M_{\rm ig}$ which is the envelope mass at the epoch of maximum nuclear luminosity $L_{\rm nuc}^{\rm max}$, maximum temperature at the epoch of maximum nuclear luminosity $T_{\rm nuc}^{\rm max}$, nova outburst duration $t_{\rm nova}$ for $L_{\rm ph}>10^{4}~{}L_{\sun}$, duration of the wind phase $t_{\rm wind}$, durations of the X-ray flash $t_{\rm X-flash}$, and SSS phase $t_{\rm SSS}$, rising time $t_{\rm rise}$ which is the time elapsed between the time of the peak nuclear luminosity and the time of the peak mass loss rate, decaying time $t_{\rm decay}=t_{\rm nova}-t_{\rm rise}$. Here, the durations of the X-ray flash and SSS phase are defined as the time span when $L_{\rm X}>10^{4}~{}L_{\sun}$ for 1.38 $M_{\sun}$ followed by Kato et al. (2016), but $L_{\rm X}>10^{3}~{}L_{\sun}$ for 1.2 $M_{\sun}$ because of faint X-ray luminosity, where $L_{\rm X}$ means the luminosity in a band of $0.3-1.0$ keV. Note that the outburst duration $t_{\rm nova}$ is different from the summation of the three phases, $t_{\rm X-flash}+t_{\rm wind}+t_{\rm SSS}$, because the definition is different.

Figure 1 shows the last cycle of our 1.2 $M_{\odot}$ WD model. We stopped the mass accretion ($\dot{M}_{\rm acc}=2.0\times 10^{-7}~{}M_{\odot}$ yr*-1*) when the luminosity increases to $\log L_{\rm ph}/L_{\sun}=4.32$ ($\log T_{\rm ph}$ (K)= 5.85), and resumed mass accretion after the flash when the luminosity decreases to $\log L_{\rm ph}/L_{\sun}=4.0$ ($\log T_{\rm ph}$ (K)$=5.83$). The period of mass accretion is depicted in Figure 1(b). The recurrence period is obtained to be 9.9 yr.

Figure 1(a) shows the photospheric luminosity, $L_{\rm ph}$, nuclear burning energy release rate integrated from the center of the WD to the surface, $L_{\rm nuc}$, and total gravitational energy release rate, $L_{\rm G}$, defined as the integration of $\epsilon_{\rm g}$ from the center of the WD to the photosphere (equation (3) in Kato et al., 2016). The supersoft X-ray (0.3 $-$ 1.0 keV) luminosity $L_{\rm X}$ is calculated assuming the blackbody emission of photospheric temperature $T_{\rm ph}$. The lower panel (b) shows the change of the photospheric temperature $T_{\rm ph}$ and radius $R_{\rm ph}$.

A closeup view in a very early phase is shown in Figure 2. The nuclear energy release rate reaches $L_{\rm nuc}=2.2\times 10^{6}L_{\odot}$ at maximum. This energy is mostly absorbed in the lower part of the burning region as shown by a large negative value of $L_{\rm G}(<0)$. Thus, only a very small part of $L_{\rm nuc}$ is transported outward up to the photosphere. As a result, the photospheric luminosity $L_{\rm ph}$ does not exceed the Eddington luminosity. These properties are the same as those reported in the $1.38~{}M_{\sun}$ model (Figure 2 in Kato et al., 2016). Comparing the two models in Table 1, we find that $L_{\rm nuc}^{\rm max}$ of the 1.2 $M_{\sun}$ model is almost a half of the $1.38~{}M_{\sun}$ model, i.e., the outburst is very weak.

The absorbed energy ($L_{\rm G}<0$) is released in the later phase of the outburst. As shown in Figure 1(a), $L_{\rm G}~{}(>0)$ reaches about 10 % of $L_{\rm ph}$ during the outburst. After the shell flash ends, $L_{\rm nuc}$ quickly decreases, but $L_{\rm ph}$ slowly decreases because the mass-accretion resumes and $L_{\rm G}$ contributes to $L_{\rm ph}$. During the quiescent phase $L_{\rm ph}\approx L_{\rm G}$ is almost constant ($\log L_{\rm ph}/L_{\sun}\gtrsim 2.4$) as shown in Figure 1(a).

3.1. H-R Diagram

Figure 3 shows the track in the H-R diagram for the same period in Figure 1. During the quiescent phase the accreting WD model stays around point A. After the unstable nuclear burning sets in, the bolometric luminosity quickly increases whereas the photospheric radius stays almost constant (from point A to B). The nuclear energy generation rate $L_{\rm nuc}$ reaches its maximum at point B. As the envelope expands, the photospheric temperature decreases and optically thick winds start to blow at point C. The winds continue until it reaches point G. At point F the photospheric temperature attains a minimum. After that the photospheric radius decreases and the temperature $T_{\rm ph}$ increases with time. The envelope mass decreases owing to wind mass loss and hydrogen burning.

Figure 3(b) shows a closeup view of the same track. There is a zigzag part (D-E-F) arising from the change of chemical composition in the surface layers. Just after the thermonuclear runaway sets in, convection develops throughout the envelope up close to the photosphere. In the uppermost layer where the convection never reaches, the chemical composition remains to be the same as that of the accreted matter, i.e., $X=0.7,~{}Y=0.28$ and $Z=0.02$. In the region where the convection reached the hydrogen mass-fraction decreased to $X\sim 0.6$. Figure 4 shows distributions of hydrogen mass-fraction in the envelope at selected stages, in which all the layer with $X=0.7$ have been blown off in the wind a little after stage E. The gradient of $X$ results in the opacity gradient (for larger $X$, opacity is larger owing to electron scattering) that causes additional acceleration of the winds. Thus, the winds blow relatively stronger until the region with $X=0.7$ is completely blown off. As the chemical composition distribution approaches uniform winds become weak and the photospheric temperature rises with time. As the shell flash goes on the envelope expands and the photospheric temperature decreases again. Thus, the track in the H-R diagram has a zigzag excursion before it reaches point F in Figure 3(b). The envelope mass continuously decreases through these stages. Note that the distribution of hydrogen abundance $X$ in Figure 4 is that of the evolution model. In wind solution we assumed X-profile mimic to the corresponding evolution model.

The dotted and dashed lines in Figure 3(b) are the steady-state/static sequences that represent the nova decay phase (Kato, 1999; Kato & Hachisu, 1994). The dashed line is for a $1.2~{}M_{\sun}$ WD with the uniform chemical composition of $X=0.7,~{}Y=0.28$ and $Z=0.02$ while the dotted is for $X=0.6,~{}Y=0.38$ and $Z=0.02$. The WD radius, which is needed in obtaining steady-state models, i.e., the radius at the bottom of the hydrogen-rich envelope, is taken from our evolution calculation as $\log R/R_{\sun}=-2.22$.

In the decay phase (stages F, G, and later) the evolution path is very close to that of the upper steady-state sequence (dotted line). In this phase, the envelope composition is almost uniform at $X=0.6$, $Y=0.38$, and $Z=0.02$, and the envelope structure is close to that of the steady-state (stage F to G) and of hydrostatic balance (after stage G) as well as thermal equilibrium. Thus, the nova evolves along with the steady-state/static sequence.

In the rising phase, on the other hand, the envelope does not yet reach thermal balance but approaches the balance with time. The envelope still has an $X=0.7$ layer on its top, so the evolution path is close to the lower dashed line.

3.2. Envelope structure

Figure 5 shows the envelope structures of the wind solution (solid lines) and Henyey code solution (dotted lines) at the maximum expansion of the photosphere (epoch F). The photosphere reaches $\log R_{\rm ph}$ (cm)=10.61 ($0.58~{}R_{\sun}$). This figure also shows the distribution of the local luminosity $L_{\rm r}$ and the local Eddington luminosity defined by

[TABLE]

where $\kappa$ is the opacity. We adopt OPAL opacity (Iglesias & Rogers, 1996). The local Eddington luminosity has a small dip corresponding to a small peak in the OPAL opacity owing to ionized O and Ne ($\log T$ (K) $\sim 6.2$–6.3), and a large decrease corresponding to a large Fe peak at $\log T$ (K) $\sim 5.2$ (see Figure 6 of Kato et al., 2016). The wind is accelerated where the opacity increases outward. The wind velocity barely exceeds the escape velocity $v_{\rm esc}=\sqrt{2GM_{\rm WD}/r}$.

The dotted lines show the interior structure of the model calculated with the Henyey code in step (1) in Section 2.4. The connecting point with the wind solution is depicted by filled circles. To demonstrate the smooth fitting, we plot the model structure beyond the connecting point. It shows a good agreement with the steady state solution not only on the connecting point but also in the wide outside region.

Figure 6 demonstrates how the optically thick winds ’quietly’ cease at stage G. In the decay phase of the nova outburst, the photospheric temperature increases with time and the dip of $L_{\rm Edd}$ becomes shallower. The solid line depicts the envelope structure just before stage G. The wind is still accelerated, but the acceleration region, where $L_{r}>L_{\rm Edd}$, is getting narrower. This region disappears at stage G (dotted line). The resemblance of the two structures demonstrates that the optically thick winds gradually weaken and stop quietly. This property was already pointed out in the steady-state sequences (Kato & Hachisu, 1994) that the optically thick winds quietly ceases when the photospheric temperature rises to beyond the Fe peak.

3.3. X-ray Light Curve and Wind Mass Loss

In Figure 7 the upper panel (a) shows the bolometric and supersoft X-ray luminosities during the outburst, while the lower panel (b) shows its closeup view with wind mass loss rates. The overall one-cycle view is already shown in Figure 1. The first bright X-ray phase is the X-ray flash that occurs just after the onset of thermonuclear runaway (Kato et al., 2016). The tiny X-ray peak shortly after the X-ray flash corresponds to the small leftward excursion toward stage E in the H-R diagram (Figure 3). The supersoft X-ray phase begins at stage G and ends shortly after stage H.

The panel (b) also shows the temporal change of the wind mass loss rate which is closely related to the change in the X-ray light curve. The X-ray flash ends when the envelope expands and the optically thick winds starts at stage C. As the photospheric radius increases (photospheric temperature decreases) toward stage F, the wind mass loss rate increases with time. After stage F the mass loss rate decreases as the star evolves back from stage F to G in the H-R diagram (Figure 3). The wind stopped at stage G when the SSS phase begins. In this way the X-ray turn on/off is closely related to the occurrence of the optically thick winds. The mass loss rate increases (decreases) as the photospheric temperature decreases (increases) which will be shown later (Figure 12 in Section 5.1).

A sequence of static and steady-state solutions is also shown in Figure 7(a) for comparison. The dotted red line depicts the X-ray light curve calculated from the model shown by the dotted line in Figure 3. The rightmost point corresponds to the epoch of extinction of hydrogen burning. This X-ray light curve shows a good agreement with the evolutional one but the termination point comes earlier. This difference is explained as follows. The evolution time of the equilibrium sequence is calculated from the mass decreasing rate owing to nuclear burning. On the other hand, in the time-dependent calculation, the evolution time is influenced by the gravitational energy release at the bottom of the envelope, i.e., the evolution time is determined by the mass decreasing rate owing to nuclear burning and additional energy source of gravitational energy release rate $L_{\rm G}$. This additional energy source, amounts about $10\%$ of the photospheric luminosity (Figure 1), lengthens the SSS lifetime. Considering this effect we can conclude that the steady-state sequence is consistent with the time-dependent calculation.

4. 1.38 $M_{\sun}$ Model

Our multicycle shell flash calculation of a 1.38 $M_{\odot}$ WD with $\dot{M}_{\rm acc}=1.6\times 10^{-7}~{}M_{\odot}$ yr*-1* was published (Kato et al., 2016), but it included only the early phase of the X-ray flash. Using the method of calculation discussed in Section 2.4 we improved the fitting in the wind mass loss phase. We adopt $\log R_{0}/R_{\sun}=-1.00$ and $\log R_{1}/R_{\sun}=-1.15$. The recurrence period is 0.95 yr. The characteristic properties of this model are summarized in Table 1. Here, we report our results focusing on the difference from the 1.2 $M_{\odot}$ case.

Figure 8 shows the last one cycle of our calculation. The outburst duration ($L_{\rm ph}>10^{4}~{}L_{\sun}$) is 29 days. In the interpulse phase, the gravitational energy release rate $L_{\rm G}$ dominates the photospheric luminosity $L_{\rm ph}$. The energy budget in the early phase, corresponding to Figure 2, is already published in Kato et al. (2016). These characteristic properties are same as those in the 1.2 $M_{\odot}$ WD model.

Figure 9(a) shows one cycle of nova outbursts in the H-R diagram. A closeup view is shown in panel (b). Similarly to the 1.2 $M_{\odot}$ model (Figure 3), stages C-D-E correspond to the period that the uppermost layer with $X=0.7$, of mass $1.6\times 10^{-9}~{}M_{\sun}$, are blown off in the wind. After that, the chemical composition in the envelope is almost uniform at $X=0.58$, $Y=0.4$ and $Z=0.02$, and the envelope evolves along with the line of steady-state/static sequence.

Figure 10 shows the internal structure at stage F. The velocity quickly rises (around $\log r$ (cm) $\sim$10.0) where the local Eddington luminosity quickly decreases corresponding to the Fe opacity peak. The envelope mass at this stage is as small as $1.1\times 10^{-7}~{}M_{\sun}$ which is not enough to expand to a giant size, so the photospheric radius is as small as $\log R_{\rm ph}$ (cm) =10.52 ($0.48~{}R_{\sun}$), and the surface temperature is rather high at $\log T_{\rm ph}$ (K)=5.10.

To demonstrate a good matching between the interior and wind solutions, the interior structure outside the fitting point is also plotted; this outer part was replaced with the wind solution in step (4) in Section 2.4. The dotted line is very close to the steady-state wind solution, indicating that the interior structure is smoothly connected to the wind solution.

Figure 11(a) shows the temporal changes in $L_{\rm ph}$, $L_{\rm nuc}$, and $L_{\rm G}$ as well as the change of $T_{\rm ph}$. The supersoft X-ray light curve shown in panel (b) is calculated from $L_{\rm ph}$ and $T_{\rm ph}$. The X-ray flash (when $L_{\rm X}>10^{4}~{}L_{\sun}$) lasts 0.77 days and the SSS phase 7.8 days. Between them there is a faint X-ray peak corresponding to stage E.

This model may be compared with the observational properties of the 1 yr recurrence period nova M31N 2008-12a. Kato et al. (2015) estimated the WD mass of M31N 2008-12a to be $1.38~{}M_{\sun}$ from the SSS duration that lasts about 8 days (Henze et al., 2014, 2015; Darnley et al., 2016). In the present paper, we do not consider the optical light curves because the relation between the optical brightness and wind mass loss rate in recurrent novae is uncertain, which remains to be a future work. If we assume that the optical magnitude peaks at stage F, i.e., when the wind mass loss rate is maximum, the X-ray turn-on time is 4 days after the optical maximum, which is roughly consistent with M31N 2008-12a ($5.9\pm 0.5$ days: Henze et al., 2015), considering the insufficient fitting quality of our $1.38~{}M_{\sun}$ model. Detailed comparison with the observation is beyond the scope of the present paper.

5. DISCUSSION

5.1. Wind Mass Loss Rate

Figure 12 shows the mass loss rate of the optically thick winds calculated in the fitting process (Section 2.4) against the photospheric temperature $T_{\rm ph}$. The mass loss rate increases as $T_{\rm ph}$ decreases. In the extended envelopes the acceleration region (where the opacity increases outward) locates deep below the photosphere where the density is large (see Figures 5 and 6), which results in a large mass loss rate. The mass loss rate also weakly depends on the WD mass. The $1.38~{}M_{\sun}$ model shows a slightly higher mass loss rates because the surface gravity is stronger in more massive WDs.

This figure also shows the two $\dot{M}_{\rm wind}$ - $T_{\rm ph}$ relations (solid lines) obtained from the steady-state sequences of 1.2 and 1.38 $M_{\sun}$, corresponding to the upper sequence in each H-R diagram (Figure 3(b) for 1.2 $M_{\sun}$ and Figure 9(b) for 1.38 $M_{\sun}$). The dependence on the chemical composition is very weak, so the solar composition sequences are omitted. The mass loss rates obtained from our time-dependent calculations agree very well with those of steady-state sequences throughout the wind phase, including the zigzag path, C-D-E-F-G.

This agreement of the $\dot{M}_{\rm wind}$ - $T_{\rm ph}$ relations is reasonable because the interior structure of our time-dependent calculation shows a good agreement with the steady state solution. The envelope is accelerated where the opacity rapidly increases outward at $T\sim 2\times 10^{5}$ K, so that similar envelope structures result in similar photospheric temperatures and wind mass-loss rates.

5.2. Details of Numerical Calculation

In Section 2.1 we described the mass loss rate $\dot{M}_{\rm evol}$ which is initially assumed. After several time-consuming iterations, we found that better results are obtained if we choose $R_{0}$ to be close to those where the optically thick wind stops in the steady-state sequence (Kato & Hachisu, 1994), i.e., when the photospheric temperature is close to that where $\dot{M}_{\rm wind}=0$ in Figure 12.

Figure 7(b) shows $\dot{M}_{\rm evol}$ (dotted purple line) and resultant wind mass-loss rate (solid red line) in the $1.2~{}M_{\sun}$ model. The resultant wind mass loss rates are somewhat larger than the assumed value ($19~{}\%$ at the time of maximum wind mass loss) but traces well the global change including early small peak in the episode of surface composition change (see Section 3.1). This agreement is satisfactorily well, considering numerical difficulties in step (1).

Figure 11(b) shows $1.38~{}M_{\sun}$ models for two sets of $\dot{M}_{\rm evol}$ (dotted lines) and corresponding resultant (solid lines) mass loss rates $\dot{M}_{\rm wind}$. Blue dotted line show $\dot{M}_{\rm evol}$ obtained with $\log R_{0}/R_{\sun}=-1.0$ and $\log R_{1}/R_{\sun}=-1.15$, which results in the wind mass loss rate shown by the solid red line. The $25\%$ matching is obtained at the maximum mass loss rate (of the red line) but there is a short period of a large discrepancy of factor 2 just before the winds stop.

Dotted green line shows $\dot{M}_{\rm evol}$ obtained adopting $\log R_{0}/R_{\sun}=-1.5$ and $\log R_{1}/R_{\sun}=-1.3$. In this case $\dot{M}_{\rm evol}$ has no large enhancement in the final wind phase, but a sharp narrow peak in the very early stage. The resultant wind mass loss rates, $\dot{M}_{\rm wind}$, denoted by the solid green line, however, are much small in the later phase. To summarize, with a smaller $\dot{M}_{\rm evol}$ (dotted blue line) we get a slightly larger $\dot{M}_{\rm wind}$ except the final peak, and with a larger $\dot{M}_{\rm evol}$ (dotted green line), we get a much smaller $\dot{M}_{\rm wind}$. We suppose that a well converged model should have wind mass loss rates between them, which may be closer to the red line model than the solid green line. Therefore, we adopted the red line model.

5.3. Comparison with other works

In the present work we calculated nova evolution models with very fine mass zonings and short timesteps. As a results, we are able to follow the temporal change of the wind mass loss, including the gradual increase and decrease, as in Figures 7 and 11. We also showed how the upper envelope with $X=0.7$ is blown off (Figure 4 and zigzag tracks in Figures 3 and 9).

Many time-dependent calculations on novae have been presented so far, but none of them showed a zigzag behavior in the H-R diagram. For example, Prialnik & Kovetz (1995) calculated nova multicycle outbursts for various WD masses, of which light curves are recently published by Hillman et al. (2014). Their light curves show no zigzag behavior like ours. The difference might arise from the difference in zonings of the envelope, because if the mass zoning did not resolve well the superficial $X=0.7$ region, the model would be unable to follow the temporal change of the chemical composition of the superficial layers. This reasoning may be supported by the description “convection had already reached the surface of the star” (Section 2.1 Prialnik & Kovetz, 1995). In addition, in their nova outbursts, the photospheric temperature suddenly jumps to the minimum value, and stay there, and jumps back. This corresponds to, in our Figure 3, that the star jumps from stage C to F instantaneously and jumps back to stage G. Such a behavior could occur if large timesteps were adopted, skipping the gradual increase/decrease periods of the mass loss rates.

Kovetz (1998) presented an algorithm to place an optically thick wind solution on top of a hydrostatic stellar configuration. The calculation shows no indication of the zigzag shape. This is because the adopted small number of mass zones (less than 50) is not sufficient to treat the very thin $X=0.7$ layer. Despite the detailed description of the algorithm, no information on the connecting point is given. In the 1.0 $M_{\sun}$ model, the wind begins at $\log T_{\rm eff}$ (K)=5.36 and ceases at $\log T_{\rm eff}$ (K)=5.5, which is roughly consistent with our model with the OPAL opacity. At the occurrence of the wind the mass loss rate suddenly jumps from zero to the maximum value followed by a gradual decrease. The gradual decrease is consistent with our calculation.

Idan et al. (2013) adopted Kovetz’ (1998) mass loss scheme with sufficient number of mass zones ($>8000$), however, the paper discusses interior structures in detail with little description of the surface region.

In the previous paper (Kato et al., 2015) we presented a multi-wavelength light curve model of $1.38~{}M_{\odot}$ WD for the 1 yr recurrence period nova M31N 2008-12a. We adopted a similar method of numerical calculation to the present work, but connected the interior structure at a much deeper place with the wind solutions. One of the reason why we chose such fitting point was that our numerical calculation with the Henyey code did not converge unless large mass loss rates were assumed, and the fitting place was chosen for the structure to be almost independent of the assumed mass loss rate. However, in such a deep envelope with a larger mass loss rate, the gravitational energy release rate per unit mass $\epsilon_{\rm g}$ is not negligible, so the fitting of steady-state solutions with the interior yielded mediocre models. The resultant mass loss rates were about 3 times larger than the present work. Thus, the wind phase was shorter by a factor of $\sim 3$. This roughly explains the difference between the previous estimate of a duration of 5 days and our present estimate of 19 days.

In the present work we have improved the numerical technique and adopt much thinner atmosphere and mass zoning, and connect to the interior structure at a place where $\epsilon_{\rm g}$ is negligible. As a result, matching of inner and outer solutions are sufficiently well as shown in Figures 5 and 10. Our $1.2~{}M_{\sun}$ model shows very good convergence of iteration process (Section 2.4) throughout the wind phase that support our fitting procedure. Our $1.38~{}M_{\sun}$ model shows a better convergence except a short period just before the SSS phase. This model must be improved further in the future.

6. CONCLUSIONS

Our main results are summarized as follows.

1. We present a new calculation procedure to follow nova outbursts in which the optically thick wind is consistently introduced. We have found that steady-state wind solutions smoothly match with the interior calculated with the Henyey-type code.

2. We calculated one cycle nova outburst on 1.2 and 1.38 $M_{\sun}$ WDs with the mass-accretion rates corresponding to the recurrence periods of $P_{\rm rec}=10$ and 1 yr, respectively. The wind mass loss rate is accurately obtained in the new calculation procedure. The X-ray flash ends when the optically thick winds begin. The SSS phase starts when the winds cease. At the maximum expansion of the photosphere, its temperature is as high as $10^{5}$ K in the both cases.

3. We obtained a zigzag track in the H-R diagram corresponding to the epoch that the very surface zone of $X=0.7$ is blown off in the wind, uncovering the region mixed by convection during the runaway, which therefore has a lower H fraction correspondingly different opacity.

4. Our $1.38~{}M_{\sun}$ model is roughly consistent with the timescale of the 1 yr recurrence period nova M31N 2008-12a.

5. The relation between the wind mass loss rate and photospheric temperature in our evolution model agrees well with those of steady-state/static sequence. The envelope structure at each stage are also very similar to those of steady-state/static sequence in the wide region except the deep interior where the gravitational energy release rate per unit mass $\epsilon_{\rm g}$ is not negligible.

MK and IH express our gratitude to T. Iijima and Astronomical Observatory of Padova (Asiago) for warm hospitality during which we initiated the present work. We also thank the anonymous referee for useful comments that improved the manuscript. This research has been supported in part by Grants-in-Aid for Scientific Research (15K05026, 16K05289) of the Japan Society for the Promotion of Science.