© ESO, 2016

1. Introduction

Space-based observations have shown that the inner cavities of many round or elliptical planetary nebulae (PNe) contain a tenuous and very hot gas that emits predominantly in the soft X-ray domain. From a gas-dynamical point of view, the existence of this gas is expected; the fast central-star wind collides with the slower and denser inner parts of the nebula and becomes shock heated. Given the density and velocity of the stellar wind, the wind shock is adiabatic and the shocked gas is expected to reach temperatures of 10⁷ K or more. The existence of this “hot bubble” is essential for the dynamics of a PN. Its pressure is not only sufficient to support the nebula shell against collapse, it also compresses and accelerates the inner parts of the nebula.

All spectral analyses of X-ray emission reveal unexpectedly low characteristic X-ray temperatures of a few 10⁶ K, or even below 10⁶ K. The emission measure is also much higher than expected. Ruiz et al. (2013) overviews properties of the diffuse X-ray emission from PNe (their Table 3). The present status and preliminary results of the Chandra Planetary Nebula Survey (ChanPlaNS) are presented by Kastner et al. (2012), Freeman et al. (2014), and Montez et al. (2015).

There are two possibilities to remedy the obvious temperature discrepancy between theoretical predictions and observations: (i) heat conduction acts across the interface between the bubble and the nebula as proposed for H ii regions by Weaver et al. (1977); and (ii) there is mixing due to hydrodynamical instabilities, which develop at the same interface. Both processes can reduce the electron temperature and increase the density within the bubble, but quantitative results − i.e., observable quantities that can be tested against real objects − can only be gained by detailed modeling of the relevant physical processes.

First attempts to model mixing processes by means of two-dimensional radiation-gasdynamics simulations are presented by Stute & Sahai (2006), and more recently in more detail by Toalá & Arthur (2014). A comparison of the numerical results with observations is still pending. For the first time, Steffen et al. (2008, hereafter SSW08) study the influence of heat conduction on the structure of hot bubbles numerically. We include heat conduction in our one-dimensional radiation-gasdynamics code, which is tailored for studies of the formation and evolution of PNe. In SSW08, we are able to show that, while the dynamics of a model nebula remains virtually unchanged, the bubble structure and properties such as the X-ray characteristic temperature and luminosity can be well explained by PN models that include heat conduction.

The heat conduction depends on the chemical composition of the plasma, as is shown by, e.g., Spitzer (1962, hereafter S62). In practical applications, a formalism based on a fully ionized pure hydrogen plasma is mostly used, which is a reasonable approximation for stellar and cosmic plasmas in general.

A small fraction of PNe harbor nuclei with hydrogen-deficient as well as helium- and carbon-enriched surfaces (such as the [WC] Wolf-Rayet objects); these objects also show hydrogen-deficient winds and hot bubbles. A hydrogen-deficient plasma − that is simultaneously helium- and carbon-rich − is expected to show heat conduction and line cooling properties that differ considerably from those of a hydrogen-rich plasma. Additionally, because of heat conduction, hydrogen-rich matter from the nebula evaporates into the hydrogen-deficient bubble, where it changes the X-ray and cooling properties of the bubble. A detailed modeling of such hot bubbles that may show a chemical stratification, together with appropriate observations, will be a harbinger for a better understanding of the formation and evolution of hydrogen-deficient central stars of PNe.

A comment is needed on the role of magnetic fields in the context of heat conduction and PNe. The two currently favored mechanisms that could explain the shaping of non-spherical PNe are central-star (CS) binarity and magnetic fields (see, e.g., the review of Balick & Frank 2002). Already very weak magnetic fields inhibit heat conduction efficiently in the direction perpendicular to the field lines (see, e.g., S62; Borkowski et al. 1990; Soker 1994). Our measurements of the surface magnetic fields in CSPNe (Steffen et al. 2014) reveal low values of B< 100 G, translating to field strengths of 1−100 mG in the central hot bubble where heat conduction potentially plays an important role. It turns out that magnetic fields much weaker than 1 μG would still completely suppress heat conduction perpendicular to the field lines for the physical conditions of hot bubbles of PNe.

In the absence of heat conduction, the diffuse X-ray emission of PNe should be much lower and the characteristic X-ray temperature of the hot bubbles should be higher than observed (e.g., SSW08). The observed X-ray properties thus indicate that either any present magnetic fields are non-perpendicular to the radial temperature gradient or some other physical mechanism of the same efficiency as heat conduction is at work (e.g., hydrodynamic instabilities, see Toalá & Arthur 2014). We restrict our present study to spherical geometry and, as in our previous work, we assume that magnetic fields are absent.

It is the goal of this paper to develop an algorithm for the calculation of heat conduction coefficients for arbitrary chemical compositions and to apply it to the calculation of bubbles with stratified chemical compositions. We overview existing theories in Sect. 2, where we also present a simplified method that suffices for all astrophysical applications. We compare the simpler approach with a more detailed and complex approach in Sect. 3 to show that differences are small. Section 4 is devoted to a detailed comparison between our new, element-dependent expressions and the pure hydrogen case. We used our one-dimensional radiation-gasdynamics code to model the evolution of a PN with normal composition around a hydrogen-deficient central star with a correspondingly hydrogen-deficient wind. These simulations are guided by properties of a particular object, namely BD+30° 3639. Calculated observables are discussed and compared with existing observations in Sect. 5. We finish the paper with conclusions in Sect. 6.

2. Brief overview of heat conduction theories

The heat conduction theory we used is based on the Fokker-Planck equation, and was first developed by Cohen et al. (1950) and Spitzer & Härm (1953), see also S62. The authors note that the accuracy of the theory is about 5−10%. Assumptions of this theory are that the gas is fully ionized and that electrons are responsible for all interactions. Devoto (1966, 1967b,a) develops an alternative description of heat conduction (and other transport coefficients) based on a fourth-order series expansion of the Chapman-Enskog-Burnett theory. This alternative description works with partially ionized gases, and is also more computationally demanding than the description of S62; we compare the two theories in Sect. 3. As in most previous studies of heat conduction, we neglected magnetic fields.

Heat conduction can be described as a diffusion process. The heat flux q is a linear function of the gradient of the electron temperature T_e, (1)where D is the diffusion coefficient. Following SSW08, we write D as a function of the mean free path λ (see below). We used the same approach as Cowie & McKee (1977) and multiplied the electron-electron equipartition time t_eq (see Eq. (5-31) in S62) with a characteristic thermal velocity, v_char, to calculate a mean free path,

which is why (2)Here, m_e is the electron mass, k_B the Stefan-Boltzmann constant, n_e the electron density, e the electron charge, lnΛ the Coulomb Logarithm, Z the effective (or mean) charge, and C₀ = 2.62479 × 10⁵ [cm^-2 K^-2]. Spitzer (1962) expresses Λ, neglecting shielding by positive ions, as (3)where h is the Debye length, p₀ the average closest impact parameter, and Z_e (=1) the electron charge. We used the electron-density weighted effective charge Z of all contributing ions with density n_i and charge Z_i, (4)It is necessary to multiply Λ by the following factor when T_e> 4.2 × 10⁵ K (S62), (5)where α is the fine structure constant, c the speed of light, and w the speed of a test electron. We write the logarithm of Λ for T_e ≤ 4.2 × 10⁵K as (6)and replace this expression with the following for higher temperatures¹, T_e> 4.2 × 10⁵K, (7)

Fig. 1 Normalized Z-dependent term ζZ of λZ in D versus effective charge Z.

Writing the expression for the coefficient of heat conduction, we follow S62, but keep the charge-dependent factors ϵ and δ_T in the expression (these are provided in Table 1),

The parameter λ_Z is the Z-dependent mean-free path. For a pure hydrogen composition, D agrees with Eq. (4) in SSW08. We show the charge-dependent term of D, ζ_Z, versus Z in Fig. 1; this term excludes the Z-dependent component of lnΛ, which typically is negligible at high temperatures (T_e ≳ 10⁶K). For given T_e and n_e, D decreases with Z, but the heat flux q also depends on the temperature gradient, which changes because of the heat conduction. The impact of the effective charge is modest in all cases. For an atmosphere with pure fully ionized carbon (Z = 6), D is 43 per cent of the value of a pure hydrogen plasma.

The heat flux cannot be larger than the heat content (per unit volume; E_th) times a characteristic electron transport velocity (v_tr; e.g., Cowie & McKee 1977). Following SSW08, we define a limiting heat flux q_sat using this property, (11)This term agrees with Eq. (7) in Cowie & McKee (1977) for pure hydrogen when ϵ = 0.4. We rewrite Eq. (11) using Eqs. (1), (8), and (10) and get an expression for the mean free path λ_sat,Z under saturated conditions, (12)where f_ϵ varies very weakly with Z through ϵ, and f_ϵ,1 = 0.255604 for the pure hydrogen case.²

The heat conduction terms are included in our nebel code following our approach in SSW08 (see Sect. 2.2.2 therein). Accounting for both regular and saturated heat flux in the same expression, we replaced the mean free path of methods 1 and 2 (see Eqs. (7) and (8) in SSW08) with We interpolated the values of ϵ and δ_T in Table 1 using a one-dimensional rational spline (Späth & Meier 1990).

Fig. 2 Diffusion-coefficient ratios [D] 4/D (gray lines) and [D] 3/D (black lines) drawn versus the electron temperature Te and the electron density ne. The used densities are ne = 100 cm-3 (no symbol), ne = 102 cm-3 (diamonds), and ne = 104 cm-3 (squares). The six panels show a) m1; b) m2; c) m3; d)m4 (solid) and m5 (dotted); e) m6; f) m7. The vertical gray lines indicate Te = 4.2 × 105K.

3. Simplified comparison with an alternative theory

We compare the diffusion coefficient D (Eq. (8)) with the corresponding fourth- and third-order accurate expressions of Devoto (1967b), The lengthy q^mp elements (m,p ∈N) are provided as functions of Coulomb cross sections in the appendix of Devoto (1967b); however, there is a typo in those expressions. A numerical check shows that the sign of the second term of q¹³ should be positive as in the other equations (Eq. (a10); “+ 8 ∑ _jn_en_j[...”). We used the Coulomb cross sections of Devoto (1967a, see Eq. (5) therein). The two theories compared here (S62 and Devoto 1967b,a) use different values for the Debye length; the ratio of Λ (see Eq. (3)) between the two theories is .

To study both hydrogen-rich and hydrogen-deficient plasmas, we used a set of seven models (see Table 2). We only considered fully ionized gases and electron-electron interactions; our tests with electron-ion and ion-ion interactions revealed a negligible contribution (this conclusion was reached by comparing the electron-electron terms with the heavy-particle contributions as outlined by Ulmschneider 1970, see his Eq. (10)). The temperature and the density intervals we considered are 10⁴ ≤ T ≤ 10⁸ K and 10⁰ ≤ n_e ≤ 10⁴ cm^-3. While this is a qualitative study, in reality none of the elements heavier than hydrogen is fully ionized at T_e = 10⁴K; for studies at such low temperatures we instead refer to Ulmschneider (1970).

The ratios [D] ₄/D and [D] ₃/D for the different chemistry sets of Table 2 are illustrated in Fig. 2. The first impression is that differences are modest, they amount to less than 14% in the shown examples. Here we compare outcome of mixtures presented in Table 2.

• For a pure hydrogen plasma (m₁; Fig. 2a), the [D] ₄/D ([D] ₃/D) ratio shows small differences 0.900−1.04% (0.980−1.12%) lower than D. Differences are slightly smaller at higher densities than at lower values.

• For a hydrogen-helium mixture, differences are larger (m₂; Fig. 2b). The diffusion coefficient [D] ₄ ([D] ₃) is 3.68−3.78% (3.58−3.68%) higher than D.

• For a helium-dominant mixture (m₃; Fig. 2c), [D] ₄ ([D] ₃) is 2.2−3.0% (2.5−3.3%) lower than D.

• For a carbon-only mixture (m₄ and m₅; Fig. 2d), [D] ₄ ([D] ₃) is 4.0−5.3% (2.2−3.4%; m₄) and 2.0−3.4% (0.3−1.9%; m₅) higher than D.

• Differences are higher with more WR-like abundances (m₆ and m₇; Figs. 2e and f). [D] ₄ ([D] ₃) is 12.4−13.8% (11.5−12.9%; m₆) and 10.5−12.0% (9.81−11.1%; m₇) higher than D.

We conclude that differences are small in general. The simpler expression in Eq. (8) can be used if there is only hydrogen, or if theory-dependent errors in D of up to 15% can be tolerated. When the chemistry is more complex and more accuracy is required, the fourth-order expression [D] ₄ should be used; using this expression instead of [D] ₃ the diffusion coefficient is up to 2% higher. Differences in calculation times are negligible with the computing power available today. Considering the many uncertainties in the other parts of this study, we used the simpler (and faster) formulation based on S62 as outlined in Sect. 2.

4. Comparison of models using the general or the pure-hydrogen heat conduction equations

Here we study how the weak abundance dependence of the heat conduction is reflected in the outcome of two models that only differ in the use of the heat conduction equation. Model HC uses the chemistry-dependent heat conduction formulation of this paper (D, Eq. (8)), and model HC_H uses the pure hydrogen formulation of SSW08 (D_SSW08, their Eq. (4)). For given T_e and n_e, the diffusion coefficient D of the two formulations differs by the Z-dependent factor ζ_Z shown in Fig. 1.

4.1. Considerations when setting up models for this study

We describe our one-dimensional radiation hydrodynamic (RHD) models of PN envelopes in detail in Perinotto et al. (1998, 2004, and references therein). We emphasize that our models calculate time dependent ionization, recombination, heating, and cooling in a region that extends out to about one parsec away from the CS. The cooling function is composed of the contribution of all considered ions (see Table 3), and up to twelve ionization stages are taken into account for every individual ion. Physical input parameters to the models include properties of the coupled CS model, element abundances, and the initial density and velocity structures of matter in the asymptotic giant branch (AGB) envelope.

The purpose of the study is to better model a full system with a Wolf-Rayet star, the preceeding AGB wind, the fast wind, a bubble, and a nebula. To do this we needed to use conditions for these objects that were as realistic as possible. This includes abundances, the AGB evolution, the fast wind evolution, the central-star evolution, and the time when the wind switched from an H-rich state to an H-deficient state. Several of these properties are unknown or poorly known.

Fig. 3 Initial AGB wind-envelope structures. The left-side axis is used with heavy particle (electron) densities drawn with a thick solid (dotted) line. The right-side axis is used with flow velocities drawn with a thin solid line. All properties are drawn as function of radius. All lines are drawn as histograms; individual gridpoints are visible in the innermost part of the envelope. The central-star mass is indicated in the panel. The large density trough at 1018cm is due to the last thermal pulse 40 000 yr before the star left the AGB.

4.2. Detailed physical model properties

In the initial model, we assumed mean Galactic-disk PNe abundances (, these were first quoted by Perinotto et al. 1998) for nine elements (see Table 3); except for values of carbon and nitrogen, is close to solar. The abundance values ϵ_i are specified in both logarithmic number fractions relative to hydrogen, i.e., ϵ_i = log  (n_i/n_H) + 12, and as mass fraction percentages β_i = A_iX_Hn_i/n_H, where A_i is the atomic mass relative to hydrogen and X_H the hydrogen mass fraction. Other properties of the initial model are the central-star mass M = 0.595 M_⊙, the effective temperature T_eff = 5865 K, the luminosity L = 5601 L_⊙, and properties of the slow AGB wind such as the expansion velocity v_∞ = 10 km s^-1 (see Steffen et al. 1998). The density and the velocity structures of the initial AGB wind-envelope model are shown in Fig. 3.

Immediately at the start of the RHD calculations, the abundances of the CS fast wind are replaced with a second set of abundances, . Assuming empirically determined abundances of the object BD+30° 3639, we based these abundances on values of Marcolino et al. (2007, see Table 2 therein). The CS fast wind adds mass at the inner boundary of the grid.

Fig. 4 Models, scaled models, and empirical values of fast winds of PNe: a) the mass-loss rate Ṁ; b) the outflow velocity v; and c) the mechanical wind luminosity Lwind. Scaled values of Ṁ are illustrated with both magenta () and blue () star symbols. All symbols are indicated with both object (Table 4) and gray interconnecting lines for easy identification. Thick and thin lines indicate models based on Pauldrach et al. (1988) and actually used relations, respectively. The shaded part where Teff ≥ 60 000 K indicates regions where both the models and the used relations fit empirical values poorly.

Pauldrach et al. (1988) present fast-wind models derived for hydrogen-rich conditions; the model mass-loss rate Ṁ decreases with time, while the outflow velocity v increases with the effective temperature (i.e., with time). More recently, Pauldrach et al. (2004) find that the mass-loss rate instead appears to increase with time for T_eff ≲ 50 000 K. Leuenhagen et al. (1996) present results of observations of CSPNe with hydrogen-deficient compositions, i.e., [WC] stars, where the mass-loss rate increases with time. Outflow velocities simultaneously increase less steeply than with a hydrogen-rich composition. Empirical values of effective temperatures, luminosities, and mass-loss rates are listed for thirteen Wolf-Rayet stars in Table 4. The table shows both mass-loss rates of the respective original reference and values that were rescaled, Ṁ_s. The rescaled values use a standard luminosity log L = 3.7 L_⊙ and a clumping density ratio . The rescaling was made using the invariance of the so-called transformed radius R_t (Schmutz et al. 1989; Hamann & Koesterke 1998), (17)where R_∗ is a stellar photospheric radius (). This expression yields when R_t and T_eff are kept fixed. Mass-loss rates (assuming both and ), outflow velocities, and mechanical wind luminosities (L_wind = Ṁv²/ 2) of both models and scaled models are shown in Fig. 4, together with the model relations of Pauldrach et al.

The mass-loss descriptions of Pauldrach et al. (1988) and Pauldrach et al. (2004) are both poor fits to the scaled empirical values of [WC] stars. Instead, we made a linear fit to the empirical mass-loss rates as function of effective temperature, where T_eff ≲ 50 000 K; we accounted for all objects but M 2-43, NGC 40, and He 2-459. For the outflow velocities, we simply scaled the expression of Pauldrach et al. (1988) with the factor 0.583; the scaled relation lies higher than the empirical values since we found the high velocities necessary to form a bubble at about T_eff ≃ 47 000 K, which is the currently measured temperature of BD+30° 3639. Despite a flattening mass-loss rate as is seen in the measurements, the mechanical luminosity increases monotonically with time. The figure illustrates poor fits where T_eff ≳ 60 000 K; we show this figure to emphasize the need for models that work with H-deficient abundances at low as well as high effective temperatures.

Heat conduction is always active during the RHD calculations. While our nebel code can use both numerical approaches HC1 and HC2 presented in SSW08 to account for heat flux below and above the saturation limit, we only used the HC2 limiter here (Eq. (14)). After the RHD calculations were started, the full model evolution was followed across the Hertzsprung-Russell diagram until the effective temperature reached T_eff ≈ 70 000 K. We are only interested in the formation and early bubble evolution, which is why we did not consider the subsequent evolution.

4.3. Comparing the HC and HC_H models

In comparison to the models that use an H-rich fast wind, models that use an H-deficient abundance composition in the fast wind experience two differences: the different abundances result in much more efficient cooling, and the mechanical luminosity of the fast wind is higher. More efficient cooling results in less available energy to form a bubble, and the PN evolution therefore changes (also see Mellema & Lundqvist 2002). Here, we merely note that in our model the bubble starts forming when T_eff ≃ 50 000 K.

Fig. 5 Physical structure of PN when Teff ≃ 50 500 K, the fast-wind speed v = 1031 km s-1, and the model age is 4246 yr: a) the electron temperature Te; b) the total ion density nA; c) the effective charge Z and the diffusion coefficient D. The electron temperature and the total ion density are shown at three scales as indicated in each panel. Values of model HC (HCH) are drawn with black (magenta) lines; for the diffusion coefficient the lines are blue (light blue). The heat-conduction flux of model HC is saturated in the region that is marked with a blue box. The chemical discontinuity between the H-rich nebula and the H-deficient inner parts is in the bubble-nebula interface. All lines are drawn as histograms to emphasize individual gridpoints. The fast wind (bubble) is drawn on a medium gray (light gray) background for model HC.

Fig. 6 Physical structure of PN when Teff ≃ 60 000 K, the fast-wind speed v = 1154 km s-1, and the model age is 4914 yr. See Fig. 5 for more details.

Fig. 7 Physical structure of PN when Teff ≃ 70 000 K, the fast-wind speed v = 1575 km s-1, and the model age is 5546 yr. See Fig. 5 for more details.

Fig. 8 For model HC at Teff ≃ 70 000 K (Fig. 7), the three panels show radial profiles of the ionization fractions of: a) helium, b) carbon, and c) oxygen. The fast wind (bubble) is drawn on a medium gray (light gray) background.

Fig. 9 Zoom in of Fig. 8 at the bubble-to-nebula interface. The histogram lines show gridpoints of individual ions as indicated in the respective panel.

The physical structure of the inner regions of the PN is shown as radial (r) plots, which include the inner parts of the nebula, the bubble, and outer parts of the fast wind at three different stages during the bubble evolution: Fig. 5) when the bubble forms at T_eff ≃ 50 500 K; Fig. 6) at T_eff ≃ 60 000 K; and Fig. 7) at T_eff ≃ 70 000 K (this case also shows all of the extent of the remaining nebula). Abscissae and ordinates use different scales. The ionization structures of the bubble at T_eff ≃ 70 000 K are shown for helium, carbon, and oxygen in Figs. 8 and 9.

The nebula already starts forming at T_eff ≃ 13 000 K. When the bubble starts forming − here its range is 1.47 ≤ r ≤ 1.525 × 10¹⁷cm − at T_eff ≃ 50 000 K, the nebula already extends 1.525 ≤ r ≤ 2.6 × 10¹⁷cm; the relative size of the bubble to the nebula is 5 per cent. The two models are very similar at this point. At this time, the outermost part of the fast wind has already reached the same temperature as the bubble. The density reaches a lower value of n_A = 0.5 cm^-3 in the fast wind. In the outer bubble, it increases to n_A ≃ 300 cm^-3. The density is the highest immediately ahead of the heat conduction front in the nebula, where n_A ≲ 2 × 10⁴ cm^-3. The temperature in the bubble (nebula) is 4.8−2.5 × 10⁶K (about 8000 K). The effective charge Z is markedly higher in the inner part of the bubble of the HC model than in the outer part. The diffusion coefficient is the property that differs the most between HC and HC_H, owing to its dependence on the effective charge and its high sensitivity to the temperature. However, heat conduction is negligible in the fast wind and the nebula when compared to the bubble (large temperature gradients are required to make it significant). Most of the bubble mass and energy are located at the front of the bubble where the two models are very similar.

Once the effective temperature has evolved to T_eff ≃ 60 000 K, the nebula extends 1.85 ≤ r ≤ 2.9 × 10¹⁷cm, and the bubble 1.60 ≤ r ≤ 1.85 × 10¹⁷cm. The models are still very similar except in the bubble. In the bubble, the HC model reaches T_e ≃ 7 × 10⁶ K at the inner boundary, while the HC_H model where D_HCH is higher than D reaches T_e ≃ 6 × 10⁶ K; the same differences are seen in the outermost part of the fast wind. The temperature differences are reflected in the diffusion coefficient. Because of the less efficient heat conduction of the HC model, the bubble is less efficiently cooled (which could also be due to lower density, and possibly to a slightly more efficient evaporation). The bubble density is slightly lower in the HC model (12−36 per cent). The hot region in the fast wind is slightly smaller in the HC_H model − exact details are unimportant as there is very little mass and energy in this region. The bubble has grown faster than the nebula; the bubble is now 30 per cent the size of the nebula.

At T_eff ≃ 70 000 K, the rim has nearly caught up with the shell, and the nebula is now a thin dense shell; the remaining nebula extends 3.6 ≤ r ≤ 3.8 × 10¹⁷cm. The bubble extends over a much larger region, 2.7 ≤ r ≤ 3.6 × 10¹⁷cm. The temperature difference is larger in this model, T_e,HC ≃ 13 × 10⁶K versus T_e,HCH ≃ 11 × 10⁶K. However, the amount of mass in the inner part of the bubble is a small fraction of the total amount of mass of the bubble. The models are indistinguishable outside of the bubble; this indicates the small influence of the physical details on the PN evolution. The outermost part of the hot fast wind is larger for the HC model, as it cools this region less efficiently.

The lower heat conduction efficiency results in higher temperatures in the new heat conduction formulation. Higher temperatures lead to slightly steeper temperature gradients, which partly compensates for the lower efficiency. The heat conduction difference is smaller in the bubble front (Z closer to 1), which also explains the very similar outcome of HC and HC_H. More nebular mass is evaporated and added to the bubble with more efficient heat conduction; the additional mass due to the evaporation increases the cooling capability of the bubble to further decrease its temperature; this explains the 10⁶K temperature difference between the two models. Finally, Figs. 6 and 7 show how the outer part of the bubble after 632 yr evolution has been filled with evaporated hydrogen-rich matter in the HC model (see Fig. 7 inside the contact discontinuity where 3.3 ≲ r ≲ 3.6 × 10¹⁷cm).

The radial ionization structures of model HC for all ions of helium, carbon, and oxygen are shown in Fig. 8; a zoom in of the interface between the bubble and the interface is shown in Fig. 9. The figures show nearly fully ionized atoms in the bubble. More atoms assume lower ionization levels in the interface to the nebula, as expected from the lower temperature there.

5. Influence of heat conduction on the diffuse X-ray emission

In the following, we investigate the observable consequences of the small differences in the physical structure of the hot bubble between models HC and HC_H. We computed the X-ray emission based on the different hydrodynamical hot bubble structures at selected positions along the evolutionary sequences for this purpose in a post-processing step.

In addition to the two previously discussed models, we have considered a third model, named HC, where the diffusion coefficient is computed exactly as in the case of HC_H, but with a global scaling factor of 1/2−D(HCHC_H)/2. The setup of this simulation is otherwise identical to that of the other two models. The scaling factor of 1/2 roughly corresponds to the ratio D/D_SSW08 at Z = 4 (ionized carbon). The heat conduction efficiency of model HC is in this sense intermediate between models HC_H and HC.

As described in more detail in SSW08, we used the well-documented chianti code (Dere et al. 1997, in version 6.0 Dere et al. 2009), which has been used extensively in various astrophysical contexts. In a first step, the X-ray emission was calculated individually for all spherical shells of a given model, where T_e> 10⁵ K (the contribution of cooler regions is negligible in the X-ray spectral range). If the shell is hydrogen rich, we used the abundance mix (see Table 3) in the calculation of the synthetic X-ray spectrum, otherwise we assumed a chemical composition similar to ³. The slight inconsistency between the chemical composition used in the hydrodynamical simulations and the a posteriori calculation of the X-ray emission is considered irrelevant to our strictly differential comparison.

For given temperature and particle densities, the synthetic spectra include contributions of emission lines of all considered elements and various continua (free-free, free-bound, and two-photon continuum). Since the X-ray emission is optically thin, the total X-ray spectrum emitted by the hot bubble can be computed by adding up the volume-weighted emission of the individual shells. In the present context, we define the X-ray luminosity L_X (erg s^-1) as the total X-ray emission of the hot bubble in the energy range 0.3−2.0 keV (6.2−41.3 Å). For the same energy band, we also computed the characteristic X-ray temperature T_X, which is the emissivity-weighted temperature of the X-ray emitting region, as defined in SSW08 (Eq. (17)).

Fig. 10 Evolution of the synthetic X-ray luminosity LX (top panel) and temperature TX (bottom panel) integrated over the energy range 0.3−2.0 keV versus Teff (which increases with time) for the three different models HC (black diamonds), HCH (red crosses), and HC (blue triangles). The observed position of the PN BD+30° 3639 with a WR-type CS is indicated by a black star for comparison.

We show the evolution of the synthetic X-ray luminosity L_X and temperature T_X versus T_eff (which increases with time) for the three different models in Fig. 10. We conclude from these results that the X-ray properties of models HC_H and HC are indistinguishable within the uncertainties of the numerical models. In addition, our experiment shows that while the provision of a Z-dependent diffusion coefficient hardly influences the X-ray properties of the hot bubble, a global scaling of the pure hydrogen diffusion coefficient by a factor 1/2 has a clearly noticeable impact. In the latter case, the bubble forms earlier in the PN evolution, namely at the time t ≈ 3580 yr when the central star has reached an effective temperature of T_eff ≈ 42 000 K and a wind speed of v ≈ 800 km s^-1. In the case of the HC_H and HC models, the hot bubble forms significantly later, at t ≈ 4250 yr, when T_eff ≈ 50 500 K and v ≈ 1000 km s^-1 (see Fig. 5). The bubble formation is associated with a steep increase of the X-ray luminosity, as is apparent in the upper panel of Fig. 10 near log T_eff = 4.7 and log T_eff = 4.6, respectively.

Fig. 11 Synthetic (unabsorbed) X-ray surface brightness integrated over the energy range 0.3−2.0 keV versus the impact parameter p for the three models HC (black solid line), HCH (red solid line), and HC (blue dotted line) at t = 5545 yr and Teff ≈ 70 000 K.

It is also seen that the reduced heat conduction of HC results in a slightly lower X-ray luminosity and a significantly higher X-ray temperature at any given time. This behavior is in qualitative agreement with our findings in Paper V that heat conduction leads to a delayed formation of the central hot bubble, while enhancing the X-ray luminosity and reducing the characteristic X-ray temperature of its diffuse X-ray emission. However, the increase of L_X and T_X with T_eff is much steeper here than found previously for similar CS masses (see Fig. 18 in SSW08, and Figs. 5 and 6 in Ruiz et al. 2013); this is attributed to the different wind model used in the present work, where the wind luminosity attains much higher values over the same time period (see Fig. 4).

We have also used our models to compute the X-ray surface-brightness distribution in the plane of the sky (in erg cm^-2 s^-1 sr^-1) in the same energy band as above. The X-ray emissivity is for this purpose integrated along lines of sight perpendicular to the plane of the sky for a series of impact parameters, p. Figure 11 shows the radial profile of the X-ray surface brightness for the three heat-conduction models at the end of the considered time evolution (t = 5545 yr) when the CS has reached T_eff ≈ 70 000 K.

The chemical discontinuity is clearly visible as a jump in the radial intensity profile near r = 3.2 × 10¹⁷ cm. The hydrogen-deficient inner bubble appears slightly more extended in the case of HC than in HC_H. Presumably, this is a consequence of the somewhat lower evaporation rate in HC. The higher bubble temperature (see Fig. 7) may also play a role. Together with the slightly lower density, this leads to a somewhat reduced overall X-ray luminosity of HC. We cannot see a difference in the radial expansion rate of the outer edge of the hot bubble between the HC and HC_H, as all physical parameters (including the diffusion coefficient) are very similar in the hydrogen-rich outer part of the hot bubble (Fig. 7). The inner wind shock at r = 2.7 × 10¹⁷ cm is not seen in X-rays owing to the markedly low density in the inner bubble.

The total radial thickness of the HC model is distinctly larger than in the other two models because the bubble of the HC model formed earlier and has had more time to grow. As a result of the general reduction of the diffusion coefficient by a factor 2, the overall temperature of the hot bubble is higher, but the total particle density is lower, such that the gas pressure inside the bubble is the same in all three cases. Remarkably, the thickness of the hydrogen-rich outer part of the bubble, which is the result of the evaporation of cool nebular material, is not significantly reduced in the case of HC despite the lower heat conduction efficiency. The higher temperatures plausibly lead to a slightly steeper temperature gradient at the evaporation front at the outer edge of the hot bubble, which partly compensates for the lower diffusion coefficient assumed in the HC scheme.

A comparison of our model predictions with the observed X-ray luminosity of BD+30° 3639 − which is a PN with a WR-type central star − indicates that our current model is an inadequate representation of this object. As is demonstrated in Fig. 10, the formation of the hot bubble obviously takes place too late (at too high T_eff and v) in the HC and HC_H models. In contrast, the model with reduced heat conduction, HC, appears to work much better. However, it is premature to conclude that the heat conduction schemes HC and HC_H are unrealistic since other factors also influence the circumstances of the hot bubble formation. These factors include the kinematics and the detailed chemical composition (cooling function) of the fast wind. The detailed modeling of BD+30° 3639 is beyond the scope of this paper.

6. Conclusions

We extended our radiation hydrodynamic PN models with a formulation of heat conduction that allows the use of any chemistry. Through this physics extension, we can now model chemically stratified plasmas such as PNe that show hydrogen-deficient hot bubbles.

The heat conduction equation is described by a diffusion coefficient multiplied by a temperature gradient where the chemistry is represented by additional terms in the diffusion coefficient that depend on the effective charge in the plasma. The higher effective charge in a hydrogen-deficient plasma results in a somewhat lower diffusion coefficient than when there is only hydrogen. In the case of PNe, a lower diffusion coefficient simultaneously results in a somewhat higher temperature and a slightly steeper temperature gradient near the interface between the hot bubble and the surrounding cooler nebula. Differences in the heat conduction efficiency between models using the new formulation and the old formulation assuming a pure hydrogen composition are therefore very small. However, the largest deviation is found in the hot bubble, where there is very little mass. We compared two different theories of the heat conduction terms: the Fokker-Planck equation (S62) and the Chapman-Enskog-Burnett theory (Devoto 1966, 1967b,a). We have found that the differences in outcome between the two theories are marginal, and so we used the simpler approach of S62.

We modeled a chemically stratified PN using both the new chemistry-dependent formulation and the old hydrogen-rich formulation of heat conduction. The model consists of a regular hydrogen-rich slow AGB wind where a hydrogen-deficient fast wind starts blowing out into the slow wind when the computations begin. The result is a PN with a nebula outside of a hot bubble. We used empirical measurements of [WC]-type PNe to describe the mass-loss and velocity evolutions of the fast wind. The model parameters were otherwise chosen to be similar to those of a real object where a hot bubble forms when the effective temperature of the central star is T_eff = 50 000 K. We found that the two models evolve nearly identically. The only differences are seen in the hot and tenuous bubble, where the model that uses the chemistry-dependent formulation of heat conduction is up to 2 × 10⁶K hotter than before.

We also calculated observational properties for the models in form of the synthetic X-ray temperature, luminosity, and surface-brightness distribution in the plane of the sky. Here, we also computed a third model that uses the old pure-hydrogen

formulation, but where the diffusion-coefficient is half as high throughout the model domain, to mimic the effect the new-formulation diffusion coefficient inside the hydrogen-poor bubble has on the entire model. The X-ray properties of the first two models are seen to be indistinguishable, while the third model results in a somewhat higher temperature, lower luminosity, and a larger bubble.

In comparison to BD+30° 3639, the hot bubble in our models forms too late at an excessive effective temperature and velocity of the fast wind. These models are therefore inadequate representations of this PN. While we improved the description of the heat conduction, there are several other factors that determine the model evolution. The physical properties of the fast wind, including the mass-loss and velocity evolutions, and the cooling function are such factors, whose descriptions need to be improved further to help the modeling of this and other PNe.

Acknowledgments

We thank W.-R. Hamann, L. Oskinova, and H. Todt (Univ. Potsdam) for fruitful and lively discussions on heat conduction in the context of hydrogen-deficient PNe. C.S. was supported by funds of DFG (SCHO-394/29-1), Land Brandenburg (SAW funds from WGL), PTDESY-05A12BA1, and the BMBF VIP program 03V0843. CHIANTI is a collaborative project involving George Mason University, the University of Michigan (USA), and the University of Cambridge (UK).

References

All Tables

All Figures

	Fig. 1 Normalized Z-dependent term ζZ of λZ in D versus effective charge Z.
In the text

	Fig. 2 Diffusion-coefficient ratios [D] 4/D (gray lines) and [D] 3/D (black lines) drawn versus the electron temperature Te and the electron density ne. The used densities are ne = 100 cm-3 (no symbol), ne = 102 cm-3 (diamonds), and ne = 104 cm-3 (squares). The six panels show a) m1; b) m2; c) m3; d)m4 (solid) and m5 (dotted); e) m6; f) m7. The vertical gray lines indicate Te = 4.2 × 105K.
In the text

Fig. 3 Initial AGB wind-envelope structures. The left-side axis is used with heavy particle (electron) densities drawn with a thick solid (dotted) line. The right-side axis is used with flow velocities drawn with a thin solid line. All properties are drawn as function of radius. All lines are drawn as histograms; individual gridpoints are visible in the innermost part of the envelope. The central-star mass is indicated in the panel. The large density trough at 1018cm is due to the last thermal pulse 40 000 yr before the star left the AGB.

In the text

Fig. 4 Models, scaled models, and empirical values of fast winds of PNe: a) the mass-loss rate Ṁ; b) the outflow velocity v; and c) the mechanical wind luminosity Lwind. Scaled values of Ṁ are illustrated with both magenta () and blue () star symbols. All symbols are indicated with both object (Table 4) and gray interconnecting lines for easy identification. Thick and thin lines indicate models based on Pauldrach et al. (1988) and actually used relations, respectively. The shaded part where Teff ≥ 60 000 K indicates regions where both the models and the used relations fit empirical values poorly.

In the text

Fig. 5 Physical structure of PN when Teff ≃ 50 500 K, the fast-wind speed v = 1031 km s-1, and the model age is 4246 yr: a) the electron temperature Te; b) the total ion density nA; c) the effective charge Z and the diffusion coefficient D. The electron temperature and the total ion density are shown at three scales as indicated in each panel. Values of model HC (HCH) are drawn with black (magenta) lines; for the diffusion coefficient the lines are blue (light blue). The heat-conduction flux of model HC is saturated in the region that is marked with a blue box. The chemical discontinuity between the H-rich nebula and the H-deficient inner parts is in the bubble-nebula interface. All lines are drawn as histograms to emphasize individual gridpoints. The fast wind (bubble) is drawn on a medium gray (light gray) background for model HC.

In the text

	Fig. 6 Physical structure of PN when Teff ≃ 60 000 K, the fast-wind speed v = 1154 km s-1, and the model age is 4914 yr. See Fig. 5 for more details.
In the text

	Fig. 7 Physical structure of PN when Teff ≃ 70 000 K, the fast-wind speed v = 1575 km s-1, and the model age is 5546 yr. See Fig. 5 for more details.
In the text

	Fig. 8 For model HC at Teff ≃ 70 000 K (Fig. 7), the three panels show radial profiles of the ionization fractions of: a) helium, b) carbon, and c) oxygen. The fast wind (bubble) is drawn on a medium gray (light gray) background.
In the text

	Fig. 9 Zoom in of Fig. 8 at the bubble-to-nebula interface. The histogram lines show gridpoints of individual ions as indicated in the respective panel.
In the text

	Fig. 10 Evolution of the synthetic X-ray luminosity LX (top panel) and temperature TX (bottom panel) integrated over the energy range 0.3−2.0 keV versus Teff (which increases with time) for the three different models HC (black diamonds), HCH (red crosses), and HC (blue triangles). The observed position of the PN BD+30° 3639 with a WR-type CS is indicated by a black star for comparison.
In the text

	Fig. 11 Synthetic (unabsorbed) X-ray surface brightness integrated over the energy range 0.3−2.0 keV versus the impact parameter p for the three models HC (black solid line), HCH (red solid line), and HC (blue dotted line) at t = 5545 yr and Teff ≈ 70 000 K.
In the text