ELEMENTAL ABUNDANCE DIFFERENCES IN THE 16 CYGNI BINARY SYSTEM: A SIGNATURE OF GAS GIANT PLANET FORMATION?

The iron lines alone can be used to determine the atmospheric parameters using the conditions of excitation and ionization equilibrium. We start by assuming that the atmospheric parameters T_eff (effective temperature), log g (surface gravity), [Fe/H] (iron abundance),⁹ and v_t (microturbulent velocity) of 16 Cyg A and B are solar (i.e., T_eff, log g, [Fe/H], v_t = 5777 K, 4.437, 0.00, 1.00 km s⁻¹).¹⁰ Then, T_eff and v_t were iteratively modified so that trends of iron abundance (from Fe i lines only) with excitation potential (EP) and reduced equivalent width (log EW/λ), respectively, approached zero. Simultaneously, log g was modified so that the mean difference in iron abundance from Fe i and Fe ii lines also approached zero. Errors in the derived parameters were computed from the uncertainty in the measured slope of abundance versus EP (for T_eff) and from the standard deviation of the mean difference of Fe i and Fe ii abundances (for log g).

Using the excitation/ionization equilibrium method we obtain: T_eff = 5818 ± 38 K, log g = 4.31 ± 0.05, [Fe/H] = +0.104 ± 0.018 for 16 Cyg A; and T_eff = 5742 ± 30 K, log g = 4.31 ± 0.05, [Fe/H] = +0.057 ± 0.017 for 16 Cyg B (see Figure 3). If we use Kurucz no-overshoot models instead of MARCS we find a shift of only +5 K in T_eff, +0.01 dex in log g, and +0.002 dex in [Fe/H], systematically consistent for both stars. Moreover, adopting a slightly different solar microturbulence, v_t = 1.0 ± 0.1 km s⁻¹, we find changes of only ±0.003 dex in [Fe/H] and, as expected, shifts of ±0.1 km s⁻¹ in the microturbulent velocities of 16 Cyg A and B. This shows that the choices of model atmosphere flavor and solar microturbulence have an almost negligible impact on our differential abundance results.

The calculation described in this section suggests that 16 Cyg A is slightly more metal-rich than 16 Cyg B. However, since there is a non-negligible difference in the effective temperature of the stars, small differential systematic effects might still be present and it is therefore important to check the parameters derived here using other independent methods.

Casagrande et al. (2010) have implemented the so-called infrared flux method (IRFM) to derive the effective temperatures of a large sample of dwarf and subgiant stars. The zero point of their temperature scale has been determined with an unprecedented precision of only 15 K. Metallicity-dependent color–T_eff transformations provided in that work allow us to estimate T_eff values for the 16 Cyg A and B pair from their observed photometry. We used the (B − V) and (b − y) colors and our derived [Fe/H] values to obtain: T_eff = 5792 K for 16 Cyg A and T_eff = 5733 K for 16 Cyg B. Photometric errors translate into an error of about 16 K for both stars, similar to the color-to-color error. The uncertainties in [Fe/H] translate into a 5 K uncertainty in T_eff. Adding these errors in quadrature we estimate a random error of 23 K, which has to be added linearly to the systematic error in the zero point of the Casagrande et al. (2010) T_eff scale to estimate the uncertainties in our derived T_eff values. Thus, the effective temperatures we derived using colors have errors of 38 K.

The color temperatures, which are based on the most reliable IRFM temperature scale to date, confirm that 16 Cyg A is warmer than 16 Cyg B, although the difference seems to be smaller when using the IRFM T_eff scale (59 K) instead of the spectroscopic method described in the previous section (76 K). Nevertheless, the T_eff values from the two methods are in good agreement within the 1σ errors.

The wings of the Balmer lines are known to be sensitive to T_eff and if log g and [Fe/H] can be estimated with high precision independently, as in our case, model fits to the Hα line profile can provide very precise T_eff values (e.g., Fuhrmann et al. 1994; Barklem et al. 2002; Ramírez et al. 2006). This is particularly the case if a good continuum normalization is made using very high quality spectra.

We used least-squares minimization to find the best model fits to our Hα data as well as the formal errors in T_eff. Careful inspection of our spectra allowed us to find several wavelength windows free from spectral lines, both stellar and telluric, other than Hα, which were then used in the model fits. The same windows were used for the three spectra. In Figure 4, we show the observed normalized Hα lines for 16 Cyg A, B, and the Sun (day sky) along with their best-fit models. The problems described in Section 2 for the solar spectrum are not as important for this broad feature, in particular because we are only interested in the wings, so we can use reliably our day-sky spectrum as solar reference in this case. For our solar data, we derive T_eff = 5732 ± 32 K, similar to the value obtained by Barklem et al. (2002).¹¹ This suggests that there is still room for improvement in the modeling of Balmer lines. Hereafter, we apply a zero-point correction of +45 K to the effective temperatures derived from Hα since this brings our derived solar T_eff up to the nominal value of 5777 K.

Zoom In Zoom Out Reset image size

We computed T_eff values from the Hα line profile fits for a small grid of log g = 4.2, 4.3, 4.4 and [Fe/H] = 0.00, 0.05, 0.10. Then, we used the grid of T_eff results to interpolate to our preferred log g and [Fe/H] parameters as well as to compute the error introduced by uncertainties in log g and [Fe/H]. With our final log g and [Fe/H] values we obtain: T_eff = 5819 ± 25 K for 16 Cyg A and T_eff = 5762 ± 26 K for 16 Cyg B. An error of 0.02 dex in log g translates into only 2 K whereas an error of 0.03 dex in [Fe/H] translates into about 1 K of T_eff. The uncertainties in our derived values of log g and [Fe/H] thus have negligible impact on this calculation. Moreover, given the very high S/Ns of our spectra in this region (S/N ≃ 500, obtained because the spectrograph setup was chosen to maximize S/N at 6562 Å), our T_eff values derived from Hα line profile fitting are the most reliable, as the smaller errors in comparison with T_eff values derived with other methods suggest.

Our three estimates of T_eff are in excellent agreement with each other for both 16 Cyg A and B. Since these techniques are independent, we can obtain a final T_eff value by averaging them, using the estimated errors as weights. We find: T_eff = 5813 ± 18 K for 16 Cyg A and T_eff = 5749 ± 17 K for 16 Cyg B.¹² The degree of precision with which we have been able to measure the effective temperatures of these stars is almost unprecedented for any other star but the Sun. The only other exception is the case of the solar twin star HIP 56948, for which Meléndez et al. (2011) report a T_eff with an error of 7 K.

The parallax of the 16 Cyg pair has been accurately measured by the Hipparcos mission (Perryman et al. 1997). The new reduction by van Leeuwen (2007) gives values of 47.44 ± 0.27 mas for 16 Cyg A and 47.14 ± 0.27 mas for 16 Cyg B. The V magnitudes of these stars are also precisely known: V(A) = 5.959 ± 0.009, V(B) = 6.228 ± 0.019. Using these numbers we can infer the absolute V magnitudes of the stars: M_V(A) = 4.340 ± 0.015, M_V(B) = 4.595 ± 0.023. Since their effective temperatures are also known with very high precision, we can reliably place the stars on the H-R diagram and compare their location with theoretical stellar evolution predictions (isochrones) to infer their masses, ages, and surface gravities.

Our particular implementation of stellar parameter determination using isochrones is described in Meléndez et al. (2011, see also Baumann et al. 2010). Our method follows the probabilistic scheme explained in R03 (their Section 3.4) and also Allende Prieto et al. (2004, their Appendix A) but the choice of isochrones is different. Basically, all isochrone points within a radius of 3σ from the observed parameters are used to compute mass, age, and log g probability distribution functions, from which a most likely value and 68% (1σ-like) and 95% (2σ-like) confidence limits can be derived. Finely spaced Yonsei–Yale isochrones (e.g., Yi et al. 2001, 2003; Kim et al. 2002) were used in these computations.

Using our isochrone method we infer masses of M(A) = 1.05 ± 0.02 M_☉ and M(B) = 1.00 ± 0.01 M_☉. For the surface gravities we find: log g(A) = 4.28 ± 0.02 and log g(B) = 4.33 ± 0.02. Since the two stars were born at the same time, it is expected for the primary to be slightly more evolved given its somewhat larger mass.

The ages that we infer for the stars in the 16 Cyg system are τ(A) = 7.15^+0.04_{− 1.03} Gyr and τ(B) = 7.26^+0.69_{− 0.33} Gyr. Error bars correspond to the 1σ-like lower and upper limits. The age probability distributions (APDs) for these stars are shown in Figure 5. The asymmetry of these APDs explains the very different upper and lower 1σ errors given above. Although they are somewhat broad, the APDs of 16 Cyg A and B peak at nearly the same age. A much more precise age can be obtained by combining the two APDs, as shown with the solid line histogram in Figure 5. This combined 16 Cyg APD is simply the product of the APDs of the components of the 16 Cyg binary system. Using the combined APD we infer an age of 7.10^+0.18_{− 0.35} Gyr for 16 Cyg. Note that this very small age error bar does not include systematic errors due to model uncertainties or errors introduced by statistical biases (see below). It represents only the internal precision of our method when applied to a binary system.

Zoom In Zoom Out Reset image size

Stellar parameters derived from isochrones are subject to a number of systematic errors as well as statistical biases due to isochrone sampling which are described in detail by, for example, Nordström et al. (2004, their Section 4.5.4). Bayesian approaches have been proposed to tackle these problems (e.g., Pont & Eyer 2004; Jørgensen & Lindegren 2005; da Silva et al. 2006; Casagrande et al. 2011). We can compare our results with the Bayesian implementation by da Silva et al. (2006) since they provide a very convenient online tool that uses their method.¹³ We find log g(A) = 4.26 ± 0.01 and log g(B) = 4.31 ± 0.01, in excellent agreement with our estimates even though the choice of isochrones is different (da Silva et al. use PADOVA isochrones; e.g., Bertelli et al. 1994; Girardi et al. 2000) and we do not correct for statistical biases. The masses derived with the da Silva et al. (2006) online tool are about 0.02 M_☉ lower but the difference between 16 Cyg A and B is nearly identical. For the stellar ages, the da Silva et al. online tool suggests 6.7 ± 0.5 Gyr (A) and 7.6 ± 0.7 Gyr (B). These values are also consistent with our derived age for the 16 Cyg system considering the 1σ uncertainties.

Note that the excitation/ionization method of atmospheric parameter determination was not able to resolve the small difference in log g found in this section although the results are also consistent within 1σ. Clearly, for nearby stars with well-determined effective temperatures, parallaxes, and good photometry, the isochrone method is superior.

The wings of certain strong spectral features are mildly sensitive to the stellar surface gravity. A particularly important feature in this regard is the Mg i b triplet at about 5180 Å. The effective temperature and magnesium abundance must be known with high precision in order for this line to be useful in the determination of log g.

Figure 6 shows the spectral region around the strongest of the Mg i b lines and the red wing of the middle one along with synthetic spectra computed for values of log g that reproduce best the regions free from lines other than the Mg i b lines, such as those around 5174 Å and 5185 Å. We fine-tuned the transition probabilities of these lines to reproduce the solar Mg i b lines with an Mg abundance of A_Mg = 7.60, as measured by Asplund et al. (2009). For 16 Cyg A and B, we used synthetic models with our final T_eff and [Fe/H] values as well as our derived magnesium abundances. Synthetic profiles for log g ± 0.2 dex are also shown in Figure 6 to illustrate the sensitivity of these features to log g. From a least-squares minimization of the clean regions we find log g = 4.27 ± 0.04 for 16 Cyg A and log g = 4.33 ± 0.04 for 16 Cyg B, in excellent agreement with the values found in the previous section.

Zoom In Zoom Out Reset image size

Similar to the case of T_eff, the three methods used to derive log g are independent so a weighted average can be used to determine the best solution. We find: log g(A) = 4.282 ± 0.017, log g(B) = 4.328 ± 0.017. These and other fundamental parameters derived in this section are summarized in Tables 2 and 3.

High-precision elemental abundance work in solar twins and analogs has revealed a possible connection between the properties of a star and the process of planet formation that occurred around it. The work by M09 showed that, compared to the majority (≃85%) of solar twin stars, the Sun exhibits a deficiency of refractory elements relative to volatiles (of order 20%). They proposed that this solar abundance anomaly could be due to the formation of the terrestrial planets in the solar system since these objects are rich in refractory elements. Thus, the missing mass of refractory elements in the solar photosphere would be currently located in the terrestrial planets. This hypothesis assumes that the majority of solar twin stars did not form as many rocky objects as in the solar system or that if they did, they were accreted by the star at later times, possibly through orbital migration. Only about 15% of solar twin stars show a photospheric composition indistinguishable from the solar one. These objects are potentially terrestrial planet hosts. Later work by R09 and Ramírez et al. (2010) confirmed the observational results of M09 and expanded the sample toward higher metallicities, where they found that, based on the abundances measured, the fraction of terrestrial planet hosts increases toward higher [Fe/H].

Thus, it seems likely that a "planet signature" has been imprinted in the chemical composition of the host stars and, arguably, it appears to relate only to the formation of rocky objects. Indeed, M09 showed that metal-rich solar analogs with giant planets detected have abundance patterns relative to iron similar to those of the majority of solar twins, suggesting that the solar abundance anomaly could not be related to the formation of the gas giants. This type of analysis, however, should be examined more carefully because the impact of systematic errors and GCE could be important in this higher [Fe/H] regime (see also Schuler et al. 2011). Data by independent groups have confirmed the observational result by M09 and R09 (Ramírez et al. 2010; González Hernández et al. 2010; Gonzalez et al. 2010b). Although the planet formation signature interpretation has been contended by González Hernández et al. (2010), a response to their criticism has been provided by Ramírez et al. (2010). Alternative scenarios to explain the solar abundance anomaly have also been explored, without success (Ramírez et al. 2010; Meléndez et al. 2011; Kiselman et al. 2011).

The planet signature hypothesis described above has been quantitatively tested by Chambers (2010), who computed the effect of adding two Earth masses (2M_⊕) of Earth-like material and 2M_⊕ of CM chondrite material to the present-day solar convective envelope. He finds that in this experiment the solar abundance anomaly would disappear, i.e., the Sun would not show the deficiency of refractory elements anymore. Thus, his calculations support the hypothesis that the missing mass of rock in the solar photosphere is currently locked up in the rocky bodies of the solar system. A similar estimate has been recently made by Meléndez et al. (2011) who suggest an even larger amount of CM chondrite material (4M_⊕). Although such large mass is not observed in the present-day asteroid belt, it has been suggested that small rocky bodies were very abundant when the Sun was formed, but most of the material from the asteroid belt was probably removed by strong resonances (e.g., Weidenschilling 1977; Chambers & Wetherill 2001; Petit et al. 2001).

Excluding the heavy elements (i.e., for Z ⩽ 30), the abundance patterns of both 16 Cyg A and 16 Cyg B are compatible with those observed in most solar twin stars (Figure 7). Within the planet signature scenario, this implies that none of these objects host terrestrial planets or that any rocky body that once formed there has been accreted by the host star. Perhaps binary systems like 16 Cyg are not stable enough to retain any terrestrial planets that might form. The gas giant planet around 16 Cyg B is highly eccentric, a property that has been attributed by some authors to the influence of the stellar companion (e.g., Cochran et al. 1997; Mazeh et al. 1997). The lack of rocky planets around 16 Cyg B could also be related to the presence of the gas giant planet. As evidenced by the Kepler mission data, multiple systems with small planets rarely contain hot Jupiters (Latham et al. 2011). In any case, the abundances of 16 Cyg A and B are not unusual considering that only about 15% of solar twins share the exact same (solar) composition.

The heavy elements clearly do not follow the typical solar twin trend (Figure 8). With the exception of Eu, the other Z > 30 elements are found significantly below the mean solar twin trend. The Nd abundances are marginally consistent with the expected trend but the Y, Zr, and Ba abundances are irreconcilable with that pattern. These elements are produced primarily by neutron capture reactions on timescales slow (the s-process) or rapid (the r-process) relative to the average β-decay rates of radioactive isotopes along the reaction chain. The s-process is associated with low- and intermediate-mass stars (1.3  M_☉ ≲ M ≲ 8.0 M_☉) on the asymptotic giant branch. The site (or sites) of the r-process is not known but is suspected to be associated with core collapse supernovae (M ≳ 8.0 M_☉). These stellar timescales imply that r-process enrichment occurs earlier than s-process enrichment, which is confirmed by the observed heavy element abundance patterns in low-metallicity stars in the Galactic halo and disk (e.g., McWilliam 1997).

Each of the r- and s-process contributes about half of the elements with Z > 30 in the Sun, and, following pioneering work by Cameron (1973), the solar abundance of each isotope (or element) can be decomposed into r- and s-process fractions. For example, 72% of the solar Y was produced via the s-process (e.g., Simmerer et al. 2004). Similarly, 81% of Zr, 85% of Ba, 58% of Nd, and 3% of Eu in the Sun were produced by the s-process. These fractions represent the chemical inventory of one point in the Galactic interstellar medium 4.5 Gyr ago. The stars in the 16 Cyg system formed 2.5 Gyr earlier and they are slightly deficient in s-process material relative to the Sun. As can be seen in Figure 8, elements with a significant s-process component, such as Y, Zr, and Ba, are the most deficient. Nd, which has roughly equal contributions from the r- and s-processes, is moderately deficient, and Eu, which has a minimal s-process contribution, is not deficient at all. The fact that the heavy elements in both 16 Cyg A and B show the same abundance pattern as a function of T_C suggests that this is a primordial effect—the result of GCE. The slightly super-solar metallicity of the system is not in contradiction, as the slope of [Ba/H] versus time is steeper than that for [Fe/H] (e.g., Edvardsson et al. 1993).

The works by M09 and R09 concentrated only on solar twin stars, therefore minimizing these relatively small GCE effects. Nevertheless, it is interesting to note that the solar twin abundances of Y, Zr, and Ba by R09 show larger star-to-star error than those of lighter elements (see their Figures 1 and 3), a result that could be associated with the age span of their sample. The opposite effect of that observed in the 16 Cyg system (i.e., neutron capture element abundances higher than the mean solar twin trend) has been recently observed in the younger solar twin star 18 Sco (J. Meléndez et al. 2011, in preparation), reinforcing the idea that this is related to GCE. Thus, our results for the heavy elements tell us that extreme caution must be practiced when interpretations of abundance ratios versus condensation temperature are made for stars that are not so similar to the Sun, in particular stars that are significantly older or younger, because GCE effects are not entirely negligible.

Based on the discussion given in the previous sections, we are led to conclude that the metallicity difference between 16 Cyg A and B is not primordial but acquired. Either the photosphere of 16 Cyg A became more metal-rich or that of 16 Cyg B became more metal-poor after they were formed (or while they were being formed) from the same molecular cloud with the same initial chemical composition. Laws & Gonzalez (2001) suggested that 16 Cyg A accreted either 2.5 Earth masses of chondritic material or 0.3 M_J of gas giant material to explain its higher [Fe/H] (they derived Δ[Fe/H] = +0.025 ± 0.009). They also argued that the late accretion of planetary material could increase the lithium abundance on the star's convective envelope. Indeed, inspection of the 6708 Å region clearly shows that 16 Cyg B is lithium-poor compared to 16 Cyg A (Figure 10). Our derived Li abundances are A(Li) = 1.34 ± 0.04 for 16 Cyg A and A(Li) = 0.73 ± 0.06 for 16 Cyg B, in good agreement with previously published estimates (King et al. 1997; Luck & Heiter 2006; Takeda et al. 2007; Israelian et al. 2009; Baumann et al. 2010). Thus, 16 Cyg A has about 4.1 times (≃ 0.6 dex) more lithium than 16 Cyg B in its photosphere.

Zoom In Zoom Out Reset image size

The behavior of Li abundances in solar-type stars is very complex and not yet fully understood (e.g., Chaboyer 1998; Pinsonneault 2010). In general, however, inspection of large samples clearly shows that at around 1 M_☉ and T_eff = 5780 K a very steep slope of increasing Li abundances with stellar mass (or T_eff) is present (e.g., Takeda et al. 2007; Baumann et al. 2010). Fitting a straight line to the Li abundance versus T_eff relation of stars in the 5700–5900 K range we find that the lithium abundances of stars with T_eff values similar to those of 16 Cyg A and B should differ by a factor of 3.5 (≃ 0.55 dex), i.e., very close to the observed value. J. R. Fish et al. (in preparation) have recently derived lithium abundances of more than 700 solar-type stars and compiled and homogenized lithium abundance data for more than 1400 FGK stars. They also derive masses and ages in a consistent manner employing the exact same method used in this work. The lithium versus age trends for stars of mass and metallicity similar to those of 16 Cyg A and B, taken from the Fish et al. catalog, are shown in Figure 11. The lithium abundances of both 16 Cyg A and B are perfectly normal for stars of their mass, age, and metallicity, despite the somewhat large scatter, in particular for stars like 16 Cyg B.

Zoom In Zoom Out Reset image size

If the convective envelope of 16 Cyg A, after it reached its present size, accreted 0.3 M_J of planetary material that has not experienced lithium depletion, as suggested by Laws & Gonzalez (2001), the photospheric lithium abundance would have increased by about 0.4 dex (a factor of 2.5). Thus, if 16 Cyg A was polluted by planetary material, considering the mass effect on lithium depletion, the lithium abundance of 16 Cyg A should be about one order of magnitude larger than that of 16 Cyg B, which is significantly larger than the observed. Such large lithium abundance for 16 Cyg A would also make the star somewhat lithium-rich compared to stars of similar mass and metallicity (the top panel in Figure 11). Already the mass effect explains most, if not all, of the observed difference. The work by Deliyannis et al. (2000) on beryllium abundances also provides no evidence for accretion of planetary material by 16 Cyg A, although the authors warn that the errors are too large to make a firm conclusion.

A different possibility to explain the elemental abundance differences in the 16 Cyg system is that it was produced by the formation of the planet around 16 Cyg B. Similar to M09, we propose that the convective zone of 16 Cyg B has been depleted of metals due to the formation of a planet, i.e., that the missing mass of metals in 16 Cyg B is currently locked up in its planet.

Before exploring this hypothesis in detail, we should re-address the problem of lithium abundance. A number of authors have suggested that lower lithium abundances could be correlated with the presence of gas giants (e.g., Israelian et al. 2009; Gonzalez et al. 2010a), although their preferred explanation is that it relates to physical processes that affect the rotational evolution of stars with planets and hence their lithium depletion history. A few planet hosts and single stars are plotted in Figure 11. Age, mass, and metallicity effects fully explain the observed abundances of both types of stars, without invoking the presence or absence of gas giant planets, consistent with the works by, for example, Baumann et al. (2010) and Ghezzi et al. (2010). This is the case also for the 16 Cyg stars, i.e., the significantly lower lithium abundance of 16 Cyg B (≃ 0.6 dex relative to 16 Cyg A) is not due to the presence of its planet. In our hypothesis, however, a very small fraction of this difference (≃ 0.04 dex) could be due to the formation of the gas giant planet around 16 Cyg B. However, the effect is not related to enhanced lithium depletion in the photosphere of 16 Cyg B but to the fact that the planet around 16 Cyg B retained a certain amount of lithium that would have otherwise been accreted by the star.

As shown in Figure 9, the 16 Cyg A–B abundance differences do not correlate with dust condensation temperature as in M09. If we force the volatiles to match the mean abundance pattern of solar twins, the refractories are too low by 0.06 dex and vice versa (see solid lines in Figure 9). If the metallicity difference between 16 Cyg A and B is due to the formation of the planet around 16 Cyg B, we must assume that the abundance ratios of metals in this planet are very similar to those in the star 16 Cyg B. In particular, we have to assume that the volatile to refractory abundance ratios are equal to those in the photosphere of 16 Cyg B (excluding H and He, of course). Thus, the planet must be a gas giant, as indeed observed.

We explore our hypothesis quantitatively with a toy model. If the planet around 16 Cyg B were to be added to the convective zone of the star, the stellar metallicity would change by

$$\begin{equation} \Delta \mathrm{[M/H]}=\log \left[\frac{(Z/X)_\mathrm{CZ}M_\mathrm{CZ}+(Z/X)_\mathrm{p}M_\mathrm{p}}{(Z/X)_\mathrm{CZ}M_\mathrm{CZ}}\right], \end{equation} \tag{ 1 }$$

where (Z/X)_CZ is the ratio of the fractional abundance of metals (Z) relative to hydrogen (X) in the unperturbed convective zone (loosely called "metallicity" in this context), M_CZ is the mass of the convective zone, (Z/X)_p is the metallicity of the planet, and M_p is the mass of the planet. According to Asplund et al. (2009), (Z/X)^☉_CZ = 0.018. Thus, for 16 Cyg B ([M/H] = 0.061) we adopt (Z/X)_CZ = 10^0.061 × 0.018 = 0.021.

We can use Equation (1) to calculate the mass of the convective envelope necessary to explain the observed metallicity difference of Δ[M/H] = +0.04 given the planetary parameters (Z/X)_p and M_p. Since a minimum mass estimate of 1.5M_J is available from the work by Cochran et al. (1997), we explore Equation (1) for a range of planet masses from 1.5 to 9.5 M_J. The results from our toy model are illustrated in Figure 12.¹⁴

Zoom In Zoom Out Reset image size

This toy model is intended to counteract the immediate aftermath of the planet formation process. Thus, the mass of the convective envelope that appears in Equation (1) corresponds to the convective zone at the time of planet formation. Standard models of early stellar evolution suggest that solar analog stars begin fully convective and only after about 30 Myr do their convective envelopes reach their present-day masses (e.g., D'Antona & Mazzitelli 1994; Serenelli et al. 2011). On the other hand, observations of disks around young stars suggest a shorter timescale for the process of planet formation, of order 10 Myr (e.g., Mamajek 2009; Meyer 2009; Fedele et al. 2010). Thus, we investigate our simple model keeping M_CZ as a variable quantity but constraining its temporal evolution as computed by Serenelli et al. (2011, their Figure 1). The latter corresponds to a 1 M_☉ star of solar metallicity. To take into account the small differences in mass and [Fe/H] for 16 Cyg B, we scaled the Serenelli et al. (2011) result using the present-day convective envelope masses as a function of T_eff, [Fe/H] by Pinsonneault et al. (2001). The dashed lines plotted in Figure 12 correspond to the mass of the convective zone at different stellar ages according to standard models. In this figure we also show, with dotted lines, the time evolution of M_CZ for one of the non-standard models with episodic accretion by Baraffe & Chabrier (2010). We discuss this case later in this section.

According to our pollution model with standard M_CZ evolution, if the planet has (Z/X)_p = 0.1 and M_p = 1.5 M_J, it should have formed when the star was about 25 Myr old but if the planet has M_p = 9.5 M_J, it should have formed 10 Myr after the star was born, when the convective envelope was more massive, therefore diluting more the larger amount of metal-depleted gas accreted. The planet around 16 Cyg B could not have formed 25 Myr after the stellar birth since the convective envelope at that point was too small. Only if the metallicity of the planet is significantly lower (<0.03) could it have formed at these late times. On the other hand, if the planet is more metal-rich, for example (Z/X)_p = 0.7, it must have formed in the first 10 Myr since the birth of the star, and earlier for larger planet masses. Such large planet metallicity appears, however, unphysical for gas giants. The amount of metals in Jupiter (core and envelope) is about 12–37 M_⊕ (Guillot 2005; Fortney & Nettelmann 2010) and it seems that all giant planets have at least 10–15 M_⊕ of metals (Miller & Fortney 2011). Since Jupiter's mass is 318 M_⊕, a total amount of metals of 12–37 M_⊕ would imply (Z/X)_p = 0.04–0.12.

Theoretical models of planet formation predict that gas giant planets form relatively quickly, on timescales shorter than 10 Myr (e.g., Pollack et al. 1996; Guilera et al. 2011). Yet, for reasonable values of the planet's metallicity and mass, planet formation timescales from our pollution model with standard M_CZ evolution are in the 10–30 Myr range, i.e., we require the process of gas giant planet formation to be significantly longer. The exact same problem is faced by the M09 scenario for the formation of terrestrial planets. They have to assume that while the terrestrial planets were being formed, the mass of the solar convective envelope had already reached its present size. A possible solution to this problem may lie in the new developments in the field of early stellar evolution. In particular, the impact of episodic accretion on the internal structure of stars may help to solve this discrepancy.

Models for the early evolution of stars with non-steady accretion have been computed by Baraffe & Chabrier (2010). They demonstrate that the internal evolution of stars is highly dependent on the accretion history. Of particular interest for our work is the fact that their accretion models predict higher central temperatures and hotter overall structures leading to a more rapid formation of the radiative core in solar-type stars. In one of their models, the mass of the convective envelope of a 1 M_☉ star reaches its present-day size after only about 10 Myr from the stellar birth line. The evolution of M_CZ for this model in particular is represented with the dotted lines in Figure 12.

One possibility to explain the metallicity difference between 16 Cyg A and B in the context of gas giant planet formation, and at the same time remain compatible with the terrestrial planet signature hypothesis by M09, consists in assuming that the gas giants formed during the first ≃ 5 Myr after the birth of the star whereas the terrestrial planets formed later (≳ 10 Myr). More specifically, the timescales of accretion of metal-deficient and fractionated gas should be consistent with this assumption. In addition, the evolution of M_CZ should be quicker than predicted by standard models, i.e., M_CZ should reach its present-day mass in about 10 Myr, as is the case of some of the Baraffe & Chabrier (2010) models.

An interesting consequence of the proposed scenario is that the solar nebula would have been more metal-rich than the present-day solar photosphere. In fact, we could re-construct the original solar abundances if we had a good knowledge of the composition of the gas giants. The combined mass of the gas giants in the solar system is about 1.4 M_J, slightly lower than the minimum mass of the planet around 16 Cyg B. On the other hand, the convective envelope of the Sun is somewhat less massive than that of 16 Cyg B. Given the uncertainties in the many parameters involved to calculate the metallicity shift due to the formation of the gas giants in the solar system, we will assume for simplicity that their impact was similar to the case of 16 Cyg B. Thus, the formation of the gas giants in the solar system could have shifted the solar abundances of all elements down by about 0.04 dex. Later, the formation of rocky planets lowered even more the abundances of refractory elements by up to 0.08 dex while keeping the volatile abundances nearly constant (cf. M09, R09). Thus, for example, the abundances of some important elements in the solar nebula would have been [C, N, O/H] ≃ +0.05, [Fe/H] ≃ +0.09, and [Al, Sc, Ti/H] ≃ 0.12. Note that in this calculation we have ignored the impact of diffusion, which would further decrease the present-day-observed solar abundances by about 0.05 dex (Turcotte & Wimmer-Schweingruber 2002). Naturally, if the mass of the planet around 16 Cyg B is significantly larger than 1.5 M_J, the impact of the formation of Jupiter on the solar abundances would be lower than 0.04 dex.

The elemental abundances measured in CM chondrite meteorites, which are often assumed to have retained the solar nebula composition, agree very well with those measured in the present-day solar photosphere (e.g., Lodders 2003; Asplund et al. 2009). The impact of planet formation and diffusion on the solar abundances described above is not necessarily in contradiction with this observation because when these comparisons are made the meteoritic abundances have to be normalized using the solar abundance of Si (or a mixture of several metals with low abundance uncertainties), therefore removing any overall metallicity offsets. The trend with condensation temperature, however, would not disappear in this way but it has not yet been detected, probably owing to the still relatively large errors in the determination of absolute solar elemental abundances.

According to our proposed scenario, the differences between the initial and present-day solar abundances are not negligible and they could be relevant, for example, in studies of GCE, chemical tagging of dynamical groups, and stellar interior modeling. Further work on the impact of planet formation on stellar abundances, both from the observational and theoretical sides, is therefore urgently needed.