A&A 484, L9-L12 (2008)
DOI: 10.1051/0004-6361:200809703

LETTER TO THE EDITOR

3D simulations of RS Ophiuchi: from accretion to nova blast^,

R. Walder^1,2,3 - D. Folini² - S. N. Shore^3,4

1 - Observatoire de Genève, Université de Genève, 51 ch. des Maillettes, 1290 Sauvernay, Switzerland
2 - Institute for Astronomy, ETH Zürich, ETH-Zentrum, SEC D2, 8092 Zürich, Switzerland
3 - Dipartimento di Fisica ``Enrico Fermi'', Università di Pisa, largo Pontecorvo 3, Pisa 56127, Italy
4 - INFN/Pisa, largo Pontecorvo 3, Pisa 56127, Italy

Received 3 March 2008 / Accepted 15 April 2008

Abstract
Context. The binary star system RS Ophiuchi is a recurrent nova, with outbursts occurring about every 22 years. It consists of a red giant star (RG) and a wind accreting white dwarf close to the Chandrasekhar limit. This system is considered a prime candidate for evolving into an SNIa. For its most recent outbursts in 1985 and 2006, exquisite multiwavelength observational data are available.
Aims. Deeper physical insight is needed regarding the inter-outburst accretion phase and the dynamical effects of the subsequent nova explosion in order to improve the interpretation of the observed data and to shed light on whether the system is an SNIa progenitor.
Methods. We present a 3D hydrodynamic simulation of the quiescent accretion with the subsequent explosive phase.
Results. The computed circumstellar mass distribution in the quiescent phase is highly structured with a mass enhancement in the orbital plane of about a factor of 2 as compared to the poleward directions. The simulated nova remnant evolves aspherically, propagating faster toward the poles. The shock velocities derived from the simulations agree with those derived from observations. For $v_{\rm RG} = 20$ km s^-1 and for nearly isothermal flows, we find that 10% of the mass lost by the RG is transfered to the WD. For an RG mass loss of $10^{-7}~M_{\odot}$ yr^-1, the orbit of the system decays by 3% per million years. With the derived mass transfer rate, multi-cycle nova models provide a qualitatively correct recurrence time, amplitude, and fastness of the nova.
Conclusions. Our 3D hydrodynamic simulations provide, along with the observations and nova models, the third ingredient for a deeper understanding of the recurrent novae of the RS Oph type. In combination with recent multi-cycle nova models, our results suggest that the WD in RS Oph will increase in mass. Several speculative outcomes then seem plausible. The WD may reach the Chandrasekhar limit and explode as an SN Ia. Alternatively, the mass loss of the RG could result in a smaller Roch volume, a common envelope phase, and a narrow WD + WD system. Angular momentum loss due to gravitational wave emission could trigger the merger of the two WDs and - perhaps - an SN Ia via the double degenerate scenario.

Key words: stars: binaries: symbiotic - stars: novae, cataclysmic variables - accretion, accretion disks - hydrodynamics - methods: numerical - stars: individual: RS Oph

1 Introduction

Type Ia supernovae (SNe Ia) are cornerstones of modern cosmology as a measure for cosmological distances, and they are crucial building blocks of the universe, as the production sites of a large part of the iron group elements. The recurrent nova RS Ophiuchi (RS Oph) (Hachisu & Kato 2001; Hachisu et al. 1999) is a candidate for the still mysterious progenitors of SN Ia. It is a symbiotic-like binary star system with a red giant (RG) and a white dwarf (WD) that undergoes a nova outburst about every 22 years (Fekel et al. 2000; Dobrzycka et al. 1996; Anupama 2002). For the most recent outbursts in 1985 and 2006, exquisite panchromatic observational data are available (Worters et al. 2007; O'Brien et al. 2006; Hachisu et al. 2007; Bode et al. 2006; Sokoloski et al. 2006; Bode et al. 2007; Shore et al. 1996; Evans et al. 2006; Das et al. 2006). Here we present first 3D hydrodynamical simulations of RS Oph, from pre-outburst accretion through nova explosion, with the goal of improving our insight into the physics of the system and interpreting the unique observational data.

RS Oph has an orbital period of 455 days (Fekel et al. 2000), RG and WD masses of 2.3 $M_{\odot}$ and close to 1.4 $M_{\odot}$ (Worters et al. 2007; Anupama 2002), respectively, and a separation between the components of a = 2.68 $\times$ 10¹³ cm. The RG mass is, however, still uncertain, much lower values (Fekel et al. 2000) have also been suggested. We supplement the system parameters by assuming a terminal wind velocity of $v_{\rm RG} = 20$ km s^-1 in the rest frame of the RG and a mass loss rate of $\dot{M}_{\rm RG} = 10^{-7}$ $M_{\odot}$ yr^-1. Values for $\dot{M}_{\rm RG}$ are still not well known. While Seaquist & Taylor (1990) find $\dot{M}_{\rm RG} \sim 10^{-7}~M_{\odot}$ yr^-1 for the majority of RG in symbiotic systems, the scatter is significant and values of less than $10^{-8}~M_{\odot}$ yr^-1 have been found. For the range of spectral types suggested (Worters et al. 2007) for RGs in the RS Oph system, its radius is always smaller than its Roche lobe, and accretion by the WD occurs only from the RG wind.

	Figure 1: The grid structure of the nine levels of refinement, shown for the entire computational domain ( left) and around the accreting WD ( right). While the coarser grids of levels 1 through 6 are fixed in space, the meshes of levels 7 through 9 follow the WD. The RG is shown in red, the accreting WD in blue.
Open with DEXTER

2 3D simulations of RS Oph

2.1 Computational method

Our numerical simulation was performed with the A-MAZE-code (Walder & Folini 2000), a parallel, block-structured, adaptive mesh refinement (AMR) hydrodynamical code using Cartesian meshes and a multidimensional high-resolution finite-volume integration scheme. This code has been used for accretion studies before (Walder 1997; Zarinelli et al. 1995; Harper et al. 2005; Dumm et al. 2000). The Euler equations are solved in three dimensions with a nearly isothermal polytropic equation of state with $\gamma =1.1$, resulting in a thermal structure that comes close to what is obtained from photoionization models of related symbiotic binary systems (Nussbaumer & Vogel 1987; Nussbaumer & Walder 1993).

The simulation was carried out in a Eulerian frame of reference with the stars moving within the computational domain which measured 10¹⁵ cm a side. The computational grid consisted of nine nested levels of refinement (Fig. 1). The coarsest level consisted of 50 cells cubed. From one level to the next, grid cells were refined by a factor of two. Each level comprised between 8 and 64 individual grids, and the entire mesh consisted of 233 grids and 2 $\times$ 10⁷ cells. The decomposed grid structure was exploited for parallelization under OpenMP.

We imposed free outflow at the outer domain boundary. The WD was represented as a spherical low-pressure, low-density, zero-velocity region of radius R= 2.2 $\times$ 10¹¹ cm. The accretion rates of mass, momentum, energy, and angular momentum onto the WD were derived from the flow quantities entering the region. The rates were not sensitive to the size, pressure, and density of the region representing the WD. The RG was represented as a spherical region of 150 $R_{\odot}$. In the rest frame of the RG star, which rotated and orbited with respect to the computational domain, we chosen a terminal velocity of the RG wind of 20 km s^-1. The absolute value of the density in the cells representing the RG was adjusted to achieve the desired mass loss of $10^{-7}~M_{\odot}$ yr^-1. The wind temperature was set to 8000 K, a value between the RG photospheric temperature and the wind temperature closer to the WD as computed by elaborate photo-ionization models (Nussbaumer & Vogel 1987).

At the beginning, the entire computational domain was filled with the RG wind. The accreting system was then relaxed over seven orbital periods, about the crossing time from the RG to the domain boundary at the RG wind speed. Mass and angular momentum losses out of the system were then computed through spherical shells around the center of mass.

The losses became constant for shell radii larger than 10¹⁴ cm. The nova explosion was then launched into the relaxed density structure.

Figure 2: Density structure during the accretion phase. Shown is density (logarithmic scale, g/cm3) in the orbital plane (xy-plane, left), and in a plane perpendicular to the orbit (yz-plane, right) for the entire computational domain (1015 cm a side, top) and a zoom around the accreting WD ($\approx$1013 cm a side, bottom). Self-interacting high-density spirals dominate the inner region up to a few times the binary separation. The RG is shown in red, the WD in blue.

Open with DEXTER

2.2 Quiescent phase: density spirals and wind accretion

The inter-outburst accretion phase leads to a circumstellar density distribution that deviates substantially from the 1/r² density distribution of a single RG wind out to several system separations. The deviations, usually neglected when interpreting RS Oph observations, arise from the orbital motion of the two components and the accretion onto the WD, and are subsequently transported outward in the flow. The resulting density structure, Fig. 2, depends critically on the ratio $R = v_{\rm RG}/v_{\rm orb}$, where $v_{\rm orb} =2 \pi a / P$, with P the orbital period. For $R \gg 1$, an Archimedian spiral occurs in the orbital plane, as observed in the colliding-wind Wolf Rayet binary star system WR 104 (Tuthill et al. 1999). For RS Oph, with $R \le 1$, the spiral pattern becomes self-interacting.

In the vicinity of the WD, a high-density accretion disk forms that is clearly visible in Fig. 2. Its diameter is roughly 1/10 of the system separation, and velocities are non-Keplerian. Strong spiral shock waves are responsible for the transport of angular momentum and mass, and stable accretion occurs. The disk is geometrically thick, and the opening angle perpendicular to the orbital plane (Fig. 2, lower right panel) measures roughly $\pm$ $60^{\circ }$ in the direction of the RG and $\pm$ $30^{\circ}$ away from it.

On length scales less than the system separation, radial density profiles as seen from the WD (Fig. 3) show density contrasts between the orbital plane and poleward directions of one to two orders of magnitude. Most prominent here, and reaching farthest out, is the trailing accretion wake of the WD. This density structure affects the early evolution of the nova remnant. At greater distances, densities in the orbital plane remain higher than toward the poles, but only by a factor of 2-3. This is, however, enough to cause an asymmetric evolution of the nova remnant. Although the density decreases as 1/r² on average, the relative amplitude of local variations due to the self-interacting spiral pattern are up to a factor of 10, especially in the orbital plane. We expect these local variations to vanish at even larger scales that are beyond our computational domain.

The velocity of the systemic outflow is a superposition of the RG wind velocity, its rotational velocity, the orbital velocity, and hydrodynamical effects. On larger scales than $\approx$5 $\times$ 10¹³ cm, the flow field is essentially directed radially outward. Depending on the exact angle of an observer, any radial speed between 20 km s^-1 (polewards) and $\approx$50 km s^-1 (mostly in the orbital plane) can be found in the computational data. This range is consistent with the different values of the wind speeds derived from observations, e.g. 36 km s^-1 by Iijima (2008) and 40-60 km s^-1by Shore et al. (1996).

Figure 3: Density profiles along radial rays originating from the WD. Bold colored lines denote rays in the orbital plane along the x- and y-coordinate axes, the RG being located in the negative x-direction (-x green; +x magenta; -y blue; +y red). Perpendicular, green bold lines denote the position of the RG. The bold black line is in a poleward direction ( $+90^{\circ }$), thin grey lines belong to an inclination of $60^{\circ }$ and positive and negative x- and y-directions, respectively.

Open with DEXTER

The computed accretion rate is 10% of the mass loss rate of the RG, independent of the absolute mass loss of the RG. To our knowledge, this is the first quantitative self-consistent prediction of the accretion rate of a symbiotic-like recurrent nova. A higher accretion rate is expected for a less massive RG (Fekel et al. 2000; Dobrzycka & Kenyon 1994), since a smaller system separation would result. The computed value is in line with the values found in 3D simulations of related symbiotic binary star systems, for example, 6% in RW Hydrae (Dumm et al. 2000) and 10%-25% in Z Andromeda (Mitsumoto et al. 2005). In absolute terms, the WD in our simulation accretes $10^{-8}~M_{\odot}$ yr^-1. This value also agrees well with the accretion rates required by multi-cycle nova evolution models (Yaron et al. 2005) for RS Oph like binary systems, with $M_{\rm WD} = 1.4~M_{\odot}$ and a recurrence period of 22 years.

Conversely, the nova models, together with the above accretion rate of 10% of the RG mass loss, quantitatively constrain the mass loss of the RG to values around $10^{-7}~M_{\odot}$ yr^-1. This is yet another, completely independent, estimate of the mass loss rate of the RG in RS Oph. It is more at the upper limit of the values given by Seaquist & Taylor (1990) for RG in classical symbiotics. A much higher estimate of $10^{-5}~M_{\odot}$ yr^-1 was derived by Shore et al. (1996). A more recent analysis indicates the rate could be a factor of about 10 lower (Shore 2008) and should only be taken as an upper limit. However, evolutionary models of single stars predict values greater than $10^{-7}~M_{\odot}$ yr^-1 for only a very short time on the RG branch (Maeder & Meynet 1988). At the moment, we are not able to resolve this spread in the estimates.

	Figure 4: Density structure of the nova remnant. Shown are density (logarithmic scale, g/cm3) and velocity (cm/s) in the orbital plane (xy-plane, left) and in a plane perpendicular to the orbit (yz-plane, right) at 29 h ( top panels) and 21 days ( bottom panels) after explosion. The RG is shown in red, the WD in blue.
Open with DEXTER

2.3 Outburst: bipolar blast due to density stratification

Launching the ejecta of a nova explosion into this complex circumstellar environment naturally results in a strongly asymmetric shock front, visible as a high-density shell in Fig. 4. The time evolution of the nova remnant is shown in a supplementary movie. The explosion was simulated by instant injection with energy and mass of 2.2 $\times$ 10⁴³ erg and 2 $\times$ $10^{-7}~M_{\odot}$, respectively, at the position of the WD. The assumed explosion mass is on the lower end of published estimates for the ejected mass (O'Brien et al. 1992; Sokoloski et al. 2006), which range from $10^{-7}~M_{\odot}$ to $10^{-6}~M_{\odot}$, but is comparable to the 2.2 $\times$ $10^{-7}~M_{\odot}$ obtained from the self-consistently computed accretion phase. An isothermal equation of state with $\gamma =1.1$ is used, because inclusion of detailed cooling as in recent 1D models (Vaytet et al. 2007) is currently beyond the reach in 3D simulations.

The initial velocities of the modeled ejecta are $\sim$3500 km s^-1, within the range of observed values (Evans et al. 2006; Bode et al. 2006; Das et al. 2006). A larger extension of the shock front in poleward directions, where pre-outburst densities are lowest, is consistent with observations (O'Brien et al. 2006; Chesneau et al. 2007). In the orbital plane, the shock front displays a roughly circular shape for about 2/3 of its arc. Toward the RG and the former accretion wake, where densities are particularly high, the shock front advances more slowly.

Quantitative evaluation of the shock position as a function of time, Fig. 5 yields shock velocities scaling with time as $v_{\rm s} \propto t^{-\alpha}$, with $0.2 \leq \alpha \leq0.5$for most times and directions. Poleward values of $\alpha$ are slightly lower after about day 10, while $\alpha = 0.5$ for directions toward the former accretion wake of the WD and times before day 6 after the explosion. Values for $\alpha$ derived from observational data of the 2006 nova outburst (Das et al. 2006; Bode et al. 2006; Sokoloski et al. 2006) are in the same range up to day 6, and are somewhat higher ( $0.45 \le \alpha \le 0.79$) later on. The present simulation suggests that the observed range of values is real and a consequence of different observational diagnostics probing different regions of the expanding shock front. A clear change in the velocity of the shock as a function of time is seen only toward the former accretion wake.

There is no evident change in the expansion when the accumulated mass equals the ejected mass of 2 $\times$ $10^{-7}~M_{\odot}$, which occurs after about eight to ten days. Observations taken much later will show only mixed ejecta with the nova processed material significantly diluted.

Figure 5: Time evolution of outer shock position relative to the WD. Bold colored lines belong to directions within the orbital plane (green: toward the RG or $+90^{\circ }$; blue: toward the trailing spiral of the WD or $0^{\circ }$; magenta: away from the RG or $-90^{\circ }$; red: in direction of the orbital motion or $-180^{\circ }$). The bold black line denotes the poleward direction. Thin lines belong to directions within the orbital plane, into the trailing accretion wake of the WD ( $50^{\circ }$ to $60^{\circ }$, solid black to gray lines) and out of the system ( $-35^{\circ }$, dashed orange line).

Open with DEXTER

2.4 Shrinking of the binary orbit

The orbit shrinks at a rate of $\dot{a}/a \approx -2.8$ $\times$ 10^-8 per year, or 3% per million years. This finding is in line with previous suggestions (Hachisu et al. 1999; Jahanara et al. 2005). The absolute losses of mass and angular momentum from the system, $\dot{M}_{\rm S}$ and $\dot{J}_{\rm S}$, scale linearly with the RG mass loss. For an RG mass loss rate of $10^{-7}~M_{\odot}$per year, the fractional losses per year are $\dot{M}_{\rm S} / M_{\rm S} = -2.2$ $\times$ 10^-8 and $\dot{J}_{\rm S} / J_{\rm S} = -4.4$ $\times$ 10^-8, respectively, implying a dimensionless specific systemic angular momentum loss of $(\dot{J}_{\rm S} / J_{\rm S}) / (\dot{M}_{\rm S} / M_{\rm S}) = 2$. This time scale does not become significantly larger if the effect of the nova blast on orbital parameters is taken into account, provided not much more mass is ejected than was accreted during the 22 years.

3 Summary and conclusions

Our 3D simulations of the quiescent phase result in a density distribution of the circumstellar matter that is strongly stratified in poleward directions, leading to an expansion of the nova remnant that is roughly twice as fast in the poleward directions than within the orbital plane. With a reasonable RG mass loss rate of about $10^{-7}~M_{\odot}$ yr^-1, the computed accretion rate of about 10% of the RG mass loss rate results in an accretion rate of $10^{-8}~M_{\odot}$ yr^-1 for the WD. Inserting this value into the multi-cycle nova models by Yaron et al. (2005) correctly predicts on a qualitative level - though not in all quantitative details - the recurrence time, amplitude, and fastness of the RS Oph nova. Our simulations further predict shrinking of the binary orbit at a rate of about 3% per million years, which would improve conditions for enhanced mass transfer. According to nova models by Yaron et al. (2005), enhanced mass loss should lead to shorter nova cycles and favor net mass gain over multiple nova cycles, a necessary condition for the evolution toward an SN Ia. On the other hand, the above RG mass loss suggests that the mass and thus the Roche volume of the RG will drastically shrink on a time scale of $\approx$10⁶ years, depending on the still not fixed mass of the star. The system then would likely enter a common envelope phase and produce a narrow WD-WD system. Further angular momentum losses by gravitational wave emission will drive the merger of the two WDs, an event that should be detectable with present GW detectors, e.g. Virgo and LIGO. Depending on what the result of the merger is, an accretion induced collapse or an SN Ia will end the life of the system RS Oph.

Acknowledgements

We thank Dr. J. Favre from the Swiss National Supercomputing Center (CSCS) for graphics support and production of the movie accompanying this paper. Computations were performed at CSCS, at CINECA (Italy), and at the computing center of ETH Zurich. The authors would like thank the staff of all supercomputing centers for substantial support. R.W. acknowledges financial support from the EGO Fellowship (VESF) and the INFN-Sezione Pisa, as well as the hospitality of the Max Planck Institute for Astrophysics, Garching, Germany, where part of the work was done. S.N.S. thanks J. José, M. Bode, A. Evans, and S. Starrfield for valuable discussions.

References

• Anupama, G. C. 2002, in Classical Nova Explosions, ed. M. Hernanz, & J. José, in Amer. Inst. Phys. Conf. Ser., 637, 32
• Bode, M. F., O'Brien, T. J., Osborne, J. P., et al. 2006, ApJ, 652, 629 [NASA ADS] [CrossRef]
• Bode, M. F., Harman, D. J., O'Brien, T. J., et al. 2007, ApJ, 665, L63 [NASA ADS] [CrossRef]
• Chesneau, O., Nardetto, N., Millour, F., et al. 2007, A&A, 464, 119 [NASA ADS] [CrossRef] [EDP Sciences]
• Das, R., Banerjee, D. P. K., & Ashok, N. M. 2006, ApJ, 653, L141 [NASA ADS] [CrossRef]
• Dobrzycka, D., & Kenyon, S. J. 1994, AJ, 108, 2259 [NASA ADS] [CrossRef]
• Dobrzycka, D., Kenyon, S. J., Proga, D., Mikolajewska, J., & Wade, R. A. 1996, AJ, 111, 2090 [NASA ADS] [CrossRef]
• Dumm, T., Folini, D., Nussbaumer, H., et al. 2000, A&A, 354, 1014 [NASA ADS]
• Evans, A., Tyne, V. H., van Loon, J. T., et al. 2006, MNRAS, 373, L75 [NASA ADS]
• Fekel, F. C., Joyce, R. R., Hinkle, K. H., & Skrutskie, M. F. 2000, AJ, 119, 1375 [NASA ADS] [CrossRef]
• Hachisu, I., & Kato, M. 2001, ApJ, 558, 323 [NASA ADS] [CrossRef]
• Hachisu, I., Kato, M., & Luna, G. J. M. 2007, ApJ, 659, L153 [NASA ADS] [CrossRef]
• Hachisu, I., Kato, M., & Nomoto, K. 1999, ApJ, 522, 487 [NASA ADS] [CrossRef]
• Harper, G. M., Brown, A., Bennett, P. D., et al. 2005, AJ, 129, 1018 [NASA ADS] [CrossRef]
• Iijima, T. 2008, in Proceedings of the workshop RS Ophiuchi (2006), ed. N. Evans, M. Bode, & T. O'Brien, preprint, http://www.astro.keele.ac.uk/rsoph/index.html (In the text)
• Jahanara, B., Mitsumoto, M., Oka, K., et al. 2005, A&A, 441, 589 [NASA ADS] [CrossRef] [EDP Sciences]
• Maeder, A., & Meynet, G. 1988, A&AS, 76, 411 [NASA ADS] (In the text)
• Mitsumoto, M., Jahanara, B., Matsuda, T., et al. 2005, Astron. Rep., 49, 884 [NASA ADS] [CrossRef] (In the text)
• Nussbaumer, H., & Vogel, M. 1987, A&A, 182, 51 [NASA ADS]
• Nussbaumer, H., & Walder, R. 1993, A&A, 278, 209 [NASA ADS]
• O'Brien, T. J., Bode, M. F., & Kahn, F. D. 1992, MNRAS, 255, 683 [NASA ADS]
• O'Brien, T. J., Bode, M. F., Porcas, R. W., et al. 2006, Nature, 442, 279 [NASA ADS] [CrossRef]
• Seaquist, E. R., & Taylor, A. R. 1990, ApJ, 349, 313 [NASA ADS] [CrossRef] (In the text)
• Shore, S. N. 2008, in RS Ophiuchi (2006), ed. N. Evans, M. Bode, & T. O'Brien, ASP Conf. Ser., in press (In the text)
• Shore, S. N., Kenyon, S. J., Starrfield, S., & Sonneborn, G. 1996, ApJ, 456, 717 [NASA ADS] [CrossRef]
• Sokoloski, J. L., Luna, G. J. M., Mukai, K., & Kenyon, S. J. 2006, Nature, 442, 276 [NASA ADS] [CrossRef]
• Tuthill, P. G., Monnier, J. D., & Danchi, W. C. 1999, Nature, 398, 487 [NASA ADS] [CrossRef] (In the text)
• Vaytet, N. M. H., O'Brien, T. J., & Bode, M. F. 2007, ApJ, 665, 654 [NASA ADS] [CrossRef] (In the text)
• Walder, R. 1997, in IAU Colloq., 163, Accretion Phenomena and Related Outflows, ed. D. T. Wickramasinghe, G. V. Bicknell, & L. Ferrario, ASP Conf. Ser., 121, 822
• Walder, R., & Folini, D. 2000, in Thermal and Ionization Aspects of Flows from Hot Stars, ed. H. Lamers, & A. Sapar, ASP Conf. Ser., 204, 281 (In the text)
• Worters, H. L., Eyres, S. P. S., Bromage, G. E., & Osborne, J. P. 2007, MNRAS, 379, 1557 [NASA ADS] [CrossRef]
• Yaron, O., Prialnik, D., Shara, M. M., & Kovetz, A. 2005, ApJ, 623, 398 [NASA ADS] [CrossRef] (In the text)
• Zarinelli, A., Walder, R., & Nussbaumer, H. 1995, A&A, 301, 922 [NASA ADS]

 

  
Online Material

RS_Oph.Explosion_Density.mov