Free Access
Issue
A&A
Volume 527, March 2011
Article Number A5
Number of page(s) 7
Section Stellar structure and evolution
DOI https://doi.org/10.1051/0004-6361/201015895
Published online 18 January 2011

© ESO, 2011

1. Introduction

Classical novae are cataclysmic stellar events. Their thermonuclear origin, theorized by Schatzmann (1949, 1951) and Cameron (1959) – see also Gurevitch & Lebedinsky (1957) and references therein – has been established through multiwavelength observations and numerical simulations pioneered by Sparks (1969), who performed the first 1-D, hydrodynamic nova simulation. These efforts helped to establish a basic picture, usually referred to as the thermonuclear runaway model (TNR), which has been successful in reproducing the gross observational properties of novae, namely the peak luminosities achieved, the abundance pattern, and the overall duration of the event; see Starrfield et al. (2008), José & Shore (2008), José & Hernanz (2007) for recent reviews.

Many details of the dynamics of nova explosions remain to be explored. In particular, there are many observed cases of nonspherical ejecta, inferred from line profiles during the early stages of the outburst and from imaging of the resolved ejecta, including multiple shells, emission knots, and chemical inhomogeneities. Although the broad phenomenology of the outburst can be captured by 1-D calculations, it is increasingly clear that the full description requires a multidimensional hydrodynamical simulation of such outbursts. To match the energetics, peak luminosities, and the abundance pattern, models of these explosions require mixing of the material accreted from the low-mass stellar companion with the outer layers of the underlying white dwarf. In fact, because of the moderate temperatures achieved during the TNR, a very limited production of elements beyond those from the CNO-cycle is expected (Starrfield et al. 1998, 2009; José & Hernanz 1998; Kovetz & Prialnik 1997; Yaron et al. 2005), and the specific chemical abundances derived from observations (with a suite of elements ranging from H to Ca) cannot be explained by thermonuclear processing of solar-like material. Mixing at the core-envelope interface represents a likely mechanism.

The details of the mixing episodes by which the envelope is enriched in metals have challenged theoreticians for nearly 40 years. Several mechanisms have been proposed, including diffusion-induced mixing (Prialnik & Kovetz 1984; Kovetz & Prialnik 1985; Iben et al. 1991, 1992; Fujimoto & Iben 1992), shear mixing at the disk-envelope interface (Durisen 1977; Kippenhahn & Thomas 1978; MacDonald 1983; Livio & Truran 1987; Kutter & Sparks 1987; Sparks & Kutter 1987), convective overshoot-induced flame propagation (Woosley 1986), and mixing by gravity wave breaking on the white dwarf surface (Rosner et al. 2001; Alexakis et al. 2004). The multidimensional nature of mixing has been addressed by Glasner & Livne (1995) and Glasner et al. (1997, 2005, 2007) with 2-D simulations of CO-novae performed with VULCAN, an arbitrarily Lagrangian Eulerian (ALE) hydrocode capable of handling both explicit and implicit steps. They report an effective mixing triggered by Kelvin-Helmholtz instabilities that produced metallicity enhancements to levels in agreement with observations. Similar studies (using the same initial model as Glasner et al. 1997) were conducted by Kercek et al. (1998, 1999) in 2-D and 3-D, respectively. Their results, computed with the Eulerian code PROMETHEUS, displayed mild TNRs with lower peak temperatures and velocities than Glasner et al. (1997) and insufficient mixing. While Glasner et al. (1997) argue that substantial mixing can naturally occur close to peak temperature, when the envelope becomes fully convective and drives a powerful TNR, Kercek et al. (1998) conclude instead that mixing must take place much earlier: if it occurs around peak temperature, it leads to mild explosions or to events that do not resemble a nova.

The differences between these studies have been carefully analyzed by Glasner et al. (2005), who conclude that the early stages of the explosion, before TNR ignition when the evolution is quasi-static, are extremely sensitive to the adopted outer boundary conditions. They show that Lagrangian simulations, in which the envelope is allowed to expand and mass is conserved, lead to consistent explosions. In contrast, in Eulerian schemes with a “free outflow” outer boundary condition, the choice adopted in Kercek et al. (1998), the outburst can be artificially quenched. The scenario was revisited by Casanova et al. (2010), who show that simulations with an Eulerian scheme – the FLASH code – and a proper choice of the outer boundary conditions can produce deep-mixing of the solar-like accreted envelopes with core material. The puzzling results reported by Kercek et al. (1998) stress the need for a systematic evaluation of the effect that different choices of model parameters (e.g. the intensity and location of the initial temperature perturbation, resolution, or size of the computational domain) may have on the results. To this end, we performed a series of 9 numerical simulations in 2-D aimed at testing the influence of these parameters on the level of metal enhancement of the envelope. Here we report the results of these simulations.

Our paper is organized as follows. In Sect. 2 we explain our input physics and initial conditions. Then Sect. 3 is devoted to studying the mixing at the core-envelope interface for our fiducial model. In Sect. 4 the effects of the size of the initial perturbation are analyzed, while in Sect. 5 we discuss the effects of the size of the computational domain. In Sect. 6 we quantify the influence of the grid resolution. Finally, in Sect. 7 we discuss the significance of our results and draw our conclusions.

2. Input physics and initial conditions

The simulations reported here were performed with FLASH, a parallelized explicit Eulerian code, based on the piecewise parabolic interpolation of physical quantities for solving the hydrodynamical equations, and with adaptive mesh refinement (see Fryxell et al. 2000). As in Casanova et al. (2010), we used the same initial model as Glasner et al. (1997) and Kercek et al. (1998): a 1M CO white dwarf that accretes solar composition matter (Z = 0.02) at a rate of 5 × 10-9M yr-1. The model was evolved spherically (1-D) and mapped onto a 2-D cartesian grid, when the temperature at the base of the envelope reached  ≈ 108 K. It initially comprised 112 radial layers – including the outermost part of the CO core – and 512 horizontal layers. The mass of the accreted envelope was about 2 × 10-5M. Nuclear energy generation is handled through a network of 13 species (1H, 4He, 12,13C, 13,14,15N, 14,15,16,17O, and 17,18F), and connected through 18 nuclear reactions. We adopted the conductive and radiative opacities from Timmes (2000) and an equation of state based on Timmes & Swesty (2000). Periodic boundary conditions were imposed on both sides of the computational domain with vertical hydrostatic equilibrium with an outflow constraint at the top and a reflecting constraint at the bottom on the velocity (see Zingale et al. 2002). A summary of the main characteristics of the 9 models computed in this work is given in Table 1, where H is the distance from the perturbation to the initial core-envelope interface, Rx and Ry1, δT, and δt are the size, strength, and duration of the temperature perturbation, and Z the mass-averaged metallicity of the envelope at the end of the calculations.

3. 2-D simulations of mixing at the core-envelope interface

In this section, we describe the basic features of our fiducial model A, as a framework for further discussion of the effect of the parameter choices on our results. A movie, showing the development of Kelvin-Helmholtz instabilities, in terms of the 12C content, up to the time when the convective front hits the upper computational boundary, ModelA-2D.wmv, is available online or at http://www.fen.upc.edu/users/jjose/Downloads.html. The simulation was performed for the conditions of model A, as summarized in Table 1.

For all sequences reported in this work, the relaxation of the initial model to guarantee hydrostatic equilibrium, together with the small amount of numerical viscosity – in contrast with the simulations performed by Glasner et al. (1997) – requires an initial perturbation close to the core-envelope interface to trigger the onset of instabilities early in the calculations. The initial perturbation is applied to the temperature using four parameters: strength, location, size and duration. Model A assumes a top-hat temperature perturbation wherever ((x − x0)/Rx)2 + ((y − y0)/Ry)2 ≤ 1, where x and y are the space coordinates measured from the center of the perturbation, (x0, y0), and Rx and Ry indicate its spatial extent. We fixed x0 = 5 × 107 cm in all sequences. The strength of the perturbation is 5% in temperature in all cases but one (see Table 1). It is 1 km wide, applied only during the initial timestep (that is, the temperature is fixed only during 10-10 s), and imposed on the core-envelope interface (y0 = 5.51 × 108 cm). The resolution adopted in model A is 1.56 × 1.56 km, and the size of the computational domain is 800 × 800 km.

thumbnail Fig. 1

Snapshot of the development of early instabilities, which later spawn Kelvin-Helmholtz instabilities, shown in terms of the 12C mass fraction (in logarithmic scale) for model A, 158 s from the start of the simulation when the core-envelope interface temperature is Tbase ~ 1.36 × 108 K.

Open with DEXTER

The initial perturbation drives a shear flow that triggers the formation of instabilities (Fig. 1), about 150 s after the start of the simulation. As soon as material from the core is mixed into the envelope, small convective cells develop. At this early stage, the fluid has a large Reynolds number, with a characteristic eddy length of 50 km, fluid velocities ranging between v = 105 − 106 cm s-1, and a dynamic viscosity2 of 2 × 104 P. The fluid velocity v remains below the speed of sound cs (that is, the Mach number Ma = v/cs is always less than unity, see Fig. 2), hence, the fluid displays a deflagration rather than a detonation – see Williams (1985), for a thorough analysis of the differences between flame propagation under detonation and deflagration conditions. At t = 235 s, the characteristic eddy turnover time is lconv/vconv ~ 10 s. The merging of the small convective cells into large eddies, characteristic of 2-D simulations, with a size comparable to the height of the envelope, reinforces the injection of CO-rich material into the envelope. Convection becomes more turbulent. At this stage (t = 450 s), the nuclear energy generation rate exceeds 1015 erg g-1 s-1, while the characteristic convective timescale decreases to  ~5 s. The convection front propagates progressively upwards (Fig. 3, left panel), with a velocity of  ~ 10 km s-1, and eventually reaches the top of our computational domain. The envelope base reaches a peak temperature of 1.64 × 108 K. At this time (t = 496 s), when matter starts to cross the outer boundary of the computational domain, we stop the calculations because of the Eulerian nature of the FLASH code. At this final stage, the mean mass-averaged metallicity in the envelope reaches Z ~ 0.21. It is worth noting, however, that the convective eddies are still pumping metal-rich matter through the core-envelope interface. Hence, it is likely that the final metallicity in the envelope will be larger. The simulation shows that the induced Kelvin-Helmholtz vortices can naturally lead to self-enrichment of the accreted envelope with core material to levels that agree with observations and that the expansion and progress of the runaway is almost spherically symmetric for nova conditions even for a point-like TNR ignition.

4. Effect of the initial perturbation

To quantify the influence of the initial perturbation on our results, we have performed a series of 2-D hydrodynamic tests for a set of different durations, strengths (intensities), locations and sizes of the perturbation. For simplicity, a top-hat perturbation, centered at x0 = 5 × 107 cm, has been adopted in all models reported in this work.

The effect of the duration of the perturbation was checked by means of a test case (model B), identical to model A but with a perturbation lasting for 10 s. As shown in Table 1, the characteristic timescales for model B, such as the time required for the first instabilities to show up, TKH, or the time needed by the convective front to hit the outer boundary, tY, become shorter. The role played by a temperature perturbation can be understood in terms of the energy injected into the envelope: the longer the duration of the perturbation, the larger the energy injected, and thus, the shorter the characteristic timescales of the TNR. This has little effect, however, on the overall metallicity enhancement in the envelope since a final CNO mass fraction of  ~0.212 was found in model B, whereas  ~0.224 resulted in model A.

Both models A and B assumed temperature perturbations of δT ~ 5% during the initial timestep (~10-10 s). To test the possible influence of the strength of the perturbation, a test case with δT ~ 0.5% (model C) has also been computed. As shown in Table 1 and Fig. 4, the time evolution of models A and C is very similar, and hence, similar final mean CNO mass fractions at the end of the simulations were found (with Z = 0.209 in model C).

Table 1

Models computed.

thumbnail Fig. 2

Mach number at two different moments of the simulation, t = 230 s (left panel) and 496 s (right panel), for model A.

Open with DEXTER
thumbnail Fig. 3

Left panel: propagation of the convective front as a function of time, for models A, H, and I. Right panel: temperature profile versus radius at two different times, t = 0 s (solid line; Tbase = 9.84 × 107 K) and t = 496 s (dashed line; Tbase = 1.64 × 108 K), for model A.

Open with DEXTER

The effect of the location of the perturbation along the vertical axis has also been studied: whereas model A assumed a temperature perturbation of  ~5%, applied at the innermost envelope shell (y0 = 5.51 × 108 cm), in model D, a similar perturbation was placed  ~5 km above the core-envelope interface (y0 = 5.515 × 108 cm). Both models exhibit a very similar temporal evolution, with almost identical times for the appearance of the first instabilities and for the time required to reach the outer boundary. Similar envelope mean CNO mass fractions (0.224 and 0.235, respectively) were also found.

Finally, the influence of the size of the perturbation has also been analyzed. Whereas model D was evolved with an initial temperature perturbation of size Rx = 1 km and Ry = 1 km, model E assumed Rx = 5 km and Ry = 5 km. As before, very similar characteristic timescales (see Table 1) and final mean CNO mass fractions (0.235 and 0.209, respectively) were found.

To summarize, the specific choice of the parameters that define the initial temperature perturbation has a negligible effect on metallicity enhancement of the envelope.

5. Effect of the size of the computational domain

The choice of the computational domain represents a compromise between computational time requirements and numerical accuracy. Several considerations constrain its minimum size. On one hand, the merger of large convective eddies often found in 2-D simulations may be severely affected by the adoption of a small computational domain. On the other hand, nova outbursts eventually result in mass ejection. With an Eulerian code such as FLASH, it is not possible to track the material that flows off the grid, and hence, it is important to use domains that are as large as possible along the radial direction (while being sufficiently wide along the horizontal axis). Unfortunately, when the initial 1-D model is mapped into the 2-D grid, and relaxed to guarantee hydrostatic equilibrium, densities quickly underflow values for large heights (Zingale et al. 2002).

The specific size adopted for most of the models computed in this work, i.e. 800 × 800 km, is a bit smaller than those used in Glasner et al. (1997) – 0.1πrad, in spherical-polar coordinates – and in Kercek et al. (1998) – 1800 × 1100 km, in a cartesian, plane-parallel geometry. In this section, we analyze possible dependences of the results on the adopted size of the computational domain. To this end, two additional simulations were performed. In the first one (model F), a wider computational domain has been adopted (i.e., 1600 × 800 km). In the second (model G), aimed at testing the influence of the vertical (radial) length, a domain of 800 × 1000 km has been used (where the choice of 1000 km results from numerical restrictions that limit the vertical extent of our computational domain).

As shown in Table 1, the horizontal width (model F) has no noticeable effect on the timescales of the simulations, either for the time required for the onset of the first instabilities or for the time required for the convective front to reach the outer boundary. The mass-averaged CNO abundance in the envelope reached  ~0.206 at the end of this simulation, close to the value found for model A. These results confirm that 800 km is an appropriate choice for the width of the computational domain, stressing that above a threshold value the course of the TNR is insensitive to the adopted width, in agreement with the sensitivity study performed by Glasner et al. (2007).

thumbnail Fig. 4

Left panel: time evolution of the nuclear energy generation rate (in erg s-1) for the 9 models computed in this work. Right panel: final CNO mass fraction versus radius.

Open with DEXTER

The specific length adopted along the vertical direction (see model G), while unimportant for the time of appearance of the instabilities (around 155 s after the start of the simulation, as in model A), affects the time required to reach the outer boundary, located 200 km above the value adopted for model A. Moreover, the larger extension of the computational domain along the radial (vertical) direction allows the convective eddies to pump additional metal-rich core material into the envelope compared with all the simulations reported previously in this paper. Indeed, the mean, mass-averaged metallicity in model G achieves the largest value of all the simulations reported,  ~0.291. This result suggests that the likely mean mass-averaged metallicity driven by Kelvin-Helmholtz instabilities should be Z ≈ 0.3. In summary, we conclude that the size of the computational domain, above a certain threshold value, has little influence on the physical quantities that are more directly related with the mixing process at the core-envelope interface.

6. Effect of the grid resolution

thumbnail Fig. 5

Snapshots of the 1H (upper panels) and 12C (middle panels) mass fractions at t ~ 395 s (model A; left panels), and 688 s (model I; right panels). Lower panels: the number of blocks administered, at this stage, is 3184 for model A, and 43 800 for model I. In both simulations, FLASH divides each block in 8 cells. Structures such as vortexs are better resolved in the finer resolution model I.

Open with DEXTER

All simulations discussed so far (e.g., models A to G) were performed with a resolution of 1.56 × 1.56 km, a value similar to the minimum resolutions adopted in Glasner et al. (2007) which is roughly  ~1.4 × 1.4 km, and in Kercek et al. (1998), 1 × 2 km. To quantitatively assess the possible effect of the resolution, two additional test cases were computed with exactly the same input parameters as in model A but with two different resolutions: 1 × 1 km (model H) and 0.39 × 0.39 km (model I)3.

As shown in Table 1, the increase in resolution produces a delay in the time required for the first instabilities to develop, tKH. This seems to be a numerical artifact. In models with a coarser resolution, the larger size of the blocks artificially generates a larger numerical diffusion compared to models with a finer resolution (a similar resolution dependence is clearly seen as well in the Kercek et al. 1998, simulations). Actually, the ratio of differences in the initial build up times (i.e., (model I-model A)/(model H-model A)) scales approximately as the zone size dimensions to the power of two. This is a purely numerical perturbation that forces the development of instabilities. To test this hypothesis, we computed an additional test case (not included in Table 1), identical to model A but without any initial perturbation. The onset of the instabilities in such an extremely low numerical diffusion regime is substantially delayed. The simulations reported by Glasner et al. (1997) also show the early appearance of instabilities in a model with substantial numerical noise: within a very short time (about 10 s), thenumerical noise (round-off) seeds an intense convective flow in the envelope without any artificial perturbations.

A similar behavior is also found for the time required for the convective front to reach the outer boundary, tY, and for the history of the nuclear energy generation rate (Fig. 4). As expected, filamentary structures and convective cells are better resolved in the finer resolution model I, compared to those computed with somewhat coarser grids (models A and H; see Fig. 5). These minor differences do not, however, show significant variations in the final, mean CNO abundances achieved in the envelope: while Z ~ 0.224 in model A, models H and I yield 0.201 and 0.205, by mass, respectively. Similar agreement is found in the peak temperatures achieved and in the overall nuclear energy generation rates (Fig. 4).

Thus, the adopted resolution has not a critical effect for the mixing models presented in this work. The variation in the final mean CNO abundance in the envelope, under the range of resolutions adopted, is only about 12% (when comparing results for models A, H, and I), a quite reasonable value.

7. Discussion and conclusions

In this paper we have reported results for a series of nine 2-D numerical simulations that test the influence of the initial perturbation (duration, strength, location, and size), the resolution of the grid, and the size of the computational domain on the results. We have shown that mixing at the core-envelope interface proceeds almost independently of the specific choice of such initial parameters, above threshold values.

The study confirms that the metallicity enhancement inferred from observations of the ejecta of classical novae can be explained by Kelvin-Helmholtz instabilities, powered by an effective mesoscopic shearing resulting from the initial buoyancy. Fresh core material is efficiently transported from the outermost layers of the white dwarf core and mixed with the approximately solar composition material of the accreted envelope. As soon as 12C and 16O are dredged up, convection sets in and small convective cells appear, accompanied by an increased nuclear energy generation rate. The size of these convective cells increases in time. Eventually, these cells merge into large convective eddies with a size comparable to the envelope height. The range of mean mass-averaged envelope metallicities obtained in our simulations at the time when the convective front hits the outer boundary, 0.21 − 0.29, matches the values obtained for classical novae hosting CO white dwarfs.

It is, however, worth noting that the convective pattern is actually produced by the adopted geometry (e.g., 2-D), forcing the fluid motion to behave very differently than 3-D convection (Shore 2007; Meakin & Arnett 2007). Nevertheless, the levels of metallicity enhancement found in our 2-D simulations will likely remain unaffected by the limitations imposed by the 2-D geometry (Arnett, private communication). Fully 3-D simulations aimed at testing this hypothesis are currently underway.


1

The different values adopted for Rx and Ry in models F and G are used to account for the assumption of a rectangular (rather than square) computational domain.

2

The dynamic viscosity evaluates the resistance to flow of a fluid under an applied force. More precisely, it is defined as the tangential force per unit area required to move one horizontal plane with respect to the other at unit velocity when maintaining a unit distance apart by the fluid.

3

For comparison, whereas a maximum number of 5300 blocks are administered in model A, the number of blocks increases up to 83 000 in model I. The total CPU time spent in both simulations, using 256 processors of the MareNostrum supercomputer, has been 3 and 110 kh, respectively.

Acknowledgments

We thank the referee for a careful reading of the manuscript and for constructive comments. The software used in this work was in part developed by the DOE-supported ASC/Alliances Center for Astrophysical Thermonuclear Flashes at the University of Chicago. This work has been partially supported by the Spanish MEC grants AYA2010-15685 and AYA2008-04211-C02-01, by the E.U. FEDER funds, and by the ESF EUROCORES Program EuroGENESIS. We also acknowledge the Barcelona Supercomputing Center for a generous allocation of time at the MareNostrum supercomputer.

References

Online material

Movie (Access here)

All Tables

Table 1

Models computed.

All Figures

thumbnail Fig. 1

Snapshot of the development of early instabilities, which later spawn Kelvin-Helmholtz instabilities, shown in terms of the 12C mass fraction (in logarithmic scale) for model A, 158 s from the start of the simulation when the core-envelope interface temperature is Tbase ~ 1.36 × 108 K.

Open with DEXTER
In the text
thumbnail Fig. 2

Mach number at two different moments of the simulation, t = 230 s (left panel) and 496 s (right panel), for model A.

Open with DEXTER
In the text
thumbnail Fig. 3

Left panel: propagation of the convective front as a function of time, for models A, H, and I. Right panel: temperature profile versus radius at two different times, t = 0 s (solid line; Tbase = 9.84 × 107 K) and t = 496 s (dashed line; Tbase = 1.64 × 108 K), for model A.

Open with DEXTER
In the text
thumbnail Fig. 4

Left panel: time evolution of the nuclear energy generation rate (in erg s-1) for the 9 models computed in this work. Right panel: final CNO mass fraction versus radius.

Open with DEXTER
In the text
thumbnail Fig. 5

Snapshots of the 1H (upper panels) and 12C (middle panels) mass fractions at t ~ 395 s (model A; left panels), and 688 s (model I; right panels). Lower panels: the number of blocks administered, at this stage, is 3184 for model A, and 43 800 for model I. In both simulations, FLASH divides each block in 8 cells. Structures such as vortexs are better resolved in the finer resolution model I.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.