Slowing down atomic diffusion in subdwarf B stars

Slowing down atomic diffusion in subdwarf B stars:
mass loss or turbulence?

Haili Hu, C. A. Tout, E. Glebbeek and M.-A. Dupret
Institute of Astronomy, The Observatories, Madingley Road, Cambridge CB3 0HA
Department of Physics & Astronomy, McMaster University, 1280 Main Street West, Hamilton, Ontario, Canada L8S 4M1
Institute d’Astrophysique et Géophysique de l’Université de Liège, Allée du 6 Août 17, 4000 Liège, Belgium
Accepted 2011 July 20. Received 2011 July 19; in original form 2011 June 14

Subdwarf B stars show chemical peculiarities that cannot be explained by diffusion theory alone. Both mass loss and turbulence have been invoked to slow down atomic diffusion in order to match observed abundances. The fact that some sdB stars show pulsations gives upper limits on the amount of mass loss and turbulent mixing allowed. Consequently, non-adiabatic asteroseismology has the potential to decide which process is responsible for the abundance anomalies. We compute for the first time seismic properties of sdB models with atomic diffusion included consistently during the stellar evolution. The diffusion equations with radiative forces are solved for H, He, C, N, O, Ne, Mg, Fe and Ni. We examine the effects of various mass-loss rates and mixed surface masses on the abundances and mode stability. It is shown that the mass-loss rates needed to simulate the observed He abundances () are not consistent with observed pulsations. We find that for pulsations to be driven the rates should be . On the other hand, weak turbulent mixing of the outer can explain the He abundance anomalies while still allowing pulsations to be driven. The origin of the turbulence remains unknown but the presence of pulsations gives tight constraints on the underlying turbulence model.

asteroseismology – diffusion – methods: numerical – stars: chemically peculiar – stars: evolution – stars: mass loss.
pagerange: LABEL:firstpageLABEL:lastpagepubyear: 2011

1 Introduction

Subdwarf B (sdB) stars are subluminous B-type stars that show strong broad Balmer lines and weak or no He I absorption lines. They are ubiquitous not only in our own Galaxy, where they are found in all stellar populations, but also in old galaxies where they are believed to cause the UV-upturn (Brown et al., 1997). Heber (1986) linked the sdB evolutionary stage to that of Extreme Horizontal Branch (EHB) stars which are core He burning stars with an unusual thin H envelope ( M). Such a star must have formed either from a red giant that lost nearly its entire H envelope or from a He white dwarf merger. Although the formation channels are well studied (Han et al. 2003 and references therein), they cannot explain all observed features such as binary period distribution and slow rotation. A promising way to discriminate between different formation channels is with asteroseismology (Hu et al., 2008) but, before this can be achieved, the theoretical models need improving. Currently, the uncertainties in the models are larger than the observational errors obtained with space missions such as CoRoT and Kepler. This paper is part of our ongoing work (Hu et al., 2009; Hu et al., 2010) to improve the input physics used in seismic modelling of sdB stars.

Interestingly, two subgroups of sdB stars are pulsating with variable class names V361 Hya and V1093 Her. The discovery of short-period ( s) pulsations in V361 Hya stars by Kilkenny et al. (1997) coincided with the prediction of unstable pressure()-mode pulsations in sdB stars by Charpinet et al. (1996). The excitation was attributed to the opacity mechanism enabled by Fe accumulating diffusively in the stellar envelope (Charpinet et al., 1997). A few years later, Green et al. (2003) discovered long-period pulsations ( min) in the cooler V1093 Her stars. Fontaine et al. (2003) showed that the same opacity mechanism could also excite long-period gravity()-modes. However, they found too cool a theoretical blue-edge of the -mode’s instability strip compared with observations. The problem can be partly solved by using OP opacities (Badnell et al., 2005) instead of OPAL’s (Iglesias & Rogers, 1996) ànd by assuming that Ni accumulates as well as Fe (Jeffery & Saio, 2006). Furthermore, inclusion of H-He diffusion shifts the theoretical blue-edge even closer to the observed value of (Hu et al., 2009). It is evident that correct modelling of the V1093 Her instability strip awaits seismic models with time-dependent diffusion such as we produce here.

Also interesting is that sdB stars are chemically peculiar without significant spectroscopic differences between variable and constant stars (O’Toole & Heber, 2006; Blanchette et al., 2008). Typically He is deficient with number ratios and there is a trend for He abundance to increase with (Edelmann et al., 2003). The Fe abundance is close to solar irrespective of stellar population and whereas other iron-group elements, such as Ni, can be strongly enriched. Light metals (C, N, O, Ne, Mg) are usually depleted, although there is some scatter between stars (Edelmann et al., 2006; Geier et al., 2010). The observed abundance anomalies can be explained in part by atomic diffusion acting in the radiative envelope. Atomic diffusion alone, however, would cause all surface He to sink on a very short time-scale (about  yr) compared to the EHB evolutionary time-scale (about  yr). Hence there must be a competing process which is generally accepted to be mass loss. Recently turbulence has been proposed to play this role instead (Michaud et al., 2011).

The observed He abundances have been reproduced with models that include atomic diffusion and mass-loss rates of (Fontaine & Chayer, 1997; Unglaub & Bues, 2001). For higher rates the effects of diffusion become unnoticeable, whereas for lower rates He would sink too quickly. Vink & Cassisi (2002) presented a mass-loss recipe for (E)HB stars derived from line-driven wind models. Their results give upper limits for the mass loss rates () and might apply to the most luminous sdB stars that show anomalous H line profiles (Heber et al., 2003; Vink, 2004). However, for most sdBs no observational evidence of mass loss has been found so far. An independent wind model by Unglaub (2008) showed that for weak winds, , the metals decouple from H and He and, for rates below , the winds become purely metallic. Such selective winds could be responsible for some of the observed metal abundance anomalies but they cannot explain the observed He abundances.

Alternatively, Michaud et al. (2011) showed that turbulent mixing of the outer could also produce most of the observed abundance anomalies. Unfortunately, their computational setup restricted their calculations to the first of the E(HB) and to He abundances above . They do not give a physical origin for the turbulence but simply mix the outer layer in order to reproduce the solar Fe abundance observed in most sdB stars. In their models Fe reaches solar abundance in the entire mixed region. This could prevent the driving of pulsation modes.

Also mass loss works against the driving mechanism. According to Chayer et al. (2004) and Fontaine et al. (2006) mass-loss rates above would gradually destroy the Fe reservoir and an sdB pulsator becomes constant within as it evolves. Their results are based on a non-diffusive stellar model assuming an initial Fe profile from an equilibrium between gravitational settling and radiative levitation. The question is, however, can Fe build up in the first place in the presence of stellar winds? If so, shouldn’t diffusion take place at the same time that mass is lost and help build up the Fe reservoir again? To answer these questions we perform a consistent analysis where we compute the effects of mass loss and atomic diffusion simultaneously. Furthermore, we determine which process (mass loss or turbulence) is responsible for retarding atomic diffusion in sdB stars by performing a pulsational stability analysis. Since non-adiabatic effects are responsible for mode driving and damping, such a study is termed non-adiabatic asteroseismology.

Beware that the term ‘atomic diffusion’ is not always used consistently in the literature. For clarification, in this work we include all physical processes that contribute to atomic diffusion. These are gravitational settling, thermal diffusion, concentration diffusion and radiative levitation. The corresponding abundance changes are caused by pressure gradients, temperature gradients, concentration gradients and radiative forces, respectively. Indeed, we treat atomic diffusion in a multi-component fluid (Burgers, 1969), using diffusion coefficients derived from a screened Coulomb potential (Paquette et al., 1986). We account for partial ionization by using a mean ion charge per element. We calculate radiative forces from first principles using atomic data. Our approach gives a high accuracy treatment of atomic diffusion because it is free of many common assumptions, such as full ionization, simplified treatment of trace elements and approximate radiative forces.

We explain our method in Sect. 2. In particular the computations of stellar evolution (Sect. 2.1), radiative accelerations (Sect. 2.2) and atomic diffusion (Sect. 2.3) are described. In Sect. 3, we present the results for typical sdB models. The effects of atomic diffusion (Sect. 3.1), mass loss (Sect.  3.2) and turbulence (Sect. 3.3) on the abundances and mode stability are evaluated. We summarize and discuss our main results in Sect. 4. We also give a brief discussion of the possible physical origin for turbulence.

2 Computational method

This work presents the first seismic computations of stellar models with a self-consistent treatment of atomic diffusion. The changes in the chemical abundances are accounted for by continuously updating the radiative accelerations and the Rosseland mean opacity, not only during the stellar evolution but also in the pulsation calculations. We describe here the computational tools and methods we use to accomplish this.

2.1 The stellar evolution code

We compute stellar models with the stellar evolution code STARS (Eggleton, 1971). This code has been frequently updated (e.g. Pols et al. 1995) and many derivatives are circulating. In our version the necessary modifications have been made to construct stellar models suitable for seismic studies with the non-adiabatic pulsation code MAD written by Dupret (2001) (see Hu et al. 2008). Furthermore, the implementation of gravitational settling, thermal diffusion and concentration diffusion is described in Hu et al. (2010). In the present work we add the process of radiative levitation. We follow the abundance changes of H, He, C, N, O, Ne, Mg, Fe and Ni while the remaining minor species are kept constant.

2.1.1 Input physics

We take nuclear reaction rates from Angulo & et al. (1999), except for the rate for which we use the recommended value by Herwig et al. (2006) and Formicola & LUNA Collaboration (2002). Neutrino loss rates are according to Itoh et al. (1989); Itoh et al. (1992). Convection is treated with a standard mixing-length prescription (Böhm-Vitense, 1958) with a mixing length to pressure scale height ratio of . Convective mixing is treated as a diffusive process in the framework of mixing length theory.

The Rosseland mean opacity is computed in-line during the evolution to account for composition changes consistently. For this we use data and codes from the Opacity Project (OP, Seaton et al. 1994 and Badnell et al. 2005) which also allows computation of radiative accelerations (see Section 2.2). This is in contrast to all previous work with the STARS code that made use of interpolation in pre-built opacity tables. For high temperatures outside the OP range ( K), we still use OPAL tables (Iglesias & Rogers, 1996) combined with conductive opacities (Cassisi et al., 2007).

2.1.2 Simultaneous solution

The STARS code distinguishes itself from other evolution codes by its fully simultaneous, implicit and adaptive approach. In other words, the equations of stellar structure, composition and mesh are solved together, at each timestep and each iteration, and variables are from the current timestep. Within the scheme of an adaptive mesh, mass loss is simply included by adapting an outer boundary condition. More importantly, we find that such an approach is free from some of the notorious numerical instabilities that accompany diffusion calculations. Such instabilities occur because atomic diffusion, and in particular radiative levitation, can happen on a much shorter time-scale than the evolution of the stellar structure, at least in the outermost layers. In non-simultaneous approaches, where the composition is solved separately from the structure, this can be dealt with by allowing multiple composition timesteps within one evolution timestep (see e.g. Turcotte et al. 1998). In addition, evolution timesteps should remain short if diffusion velocities are kept constant during a timestep. This requires lengthy calculations because the radiative accelerations are computationally expensive to evaluate.

In our calculations, the timesteps are short ( yr) at the start of the evolution but can grow to yr, which is of the same order as for models without diffusion. The timestep size is determined by the change in variables from one timestep to the next. Normally, the temperature and degeneracy have the largest influence but if diffusion acts rapidly the abundance variations limit the timestep size. One might wonder whether the results are sensitive to the timestep size, particularly if mass loss is included. We checked, for one simulation with mass loss, that the results do not change noticeably if timesteps are kept below  yr.

Unfortunately, our numerical scheme requires some effort as well. Like most variables, the radiative accelerations are calculated implicitly; they are not kept constant during a timestep but can change from one iteration to the next. In order to find the correct stellar model, we need to know how a change in affects the stellar structure. In other words, we need the derivatives of each with respect to the variables that represent the stellar model. STARS normally computes the derivatives numerically but this gives instabilities in the case of , probably because of their irregularity. Hence their derivatives (or at least some of them) must be obtained analytically. This is done by the OP routines with some modifications, see Sect. 2.2.

2.2 Radiative accelerations

We compute radiative accelerations and Rosseland mean opacities using OP atomic data. The OP project also provides a set of codes called OPserver (Seaton, 2005; Mendoza et al., 2007) which we have rewritten to suit our specific needs. The most relevant modifications are summarized here.

2.2.1 Radiative accelerations in stellar interiors

For Rosseland mean optical depths , the radiative acceleration for element can be expressed as follows


where is the mean atomic weight, is the atomic weight of element , is the speed of light, is the Rosseland mean opacity, is the luminosity and is the radius. The dimensionless parameter depends on the monochromatic opacity data by

where is the cross-section for absorption or scattering of radiation by element , corrects for electron scattering and momentum transfer to electrons and is the number fraction. Eq. (1) is equivalent to equation (6) of Seaton (2005) but makes use of the luminosity rather than the effective temperature and stellar radius. The latter two surface quantities are not known during the iterations because the evolution code solves the stellar structure from the centre to the surface.

2.2.2 Radiative accelerations for multiple elements

For each call OPserver allows the computation of for one selected element in the mixture and for multiple mixtures with varying abundances of . Because we are interested in the diffusion of nine elements it saves CPU time to compute s for multiple elements but only for one mixture in each call.

2.2.3 Mean ion charges

As mentioned before, elements are treated as if they have an mean ion charge . It has been said by Michaud & Richer (2008) that a disadvantage of OP is that the database does not contain mean charges which are needed to determine the diffusion velocities. We found, however, that can easily be calculated from their data. That should not be too surprising because the level populations of the ionization states are also needed to compute .

2.2.4 Interpolation method

To obtain and for a specified temperature and density (converted to electron density ) OPserver interpolates within a tabulated mesh of () values. OPserver uses two different interpolation schemes. The first is a refined bi-cubic spline interpolation (Seaton, 1993) within the entire () mesh. This scheme is applicable for a single chemical mixture only. Secondly, when the chemical abundances vary with stellar depth, a less refined bi-cubic interpolation occurs within a two-dimensional array of points around the specified values. It should be clear that the second scheme applies if we want to take into account the variation of abundances owing to atomic diffusion. Unfortunately, this scheme is too crude for pulsation studies for which it is required to have smooth opacity derivatives

Therefore, we use the general idea of the second scheme but replace the interpolation method by that developed by Dupret (2003) for the purpose of non-adiabatic pulsation computations. This is a local interpolation method by a polynomial of degree 3. It ensures continuity of and and their derivatives. Note that, in our scheme, interpolation is made within an array of dimension rather than . This increases CPU time, so we only use the third interpolation scheme for the pulsation calculations. For the evolution calculations OPserver’s second scheme suffices. To demonstrate the differences between the three interpolation methods, we perform test runs on the same stellar model as used by Seaton (2005). While (not shown) itself is not visibly affected by the interpolation method, its derivatives (shown in Fig. 1) are.

2.2.5 Computation of derivatives

As explained in Sect. 2.1, analytical derivatives of each are needed in order for the evolution code to converge towards the correct stellar model. It improves the numerical stability to do this for and as well. The derivatives we compute are

The opacity derivatives are already output by OPserver (see Sect. 2.2.4). The others we added. Strictly speaking, the derivatives with respect to and are not analytic. They are however obtained from the interpolaton curve and so have an analytic form.

Figure 1: The opacity derivatives (left y-axis) and (right y-axis) throughout the star as a function of temperature. OP1 and OP2 indicate OPserver’s first and second interpolation scheme, respectively.

2.3 The diffusion equations

Having obtained the radiative accelerations , we put them in Burgers’ diffusion equations (Burgers, 1969),


including the heat flow equations,


In addition, we have two constraints, current neutrality,


and local mass conservation,


In the above equations the quantities , , , and are the partial pressure, mass density, number density, mean charge and mass for species , respectively. The total number of species (including electrons) is . The unknown variables are the diffusion velocities , the heat fluxes , the gravitational acceleration and the electric field . The resistance coefficients , , and are taken from Paquette et al. (1986). We solve Eqs (2.3)–(5) with an adapted version of the routine by Thoul et al. (1994). In Appendix A, we describe how we modify Thoul et al.’s routine to include radiative forces and Paquette et al.’s resistance coefficients. The output of the routine are the diffusion velocities . These are then inserted as an advection term in the diffusion equation governing the time evolution of abundance ,

where is the change owing to nuclear reactions, is the combined diffusion coefficient for convection, overshooting and semi-convection (Eggleton, 1972) and is the diffusion coefficient for turbulent mixing. We take from Michaud et al. (2011),



is an approximation to the He diffusion coefficient at a certain reference depth. The subscript 0 indicates the reference depth in terms of the outer mass coordinate . The assumed is times the He diffusion coefficient at and varying as . Michaud et al. (2011) used a turbulence model with , and but we shall vary these parameters to examine the influence of the efficiency of turbulent mixing on the model (see Sect. 3.3).

3 Results

We investigate the effects of atomic diffusion, mass loss, and turbulence in a typical sdB star with total mass and envelope mass . is the amount of mass above the core boundary and is defined at ZAEHB only, since the H profile can change with time. At the core boundary starts to increase from 0 to 0.7 at the surface. The initial H profile is fitted to the profile on the RGB according to Hu et al. (2009). The ZAEHB model has a and , and it has a metallicity at all depths of with metal mixture of Grevesse & Noels (1993). Although sdB stars could have a range of original metallicities, Michaud et al. (2011)’s diffusion calculations showed that the surface abundances are poorly sensitive to the initial metallicity. Thus we do not expect our results to be much affected by the choice of initial abundances.

Starting at the ZAEHB we compute the stellar evolution until the end of core He burning under different assumptions for mass loss and turbulent mixing. Thereafter, we compute non-adiabatic oscillations of the sdB models. We restrict ourselves to modes of spherical degree , because higher values are geometrically disfavoured for observation. We search for eigenfrequencies in the range of where is the dimensionless angular eigenfrequency, , is the frequency and is the dynamical time-scale.

3.1 Atomic diffusion, no mass loss/turbulence

In Fig. 2 we show radiative accelerations for C, N, O, Ne, Mg, Fe and Ni in the ZAEHB model. H and He are not shown because their are so low that they lie outside the range of the plot. Because gravitational settling and radiative levitation are the dominant diffusion processes, gravity is also plotted for comparison.

Fig. 2 tells us that Fe and Ni are levitated, whereas the other elements sink sooner or later. This is confirmed in Fig. 3 which illustrates the time evolution of the interior abundances from ZAEHB to yr. Notice that Ni accumulates to mass fractions comparable to (or even higher than) Fe in the outer layer, while it starts out with a much smaller mass fraction. This implies that Ni is (at least) as important for mode driving as Fe. Notice also the development of an iron-group convection zone within yr around . The abundance profiles in this region are flattened and do not have a detailed shape as predicted by equilibrium diffusion theory.

Figure 2: Radiative accelerations in the stellar envelope as a function of mass coordinate (bottom axis) and temperature (top axis). Gravitational acceleration is represented by the solid line. This model is at the ZAEHB with and at the surface.
Figure 3: Interior abundance profiles at different stellar ages. The horizontal dashed lines give initial abundances at ZAEHB, and the solid curves are at EHB ages , 5, 6, 7 and 8 as indicated in the legend. Note that in the case of hydrogen we plot for visibility. These simulations are with neither mass loss nor turbulence.
Figure 4: An evolutionary sequence of a model with neither mass loss nor turbulence. Top panel: time evolution of pulsation periods. Red dots indicate unstable pulsation modes while green crosses are stable modes. The shaded areas are the period ranges for the observed - (lower) and -modes (upper). Bottom panel: time evolution of surface abundances. The top x-axis displays the corresponding to the stellar age on the bottom x-axis.

In the upper panel of Fig. 4, we show the computed pulsation periods as a function of stellar age (bottom axis) and (top axis). We also plot the ranges of observed periodicities. We see that this particular simulation represents a hybrid pulsator, showing both unstable - and -modes. The total period range of unstable modes matches the observed range (100 s – 2 hr) well. In the transition region between the - and -modes (400 s – 30 min) the agreement is not as good, but unconfirmed low-amplitude peaks have been observed in that region (Baran et al. 2011 and references therein).

It is important to realize that the low-degree () -modes would not be driven at such high without Ni diffusion. We find that Ni plays an even bigger role in mode driving than first suggested by Jeffery & Saio (2006). In their models only modes, which are unlikely to be observed, are driven at such high s. This is because Ni accumulates more than Fe in the driving region which was not realized before without diffusion calculations. We shall present a detailed investigation of the sdB instability strips elsewhere (Hu et al. in prep.).

The lower panel of Fig. 4 shows the time evolution of the surface abundances. The problem of the He abundance is clearly illustrated because decreases to below within yr, or 0.01 per cent of the EHB lifetime, whereas the average observed value is around (see e.g. Edelmann et al. 2003). The other surface abundances are qualitatively in agreement with observations. C, N, O and Fe are within the observed ranges given by Geier et al. (2010). Ne and Mg are almost completely depleted in our models but this could be consistent with measurements that give only upper limits. Ni is 1.5 dex above solar. This is on the high side but still in agreement with observations that go up to dex above solar (O’Toole & Heber, 2006). For the remainder of this work, we shall discuss only the abundances of He, Fe and Ni, since only these are relevant for our argumentation for constraining the mass loss and turbulent mixing.

3.2 Atomic diffusion and mass loss

The fact that we observe sdB pulsations puts an upper limit on the mass-loss rates, at least for the pulsators. After all, Fig. 4 shows that it takes over yr before enough iron-group elements have accumulated to drive pulsations. The accumulation at this age goes down to in our models (see Fig. 3). So if this outer mass is removed within yr, or equivalently , then there is no chance for the metals to build up.

To verify this estimate we perform simulations with , , and . We assume that the mass loss is constant and chemical homogeneous. The possibilities of a variable and selective mass loss are discussed later. We found that numerical instabilities occur for models that develop an iron-group convection zone and experience mass loss at the same time (see also Charbonneau 1993). This happens for in our models. To reduce instabilities we assume, only in this case, overshooting of the iron-group convection zone to the atmosphere. The extra mixing occurs only in the outer and has a negligible effect. Still, this simulation is terminated after yr because of poor convergence and hence long computing time.

In the upper panels of Fig. 5 we show the effect of these mass-loss rates on the stability of the pulsation modes. We find that, as expected, for pulsations cannot be driven. For pulsations are driven only during the first  yr of evolution. This time-scale is smaller than that found by Fontaine et al. (2006) for a comparable mass-loss rate. It is too short compared to the evolutionary timescale of yr to explain the fraction of sdBs that are variable (about 10 per cent) as suggested by Fontaine et al. (2006). The difference is caused by their assumption that an equilibrium between gravitational settling and radiative levitation can be reached before the onset of mass loss. Furthermore, they did not include diffusion during the evolution while mass is lost. In our models, however, diffusion operates at the same time as mass loss and both processes start at the ZAHB. We find that, under these circumstances, the rates should be if pulsations are to be driven for a significant amount of time.

Remember that, if mass loss were to explain the observed He abundances, then the rates should be in the range (Fontaine & Chayer, 1997; Unglaub & Bues, 2001). This can also be seen in Fig. 5 (lower panels) where the surface abundances of He, Fe and Ni are plotted. Our results show that these mass-loss rates are not consistent with the observed pulsations in sdB stars. One could argue that the constant sdBs are experiencing mass loss while the pulsators are not (or much less). But then there should be a negative correlation between the He abundance and the presence of pulsations. This has not been observed despite many efforts (e.g. O’Toole & Heber 2006; Blanchette et al. 2008).

Another possibility is a variable mass loss during stellar evolution. When the mass-loss rate is low, diffusion can build up enough Fe and Ni to drive pulsations. When the rate increases, the Fe and Ni reservoir is emptied and the star becomes constant. We use the mass loss recipe by Vink & Cassisi (2002) to estimate how much the mass-loss rate could vary. Applying their equation (5) to our models, we find that the rate increases gradually from to during the sdB lifetime. Although the rates themselves could be overestimated due to an underestimation of the terminal velocity (Unglaub, 2008), the dependence on the stellar parameters () should hold (Vink et al., 2000). Thus we expect the increase in the mass-loss rate, which is mostly due to the increasing luminosity as the star evolves111The dependence of on has a small negative effect during the evolution because the metallicity decreases as the lighter metals sink. Although one might expect that the mass-loss rate increases with [Fe/H], one should not forget that the relative role of Fe, compared to lighter elements, diminishes for weaker winds (Vink et al., 2001). A detailed study that evaluates the contribution of different metals to the mass loss is needed but beyond the scope of this work., to remain approximately valid. This factor 5 over the sdB lifetime is not sufficient to obtain the rates needed to keep He from settling () if starting with a rate low enough to drive pulsations (). Hence we find it unlikely that mass loss is the dominant process for slowing down atomic diffusion in sdB stars, although a more sophisticated understanding of mass loss is required to draw any definite conclusions.

Fig. 5 also shows that, in the presence of mass loss, the surface abundances are poor indicators of whether a star is pulsating. For example, during the first  yr of the simulation with , both Fe and Ni reach overabundances above 1 dex while no modes are excited. This is because, although the surface is enriched, Fe en Ni are actually depleted in the driving region.222The driving is caused by the iron-group opacity bump at . This corresponds to in our stellar models. In Fig. 6 we plot the time evolution of the interior abundances of Fe and Ni. Basically what happens is that, as the outermost layers of the star are removed, the regions underneath become exposed. Thus as time progresses and more mass is peeled away, the surface abundances are determined by processes that occurred deeper in the star. After yr the surface is enriched in Fe and Ni because the region where accumulation took place earlier () is now exposed. However, the driving region at this time corresponds to a region that was previously deeper inside the star and from which Fe and Ni have been transported away. After yr so much mass is stripped away that the depleted regions are now at the surface, causing underabundances of Fe and Ni. The actual situation is more complicated because diffusion operates at the same time as mass is lost and the diffusion velocities change as the star evolves.

Figure 5: Models with mass loss. Upper panels show pulsations as a function of time for mass-loss rates (a) , (b) , (c) and (d) . The red dots are unstable modes while the green crosses are stable. Lower panels show surface abundances of He, Fe, and Ni corresponding to these mass-loss rates. The kink in the surface abundances for at yr is caused by the overshoot of the iron-group convection zone assumed in this particular simulation.
Figure 6: Interior abundance profiles of Fe and Ni at different stellar ages for the simulation with . Lines are as in Fig. 3.

3.3 Atomic diffusion and turbulent mixing

If mass loss is not responsible for retarding atomic diffusion then perhaps turbulence is. Michaud et al. (2011) showed that their turbulence model (Eq. 6) with , and M333This results in the outer M of the star being completely mixed. The mixed mass varies slightly with atomic species. can explain most observed abundance anomalies in field and globular cluster sdBs. In their models Fe is near solar throughout the entire mixed region, including the driving region. This raises doubt as to whether pulsations can be excited. As they conclude themselves, their work has yet to be tested with asteroseismology. For this purpose, we run simulations with , and , corresponding to efficient mixing of the outer , and , respectively. The other parameters are kept at and .

In Fig. 7a–c we show the effect of these three turbulence models on the stability of the pulsation modes. We see that mixing the outer mass of M allows some modes to be driven. Although Fe stays near solar, Ni can still accumulate and facilitate mode excitation. However, the He surface abundance still drops too quickly to explain the observations. Mixing more mass ( M) slows down He settling further but it also prevents the accumulation of both Fe and Ni and hence mode excitation. Interestingly, if mixing occurs down to M a few pulsation modes can be excited after yr. This can be understood from Fig. 3 where we see a second, very small accumulation of Fe and Ni appearing around M) at late stellar ages. However, the number of unstable modes is very small and cannot explain the many (dozens) frequencies observed in sdB stars.

Figure 7: Models with turbulent mixing. Upper panels show pulsations as a function of time for turbulence models (a) , , (b) , , , (c) , , and (d) , , . Lower panels show surface abundances of He, Fe, and Ni corresponding to these turbulence models. Symbols and lines are as in Fig. 5.

It would at first sight seem that turbulence cannot be reconciled with mode excitation either. However, this can be solved by adjusting the efficiency of turbulence because the amount of mixing can vary with atomic species. After all, the diffusion velocities of Fe and Ni are much larger than for He, so a weak amount of turbulence could retard He settling while having only a small effect on Fe and Ni levitation. This possibility is investigated by varying the parameters , and in Eq. (6). Indeed, we find for , and that He settles from to during the sdB lifetime while Fe and Ni can accumulate enough to excite pulsations (see Fig. 7d). Efficient mixing only occurs in the outer , but weak mixing occurs up to . This can be seen by the flattening of the abundance profiles plotted in Fig. 8. Of course, it cannot be guaranteed that the consistency between abundance predictions and observations found by Michaud et al. (2011) can be recovered. However, their turbulence model was fine-tuned to keep the Fe abundance near solar, and this remains true in our model.

Figure 8: Interior abundance profiles of Fe and Ni at different stellar ages for the simulation with , , . Lines are as in Fig. 3.

In Table 1 we summarize our results in terms of the fraction of its lifetime the sdB star has (i) He surface abundances in the typically observed range of , (ii) unstable pulsation modes or (iii) both. We see that, of all the scenarios we examined, only the sdB models with mixing in the outer spend a significant time fulfilling both conditions.

Simulation M0T0 M1 M2 M3 M4 T1 T2 T3 T4
(i) He condition 4 5 6 100 7 50 100 65
(ii) pulsation condition 0 90 0 40
(iii) both conditions 0 4 0 5 0 40 64

Simulation (M0T0) is with neither mass loss nor turbulence. The simulations with mass loss are (M1) , (M2) , (M3) and (M4) . The simulations with turbulence are (T1) , , , (T2) , , , (T3) , , and (T4) , , . The numbers give the percentage of its life the sdB star spends with (i) He surface abundance in the range , (ii) unstable pulsation modes or (iii) both conditions fulfilled. Note that simulation (M1) was terminated after yr owing to poor convergence, so we only provide an estimate here.

Table 1: Percentage of EHB lifetime ( yr) spent with He surface abundance in the observed range and/or with pulsation modes excited.

4 Summary and discussion

We have constructed sdB evolutionary and seismic models with atomic diffusion, including the often neglected radiative forces. Our aim is to discriminate between mass loss and turbulence in sdB stars using non-adiabatic asteroseismology. We performed a stability analysis on stellar models with various mass-loss rates and turbulence models. We found that the mass-loss rates required to match the observed He abundances are not consistent with observed pulsations. However, weak turbulent mixing of the outer can also explain the He abundances while still allowing pulsation modes to be driven. The presence of unstable modes is very sensitive to the turbulence model. This could explain why some sdBs are pulsating while others, with similar and , are constant.

Although our results favour turbulence over mass loss, we should remain cautious because the turbulence is not yet supported by a physical model. Possibly, thermohaline mixing in the presence of a mean molecular weight inversion could play a role (Théado et al., 2009). But a detailed investigation including the effect of radiative accelerations on the -inversion instability is still missing. Also, the possibility of overshooting beyond the surface convection zones of sdB stars should be examined with convection models that are more realistic than mixing-length theory. Such studies for A stars show that overshoot can cause the regions between surface convection zones to completely mix (Kupka & Montgomery 2002; Freytag & Steffen 2004). Our models based on the mixing-length theory develop a relatively broad iron-group convection zone on a very short time-scale. Another candidate for causing turbulent mixing could be rotation. Although sdB stars typically have slow surface rotation (), they might hide a rapidly rotating core, a relic from their previous evolution (Kawaler & Hostler, 2005). The differential rotation could cause shear turbulence but this has never been examined in sdB stars as far as we know. We intend to study the influence of these additional mixing processes in a forthcoming paper.

In this study we used observed He abundances to constrain mass-loss rates and turbulent mixing. Metals are not trustworthy in this respect because of the possibility of weak metallic winds () in sdB stars as suggested by Unglaub (2008). Such low mass-loss rates could still act to change the atmospheric metal abundances without interfering with mode driving. This complicates or even prevents the discrimination between variable and constant sdBs on the basis of atmospheric abundances, in agreement with hitherto unsuccessful attempts (O’Toole et al., 2004; O’Toole & Heber, 2006; Blanchette et al., 2008).

It is noteworthy that Ni is as important for mode driving as Fe. Indeed, the diffusive equilibrium abundances in the driving region are comparable despite nickel’s lower initial abundance. Thus Ni plays an even more important role than previous studies of Jeffery & Saio (2006) and Hu et al. (2009) indicated because there it was assumed that Ni was enhanced by the same factor as Fe. Our models show excitation of low-degree () -modes at relatively high ( K), so it is likely that the blue-edge problem can be resolved. A detailed study of the instability strips is underway (Hu et al. in preparation).


We thank J. J. Eldridge and S. Théado for helpful discussions. We also like to thank the referee for useful comments. HH is supported by the Netherlands Organisation for Scientific Research (NWO).


  • Angulo & et al. (1999) Angulo C., et al. 1999, Nuclear Physics A, 656, 3
  • Badnell et al. (2005) Badnell N. R., Bautista M. A., Butler K., Delahaye F., Mendoza C., Palmeri P., Zeippen C. J., Seaton M. J., 2005, \mnras, 360, 458
  • Baran et al. (2011) Baran A. S., Kawaler S. D., Reed M. D., Quint A. C., O’Toole S. J., Østensen R. H., Telting J. H., Silvotti R., Charpinet S., Christensen-Dalsgaard J., Still M., Hall J. R., Uddin K., 2011, \mnras, p. 851
  • Blanchette et al. (2008) Blanchette J., Chayer P., Wesemael F., Fontaine G., Fontaine M., Dupuis J., Kruk J. W., Green E. M., 2008, \apj, 678, 1329
  • Böhm-Vitense (1958) Böhm-Vitense E., 1958, \zap, 46, 108
  • Brown et al. (1997) Brown T. M., Ferguson H. C., Davidsen A. F., Dorman B., 1997, \apj, 482, 685
  • Burgers (1969) Burgers J. M., 1969, Flow Equations for Composite Gases. New York: Academic Press, 1969
  • Cassisi et al. (2007) Cassisi S., Potekhin A. Y., Pietrinferni A., Catelan M., Salaris M., 2007, \apj, 661, 1094
  • Charbonneau (1993) Charbonneau P., 1993, \apj, 405, 720
  • Charpinet et al. (1997) Charpinet S., Fontaine G., Brassard P., Chayer P., Rogers F. J., Iglesias C. A., Dorman B., 1997, \apjl, 483, L123
  • Charpinet et al. (1996) Charpinet S., Fontaine G., Brassard P., Dorman B., 1996, \apjl, 471, L103
  • Chayer et al. (2004) Chayer P., Fontaine G., Fontaine M., Lamontagne R., Wesemael F., Dupuis J., Heber U., Napiwotzki R., Moehler S., 2004, \apss, 291, 359
  • Dupret (2001) Dupret M. A., 2001, \aap, 366, 166
  • Dupret (2003) Dupret M.-A., 2003, PhD thesis, Université de Liège
  • Edelmann et al. (2003) Edelmann H., Heber U., Hagen H., Lemke M., Dreizler S., Napiwotzki R., Engels D., 2003, \aap, 400, 939
  • Edelmann et al. (2006) Edelmann H., Heber U., Napiwotzki R., 2006, Baltic Astronomy, 15, 103
  • Eggleton (1971) Eggleton P. P., 1971, \mnras, 151, 351
  • Eggleton (1972) Eggleton P. P., 1972, \mnras, 156, 361
  • Fontaine et al. (2003) Fontaine G., Brassard P., Charpinet S., Green E. M., Chayer P., Billères M., Randall S. K., 2003, \apj, 597, 518
  • Fontaine & Chayer (1997) Fontaine G., Chayer P., 1997, in Philip A. G. D., Liebert J., Saffer R., Hayes D. S., eds, The Third Conference on Faint Blue Stars: The Helium Abundance Puzzle in Subdwarf B Stars: a Theoretical Interpretation Through Mass Loss. L. Davis Press, p. 169
  • Fontaine et al. (2006) Fontaine G., Green E. M., Chayer P., Brassard P., Charpinet S., Randall S. K., 2006, Baltic Astronomy, 15, 211
  • Formicola & LUNA Collaboration (2002) Formicola A., LUNA Collaboration 2002, in W. Hillebrandt & E. Müller ed., Nuclear Astrophysics Re-investigation of the N(p,)O reaction within LUNA collaboration. p. 111
  • Freytag & Steffen (2004) Freytag B., Steffen M., 2004, in J. Zverko, J. Ziznovsky, S. J. Adelman, & W. W. Weiss ed., The A-Star Puzzle Vol. 224 of IAU Symposium, Numerical simulations of convection in A-stars. p. 139
  • Geier et al. (2010) Geier S., Heber U., Edelmann H., Morales-Rueda L., Napiwotzki R., 2010, \apss, 329, 127
  • Green et al. (2003) Green E. M., Fontaine G., Reed M. D., Callerame K., Seitenzahl I. R., White B. A., Hyde E. A., Østensen R., Cordes O., Brassard P., Falter S., Jeffery E. J., Dreizler S., Schuh S. L., Giovanni M., Edelmann H., Rigby J., Bronowska A., 2003, \apjl, 583, L31
  • Grevesse & Noels (1993) Grevesse N., Noels A., 1993, in N. Prantzos, E. Vangioni-Flam, & M. Casse ed., Origin and Evolution of the Elements: Cosmic abundances of the elements. p. 15
  • Han et al. (2003) Han Z., Podsiadlowski P., Maxted P. F. L., Marsh T. R., 2003, \mnras, 341, 669
  • Heber (1986) Heber U., 1986, \aap, 155, 33
  • Heber et al. (2003) Heber U., Maxted P. F. L., Marsh T. R., Knigge C., Drew J. E., 2003, in I. Hubeny, D. Mihalas, & K. Werner ed., Stellar Atmosphere Modeling Vol. 288 of Astronomical Society of the Pacific Conference Series, Stellar Wind Signatures in sdB Stars?. p. 251
  • Herwig et al. (2006) Herwig F., Austin S. M., Lattanzio J. C., 2006, \prc, 73, 025802
  • Hu et al. (2008) Hu H., Dupret M., Aerts C., Nelemans G., Kawaler S. D., Miglio A., Montalban J., Scuflaire R., 2008, \aap, 490, 243
  • Hu et al. (2010) Hu H., Glebbeek E., Thoul A. A., Dupret M., Stancliffe R. J., Nelemans G., Aerts C., 2010, \aap, 511, A87
  • Hu et al. (2009) Hu H., Nelemans G., Aerts C., Dupret M., 2009, \aap, 508, 869
  • Iglesias & Rogers (1996) Iglesias C. A., Rogers F. J., 1996, \apj, 464, 943
  • Itoh et al. (1989) Itoh N., Adachi T., Nakagawa M., Kohyama Y., Munakata H., 1989, \apj, 339, 354
  • Itoh et al. (1992) Itoh N., Mutoh H., Hikita A., Kohyama Y., 1992, \apj, 395, 622
  • Jeffery & Saio (2006) Jeffery C. S., Saio H., 2006, \mnras, 372, L48
  • Kawaler & Hostler (2005) Kawaler S. D., Hostler S. R., 2005, \apj, 621, 432
  • Kilkenny et al. (1997) Kilkenny D., Koen C., O’Donoghue D., Stobie R. S., 1997, \mnras, 285, 640
  • Kupka & Montgomery (2002) Kupka F., Montgomery M. H., 2002, \mnras, 330, L6
  • Mendoza et al. (2007) Mendoza C., Seaton M. J., Buerger P., Bellorín A., Meléndez M., González J., Rodríguez L. S., Delahaye F., Palacios E., Pradhan A. K., Zeippen C. J., 2007, \mnras, 378, 1031
  • Michaud & Richer (2008) Michaud G., Richer J., 2008, \memsai, 79, 592
  • Michaud et al. (2011) Michaud G., Richer J., Richard O., 2011, \aap, 529, A60
  • O’Toole & Heber (2006) O’Toole S. J., Heber U., 2006, \aap, 452, 579
  • O’Toole et al. (2004) O’Toole S. J., Heber U., Chayer P., Fontaine G., O’Donoghue D., Charpinet S., 2004, \apss, 291, 427
  • Paquette et al. (1986) Paquette C., Pelletier C., Fontaine G., Michaud G., 1986, \apjs, 61, 177
  • Pols et al. (1995) Pols O. R., Tout C. A., Eggleton P. P., Han Z., 1995, \mnras, 274, 964
  • Seaton (1993) Seaton M. J., 1993, \mnras, 265, L25
  • Seaton (2005) Seaton M. J., 2005, \mnras, 362, L1
  • Seaton et al. (1994) Seaton M. J., Yan Y., Mihalas D., Pradhan A. K., 1994, \mnras, 266, 805
  • Théado et al. (2009) Théado S., Vauclair S., Alecian G., Le Blanc F., 2009, \apj, 704, 1262
  • Thoul et al. (1994) Thoul A. A., Bahcall J. N., Loeb A., 1994, \apj, 421, 828
  • Turcotte et al. (1998) Turcotte S., Richer J., Michaud G., Iglesias C. A., Rogers F. J., 1998, \apj, 504, 539
  • Unglaub (2008) Unglaub K., 2008, \aap, 486, 923
  • Unglaub & Bues (2001) Unglaub K., Bues I., 2001, \aap, 374, 570
  • Vink (2004) Vink J. S., 2004, \apss, 291, 239
  • Vink & Cassisi (2002) Vink J. S., Cassisi S., 2002, \aap, 392, 553
  • Vink et al. (2000) Vink J. S., de Koter A., Lamers H. J. G. L. M., 2000, \aap, 362, 295
  • Vink et al. (2001) Vink J. S., de Koter A., Lamers H. J. G. L. M., 2001, \aap, 369, 574

Appendix A Solving the Burgers’ equations

With the approach of Thoul et al. (1994), the Burgers’ equations without radiative forces, together with the constraints of current neutrality and mass conservation, can be expressed as a matrix equation

Most quantities are defined as by Thoul et al. (1994) except for the coefficients and which we define as

Our now takes into account the He concentration gradient which was unnecessarily eliminated by Thoul et al. (1994). The coefficients are modified to make use of Paquette et al. (1986)’s resistance coefficients (, , and ) derived from a screened Coulomb potential. For we have

For ,

with , , and . Note that the approximations of Thoul et al. (1994) are retrieved with , and as estimated from a pure Coulomb potential.

The above modifications were already made by Hu et al. (2010) but not explicitly described then. In this work we adapted Thoul et al. (1994)’s routine to include radiative forces. After some algebra we find that Eqs (2.3)–(5) can be written as the matrix equation


The terms on the left-hand side represent the contributions to the diffusion velocity by radiative levitation, gravitational settling, thermal diffusion and concentration diffusion. The matrix equation is solved by LU decomposition for each of these terms separately and the solutions are combined to give to the total diffusion velocity. In principle, we could also add the terms on the left-hand side of Eq. (7) beforehand and solve the matrix equation in one step. However, it can be insightful to quantify the different contributions to the diffusion velocities.

Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description