A&A 484, L9-L12 (2008)
DOI: 10.1051/0004-6361:200809703
LETTER TO THE EDITOR
R. Walder1,2,3 - D. Folini2 - S. N. Shore3,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
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
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
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
and close to 1.4
(Worters et al. 2007; Anupama 2002), respectively, and a separation between the components of a = 2.68
1013 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
km s-1 in the rest frame of the RG and a mass loss
rate of
yr-1.
Values for
are still not well known.
While Seaquist & Taylor (1990) find
yr-1 for the majority of RG in symbiotic
systems, the scatter is significant and values of less than
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 |
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
,
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
1015 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
107 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
1011 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
.
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
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 1014 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 (![]() |
Open with DEXTER |
The inter-outburst accretion phase leads to a circumstellar density
distribution that deviates substantially from the 1/r2 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
,
where
,
with P the orbital period. For
,
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
,
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
in the direction of the RG and
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/r2 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 5
1013 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
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 (
![]() ![]() |
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
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
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
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
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
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 |
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
1043 erg
and 2
,
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
to
,
but is comparable to the 2.2
obtained from the self-consistently computed accretion phase. An isothermal equation of state with
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 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
,
with
for most times and directions. Poleward values of
are
slightly lower after about day 10, while
for directions
toward the former accretion wake of the WD and times before day 6
after the explosion. Values for
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 (
)
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
,
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
![]() ![]() ![]() ![]() ![]() ![]() ![]() |
Open with DEXTER |
The orbit shrinks at a rate of
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,
and
,
scale linearly with the
RG mass loss. For an RG mass loss rate of
per year, the fractional losses per year are
10-8 and
10-8, respectively, implying a
dimensionless specific systemic angular momentum loss of
.
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.
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
yr-1, the computed accretion rate of
about 10% of the RG mass loss rate results in an accretion rate of
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
106 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.