© ESO, 2017

1. Introduction

The asymptotic giant brach (AGB; Herwig 2005; Karakas & Lattanzio 2014) is occupied by low- and intermediate-mass (0.8 ≤ M ≤ 8 M_⊙) stars in the last nuclear-burning phase. At the end of the AGB phase, these stars develop thermal pulses (TP) and suffer extreme mass loss. AGB stars are thus one of the main contributors to the enrichment of the interstellar medium (ISM) of light elements (e.g. Li, C, N, F) and heavy (slow neutron capture, s-process) elements and thus to the chemical evolution of galaxies (Busso et al. 1999). AGB stars are also one of the most prominent sources of dust in galaxies and the site of origin of the vast majority of meteoritic stardust grains (e.g. Hoppe & Ott 1997; Nittler et al. 1997; Lugaro et al. 2017). In low-mass AGB stars (M< 4M_⊙) ¹²C is produced during the TP-AGB phase, and carried to the stellar surface via the third dredge-up (TDU) that can occur after each TP, transforming originally O-rich stars into C-rich stars (C/O > 1) (e.g. Herwig 2005; Karakas & Lattanzio 2007; Lugaro & Chieffi 2011). However, the more massive AGB stars (M> 4–5 M_⊙) are O-rich (C/O < 1) because the so-called hot bottom burning (HBB) process is activated. Through the CN cycle via proton captures at the base of the convective envelope, HBB converts ¹²C into ¹³C and ¹⁴N, thus preventing the formation of a carbon star (Sackmann & Boothroyd 1992; Mazzitelli et al. 1999).

The s-process allows the production of neutron-rich elements heavier than iron (s-elements such as Sr, Y, Zr, Ba, La, Nd, Tc, etc.) by the s-process. In the low-mass AGB stars (roughly < 4 M_⊙), the ¹³C(α, n)¹⁶O reaction is the dominant neutron source (e.g. Abia et al. 2001). In the more massive AGB stars, instead, neutrons are mainly released by the ²²Ne(α, n)²⁵Mg reaction, resulting in a higher neutron density (up to 10¹³ n/cm³) and temperature environment than in lower-mass AGB stars (García-Hernández et al. 2006). The Rb produced depends on the probability of the ⁸⁵Kr and ⁸⁶Rb capturing a neutron before decaying and acting as “branching points” (see van Raai et al. 2012, for more details). The probability of this happening depends on the local neutron density (Beer & Macklin 1989). The ⁸⁷Rb/⁸⁵Rb isotopic ratio is a direct indicator of the neutron density at the production site but it is not possible to distinguish individual ⁸⁷Rb and ⁸⁵Rb from stellar spectra (García-Hernández et al. 2006). However, the relative abundance of Rb to other nearby s-process elements such as Zr is very sensitive to the neutron density, and so a good discriminant of the stellar mass and the neutron source at the s-process site (Lambert et al. 1995; Abia et al. 2001; García-Hernández et al. 2006; van Raai et al. 2012). In other words, [Rb/Zr] < 0 is observed in low-mass AGB stars where the main neutron source is the ¹³C(α, n)¹⁶O reaction (Plez et al. 1993; Lambert et al. 1995; Abia et al. 2001), while [Rb/Zr] > 0 is observed in more massive AGB stars, where the neutrons are mainly released through the ²²Ne(α, n)²⁵Mg reaction (García-Hernández et al. 2006, 2007, 2009).

Chemical abundance analyses using classical MARCS hydrostatic atmospheres (Gustafsson et al. 2008) revealed strong Rb overabundances (~10³–10⁵ times solar) and high [Rb/Zr] ratios (3–4 dex) in massive AGB stars (generally very luminous OH/IR stars) of our own Galaxy and the Magellanic Clouds (MC; García-Hernández et al. 2006, 2007, 2009). This observationally confirmed for the first time that the ²²Ne neutron source dominates the production of s-process elements in these stars. However, the extremely high Rb abundances and [Rb/Zr] ratios observed in most the massive stars (and especially in the lower-metallicity MC AGB stars) have posed a “Rb problem”; such extreme [Rb/Fe] and [Rb/Zr] values are not predicted by the s-process AGB models, (van Raai et al. 2012; Karakas et al. 2012), suggesting fundamental problems in our present understanding of AGB nucleosynthesis and/or of the complex extended dynamical atmospheres of these stars (García-Hernández et al. 2009).

Zamora et al. (2014) constructed new pseudo-dynamical MARCS model atmospheres by considering the presence of a gaseous circumstellar envelope with a radial wind and applied them to a small sample of five O-rich AGB stars with different expansion velocities and metallicities. The Rb abundances and [Rb/Zr] ratios obtained were much lower than those obtained with classical hydrostatic models; in better agreement with the AGB nucleosynthesis theoretical predictions. In this paper, we use the Zamora et al. (2014) pseudo-dynamical model atmospheres to obtain the abundances of Rb and Zr in the full sample of massive Galactic AGB stars previously analysed with hydrostatic models (García-Hernández et al. 2006, 2007). These Rb and Zr abundances are then compared with the more recent AGB nucleosynthesis theoretical predictions available in the literature.

2. Sample and observational data

Our sample is composed of 21 massive Galactic AGB stars (most of them very luminous OH/IR stars) previously analysed by García-Hernández et al. (2006, 2007); we use their high-resolution (R ~ 40 000–50 000) optical echelle spectra (we refer to García-Hernández et al. 2006, 2007, for further observational details)¹. The signal-to-noise ratios (S/N) achieved in the reduced spectra vary greatly from the blue to the red (typically ~10–20 at 6000 Å and >100 at 8000 Å). The Rb and Zr abundances were determined from the resonant 7800 Å Rb I line and the 6474 Å ZrO bandhead, respectively, by using classical MARCS hydrostatic model atmospheres (García-Hernández et al. 2006, 2007). The Rb abundances and [Rb/Zr] ratios obtained from this chemical analysis are mostly in the range [Rb/Fe] ~ 0.6–2.6 dex and [Rb/Zr] ~ 0.1–2.1 dex. The atmospheric parameters and Rb abundances derived with the hydrostatic models as well as other useful observational information like the OH expansion velocity, variability period, and the presence of Li are listed in Table 1.

3. Chemical abundance analysis using pseudo-dynamical models

3.1. Modified version of the Turbospectrum spectral synthesis code

We used the v12.2 version of the spectral synthesis code Turbospectrum (Alvarez & Plez 1998; Plez 2012), which considers the presence of a circumstellar gas envelope and a radial wind, as modified by Zamora et al. (2014). The main modifications are the following: (i) The Doppler effect due to the extended atmosphere and velocity field is introduced in the routines that compute the line intensities at the stellar surface; (ii) the source function of the radiative transfer is assumed to be the same as computed in the static case (Gustafsson et al. 2008). The validity of this approximation was tested by comparing with Monte Carlo simulations (see Zamora et al. 2014); (iii) The scattering term of the source function (∝σ_λJ_λ) is not shifted to save computing time and is only incorporated for the continuum. This scattering term is computed as in the static case using the Feautrier method (Nordlund 1984; Gustafsson et al. 2008); and (iv) the velocity field is taken into account through a shift of the absorption coefficient κ_λ; the source function is built using the static σ_λJ_λ and the shifted κ_λB_λ. The emerging intensity is then computed in the observer frame by a direct quadrature of the source function.

3.2. Extended atmosphere models

For the analysis of each star in our sample, we adopted the atmosphere parameters from García-Hernández et al. (2006, 2007) and the solar reference abundances by Grevesse et al. (2007). We constructed our pseudo-dynamical models from the original MARCS hydrostatic atmosphere model structure. We expanded the atmosphere radius by a wind out to ~5 stellar radii and a radial velocity field (Zamora et al. 2014). In the MARCS hydrostatic model, the R_∗ is the radius corresponding to r(τ_Ross = 1), where r is the distance from the centre of the star and τ_Ross is the Rosseland optical depth. We computed the stellar wind following the mass conservation (Eq. (1)), radiative thermal equilibrium (Eq. (2)) and a classical β-velocity law (Eq. (3)), $\begin{array}{ccc}\mathit{\rho }\mathrm{\left(}\mathit{r}\mathrm{\right)}& \mathrm{=}& \frac{\mathit{\dot{}}\mathit{M}}{\mathrm{4}\mathit{\pi }{\mathit{r}}^{\mathrm{2}}\mathit{v}\mathrm{\left(}\mathit{r}\mathrm{\right)}}\mathit{,}\\ \mathit{r}{\mathit{T}}^{\mathrm{2}}& \mathrm{=}& \mathrm{constant}\mathrm{=}{\mathit{r}}_{\mathrm{out}}{\mathit{T}}_{\mathrm{out}}^{\mathrm{2}}\mathit{,}\\ \mathit{v}\mathrm{\left(}\mathit{r}\mathrm{\right)}& \mathrm{=}& {\mathit{v}}_{\mathrm{0}}\mathrm{+}\mathrm{(}{\mathit{v}}_{\mathrm{\infty }}\mathrm{-}{\mathit{v}}_{\mathrm{0}}\mathrm{)}{\left(\mathrm{1}\mathrm{-}\frac{{\mathit{R}}_{\mathrm{\ast }}}{\mathit{r}}\right)}^{\mathit{\beta }}\mathit{,}\end{array}$where ρ(r) is the density of the envelope radius r, Ṁ is the mass-loss rate, and v(r) is the velocity of the envelope, which is calculated by means of Eq. (3). In Eq. (3), v₀ is a reference velocity for the beginning of the wind and β is an arbitrary free parameter. We take v₀ = v(R_∗) for the onset of the wind and the extension of the envelope begins from the outer radius of the hydrostatic model. Using Eq. (2) the envelope is extended, layer by layer, out to the distance r_max, which corresponds to the maximum radius in our calculations, with T_min = 1000 K. Turbospectrum cannot compute lower temperatures due to numerical reasons (Zamora et al. 2014).

3.3. Resulting grids of synthetic spectra

The synthetic spectra are generated with the modified version of Turbospectrum by using the extended pseudo-dynamical model atmospheres as input. We constructed a mini-grid of synthetic spectra for each sample star by adopting the atmospheric parameters (e.g. effective temperature, macroturbulence²) from García-Hernández et al. (2006, 2007). Basically, the stellar mass, gravity log g, microturbulent velocity ξ, metallicity [Fe/H], and C/O ratio are fixed to 2 M_⊙, −0.5 dex, 3 km s^-1, 0.0, and 0.5 dex, respectively (we refer to García-Hernández et al. 2007, for more details). On the other hand, for the mass-loss rate Ṁ and the exponent β, we use values in the range Ṁ ~ 10^-9–10^-6M_⊙ yr^-1 in steps of 0.5 × 10^-1M_⊙ yr^-1 and the range β ~ 0.2–1.6 in steps of 0.2. We have not considered the case where β = 0.0 because the expansion velocity would be constant at any r. We assume the OH expansion velocity (v_exp(OH); see Table 1) as the terminal velocity because the OH maser emission is found at very large distances from the central star (e.g. Decin et al. 2010). Figure 1 shows examples of the β-velocity laws used in our pseudo-dynamical models based on the MARCS hydrostatic models. Finally, for the Rb and Zr abundances, we used [Rb/Fe] ~−2.6 to +3.0 dex, and [Zr/Fe] ~−1.0 to +1.0 in steps of 0.1 and 0.25 dex, respectively.

Fig. 1 Velocity vs. distance from the star in four of our AGB wind models. These velocity laws present different expansion velocities vexp(OH), mass-loss rates and β exponents. The effective temperature Teff = 3000 K, gravity log g = −0.5 and the solar chemical composition are the same in all models.

The resulting mini-grid (~4400 models) is compared to the observed spectrum in order to find the synthetic spectrum that best fits the 7800 Å Rb I line and the 6474 Å ZrO bandhead profiles and their adjacent pseudocontinua. In order to obtain the best fits, we made use of a procedure based on the comparison between synthetic and observed spectra, while in Zamora et al. (2014) the observed spectra were fitted by eye. The method is a modified version of the standard χ² test, ${\mathit{\chi }}^{\mathrm{2}\mathrm{\ast }}\mathrm{=}{\mathit{\chi }}^{\mathrm{2}}\mathrm{\times }\mathit{w}\mathrm{=}\left(\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{N}}\frac{\mathrm{[}\mathit{Yob}{\mathit{s}}_{\mathit{i}}\mathrm{-}\mathit{Ysynt}{\mathit{h}}_{\mathit{i}}\mathrm{(}{\mathit{x}}_{\mathrm{1}}\mathit{...}{\mathit{x}}_{\mathit{M}}\mathrm{)}{\mathrm{]}}^{\mathrm{2}}}{\mathit{Yob}{\mathit{s}}_{\mathit{i}}}\right)\mathrm{\times }\mathit{w,}$(4)where Yobs_i and Ysynth_i are the observed and synthetic data points, respectively, with N the number of data points, and M the number of free parameters. On the other hand, w is a vector that gives a stronger weight to the detailed spectral profiles of the Rb I line and the ZrO bandhead. This way, the lowest value of χ^{2 ∗} gives us the best fitting synthetic spectrum from the mini-grid for each sample star.

The use of the χ^{2 ∗} test to find the best fits to the observed spectra reveals the presence of important degeneracies in the resulting grids of pseudo-dynamical synthetic spectra; that is, very similar synthetic spectra are obtained from different sets of wind parameters (more details below). Moreover, in a some cases (IRAS 04404−7427, IRAS 05027−2158, IRAS 05098−6422, IRAS 06300+6058, IRAS 10261−5055, IRAS 18429−1721, IRAS 19059−2219 and IRAS 20343−3020) the use of the χ^{2 ∗} test is not enough for obtaining the synthetic spectrum that best reproduces the observed one and the best fits have to be found by eye. Unfortunately, the wind model parameters Ṁ and β are generally not known for stars in our sample (see below), complicating the abundance analysis. Thus, here we study the sensitivity of the synthetic spectra and the abundance results to variations of the stellar and wind parameters.

3.4. Sensitivity of the synthetic spectra to variations of the model parameters

Fig. 2 Illustrative examples of synthetic spectra for different stellar (Teff) and wind (Ṁ, β, and vexp(OH)) parameters in the spectral region around the 7800 Å Rb I line. The black vertical line indicates the position of the 7800 Å Rb I line.

Fig. 3 Illustrative examples of synthetic spectra for different stellar (Teff) and wind (Ṁ, β, and vexp(OH)) parameters in the spectral region around the 6474 Å ZrO bandhead. The black vertical line indicates the position of the 6474 Å ZrO bandhead.

Here, we analyse how the variations in stellar (T_eff) and wind (Ṁ, β and v_exp(OH)) parameters influence the output synthetic spectra. Figures 2 and 3 show examples of synthetic spectra for different stellar and wind parameters in the spectral regions around the 7800 Å Rb I line and 6474 Å ZrO bandhead, respectively. We note that the fraction of the absortion at 7800 Å due to other species (e.g. TiO) is typically around 20%.

The Rb I line profile is very sensitive to the wind mass-loss rate Ṁ (especially for Ṁ ≥ 10^-8M_⊙ yr^-1); the Rb I line is significantly deeper and blue-shifted with increasing Ṁ (Fig. 2, top-left panel). However, the Rb I line profile is much less sensitive to changes of the wind velocity-law (β parameter); being only slightly deeper with increasing β (Fig. 2, bottom-left panel). In addition, for β values higher than ~1.2 (shallower velocity profiles), the Rb I line profile is not sensitive to variations of the expansion velocity v_exp(OH) because the velocity profiles are very similar in our extended model atmosphere (up to ~10¹⁴ cm; see Fig. 1). Variations in the expansion velocity v_exp(OH) mainly affect the blue-shift of the Rb I line and, in addition, for large v_exp(OH) values the core of the Rb I line is less deep (Fig. 2, bottom-right panel). Finally, the Rb I absorption line is stronger with decreasing effective temperature T_eff (as expected; Fig. 2, top-right panel) but this time the wealth of TiO molecular lines and the pseudo-continua are also affected. We note that all these effects (variations in the Rb I profile in terms of depth and blue-shift) are more evident for extreme mass-loss rates (Ṁ ≥ 10^-7M_⊙ yr^-1) and higher Rb abundances.

On the other hand, the ZrO bandhead profile is not sensitive to the wind parameters Ṁ, β, and v_exp(OH) (see Fig. 3; top-left, bottom-left, and bottom-right panels, respectively). The ZrO bandhead profile (as well as the adjacent TiO lines and pseudo-continuum) are, again as expected, stronger with decreasing T_eff (Fig. 3, top-right panel). This is because ZrO is formed deeper than Rb I in the atmosphere, being much less affected by the circumstellar envelope and radial wind.

4. Results

As we have mentioned above, there are important degeneracies in the resulting mini-grids of synthetic spectra for each star. Two synthetic spectra with the same T_eff, log g and v_exp(OH), but different β, Ṁ, and [Rb/Fe] abundances could be practically identical in spite of the different wind parameters. This complicates the abundance analysis because the wind model parameters Ṁ and β are generally not known for stars in our sample. In any case, we can use some observational constraints and previous results on a few similar OH/IR stars to limit the possible variation range of these wind parameters (in particular for the mass-loss rates Ṁ, see below).

By using multiple rotationally excited lines of both ¹²CO and ¹³CO, De Beck et al. (2010) provide accurate mass-loss rates Ṁ for a large sample of Galactic AGB stars. Unfortunately, only one star (IRAS 20077−0625) from our present sample of Rb-rich OH/IR massive AGB stars is included in their work and we cannot fit this star with our pseudo-dynamical models (see below). There are seven massive AGB stars of OH/IR type (WX Psc, V669 Cas, NV Aur, OH 26.5+0.6, OH 44.8−2.3, IRC −10529 and OH 104.9+2.4) previously studied in the optical by García-Hernández et al. (2007). Their variability periods and mass-loss rates range from 552 to 1620 days³ and from 1.8 × 10^-5M_⊙ yr^-1 to 9.7 × 10^-6M_⊙ yr^-1, respectively. Interestingly, all these stars are extremely obscured in the optical, being too red or without optical counterpart⁴; they likely already have entered the superwind phase. Thus, the Ṁ values in optically obscured OH/IR AGB stars can be taken as upper limits (i.e. <10^-6M_⊙ yr^-1) for our sample of OH/IR massive AGB stars with optical counterparts; that is, with useful spectra around the 7800 Å Rb I line. Indeed, we generally find that lower mass-loss rates (~10^-7–10^-8M_⊙ yr^-1) give superior fits to the observed Rb I line profiles. Mass-loss rates of ~10^-6M_⊙ yr^-1 (or higher) give strong Rb I absorption lines for solar Rb abundances (see also Zamora et al. 2014) with the consequence that all stars in our sample of OH/IR massive AGB stars would be Rb-poor. By combining the variability periods from Table 1 and the mass-loss rates estimated from the Rb I line profiles (mainly in the range ~10^-7–10^-8M_⊙ yr^-1; Table 2) into the AGB mass-loss formula by Vassiliadis and Wood (1993, their Eq. (5)), we obtain reasonable current stellar masses in the range ~2.5–6 M_⊙. In Table 2 we show the mass-loss rates obtained from the best spectral fits (Ṁ_fit) and the current stellar masses by using the mass-loss expression from Vassiliadis & Wood (1993).

The β parameter in our models (only up to ~10¹⁴ cm from the photosphere; Fig. 1) cannot be directly compared with other estimations of this parameter in the literature (e.g. Decin et al. 2010; Danilovich et al. 2015), which map many outer regions in the circumstellar envelope and that usually obtain relatively high and uncertain values (0 ≤ β ≤ 5.0). However, the effect of the β parameter on our synthetic spectra is minor compared to the mass-loss rate Ṁ and we keep it as a free parameter in our abundance analysis. We note also that the velocity profiles are very similar in our extended model atmosphere for β ≥ 1.2; that is, the Rb I line profile is not more sensitive to variations of the expansion velocity and the abundance results are very similar for β ≥ 1.2. We generally find better fits with low β values (or steeper velocity profiles; see Table 3).

As mentioned above, the parameters of the hydrostatic models providing the best fit to the observations and the Rb abundances derived are shown in Table 1. The static models use the solar abundances from Grevesse & Sauval (1998) for computing the Rb abundances (García-Hernández et al. 2006, 2007), while our pseudo-dynamical models use the more recent solar abundances from Grevesse et al. (2007). In Zamora et al. (2014) we compared the Rb abundances from static models using Grevesse & Sauval (1998) and Grevesse et al. (2007), and the Rb abundances obtained agree within ~0.2 dex in most cases.

Figure 4 shows that our pseudo-dynamical atmosphere models reproduce the observed 7800 Å Rb I line profile much better than the classical hydrostatic models in four sample stars (we refer to Appendix A for the remainder of the sample stars). On the other hand, the Zr abundances derived from the extended models are similar to those obtained with the hydrostatic models because the 6474 Å ZrO bandhead is formed deeper in the atmosphere and is less affected by the radial velocity field (Zamora et al. 2014). We could obtain the Rb and Zr abundances (or upper limits) for 17 sample stars. The remainder of our sample stars (IRAS 07222−2005, IRAS 09194−4518, IRAS 19426+4342 and IRAS 20077−0625) seem to display different Rb I line profiles (e.g. with more than one circumstellar contribution or anomalously broad profiles with red-extended wings; e.g. Fig. 5) that cannot be completely reproduced by our present version of the spectral synthesis code. In the two stars (IRAS 07222−2005 and IRAS 09194−4518) shown in Fig. 5 we cannot fit the two Rb I components (circumstellar and photospheric) at the same time; for example, we could only partially fit the blue-shifted circumstellar component using larger mass-loss rates (>10^-6M_⊙ yr^-1). Curiously, these two stars present the largest periods (see Table 1) and they may be the most extreme and evolved stars in our sample, where our extended models do not work so well (e.g. due to even more extended atmospheres). It is not completely clear, however, if the observed profiles are real because these four sample stars have the lowest-quality spectra (S/N< 30 at 7800 Å; Table 1).

For the two sample stars with unknown OH expansion velocity, IRAS 16030−5156 and IRAS 17034−1024, we explore the velocity range displayed by other stars with similar variability periods (see Table 1). Similar fits can be obtained for v_exp(OH) ~7−12 and 7−9 km s^-1 (in combination with sligthly different wind parameters) for IRAS 16030−5156 and IRAS 17034−1024, respectively, and we thus adopt average velocities of 10 and 8 km s^-1, respectively, in the abundance analysis (Table 3).

Fig. 4 Rb I 7800 Å (left panel) and ZrO 6474 Å (right panel) spectral regions in massive Galactic AGB stars. The pseudo-dynamical models (red lines) that best fit the observations (black dots) are shown in four sample stars. For comparison, the hydrostatic models are also displayed (blue lines). The location of the Rb I line and the ZrO bandhead are indicated by black vertical lines.

Fig. 5 IRAS 07222−2005 (left panel) and IRAS 09194−4518 (right panel) display two components (circumstellar and photospheric) in the 7800 Å Rb I line (black line) and cannot be reproduced by our present version of the spectral synthesis code. The red line shows the pseudo-dynamical synthetic spectrum that gives a good fit to the photospheric component only.

Table 3 shows the atmospheric and wind parameters as well as the Rb and Zr abundances (or upper limits) from the best fits to the observed spectra when the wind parameters Ṁ and β are not fixed. In most cases, the best fit is obtained for both low β (~0.2) and Ṁ (~10^-9–10^-7M_⊙ yr^-1) values. The new Rb abundances obtained from extended models are lower than those obtained using the hydrostatic models, and the difference is larger for stars with higher hydrostatic Rb abundances. In addition, this difference is smaller for lower v_exp(OH) and increases with increasing v_exp(OH), as expected. On the other hand, in the case of Zr we obtain upper limits mostly between 0.0 and +0.25 dex, as derived from the hydrostatic models. Figure 6 displays the hydrostatic and pseudo-dynamical Rb abundances versus the OH expansion velocity for the wind parameters that provide the best fits (Table 3). We plot the Rb abundances versus the expansion velocity because the v_exp(OH) can be used as a mass indicator independent of the distance in OH/IR stars (García-Hernández et al. 2007). In addition, in Fig. 6 we have marked the Li-rich stars (García-Hernández et al. 2007) by squares. About half of the stars with v_exp(OH)> 6 km s^-1 are Li-rich and most of these stars are also the more Rb-rich ones. We get pseudo-dynamical abundances lower than the hydrostatic ones and a worse correlation between⁵ the Rb abundances and v_exp(OH); the Rb-v_exp(OH) relationship is flatter (with a higher degeneration) for the pseudo-dynamical case (see also Sect. 5).

Fig. 6 Rb abundances derived both with hydrostatic (blue dots) and pseudo-dynamical model atmospheres (red triangles) with best-fit parameters plotted against the OH expansion velocity. The Li-rich stars are indicated by squares. A typical error bar of ±0.7 dex is also displayed.

Also, we carried out several tests with different β and Ṁ values in order to check the sensitivity of the derived abundances to variations of the wind parameters. In Table 4 we present the wind parameters and Rb abundances obtained when fixing β to 0.2 and 1.2. Basically, the Rb abundances are lower in the β = 1.2 case because a higher β deepens the Rb I 7800 Å line for the same Rb abundance (see Fig. 2); in a few cases, however, the Ṁ of the best fit also changes, further affecting the determination of the Rb abundance. The Zr abundances (or upper limits) are similar in most cases; the upper limits only change when the Ṁ is not the same for the β = 0.2 and 1.2 cases. Figure 7 represents the Rb abundances obtained versus v_exp(OH) for β = 0.2 and β = 1.2. By comparing the Rb abundances from Figs. 6 and 7, it is clear that the Rb abundances are slightly lower in the β = 1.2 case. Moreover, the correlation between the pseudo-dynamical Rb abundances and v_exp(OH) for different β values is worse (e.g. flatter) than the hydrostatic case. In addition, the dispersion seems to be larger for the β = 1.2 case.

Fig. 7 Rb abundances derived both with hydrostatic (blue dots) and extended model atmospheres with β = 0.2 (magenta triangles) and β = 1.2 (green squares) plotted against the OH expansion velocity.

On the other hand, Fig. 8 displays the Rb results when Ṁ is fixed to 10^-8, 10^-7 and 10^-6M_⊙ yr^-1. This could be equivalent to considering that our AGB sample stars are at a similar evolutionary stage in terms of mass loss; of course we have a strong degeneracy between the progenitor masses and mass loss/evolutionary stage. In the particular case of the Li-rich AGB stars, statistical arguments suggest that these stars should have a narrow initial mass range (Di Criscienzo et al. 2016); 4–5 or 5–6 M_⊙ according to the most recent ATON (Di Criscienzo et al. 2016) or Monash (Karakas & Lugaro 2016) AGB nucleosynthesis models, respectively. The current stellar masses from Table 2 show, however, that there is a complicated interplay (degeneracy) between Li enhancement, progenitor mass, and mass-loss rate and that the progenitor mass range of these stars may actually be broader; for example, their current stellar mass and Li abundance ranges are ~2.7–5.6 M_⊙ and ~logε(Li)~0.7−2.6 dex. Figure 8 shows that the Rb abundances decrease with increasing Ṁ and the dispersion of the Rb abundances is much lower when fixing Ṁ. The slopes (and correlation coefficients) of the Rb-v_exp(OH) correlations are more similar to the ones obtained with hydrostatic models. The Rb abundances from extended models approach the hydrostatic ones with decreasing Ṁ (both sets of Rb abundances are identical for Ṁ ≤ 10^-9M_⊙ yr^-1) because the atmosphere is less extended with decreasing Ṁ, as expected.

Finally, we fixed Ṁ and β, which could be equivalent to considering that our AGB sample stars have the same mass-loss stage and velocity profile. Figure 9 displays the pseudo-dynamical Rb abundances versus v_exp(OH) for β = 0.2 and Ṁ values of 10^-8, 10^-7 and 10^-6M_⊙ yr^-1. The Rb results, when fixing both Ṁ and β, are very similar (with a slightly tighter correlation with v_exp(OH)) to those obtained when only fixing Ṁ (see Fig. 8) because β = 0.2 is the most common value obtained from the best spectral fits (all wind parameters free); an exception is the AGB star IRAS 15193+3132 (with the lowest v_exp(OH) and high β) for which only an upper limit to the Rb abundance ([Rb/Fe]≤0.7) could be obtained because the pseudo-dynamical model does not converge for such an unusual combination of wind parameters (Ṁ = 10^-6M_⊙ yr^-1; β = 0.2; v_exp(OH)= 3 km s^-1) coupled with [Rb/Fe] < 0.7 dex (see Fig. 9).

Fig. 8 Rb abundances vs. the expansion velocity (vexp(OH)) for extended model atmospheres with Ṁ = 10-8, 10-7 and 10-6M⊙ yr-1 (green triangles, yellow squares and cyan diamonds, respectively) in comparison with those obtained from hydrostatic models (blue dots).

Fig. 9 Rb abundances vs. the expansion velocity (vexp(OH)) for extended model atmospheres with β = 0.2 and Ṁ = 10-8, 10-7 and 10-6M⊙ yr-1 (green triangles, yellow squares and cyan diamonds, respectively) in comparison with those obtained from hydrostatic models (blue dots).

Figure 10 displays the [Rb/Fe] abundances from the best spectral fits versus the variability periods P. As mentioned above, our sample is composed of AGB stars of different progenitor masses and evolutionary stages. Most stars with P > 400 days are Li-rich and present some Rb enhancement⁶, which suggests that, on average, these stars are more massive stars experiencing HBB and/or more evolved stars (because of the longer periods) than the group of non Li-rich (and generally Rb-poor) stars with P< 400 days. The stars IRAS 05027−2158 (P = 368 days) and IRAS 20343−3020 (P = 349 days) are exceptions in the latter group. IRAS 05027−2158 is slightly Li-rich and Rb-poor, suggesting that it is a relatively massive AGB star (say ~3.5−4.5 M_⊙⁷) at the beginning of the TP phase (e.g. in an inter-pulse period immediately before or after the super Li-rich phase) but is not evolved enough for efficient Rb production (García-Hernández et al. 2013). On the other hand, IRAS 20343−3020 is slightly Rb-rich and Li-poor, which suggests a more advanced evolutionary stage and a slightly higher initial mass (say ~4.0−5 M_⊙) than IRAS 05027−2158 (we refer to Fig. 1 in García-Hernández et al. 2013).

Fig. 10 [Rb/Fe] pseudo-dynamical abundances versus variability period (P). The Li-rich and Li-poor stars are marked with magenta dots and green triangles, respectively. The two stars where Li could not be estimated are marked with cyan dots (see text).

5. Comparison with AGB nucleosynthesis models

In Fig. 11 we compare our new [Rb/Fe] abundances and [Rb/Zr] ratios with solar metallicity massive (3−9 M_⊙) AGB predictions from several nucleosynthesis models: van Raai et al. (2012), Karakas et al. (2012), Karakas & Lugaro (2016) (Monash), Pignatari et al. (2016) (NuGrid/MESA) and Cristallo et al. (2015) (FRUITY⁸). The predicted [Rb/Fe] abundances and [Rb/Zr] ratio ranges are 0.00−1.35 and −0.45−0.52 dex, respectively.

The Monash models (van Raai et al. 2012; Karakas et al. 2012; Karakas & Lugaro 2016) use the stellar evolutionary sequences calculated with the Monash version of the Mount Stromlo Stellar Structure Program (Frost & Lattanzio 1996), which uses the Vassiliadis & Wood (1993) mass-loss prescription on the AGB. A post-processing code is used to obtain in detail the nucleosynthesis of a large number of species, including the s-process abundances. Due to convergence difficulties, the stellar evolution models used in the calculations are not always evolved until the end of the superwind phase and synthetic models have been used to estimate the effect of remaining TPs and to completely remove the envelope. We refer the reader to van Raai et al. (2012), Karakas et al. (2012) and Karakas & Lugaro (2016) for more details about the theoretical models. Here, we only report on the main differences between these models: i) The use of different nuclear networks; that is, the total number of nuclear species considered and the values of some reaction rates and neutron-capture cross-sections (see below); and ii) the use by Karakas et al. (2012) of a modified Vassiliadis & Wood (1993) mass-loss prescription, which delays the beginning of the superwind phase until the pulsation period reaches values of 700-800 days (instead of the value of 500 days used in the other models), resulting in a higher Rb production.

The NuGrid/MESA and FRUITY models assume AGB mass-loss prescriptions, nuclear physics inputs and treatments of convection different from the Monash models. In particular, the Blöcker (1995) and Straniero, Gallino & Cristallo (2006) mass-loss formulae for the AGB phase are assumed by the NuGrid/MESA and FRUITY models, respectively. Furthermore, these models produce self-consistently the ¹³C neutron source as a result of the different convective boundary mixing scheme and treatments of the convective borders, while in the Monash models the mixing required to produce the ¹³C neutron source is included in a parametrised way during the post processing and is typically not included in massive AGB stars, following theoretical (Goriely & Siess 2004) and observational indications (García-Hernández et al. 2013) (we refer also to Pignatari et al. 2016; Cristallo et al. 2015; Karakas & Lugaro 2016, for more details). In relation to the main results, i) the NuGrid/MESA solar metallicity massive AGB models are qualitatively similar to the Monash models in that HBB and light s-process element production (of the elements from Rb to Zr) are seen at the stellar surface, the latter due to the activation of the ²²Ne neutron source and the subsequent operation of the TDU; and ii) the FRUITY solar metallicity massive AGB models are different to the Monash and NuGrid/MESA models because these models experience very inefficient TDU, hence the signature of the nucleosynthesis due to the ²²Ne neutron source is not visible at the stellar surface.

Fig. 11 Model predictions from van Raai et al. (2012), Karakas et al. (2012), Karakas & Lugaro (2016), FRUITY database and Pignatari et al. (2016): Stellar mass vs. [Rb/Fe] (left panel) and [Rb/Zr] (right panel). The abundances from the last computed thermal pulse and from the synthetic evolution calculations are shown with dots and triangles, respectively. P represents the period for the beginning of the superwind phase and n is the number of species in the nucleosynthesis network. The shaded regions mark the range of the new Rb abundances and [Rb/Zr] ratios obtained in our sample with extended models.

Figure 11 shows that the FRUITY massive AGB model predicts final [Rb/Fe] < 0.15, which does not explain the observed range of Rb abundances and [Rb/Zr] ratios; specifically the [Rb/Zr] ratios remain negative for all masses. Another difference between the FRUITY models and the Monash and NuGrid models is that the FRUITY models do not predict HBB to occur in AGB stars unless the metallicity is very low; at least ten times lower than solar. However, spectroscopic observations of massive AGB stars demonstrate that they experience HBB; as evidenced by i) strong Li overabundances observed in massive AGB stars in the Galaxy (Fe/H = 0.0; e.g. García-Hernández et al. 2007, 2013), the Magellanic Clouds (Fe/H = −0.7–0.3; e.g. Plez et al. 1993; Smith et al. 1995; García-Hernández et al. 2009) and the dwarf galaxy IC 1613 (Fe/H = −1.6; e.g. Menzies, Whitelock & Feast 2015); and ii) N enhancements and low ¹²C/¹³C ratios in Magellanic Cloud Li-rich massive AGBs (e.g. Plez et al. 1993; McSaveney et al. 2007). The lack of HBB in the FRUITY predictions is also at odds with the observations of the so-called type-I planetary nebulae in very different metallicity environments and galaxies; which are expected to be the descendants of HBB massive AGB stars based on their strong N and He overabundances (see e.g. Stanghellini et al. 2006; Karakas et al. 2009; Leisy& Dennefeld 1996; García-Rojas et al. 2016, and references therein).

It is to be noted here that the several Monash AGB models (van Raai et al. 2012; Karakas et al. 2012; Karakas & Lugaro 2016) mentioned above notably use different rates for the ²²Ne(α,n)²⁵Mg reaction, which drives the production of s-process elements in massive AGB stars. In particular, Karakas & Lugaro (2016) use the ²²Ne(α,n)²⁵Mg reaction from Iliadis et al. (2010), neutron-capture cross-section of the Zr isotopes (Lugaro et al. 2014), and a more extended nuclear network of 328 species (from H to S, and then from Fe to Bi). The van Raai et al. (2012) models, instead, use a nuclear network of 166 species (up to Nb) and the ²²Ne(α,n)²⁵Mg reaction rate from Karakas et al. (2006), while Karakas et al. (2012) explored different networks (166, 172 and 320 species) and ²²Ne(α,n)²⁵Mg reaction rates; from Karakas et al. (2006), Iliadis et al. (2010) and Angulo et al. (1999), NACRE.

The van Raai et al. (2012) models (from 4 to 6.5 M_⊙ at Z = 0.02; Fig. 11) show that both the [Rb/Fe] abundances and [Rb/Zr] ratios increase with the initial mass of the AGB star, as the star becomes hotter and the ²²Ne(α,n)²⁵Mg reaction is more efficiently activated. However, the [Rb/Fe] abundances from the last computed TP are too low (ranging from 0.0 to 0.26 dex). The corresponding Rb abundances ([Rb/Fe] ~ 0.0–1.0 dex) from the synthetic evolution calculations cover most of the Rb abundances observed; although they cannot explain the star IRAS 05151+6312 with [Rb/Fe] = 1.3 dex. Such high Rb abundances can be reached by the synthetic calculations of the solar metallicity 6 and 7 M_⊙ AGB models with delayed superwinds of Karakas et al. (2012) when using the faster NACRE rate for the ²²Ne(α,n)²⁵Mg reaction. Finally, the Karakas & Lugaro (2016) models (from 4.5 to 8 M_⊙ at Z = 0.014 ⁹; Fig. 11) predict lower Rb abundances than the Karakas et al. (2012) models of the same mass and similar metallicity, mostly due to the implementation of the delayed superwind and the use of the NACRE rate in Karakas et al. (2012).

The NuGrid/MESA models (from 3 to 5 M_⊙ at Z = 0.02; Fig. 11) reproduce the observed [Rb/Fe] and [Rb/Zr] ranges, up to 0.9 and 0.4 dex, respectively. However, we note that only in the 5 M_⊙ case do the NuGrid/MESA models see signature of HBB and predict O-rich stars. The 3 and 4 M_⊙ cases become C-rich stars and do not experience HBB, which is at odds with our sample of O-rich stars (García-Hernández et al. 2006).

Regarding the [Rb/Zr] ratios, obviously also in this case the higher [Rb/Zr] ratios are obtained from the models with delayed superwind (P = 700–800), however, these [Rb/Zr] ratios are still lower than our observed values. The maximum value from the AGB models is [Rb/Zr] = 0.52 for M = 5 M_⊙, and the maximum value from our observations is [Rb/Zr] = 1.05. A possible explanation is that Zr could be depleted into dust (e.g. van Raai et al. 2012; Zamora et al. 2014), producing the differences between the theoretical and observational [Rb/Zr] ratios. Abundance measurements of similar s-elements such as Sr and Y would be needed in order to clarify this problem.

6. Conclusions

Here we report new Rb and Zr abundances determined from the 7800 Å Rb I line and the 6474 Å ZrO bandhead, respectively, in a complete sample of massive Galactic AGB stars, previously studied with hydrostatic models, using more realistic extended atmosphere models and a modified version of the spectral synthesis code Turbospectrum, which considers the presence of a circumstellar envelope with a radial wind. The Rb abundances are much lower (in some cases even 1–2 dex) with the pseudo-dynamical models, while the Zr abundances are close to the hydrostatic ones because the 6474 Å ZrO bandhead is formed deeper in the atmosphere and is less affected than the 7800 Å Rb I resonant line by the circumstellar effects.

We have studied the sensitivity of the determined abundances to variations in the stellar (T_eff) and wind (Ṁ, β and v_exp) parameters. The Rb abundances are very sensitive to the mass loss rate Ṁ but much less to the β parameter and v_exp(OH). The Zr abundances, instead, are not affected by variations of the stellar and wind parameters. The Rb abundances from extended models are lower than those obtained from the hydrostatic ones, and the difference is larger in the stars with the highest Rb abundances in the hydrostatic case. We have represented the hydrostatic and pseudo-dynamical Rb abundances against the v_exp(OH), which can be used as a mass indicator independent of the distance, and we have observed a flatter correlation. The difference between the hydrostatic and pseudo-dynamical Rb abundances increases with increasing the v_exp(OH), due to the fact that the presence of a circumstellar envelope affects more strongly the more massive stars. Furthermore, the dispersion of the correlation between the Rb abundance and v_exp(OH) is larger in the pseudo-dynamical case. When we fix the wind parameters Ṁ (i.e. equivalent to assuming that our AGB sample stars have a similar evolutionary stage in terms of mass loss), and/or β (the same velocity profile), the dispersion is lower.

The Monash nucleosynthesis theoretical predictions reproduce the range of new Rb and Zr abundances although [Rb/Fe] values above 1.0 can only be matched if the superwind is delayed to after the period reaches 700–800 days. We also note that the rate of the ²²Ne(α,n)²⁵Mg reaction is crucial, but still hampered by large systematic uncertainties (e.g. Bisterzo et al. 2016; Massimi et al. 2017). Underground measurements, planned, for example, at LNGS-LUNA (Laboratory for Underground Nuclear Astrophysics) will help to resolve the current issues. The FRUITY massive AGB models predict Rb abundances much lower than observed and negative [Rb/Zr] ratios, at odds with the observations. The NuGrid/MESA models of 4 and 5 M_⊙ predict [Rb/Fe] as high as 0.9 dex, however, the 4 M_⊙ model does not experience HBB and becomes C-rich, while our sample stars are clearly O-rich. The maximum observed [Rb/Zr] ratios are still more than a factor of two larger than predicted by the nucleosynthesis models. A possible explanation for this difference between the observations and the predictions is that Zr could be depleted into dust. Observations of other s-process elements Sr and Y belonging to the same first peak as Rb and Zr will help to clarify this mismatch.

In summary, the [Rb/Fe] abundances and [Rb/Zr] ratios previously derived with hydrostatic models are certainly not predicted by the most recent theoretical models of AGB nucleosynthesis. In particular, the highest [Rb/Fe] abundances and [Rb/Zr] ratios observed in massive Galactic AGBs are much larger than theoretically predicted. The new [Rb/Fe] abundances and [Rb/Zr] ratios as obtained from our simple (but more realistic) pseudo-dynamical model atmospheres are much lower in much better agreement with the theoretical predictions, significantly resolving the mismatch between the observations and the nucleosynthesis models in the more massive AGB stars. This confirms the earlier Zamora et al. preliminary results on a smaller sample of massive O-rich AGB stars but here we find that the Rb abundances are strongly dependent on the wind mass-loss Ṁ, which is basically unknown in our AGB star sample. Follow-up radio observations (e.g. the rotational lines of the several CO isotopologues) of these massive Galactic AGB stars are encouraged in order to obtain precise mass-loss rate estimates needed to break the actual models degeneracy and obtain more reliable (non-model-dependent) Rb abundances in massive AGB stars.

Acknowledgments

This work is based on observations at the 4.2 m William Herschel Telescope operated on the island of La Palma by the Isaac Newton Group in the Spanish Observatorio del Roque de Los Muchachos of the Instituto de Astrofisica de Canarias. Also based on observations with the ESO 3.6 m telescope at La Silla Observatory (Chile). We thank Marco Pignatari and Umberto Battino for providing information about the Nugrid/MESA models. V.P.M. acknowledges the financial support from the Spanish Ministry of Economy and Competitiveness (MINECO) under the 2011 Severo Ochoa Program MINECO SEV-2011-0187. D.A.G.H. was funded by the Ramón y Cajal fellowship number RYC-2013-14182. V.P.M., O.Z., D.A.G.H. and A.M. acknowledge support provided by the MINECO under grant AYA-2014-58082-P. M.L. is a Momentum (“Lendület-2014” Programme) project leader of the Hungarian Academy of Sciences. M.L. acknowledges the Instituto de Astrofísica de Canarias for inviting her as a Severo Ochoa visitor during 2015 August when part of this work was done. This paper made use of the IAC Supercomputing facility HTCondor (http://research.cs.wisc.edu/htcondor/), partly financed by the Ministry of Economy and Competitiveness with FEDER funds, code IACA13-3E-2493. This work benefited from discussions at The 12th Torino Workshop on Asymptotic Giant Branch Stars in August 2016 supported by the National Science Foundation under Grant No. PHY-1430152 (JINA Center for the Evolution of the Elements).

References

Appendix A: Complete sample

Fig. A.1 Observed spectra (black dots), best hydrostatic (blue lines) and pseudo-dynamical (red lines) fits of our sample of AGB stars in the regions of 7800 Å Rb I line (left panels) and 6474 Å ZrO bandhead (right panels). The parameters of the best fit model atmospheres are indicated in Table 3. The plots are displayed in increasing R.A. order.

All Tables

All Figures

	Fig. 1 Velocity vs. distance from the star in four of our AGB wind models. These velocity laws present different expansion velocities vexp(OH), mass-loss rates and β exponents. The effective temperature Teff = 3000 K, gravity log g = −0.5 and the solar chemical composition are the same in all models.
In the text

	Fig. 2 Illustrative examples of synthetic spectra for different stellar (Teff) and wind (Ṁ, β, and vexp(OH)) parameters in the spectral region around the 7800 Å Rb I line. The black vertical line indicates the position of the 7800 Å Rb I line.
In the text

	Fig. 3 Illustrative examples of synthetic spectra for different stellar (Teff) and wind (Ṁ, β, and vexp(OH)) parameters in the spectral region around the 6474 Å ZrO bandhead. The black vertical line indicates the position of the 6474 Å ZrO bandhead.
In the text

	Fig. 4 Rb I 7800 Å (left panel) and ZrO 6474 Å (right panel) spectral regions in massive Galactic AGB stars. The pseudo-dynamical models (red lines) that best fit the observations (black dots) are shown in four sample stars. For comparison, the hydrostatic models are also displayed (blue lines). The location of the Rb I line and the ZrO bandhead are indicated by black vertical lines.
In the text

	Fig. 5 IRAS 07222−2005 (left panel) and IRAS 09194−4518 (right panel) display two components (circumstellar and photospheric) in the 7800 Å Rb I line (black line) and cannot be reproduced by our present version of the spectral synthesis code. The red line shows the pseudo-dynamical synthetic spectrum that gives a good fit to the photospheric component only.
In the text

	Fig. 6 Rb abundances derived both with hydrostatic (blue dots) and pseudo-dynamical model atmospheres (red triangles) with best-fit parameters plotted against the OH expansion velocity. The Li-rich stars are indicated by squares. A typical error bar of ±0.7 dex is also displayed.
In the text

	Fig. 7 Rb abundances derived both with hydrostatic (blue dots) and extended model atmospheres with β = 0.2 (magenta triangles) and β = 1.2 (green squares) plotted against the OH expansion velocity.
In the text

	Fig. 8 Rb abundances vs. the expansion velocity (vexp(OH)) for extended model atmospheres with Ṁ = 10-8, 10-7 and 10-6M⊙ yr-1 (green triangles, yellow squares and cyan diamonds, respectively) in comparison with those obtained from hydrostatic models (blue dots).
In the text

	Fig. 9 Rb abundances vs. the expansion velocity (vexp(OH)) for extended model atmospheres with β = 0.2 and Ṁ = 10-8, 10-7 and 10-6M⊙ yr-1 (green triangles, yellow squares and cyan diamonds, respectively) in comparison with those obtained from hydrostatic models (blue dots).
In the text

	Fig. 10 [Rb/Fe] pseudo-dynamical abundances versus variability period (P). The Li-rich and Li-poor stars are marked with magenta dots and green triangles, respectively. The two stars where Li could not be estimated are marked with cyan dots (see text).
In the text

Fig. 11 Model predictions from van Raai et al. (2012), Karakas et al. (2012), Karakas & Lugaro (2016), FRUITY database and Pignatari et al. (2016): Stellar mass vs. [Rb/Fe] (left panel) and [Rb/Zr] (right panel). The abundances from the last computed thermal pulse and from the synthetic evolution calculations are shown with dots and triangles, respectively. P represents the period for the beginning of the superwind phase and n is the number of species in the nucleosynthesis network. The shaded regions mark the range of the new Rb abundances and [Rb/Zr] ratios obtained in our sample with extended models.

In the text

	Fig. A.1 Observed spectra (black dots), best hydrostatic (blue lines) and pseudo-dynamical (red lines) fits of our sample of AGB stars in the regions of 7800 Å Rb I line (left panels) and 6474 Å ZrO bandhead (right panels). The parameters of the best fit model atmospheres are indicated in Table 3. The plots are displayed in increasing R.A. order.
In the text